logo标签列表

使用raster包中的intersect存在内存（RAM）问题

rr-rasterr-spr-sf
4
我在R中遇到了两个大型SpatialPolygonsDataFrame的交集问题。我的多边形数据代表建筑和行政边界，我正在尝试获取它们之间的交集多边形。
我知道raster包中的intersect函数和rgeos包中的gIntersection函数可以完成此任务（有一些区别），但它们无法同时处理我所有的多边形（大约50000个多边形/实体）。
因此，我必须在循环中分割我的计算，并保存每个步骤的结果。问题是：这些函数不断填充我的物理内存，而我无法清理它。我尝试使用rm()和gc()，但没有任何改变。内存问题会导致我的R会话崩溃，我无法进行计算。
是否有一种方法在模拟过程中释放RAM，在循环中？或避免这个内存问题？
下面是一个可重复的例子，使用随机多边形。
library(raster)
library(sp)
library(rgeos)

#Generating 50000 points (for smaller polygons) and 150000 (for larger polygons) in a square of side 100000
size=100000

Nb_points1=50000
Nb_points2=150000
start_point=matrix(c(sample(x = 1:size,size = Nb_points1,replace = T),sample(x = 1:size,size = Nb_points1,replace = T)),ncol=2)
start_point2=matrix(c(sample(x = 1:size,size = Nb_points2,replace = T),sample(x = 1:size,size = Nb_points2,replace = T)),ncol=2)

#Defining different sides length
radius=sample(x = 1:50,size = Nb_points1,replace = T)
radius2=sample(x = 1:150,size = Nb_points2,replace = T)

#Generating list of polygons coordinates
coords=list()
for(y in 1:Nb_points1){
  xmin=max(0,start_point[y,1]-radius[y])
  xmax=min(size,start_point[y,1]+radius[y])
  ymin=max(0,start_point[y,2]-radius[y])
  ymax=min(size,start_point[y,2]+radius[y])
  coords[[y]]=matrix(c(xmin,xmin,xmax,xmax,ymin,ymax,ymax,ymin),ncol=2)
}

coords2=list()
for(y in 1:Nb_points2){
  xmin=max(0,start_point2[y,1]-radius2[y])
  xmax=min(size,start_point2[y,1]+radius2[y])
  ymin=max(0,start_point2[y,2]-radius2[y])
  ymax=min(size,start_point2[y,2]+radius2[y])
  coords2[[y]]=matrix(c(xmin,xmin,xmax,xmax,ymin,ymax,ymax,ymin),ncol=2)
}

#Generating 75000 polygons
Poly=SpatialPolygons(Srl = lapply(1:Nb_points1,function(y) Polygons(srl = list(Polygon(coords=coords[y],hole = F)),ID = y)),proj4string = CRS('+init=epsg:2154'))
Poly2=SpatialPolygons(Srl = lapply(1:Nb_points2,function(y)Polygons(srl =  list(Polygon(coords=coords2[y],hole = F)),ID = y)),proj4string = CRS('+init=epsg:2154'))

#Union of overlapping polygons
aaa=gUnionCascaded(Poly)
bbb=gUnionCascaded(Poly2)

aaa=disaggregate(aaa)
bbb=disaggregate(bbb)

intersection=gIntersects(spgeom1 = aaa,bbb,byid = T,returnDense = F)

#Loop on the intersect function
pb <- txtProgressBar(min = 0, max = ceiling(length(aaa)/1000), style = 3)

for(j in 1:ceiling(length(aaa)/1000)){
  tmp_aaa=aaa[((j-1)*1000+1):(j*1000),]
  tmp_bbb=bbb[unique(unlist(intersection[((j-1)*1000+1):(j*1000)])),]
  List_inter=intersect(tmp_aaa,tmp_bbb)
  gc()
  gc()
  gc()
  setTxtProgressBar(pb, j)
}

