r - g 非常大的空间对象上的交集

Question

长期读者，第一次发帖。

我正在尝试对两个非常大的 SpatialPolygonsDataFrame 对象执行 gIntersection()。第一个是所有美国县，第二个是一个 240 行 x 279 列的网格，作为一系列 66,960 个多边形。

我仅使用宾夕法尼亚州和与 PA 重叠的那块网格就成功地运行了它:
gIntersection(PA, grid, byid=TRUE)
我试图在整个美国通宵运行它，但今天早上它仍在运行，我的硬盘驱动器上有 10 GB(!)交换文件，并且没有任何进展迹象。我做错了什么，还是这是正常行为，我应该做一个逐个状态的循环？

谢谢!

Answer 1

比我希望的要晚一点，但这是我最终用于与此相关的任务的功能。它可能适用于其他应用程序。

@mdsumner 是对的，丢弃不相交的高级操作大大加快了这一过程。希望这是有用的!

library("sp")
library("rgeos")
library("plyr")

ApportionPopulation <- function(AdminBounds, poly, Admindf) { # I originally wrote this function to total the population that lies within each polygon in a SpatialPolygon object. AdminBounds is a SpatialPolygon for whatever administrative area you're working with; poly is the SpatalPolygon you want to total population (or whatever variable of your choice) across, and Admindf is a dataframe that has data for each polygon inside the AdminBounds SpatialPolygon.

  # the AdminBounds have the administrative ID code as feature IDS. I set that up using spChFID()

  # start by trimming out areas that don't intersect

  AdminBounds.sub <- gIntersects(AdminBounds, poly, byid=TRUE) # test for areas that don't intersect
  AdminBounds.sub2 <- apply(AdminBounds.sub, 2, function(x) {sum(x)}) # test across all polygons in the SpatialPolygon whether it intersects or not
  AdminBounds.sub3 <- AdminBounds[AdminBounds.sub2 > 0] # keep only the ones that actually intersect

  # perform the intersection. This takes a while since it also calculates area and other things, which is why we trimmed out irrelevant areas first

  int <- gIntersection(AdminBounds.sub3, poly, byid=TRUE) # intersect the polygon and your administrative boundaries

  intdf <- data.frame(intname=names(int)) # make a data frame for the intersected SpatialPolygon, using names from the output list from int
  intdf$intname <- as.character(intdf$intname) # convert the name to character
  splitid <- strsplit(intdf$intname, " ", fixed=TRUE) # split the names
  splitid <- do.call("rbind", splitid) # rbind those back together
  colnames(splitid) <- c("adminID", "donutshpid") # now you have the administrative area ID and the polygonID as separate variables in a dataframe that correspond to the int SpatialPolygon.
  intdf <- data.frame(intdf, splitid) # make that into a dataframe
  intdf$adminID <- as.character(intdf$adminID) # convert to character
  intdf$donutshpid <- as.character(intdf$donutshpid) # convert to character. In my application the shape I'm using is a series of half-circles

  # now you have a dataframe corresponding to the intersected SpatialPolygon object

  intdf$polyarea <- sapply(int@polygons, function(x) {x@area}) # get area from the polygon SP object and put it in the df
  intdf2 <- join(intdf, Admindf, by="adminID") # join together the two dataframes by the administrative ID
  intdf2$popinpoly <- intdf2$pop * (intdf2$polyarea / intdf2$admin_area) # calculate the proportion of the population in the intersected area that is within the bounds of the polygon (assuming the population is evenly distributed within the administrative area)
  intpop <- ddply(intdf2, .(donutshpid), summarize, popinpoly=sum(popinpoly)) # sum population lying within each polygon

  # maybe do other final processing to get the output in the form you want

  return(intpop) # done!
}