Более эффективное наложение многоугольника или извлечение() растровых данных из пространственных линий

Question

У меня есть огромный набор данных из 1,5 миллиарда пространственных линий, которые я создал, используя все комбинации из 37 000 точек. Для каждой пространственной линии я хотел бы извлечь максимальное значение многоугольника (или растра - все быстрее), к которому касается линия. По существу это очень большое «пространственное соединение» в Arc lingo. Если оверлейные линии на многоугольном слое, выход будет максимальным значением пространственной линии во всех полях атрибутов - каждый из которых составляет один месяц одного года. Я включил также набор растровых данных, который был создан только с января 1990 года из многоугольного файла с разрешением ~ 30 м. Растр представляет собой альтернативный подход, который, как я думал, может сэкономить время. Растровые слои многоугольника & представляют собой большую пространственную область: примерно 30 км x 10 км. Имеются данные here. Набор данных пространственных линий, которые я включил в .zip, имеет только 9900 строк, выборочных выборок из всего набора данных в 1,5 миллиарда строк.

Первое чтение в данных

#polygons 

poly<-readShapePoly("ls_polys_bin",proj4string=CRS("+proj=utm +zone=21 +south +datum=WGS84 +units=m +no_defs")) 
poly$SP_ID<-NULL #deleting this extra field in prep for overlay 

#raster - this represents only one month (january 1990) 
    #raster created from polygon layer but one month only 

    raster.jan90<-readGDAL("rast_jan90.tif") 
    raster.jan90<-raster(raster.jan90) #makes it into a raster 

#lines (9900 of 1.5 billion included) 

    lines<-readShapeLines("l_spatial",proj4string=CRS("+proj=utm +zone=21 +south +datum=WGS84 +units=m +no_defs")) 

Для того, чтобы данные линии более управляемым, возьмите образец 50 линий

lines.50<-lines[sample(nrow(lines),50),] 

земля все три слоя вместе

plot(raster.jan90)#where green=1 
plot(poly, axes=T,cex.axis=0.75, add=T) 
plot(lines.50, col="red", add=TRUE) 

Сначала я попробовал оверлей, но по текущему курсу весь набор данных в 1,5 миллиарда займет около 844 дней. ип на моей машине

ptm <- proc.time() #start clock 
overlays.all<-over(lines.50,poly, fn=max) 
ptm.sec.overlay<-proc.time() - ptm # stop clock 
ptm.sec.overlay #.56 sec w/ n=12 lines; 2.3 sec w/ 50 lines 

Далее я преобразовал полигоны растра (только один месяц - январь 1990), и я побежал экстракт() с пространственным линиям, но это заняло еще больше времени.

ptm <- proc.time() # Start clock 
ext.rast.jan90<-extract(raster.jan90,lines.50, fun=max, method=simple) 
ptm.sec.ext<-proc.time() - ptm # stop clock 
ptm.sec.ext #32 sec w/ n=12 lines; 191 sec w/ n=50 lines 

Мои попытки конвертировать все ячейки «0» в «NA» не помогли сэкономить время. Есть ли другой способ сделать это чудовищное наложение или экстракт() более эффективно? Обратите внимание, что эти данные в настоящее время закодированы как «1» или «0», но в конце концов я хотел бы запустить этот код для непрерывной переменной, которая работает 0: 300.

Answer 1

Здесь взломать, что должно дать хорошее приближение. Вероятно, это может быть улучшено (getCrds занимает много времени), в том числе, делая более крупные шаги (независимо от того, хорошо ли это для вас или нет, я не знаю).

library(raster) 
raster.jan90 <- raster("rast_jan90.tif") 
lines <- shapefile("l_spatial.shp", p4s="+proj=utm +zone=21 +south +datum=WGS84 +units=m +no_defs") 
lines.50<-lines[sample(nrow(lines),50),] 

test <- function(lns) { 

    getCrds <- function(i) { 
    p <- z[[i]][[1]] 
    s <- (p[2,] - p[1,])/res(raster.jan90) 
    step <- round(max(abs(s))) 
    if (step < 1) { 
     # these probably should not exist, but they do 
     return(cbind(i, cellFromXY(raster.jan90, p[1, , drop=FALSE]))) 
    } 
    x <- seq(p[1,1], p[2,1], length.out=step) 
    y <- seq(p[1,2], p[2,2], length.out=step) 
    cbind(i, unique(cellFromXY(raster.jan90, cbind(x, y)))) 
    } 

    z <- coordinates(lns) 
    crd <- sapply(1:length(z), getCrds) 
    crd <- do.call(rbind, crd) 

    e <- extract(raster.jan90, crd[, 2]) 
    tapply(e, crd[,1], max) 
} 

system.time(res <- test(lines.50)) 
# user system elapsed 
# 0.53 0.01 0.55 

system.time(res <- test(lines)) 
# user system elapsed 
# 59.72 0.85 60.58 

(684481500 * 60,58/длина (строки))/(3600 * 24) составляет около 50 дней ...

только 1 день на 50 компьютерах

Обратите внимание, что он получает относительно более эффективнее с большим количеством строк (так как относительно меньше уникальных ячеек, которые нужно запросить).

Answer 2

Я думаю, что самым быстрым способом сделать это было бы растрирование линий на тот же растре, что и ваши растровые данные.

Однако я бы не растрировал их в R. Я бы написал код С, который взял данные для растровых и 37 000 точек, а затем использовал алгоритм линейного рисования Bresenham, чтобы получить растровые местоположения линия. Попробуйте растр в этих местах и ​​сделайте все, что вам нужно, с этими данными. Быстрый код для алгоритма Брешенема должен быть легко доступен, и вы можете даже найти версии для работы на графических процессорах для массового ускорения. Что быстрее при рисовании прямых линий, чем в видеокарте?

Я сделал предположение, что ваши пространственные линии представляют собой отдельные прямые линии между двумя точками.

В качестве альтернативы, вы можете арендовать 1000 серверов у Amazon (или некоторых других поставщиков облачных услуг) на полдня.