感谢您！
1
为了避免内存问题，您可以切换到gdalUtils。 - loki
我不熟悉这个包。你能帮我了解一下吗？有什么函数可以帮助我吗？我没有看到关于内存或交集的相关内容。 - A.Rogeau
gdalUtils是一个非常好用且实用的软件包，但它在这里并没有帮助。它主要用于处理栅格数据。您使用了raster软件包，但不是用于栅格数据，因此我怀疑它是否有所帮助。 - Bastien
R对于大型GIS数据处理并不是很高效。我通常更喜欢使用R作为基础来调用其他软件。在这方面，RSAGA是我最喜欢的，其次是RQGIS，然后是更复杂的RGRASS7。所有这些都需要您安装适当的软件（可以通过OSGEO4W一次性完成）。它们应该能够成功地完成您的任务。我现在有点忙，如果以后有机会，我会发布一个示例。 - Bastien
2个回答
3
你可以考虑使用 sf 包中的 st_intersects 和 st_intersection 函数。例如：
aaa2 <- sf::st_as_sf(aaa)
bbb2 <- sf::st_as_sf(bbb)
intersections_mat <- sf::st_intersects(aaa2, bbb2)
intersections <- list()
for (int in seq_along(intersections_mat)){
  if (length(intersections_mat[[int]]) != 0){
    intersections[[int]] <- sf::st_intersection(aaa2[int,], 
    bbb2[intersections_mat[[int]],])
  }
}

将为您提供一个长度等于aaa的intersection_mat，其中包含针对aaa的每个特征与其相交的bbb元素的“索引”（如果未找到交集，则为空）：

> intersections_mat
Sparse geometry binary predicate list of length 48503, where the predicate was `intersects'
first 10 elements:
 1: 562
 2: (empty)
 3: 571
 4: 731
 5: (empty)
 6: (empty)
 7: (empty)
 8: 589
 9: 715
 10: (empty)

并且一个交集列表包含交叉多边形的列表：

>head(intersections)
[[1]]
Simple feature collection with 1 feature and 0 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 98873 ymin: 33 xmax: 98946 ymax: 98
epsg (SRID):    2154
proj4string:    +proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
                        geometry
1 POLYGON ((98873 33, 98873 9...

[[2]]
NULL

[[3]]
Simple feature collection with 1 feature and 0 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 11792 ymin: 3 xmax: 11806 ymax: 17
epsg (SRID):    2154
proj4string:    +proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
                        geometry
1 POLYGON ((11792 3, 11792 17...

(i.e., intersections[[1]] 表示 aaa 的第一个多边形和 bbb 的第571个多边形的交集)
希望对您有所帮助。
谢谢，这个包完美地运行了！我实际上直接在两个大的SpatialPolygonsDataFrame上使用了sf::intersection函数。如果交集包含不同的几何类型（例如点、线和多边形），需要通过类型来“拆分”intersection_mat。在将它们传递回作为sp包中的空间对象时，请使用sf::st_geometry_type，并用as( ,"Spatial")来表示意思。 - A.Rogeau
很高兴能帮到你。我建议使用“双重通道”，因为我无法从st_intersection的输出中找到一种简单的方法来理解不同输出多边形的“发起者”是谁。然而，我现在注意到帮助文档中提到：“返回的sfc几何列表列包含一个属性idx，它是一个n×2矩阵，每行分别是x和y对应条目的索引”。 - lbusett
因此，对interscection_mat的循环是无用的（正如您已经看到的），直接调用st_intersection即可。如果您希望，我可以修改回复。 - lbusett
1

这个例子对于我来说运行良好（8 GB RAM），在循环中进行了一些更改。请参见下文。这些更改与内存使用无关---您没有存储结果。

List_inter <- list()

for(j in 1:ceiling(length(aaa)/1000)){
    begin <- (j-1) * 1000 + 1
    end <- min((j*1000), length(aaa))
    tmp_aaa <- aaa[begin:end,]
    tmp_bbb <- bbb[unique(unlist(intersection[begin:end])),]
    List_inter[[j]] <- intersect(tmp_aaa,tmp_bbb)
    cat(j, "\n"); flush.console()
}

x <- do.call(bind, List_inter)

或者，您可以将中间结果写入磁盘，并稍后处理：

inters <- intersect(tmp_aaa,tmp_bbb)
saveRDS(inters, paste0(j, '.rds'))

或者
shapefile(inters, paste0(j, '.shp'))
谢谢你的回答。我的RAM使用量仍然会随着这个新循环而增加。实际上，我在示例中没有存储结果。 - A.Rogeau
随着您创建新对象，RAM使用应该会增加。您可以尝试将中间结果写入磁盘，并稍后再回来处理它们。我已经添加了一些代码来实现这个功能。 - Robert Hijmans
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接