如何在R中通过属性从一个shapefile中消除某些区域并创建一个新的shapefile？

Question

我已经在这方面努力了一段时间并取得了一些成功，但是当涉及到重写文件时，我却没有成功。

我正在使用一个多边形状的地图文件，其中包含都市区和小都市区，但我只对都市区感兴趣，所以我正在使用R语言来尝试从我的地图中消除小都市区。

数据来源

要下载正确的文件，请确保在“都市/小都市统计区（2010）”下选择“所有州的一个国家文件”。

以下是我目前的进展：

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

在子集之前：

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

然后：

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

子集后:

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

那么：

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

使用这段代码，我已成功地消除了被指定为“M1”的小都市地区，但当我尝试重写文件时，它不会出现在我的工作目录中。

我还尝试使用“shapefiles”包来完成此操作，但是效果更差。因此，在这个问题上任何帮助将不胜感激。

- Dan Johnson

---

你的数据链接无效。 - Spacedman

你可以通过运行list.files(".")并将结果发布来向我们证明它不在你的“wd”中。我只是不相信它不存在，我认为你把R的当前目录与你的“wd”弄混了。 - mdsumner

好的，我收回之前的话，Spacedman已经证明了它。 :) - mdsumner

我已经向地图工具维护者报告了这种行为！ - Spacedman

Roger已经添加了一个错误信息，并将很快发布新版本到CRAN。你不喜欢开源吗？ - Spacedman

Answer 1

当你执行met=met@data时，你失去了met的空间方面，只得到了一个普通的数据框。

然后，你对这个数据框进行子集操作，接着尝试使用writeSpatialShape函数。哇！submet对象没有任何多边形或坐标。此时writeSpatialShape应该会报错。但它竟然把错误信息吞噬了下去。看看这个例子：

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

正如你所观察到的，没有创建任何shapefile。哇，这种错误处理真的很糟糕。

有两个解决方案：1、直接在实际对象上工作而不是数据组件上，像这样：

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

2、使用package:rgdal和read/write|OGR，这不仅可以处理投影，而且如果您试图让它吃下一些不可口的东西，它会非常生气：

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

- Spacedman

---

感谢您的快速回复和帮助。我遵循了您的第二个建议，使用了rgdal，并获得了我想要的结果。我甚至成功在ArcGIS中打开了我的新shapefile文件，而且完美无缺。 - Dan Johnson