使用gIntersection循环多边形裁剪

debugcn 发表于 Dev

菲尔

问题

我正在尝试使用gIntersection函数和R在循环中执行多个多边形剪辑。我可以获取正确的剪辑并手动重新输入数据（这样我就可以将生成的SpatialPolygons对象变回SpatialPolygonsDataFrame对象）。我不能做的是使此循环与for()或一起工作apply()。

目前，这不是问题。我有9个英语地区（伦敦地区），因此手动设置每个剪辑并不是一个很大的挑战。但是，我最终希望在LAD中剪辑LSOA，这实际上意味着要设置> 400个剪辑。

所以，我的问题是，如何将手动剪辑变成一个工作循环？

最小的可复制示例

为了简单起见，让我们使用英语区域（n = 9）。对于9个区域中的每个区域，我将对各个县进行裁剪。以下代码加载适当的shapefile并将它们重新投影为British National Grid：

require(rgdal)
require(rgeos)
# English counties shapefile (~ 10MB zipped)
download.file(
  "https://dl.dropboxusercontent.com/s/6o0mi28pjo1kh9k/england-counties.zip",
  "ec", method = "wget")
unzip("ec")
ec <- readOGR("england-counties", "england_ct_2011")
proj4string(ec) <- CRS("+init=epsg:27700")

# English regions (~6MB zipped)
download.file(
  "https://dl.dropboxusercontent.com/s/p69m0vk2fh4xe3o/england-regions-2011.zip",
  "er", method = "wget")
unzip("er")
er <- readOGR("england-regions-2011", "England_gor_2011")
proj4string(er) <- CRS("+init=epsg:27700")

您应该剩下两个对象，er（英语区域）和ec（英语县）。两者都是SpatialPolygonsDataFrame对象。

以第一个区域（英格兰东部）为例，E12000006剪裁各县，然后将结果返回到SpatialPolygonsDataFrame对象中：

ee <- gIntersection(ec, er[er$CODE == "E12000006", ],
                    byid = T, drop_not_poly = T)
row.names(ee) <- as.character(gsub(" 0", "", row.names(ee)))
# gIntersection adds ' 0' to each row.name?
ee <- SpatialPolygonsDataFrame(ee, ec@data[row.names(ee), ])

的图ee确认此方法有效：

英格兰东部地区与县的情节

如您所见，这是仅适用于几种形状的不错的工作流程，但我确实想遍历所有区域，并最终遍历更多多边形。

我尝试过的

我对apply()循环不是很好，所以到目前为止，我尝试过的是一个for()循环（我知道这比较慢，但比输入所有内容还快！）：

regions <- as.character(er$CODE)  # length = 9 as expected
for(i in 1:length(regions)){
  as.name(paste0(regions[i], "c")) <- 
    gIntersection(ec, er[er$CODE == regions[1], ], byid = T, drop_not_poly = T)
}

而不是预期的行为，我得到以下错误：

Error in as.name(paste0(regions[1], "c")) <- gIntersection(ec, er[er$CODE ==  : 
could not find function "as.name<-"

我还尝试将对象名称包装到中，eval()但出现以下错误：

Error in eval(as.name(paste0(regions[1], "c"))) <- gIntersection(ec, er[er$CODE ==  :
could not find function "eval<-"

我想念什么？

除了gIntersection，如果可能的话，我想重新创建一个SpatialPolygonsDataFrame对象。我尝试了以下代码，手动完成了一个gIntersection，但再次无法正常工作：

for(i in 1:length(regions)){
   row.names(as.name(paste0(regions[i], "c"))) <- as.character(gsub(" 0", "",
     row.names(as.name(paste0(regions[i], "c")))))
}

我收到以下错误：

Error in `rownames<-`(x, value) : 
  attempt to set 'rownames' on an object with no dimensions

我也不能确定如何递增" 0"通过为每个新的区域，因为这将增加（" 1"，" 2"，等）

再次，手动设置第一个示例，我也无法执行最后的SpatialPolygonsDataFrame步骤：

for(i in 1:length(regions)){    
    as.name(regions[i]) <- SpatialPolygonsDataFrame(regions[i],
      ec@data[row.names(regions[i], )])
}

为此，我得到以下错误：

Error in stopifnot(length(Sr@polygons) == nrow(data)) :
trying to get slot "polygons" from an object of a basic class ("character") with no
slots

我看过的地方

以下SO示例与之关联似乎并没有帮助，或者至少我看不到如何使它们适用于此示例：

感谢您抽时间阅读。

约翰尼斯

这有帮助吗？

ee <- lapply(regions, function(x)  
  gIntersection(ec, er[er$CODE == x, ], byid = TRUE, drop_not_poly = TRUE))

这为您提供了SpatialPolygonsDataFrames的列表，每个区域一个。您可以用通常的方式访问它，例如

ee[[1]]
plot(ee[[1]])  # to plot the first region with counties

编辑

您的原始代码应稍作修改即可工作（请参见打击）。

res <- list()
for (i in 1:length(regions)) {
  ee <- gIntersection(ec, er[er$CODE == regions[i], ],
                      byid = TRUE, drop_not_poly = TRUE)
  row.names(ee) <- as.character(gsub(paste0(" ", i-1), "", row.names(ee)))
  ee <- SpatialPolygonsDataFrame(ee, ec@data[row.names(ee), ])
  res[[i]] <- ee
}

如果这样可以解决问题，那么问题就在于，行名ee总是增加一，而您没有考虑到这一点。

本文收集自互联网，转载请注明来源。

如有侵权，请联系[email protected] 删除。