78

我想根据 @data 数据框中的相应属性值从 SpatialPolygonsDataFrame 对象中简单地删除一些多边形，以便我可以绘制简化/子集的 shapefile。到目前为止，我还没有找到一种方法来做到这一点。

例如，假设我想从这个世界 shapefile中删除所有面积小于 30000 的多边形。我该怎么做呢？

或者，类似地，我怎样才能删除南极洲？

require(maptools)

getinfo.shape("TM_WORLD_BORDERS_SIMPL-0.3.shp") 
# Shapefile type: Polygon, (5), # of Shapes: 246
world.map <- readShapeSpatial("TM_WORLD_BORDERS_SIMPL-0.3.shp")

class(world.map)
# [1] "SpatialPolygonsDataFrame"
# attr(,"package")
# [1] "sp"

head(world.map@data)
#   FIPS ISO2 ISO3 UN                NAME   AREA  POP2005 REGION SUBREGION     LON     LAT
# 0   AC   AG  ATG 28 Antigua and Barbuda     44    83039     19        29 -61.783  17.078
# 1   AG   DZ  DZA 12             Algeria 238174 32854159      2        15   2.632  28.163
# 2   AJ   AZ  AZE 31          Azerbaijan   8260  8352021    142       145  47.395  40.430
# 3   AL   AL  ALB  8             Albania   2740  3153731    150        39  20.068  41.143
# 4   AM   AM  ARM 51             Armenia   2820  3017661    142       145  44.563  40.534
# 5   AO   AO  AGO 24              Angola 124670 16095214      2        17  17.544 -12.296

如果我做这样的事情，情节不会反映任何变化。

world.map@data = world.map@data[world.map@data$AREA > 30000,]
plot(world.map)

如果我这样做，结果相同：

world.map@data = world.map@data[world.map@data$NAME != "Antarctica",]
plot(world.map)

任何帮助表示赞赏！

4

5 回答 5

74

看起来您正在覆盖数据，但没有删除多边形。如果您想减少包括数据和多边形的数据集，请尝试例如

world.map <- world.map[world.map$AREA > 30000,]
plot(world.map)

[[编辑 2016 年 4 月 19 日]] 该解决方案曾经有效，但 @Bonnie 报告了较新的 R 版本（尽管数据可能也发生了变化？）： world.map <- world.map[world.map@data$AREA > 30000, ] 如果有帮助，请支持 @Bonnie 的回答。

于 2012-11-18T19:31:00.587 回答
42

当我尝试在 R 3.2.1 中执行此操作时，上述 tim riffe 的技术对我不起作用，尽管稍微修改它可以解决问题。我发现在将属性指定为子集之前，我还必须专门引用数据槽，如下所示：

world.map <- world.map[world.map@data$AREA > 30000, ]
plot(world.map)

如果其他人遇到相同的问题，请将其添加为替代答案。

于 2015-08-25T00:04:06.987 回答
16

顺便提一下，这subset也使得工作避免在条件中写入数据的名称。 

world.map <- subset(world.map, AREA > 30000)
plot(world.map)
于 2016-08-05T15:11:00.043 回答
11

我使用上述技术制作了一张澳大利亚地图：

australia.map < - world.map[world.map$NAME == "Australia",]
plot(australia.map)

事实证明，“澳大利亚”之后的逗号很重要。

这种方法的一个缺陷是它似乎保留了所有其他国家/地区的所有属性列和行，并且只是用零填充它们。我发现如果我写出一个 .shp 文件，然后使用 readOGR（rgdal 包）将其读回，它会自动删除空地理数据。然后我可以只用我想要的数据编写另一个形状文件。

writeOGR(australia.map,".","australia",driver="ESRI Shapefile")
australia.map < - readOGR(".","australia")
writeOGR(australia.map,".","australia_small",driver="ESRI Shapefile")

至少在我的系统上，它是删除空数据的“读取”函数，所以我必须在读回一次后写入文件（如果我尝试重新使用文件名，我会收到错误消息）。我确信有一种更简单的方法，但这似乎对我的目的来说已经足够好了。

于 2013-08-13T01:36:30.243 回答
1

作为第二个指针：这不适用于形状中带有“孔”的 shapefile，因为它是按索引进行子集化的。

于 2017-01-30T23:20:41.540 回答