通过R中的属性对SpatialPolygonsDataFrame(即删除多边形)进行子集化的简单方法

Question

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

例如,假设我要删除此世界shapefile中面积小于30000的所有多边形.我将如何进行此操作？

或者,同样,我如何删除Antartica？

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)

任何帮助表示赞赏!

Answer 1

看起来你正在覆盖数据,但不会删除多边形.如果要减少包括数据和多边形的数据集,请尝试例如

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, ] 如果有帮助的话,Upvote @ Bonnie的回答.

Answer 2

当我在R 3.2.1中尝试这样做时,上面的tim riffe技术对我来说不起作用,虽然修改它稍微修复了问题.我发现在指定属性到子集之前我必须专门引用数据槽,如下所示:

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

如果其他人遇到同样的问题,请添加此替代答案.

Answer 3

只是提到这subset也使得工作避免在条件中写入数据的名称.

world.map <- subset(world.map, AREA > 30000)
plot(world.map)

Answer 4

我使用上述技术制作了澳大利亚的地图:

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")

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