3

I use the package maptools to read in Shapefiles in R. I would like to know if there is a function that can built up a query/subset to select only some Polygons based on a specific attribute in the a.t. (for example one state out of a shapefile with all states in one country). And, can I save this subset as a spatial object to do further analyses in R? Best

1
  • 4
    Comment: don't use maptools for reading, use readOGR from package:rgdal - it will preserve any coordinate reference info from the shape file.
    – Spacedman
    May 12 '13 at 22:31
6

To summarize the very useful comments by @Spacedman, @JeffreyEvans and @AriB.Friedman, but also to address the questioner's second problem about saving the generated subset and to provide a reproducible example for future audience, let me add a brief example to cover the whole topic. In fact, I'm much more often working with raster data than point or polygon shapefiles, so this was also pretty useful for me.

As @Spacedman mentioned, you should rather use rgdal instead of maptools for the given reason.

library(rgdal)

Now, let's download a suitable ESRI shapefile containing all states of the US and extract it to a folder named 'data' in your current working directory. After that, you are ready to import the shapefile using readOGR().

usa.states <- readOGR(dsn = "data", layer = "states")

You can now create a subset of this SpatialPolygonsDataFrame either via accessing the @data slot (thanks @JeffreyEvans, I didn't know that...) or addressing a particular column solely via the $ operator. Let's take New York for our testing purposes:

newyork <- usa.states[usa.states$STATE_NAME == "New York", ]
# Or alternatively... :-)
newyork <- usa.states[usa.states@data$STATE_NAME == "New York", ]

Have a look at the generated subset via plot(newyork). Now that you know that everything worked, use writeOGR() to save the shapefile to the desired folder ('data' in our case). If you're not sure whether 'ESRI Shapefile' or any other driver is available on your local system, check ogrDrivers() prior to sending the writeOGR() statement to the console.

writeOGR(oregon, dsn = "data", layer = "oregon", 
         driver = "ESRI Shapefile", overwrite_layer = TRUE)

That's it. Hope that helps!

0
0

ok this was easier than i thought: spatial.subset<-states.shp[5,] if you know which line in your attribute table should be the subset. If you don't know you could do something like spatial.subset<-states.shp[states.shp@data$STATE=="California" ,] for a single state or spatial.subset<-states.shp[c(states.shp@data$POPULATION>1000000),] for several states with a population greater than 1.Mio.

6
  • 1
    This doesn't look like it's an attribute query. Perhaps you could elaborate on how it answers your question?
    – whuber
    May 12 '13 at 16:28
  • @whuber you are right. Please note the edits in the answer to make it more comprehensive.
    – Dspanes
    May 12 '13 at 20:40
  • 1
    Using "which" could have unexpected behavior and is quite unnecessary. You should also be accessing the data slot in the bracket index. This will get you what you want. states.shp[states.shp@data$STATE=="California" ,] May 13 '13 at 3:41
  • 2
    @JeffreyEvans Why is accessing the @data slot necessary? $ is overloaded for Spatial*DataFrame objects to return the named column of the @data slot, and thus reaching inside the object's structure would seem to be disfavored over using the appropriate method.... May 13 '13 at 11:46
  • 1
    Because calling the slot is much more version proof. Until very recently you could not call a column in the data slot without explicitly accessing the slot object. May 13 '13 at 17:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.