У меня проблемы с получением пересечения между двумя большими SpatialPolygonsDataFrame на R. Мои данные полигонов представляют здания и административные границы, и я пытаюсь получить полигоны пересечения между ними.
Я понимаю, что функция пересечения из растрового пакета и gIntersection из пакета rgeos может выполнять эту работу (с некоторыми отличиями), но они не могут обрабатывать все мои полигоны одновременно (около 50 000 полигонов на объект).
По этой причине мне приходится разбивать вычисления внутри цикла, сохраняя результат для каждого шага. Проблема в том, что эти функции продолжают заполнять мою физическую память, и я не могу ее очистить. Я пробовал использовать rm () и gc (), но это ничего не меняет. Проблема с памятью приводит к сбою моего сеанса R, и я не могу выполнить расчет.
Есть ли способ освободить оперативную память во время моделирования в циклах? Или чтобы избежать этой проблемы с памятью?
Вот воспроизводимый пример для случайных многоугольников.
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)
}
Спасибо !
gdalUtils
. - person loki   schedule 08.01.2018gdalUtils
- очень хороший и полезный пакет, но здесь он не поможет. В основном это игра с растрами. Вы используете растровый пакет, но не на растрах, поэтому я сомневаюсь, что это поможет. - person Bastien   schedule 09.01.2018RSAGA
- мой любимый, за ним следуетRQGIS
и более сложныйRGRASS7
. Все, что вам нужно, чтобы установить соответствующее программное обеспечение (можно сделать один снимок с OSGEO4W). Они должны преуспеть в вашей задаче. Я сейчас вроде как занят, если будет возможность позже, выложу пример. - person Bastien   schedule 09.01.2018