simple overlap feature extraction for US states

[jordansread]

Do not include territories for the time being. That could be an additional task in the future.

want to extract this automatically from spatial data:

<place>
        <placekt>Department of Commerce, 1995, Countries, Dependencies, Areas of Special Sovereignty, and 
                  Their Principal Administrative Divisions,  Federal Information Processing Standard (FIPS) 10-4, 
                  Washington, D.C., National Institute of Standards and Technology</placekt>
        <placekey>United States</placekey>
        <placekey>US</placekey>
      </place>
      <place>
        <placekt>U.S. Department of Commerce, 1987, Codes for the identification of the States, the District of Columbia and the outlying areas of the United States, and associated areas (Federal Information Processing Standard 5-2): Washington, D. C., NIST</placekt>
        <placekey>Wisconsin</placekey>
        <placekey>WI</placekey>
      </place>
      <place>
        <placekt>none</placekt>
        <placekey>Wisconsin</placekey>
      </place>

• [x] CONUS states
• [x] AK and HI

[jordansread]

An overly complex version of this exists, but I would like to be able to do it simply with the built in data from the maps package.

[berdaniera]

Here is an example using maps for the lower 48. Could be expanded for other regions. Don't know if it is faster than what you have already:

  library(rgdal)
  library(maps)
  library(maptools)
  library(rgeos)
  usa <- map("state", fill=T, plot=F) # get state bounds
  IDS <- sapply(strsplit(usa$names,":"), function(x) x[1]) # just keep state names (no sub-names)
  usa_sp <- map2SpatialPolygons(usa,IDs=IDS,proj4string=CRS(proj4string(sp))) # to sp
  overlaps <- sp_overlaps(usa_sp, sp)
  state.has.sp <- colnames(overlaps[,colSums(overlaps)>0])

[jordansread]

Awesome, thanks @berdaniera added that to upcoming PR

[berdaniera]

I can add AK and HI this afternoon. Any other locations you want?

[jordansread]

That's it. Great!

[berdaniera]

Ok, this should work to add ak and hi... seems to work on my local machine (but I haven't set up the test part yet...). Do you want to try it out?

us_48 <- map("state", fill=TRUE, plot=FALSE)
us_48$names <- paste0("USA:", us_48$names)  # because the other maps have a USA prefix
us_hi <- map("world", "USA:Hawaii", fill=TRUE, plot=FALSE)
us_ak <- map("world", "USA:Alaska", fill=TRUE, plot=FALSE)
usa <- Reduce(function(m1,m2){
 list(x=c(m1$x, NA, m2$x),
    y=c(m1$y, NA, m2$y),
    names=c(m1$names, m2$names))
}, list(us_48,us_ak,us_hi))
IDs <- sapply(strsplit(usa$names, ":"), function(x) x[2])  # now grabbing the second item
usa <- map2SpatialPolygons(usa, IDs=IDs, proj4string=CRS("+proj=longlat +datum=WGS84"))

[jordansread]

do you want to set up a PR with this change and also modify this line of test-feature_extraction.R (in tests/testthat/)

  expect_equal(state.points, c("Alaska", "Illinois", "New Mexico"))

? @berdaniera

[jordansread]

closed in #23