r - 使用来自光栅包的相交的内存 (RAM) 问题

Question

我无法在 R 上获取两个大型 SpatialPolygonsDataFrame 之间的交集。我的多边形数据表示建筑物和行政边界，我正在尝试获取它们之间的交集多边形。

我知道来自 raster 包的 intersect 函数和来自 rgeos 包的 gIntersection 可以完成这项工作(有一些差异)，但它们不能一次处理我所有的多边形(大约 50.000 个多边形/实体)。

出于这个原因，我必须在一个循环中拆分我的计算，保存每一步的结果。问题是:这些函数不断填满我的物理内存，我无法清理它。我尝试使用 rm() 和 gc()，但它没有改变任何事情。内存问题使我的 R session 崩溃，我无法进行计算。

有没有办法在模拟期间在循环内释放 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)
}

谢谢 !

Answer 1

您可以考虑使用 st_intersects和 st_intersection包的功能sf .例如:

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]],])
  }
}

会给你一个 intersection_mat长度等于 aaa ，并包含 ，对于 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)

, 和 intersection包含相交多边形列表的列表:

>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...

(即， intersections[[1]] 是 aaa 的多边形 1 与 bbb 的多边形 571 之间的交集)

哈。