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.

Answer 1

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