How can I eliminate some areas by attribute from a shapefile in R and create a new shapefile?

Question

I have been at this for a while now and have had some success, however when it comes to rewriting the file, I have had none.

The shapefile I am working with is a polygon shapefile of metro and micropolitan areas but I am not interested in the micropolitan areas so I am working with the shp in R to try and eliminate them from my map.

Data Source

To download the right file, make sure that you select "All States in one national file" under "Metropolitan/Micropolitan Statistical Area (2010)"

Here's what I have so far:

library(maptools)
met=readShapeSpatial("tl_2010_us_cbsa10.shp")
met=met@data

Before subset:

NAMELSAD10      LSAD10
Anchorage, AK Metro Area   -  M1 
Clarksdale, MS Micro Area    - M2
Richmond, VA Metro Area   -  M1
Big Spring, TX Micro Area  -   M2
Dallas-Fort Worth-Arlington, TX Metro Area  -   M1
Rio Grande City-Roma, TX Micro Area  -   M2

then:

submet=subset(met, LSAD10 == "M1")

After subset:

NAMELSAD10 LSAD10
Anchorage, AK Metro Area   -  M1
Richmond, VA Metro Area  -   M1
Dallas-Fort Worth-Arlington, TX Metro Area   -  M1
Vineland-Millville-Bridgeton, NJ Metro Area  -   M1
Casper, WY Metro Area   -  M1
Cheyenne, WY Metro Area   -  M1

then:

writeSpatialShape(submet, "tl_2010_us_ma10", factor2char = TRUE)

Using this code I have been able to successfully eliminate the Micropolitan areas designated "M1" but when I try to rewrite the file, it doesn't show up in my wd.

I have also attempted using the package "shapefiles" to get this done but had even less success. So any help in this matter would be greatly appreciated.

You can demonstrate to us that it's not in your "wd" by running list.files(".") and posting that too. I just don't believe it's not there, I reckon you've mistaken your "wd" with R's current directory. — mdsumner, Mar 07 '13 at 08:24
And I've now reported the behaviour to the maptools maintainer! — Spacedman, Mar 07 '13 at 08:31
And Roger has added an error message and will release a new version to CRAN soon. Dontcha LOVE open source? — Spacedman, Mar 07 '13 at 19:00

score 8 · Accepted Answer · answered Mar 07 '13 at 08:24

When you do met=met@data you are losing the Spatial aspect of met and just getting the plain old data frame.

Then you subset that data frame okay, and then you try and writeSpatialShape it. Whoa. The submet object doesn't have any polygons or coordinates. writeSpatialShape should barf. But it manages to swallow it whole and keep it down. Check this out:

> writeSpatialShape(1,"foo.shp")
>

and as you observed, no shapefile is created. Wow, that's pretty sucky error handling.

Two solutions: 1, work on the actual object rather than the data component, like this:

met = readOGR(dir,name)
submet = met[met$thing=="whatever",]
writeOGR(submet,dir,newname,"ESRI Shapefile")

2, use package:rgdal and read/write|OGR which not only handles projections but is violently sick if you try and make it eat something unpalatable:

> writeOGR(1,".","foo","ESRI Shapefile")
Error: inherits(obj, "Spatial") is not TRUE

Thanks for your quick response, and for your help. I followed your second recommendation and used rgdal instead and got my desired results. I even got my new shapefile opened in ArcGIS and it is perfect. — Dan Johnson, Mar 11 '13 at 00:52

How can I eliminate some areas by attribute from a shapefile in R and create a new shapefile?

1 Answers1