gpt4 book ai didi

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

转载 作者:行者123 更新时间:2023-12-03 23:16:23 25 4
gpt4 key购买 nike

长期读者,第一次发帖。

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

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

谢谢!

最佳答案

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

@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!
}

关于r - g 非常大的空间对象上的交集,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/16918767/

25 4 0
Copyright 2021 - 2024 cfsdn All Rights Reserved 蜀ICP备2022000587号
广告合作:1813099741@qq.com 6ren.com