Aumentar a velocidade de cortar, mascarar e extrair varredura por muitos polígonos no R?

29

Estou extraindo a área e a porcentagem de cobertura de diferentes tipos de uso da terra de uma varredura baseada em vários milhares de limites de polígonos. Descobri que a função extrair funciona muito mais rápido se eu percorrer cada polígono individual e cortar, mascarar a varredura até o tamanho do polígono específico. No entanto, é bem lento, e estou me perguntando se alguém tem alguma sugestão para melhorar a eficiência e a velocidade do meu código.

A única coisa que encontrei relacionada a isso é essa resposta de Roger Bivand, que sugeriu o uso de GDAL.open()e GDAL.close()como getRasterTable()e getRasterData(). Eu olhei para eles, mas tive problemas com o gdal no passado e não o conheço o suficiente para saber como implementá-lo.

Exemplo reproduzível:

library(maptools)  ## For wrld_simpl
library(raster)

## Example SpatialPolygonsDataFrame
data(wrld_simpl) #polygon of world countries
bound <- wrld_simpl[1:25,] #name it this to subset to 25 countries and because my loop is set up with that variable  

## Example RasterLayer
c <- raster(nrow=2e3, ncol=2e3, crs=proj4string(wrld_simpl), xmn=-180, xmx=180, ymn=-90, ymx=90)
c[] <- 1:length(c)

#plot, so you can see it
plot(c)    
plot(bound, add=TRUE) 

Método mais rápido até agora

result <- data.frame() #empty result dataframe 

system.time(
     for (i in 1:nrow(bound)) { #this is the number of polygons to iterate through
      single <- bound[i,] #selects a single polygon
      clip1 <- crop(c, extent(single)) #crops the raster to the extent of the polygon, I do this first because it speeds the mask up
      clip2 <- mask(clip1,single) #crops the raster to the polygon boundary

      ext<-extract(clip2,single) #extracts data from the raster based on the polygon bound
      tab<-lapply(ext,table) #makes a table of the extract output
      s<-sum(tab[[1]])  #sums the table for percentage calculation
      mat<- as.data.frame(tab) 
      mat2<- as.data.frame(tab[[1]]/s) #calculates percent
      final<-cbind(single@data$NAME,mat,mat2$Freq) #combines into single dataframe
      result<-rbind(final,result)
      })

   user  system elapsed 
  39.39    0.11   39.52 

Processamento paralelo

O processamento paralelo reduziu o tempo do usuário pela metade, mas negou o benefício dobrando o tempo do sistema. A varredura usa isso para a função de extração, mas infelizmente não para a função de corte ou máscara. Infelizmente, isso deixa uma quantidade ligeiramente maior de tempo total decorrido devido à "espera" do "IO".

beginCluster( detectCores() -1) #use all but one core

executar código em vários núcleos:

  user  system elapsed 
  23.31    0.68   42.01 

então termine o cluster

endCluster()

Método lento: o método alternativo de fazer uma extração diretamente da função raster demora muito mais e não tenho certeza sobre o gerenciamento de dados para obtê-lo da forma que desejo:

system.time(ext<-extract(c,bound))
   user  system elapsed 
1170.64   14.41 1186.14 
Luke Macaulay
fonte
Você pode tentar este criador de perfil de código R ( marcodvisser.github.io/aprof/Tutorial.html ). Pode dizer quais linhas demoram a maior parte do tempo. O link também tem diretrizes para reduzir o tempo de processamento em R.
J Kelly
Apenas meus dois centavos aqui. . . mas o método crop / getvalues ​​não funciona quando o número de pixels no corte é muito baixo. Não sei ao certo onde está o limite, mas tive problemas com as colheitas em que havia apenas 1 a 5 pixels (não determinei o motivo exato (ainda um pouco novo para os pacotes espaciais), mas aposto que a função de colheita depende limites de pixels, portanto, luta para cortar todos os pixels individuais). A extração do pacote raster não tem esse problema, mas o acordado é mais do dobro do tempo do usuário e muito mais do que duas vezes no tempo do sistema. Apenas um aviso para aqueles que têm rasters baixa resolução (e um em
Neal Barsch
2
Existe um pacote um pouco novo, o velox, que moveu a extração para o C através do pacote Rcpp. Ele está dando ~ 10 vezes aumentos de velocidade nas operações de extração usando polígonos.
Jeffrey Evans
@JeffreyEvans. Testando a resposta a esta pergunta usando o velox agora. Porém, tendo problemas ao alocar vetores extremamente grandes.
SeldomSeenSlim 3/17/17

Respostas:

23

Finalmente consegui melhorar esta função. Descobri que, para meus propósitos, era mais rápido rasterize()primeiro no polígono e em getValues()vez de usá-lo extract(). A rasterização não é muito mais rápida que o código original para tabular valores de varredura em pequenos polígonos, mas brilha quando se trata de grandes áreas de polígonos com grandes rasters a serem cortados e os valores extraídos. Também achei getValues()muito mais rápido que a extract()função.

Eu também descobri o processamento multi-core usando foreach().

Espero que isso seja útil para outras pessoas que desejam uma solução R para extrair valores raster de uma sobreposição de polígono. Isso é semelhante ao "Tabulate Intersection" do ArcGIS, que não funcionou bem para mim, retornando saídas vazias após horas de processamento, como este usuário.

#initiate multicore cluster and load packages
library(foreach)
library(doParallel)
library(tcltk)
library(sp)
library(raster)

cores<- 7
cl <- makeCluster(cores, output="") #output should make it spit errors
registerDoParallel(cl)

Aqui está a função:

multicore.tabulate.intersect<- function(cores, polygonlist, rasterlayer){ 
  foreach(i=1:cores, .packages= c("raster","tcltk","foreach"), .combine = rbind) %dopar% {

    mypb <- tkProgressBar(title = "R progress bar", label = "", min = 0, max = length(polygonlist[[i]]), initial = 0, width = 300) 

    foreach(j = 1:length(polygonlist[[i]]), .combine = rbind) %do% {
      final<-data.frame()
      tryCatch({ #not sure if this is necessary now that I'm using foreach, but it is useful for loops.

        single <- polygonlist[[i]][j,] #pull out individual polygon to be tabulated

        dir.create (file.path("c:/rtemp",i,j,single@data$OWNER), showWarnings = FALSE) #creates unique filepath for temp directory
        rasterOptions(tmpdir=file.path("c:/rtemp",i,j, single@data$OWNER))  #sets temp directory - this is important b/c it can fill up a hard drive if you're doing a lot of polygons

        clip1 <- crop(rasterlayer, extent(single)) #crop to extent of polygon
        clip2 <- rasterize(single, clip1, mask=TRUE) #crops to polygon edge & converts to raster
        ext <- getValues(clip2) #much faster than extract
        tab<-table(ext) #tabulates the values of the raster in the polygon

        mat<- as.data.frame(tab)
        final<-cbind(single@data$OWNER,mat) #combines it with the name of the polygon
        unlink(file.path("c:/rtemp",i,j,single@data$OWNER), recursive = TRUE,force = TRUE) #delete temporary files
        setTkProgressBar(mypb, j, title = "number complete", label = j)

      }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")}) #trycatch error so it doesn't kill the loop

      return(final)
    }  
    #close(mypb) #not sure why but closing the pb while operating causes it to return an empty final dataset... dunno why. 
  }
}

Portanto, para usá-lo, ajuste-o single@data$OWNERpara caber no nome da coluna do seu polígono de identificação (acho que poderia ter sido incorporado à função ...) e insira:

myoutput <- multicore.tabulate.intersect(cores, polygonlist, rasterlayer)
Luke Macaulay
fonte
3
A sugestão que getValuesfoi muito mais rápida que a extractnão parece válida, porque se você usar extract, não precisará fazer crope rasterize(ou mask). O código na pergunta original faz as duas coisas, e isso deve dobrar o tempo de processamento.
Robert Hijmans
a única maneira de saber é testando.
Dj #
Que classe é a lista poligonal aqui e o que a lista poligonal [[i]] [, j] deve fazer aqui (ELI5, por favor)? Eu sou novato em coisas paralelas, então não entendo isso muito bem. Não consegui que a função retornasse nada, até que mudei para polygonlist [[i]] [, j] para polygonlist [, j], o que parece lógico porque [, j] é o j-ésimo elemento de um SpatialPolygonsDataframe, se esse é a classe correta? Depois de mudar isso, consegui o processo em execução e algumas saídas, mas ainda há definitivamente algo errado. (Eu tento extrair o valor mediano dentro de n polígonos pequenos, então mudei um pouco do código também).
Reima #
@RobertH No meu caso, o corte (e máscara) faz com que seja executado cerca de três vezes mais rápido. Estou usando uma varredura de 100 milhões de acres e os polígonos são uma pequena fração disso. Se eu não cortar o polígono, o processo será muito mais lento. Aqui estão meus resultados: clip1 <- recorte (camada raster, extensão (única))> system.time (ext <-extract (clipe1, única)) #extraindo do sistema de usuário raster recortado decorrido 65,94 0,37 67,22> system.time (ext < -extract (rasterlayer, single)) #extraindo-se de um sistema de usuário de varredura de 100 milhões de acres decorrido 175,00 4,92 181,10 #
Luke Macaulay
4

Acelere a extração de varredura (pilha de varredura) do ponto, XY ou polígono

Ótima resposta Luke. Você deve ser um assistente de R! Aqui está um pequeno ajuste para simplificar seu código (pode melhorar um pouco o desempenho em alguns casos). Você pode evitar algumas operações usando cellFromPolygon (ou cellFromXY para pontos) e, em seguida, recorte e getValues.

Extrair dados de polígonos ou pontos de pilhas de varredura ------------------------

 library(raster)  
 library(sp)   

  # create polygon for extraction
  xys= c(76.27797,28.39791,
        76.30543,28.39761,
        76.30548,28.40236,
        76.27668,28.40489)
  pt <- matrix(xys, ncol=2, byrow=TRUE)
  pt <- SpatialPolygons(list(Polygons(list(Polygon(pt)), ID="a")));
  proj4string(pt) <-"+proj=longlat +datum=WGS84 +ellps=WGS84"
  pt <- spTransform(pt, CRS("+proj=sinu +a=6371007.181 +b=6371007.181 +units=m"))
  ## Create a matrix with random data & use image()
  xy <- matrix(rnorm(4448*4448),4448,4448)
  plot(xy)

  # Turn the matrix into a raster
  NDVI_stack_h24v06 <- raster(xy)
  # Give it lat/lon coords for 36-37°E, 3-2°S
  extent(NDVI_stack_h24v06) <- c(6671703,7783703,2223852,3335852)
  # ... and assign a projection
  projection(NDVI_stack_h24v06) <- CRS("+proj=sinu +a=6371007.181 +b=6371007.181 +units=m")
  plot(NDVI_stack_h24v06)
  # create a stack of the same raster
  NDVI_stack_h24v06 = stack( mget( rep( "NDVI_stack_h24v06" , 500 ) ) )


  # Run functions on list of points
  registerDoParallel(16)
  ptm <- proc.time()
  # grab cell number
  cell = cellFromPolygon(NDVI_stack_h24v06, pt, weights=FALSE)
  # create a raster with only those cells
  r = rasterFromCells(NDVI_stack_h24v06, cell[[1]],values=F)
  result = foreach(i = 1:dim(NDVI_stack_h24v06)[3],.packages='raster',.combine=rbind,.inorder=T) %dopar% {
     #get value and store
     getValues(crop(NDVI_stack_h24v06[[i]],r))
  }
  proc.time() - ptm
  endCluster()

sistema do usuário decorrido 16.682 2.610 2.530

  registerDoParallel(16)
  ptm <- proc.time()
  result = foreach(i = 1:dim(NDVI_stack_h24v06)[3],.packages='raster',.inorder=T,.combine=rbind) %dopar% {
        clip1 <- crop(NDVI_stack_h24v06[[i]], extent(pt)) #crop to extent of polygon
        clip2 <- rasterize(pt, clip1, mask=TRUE) #crops to polygon edge & converts to raster
         getValues(clip2) #much faster than extract
  }
  proc.time() - ptm
  endCluster()

sistema de usuário decorrido 33.038 3.511 3,288

mmann1123
fonte
Eu executei as duas abordagens e seu método saiu um pouco mais lento no meu caso de uso.
Luke Macaulay
2

Se uma perda na precisão da sobreposição não for muito importante - supondo que seja preciso começar - normalmente é possível obter velocidades de cálculo zonais muito maiores convertendo primeiro os polígonos em uma varredura. O rasterpacote contém a zonal()função, que deve funcionar bem para a tarefa pretendida. No entanto, você sempre pode subconjunto de duas matrizes da mesma dimensão usando a indexação padrão. Se você deve manter polígonos e não se importa com o software GIS, o QGIS é realmente mais rápido nas estatísticas zonais do que o ArcGIS ou o ENVI-IDL.

Adam Erickson
fonte
2

Também lutei com isso por algum tempo, tentando calcular a parcela da área das classes de cobertura do solo de um mapa de grade de ~ 300mx300m em uma grade de ~ 1kmx1km. O último era um arquivo de polígono. Tentei a solução multicore, mas ainda era muito lenta para o número de células da grade que eu tinha. Em vez disso eu:

  1. Rasterizou a grade de 1 km x 1 km, dando a todas as células da grade um número único
  2. Utilizou o allign_rasters (ou gdalwarp diretamente) do pacote gdalUtils com a opção r = "near" para aumentar a resolução da grade de 1kmx1km para 300mx300m, mesma projeção etc.
  3. Empilhe o mapa de cobertura terrestre de 300mx300m e a grade de 300mx300m da etapa 2, usando o pacote raster: stack_file <- stack (lc, grid).
  4. Crie um data.frame para combinar os mapas: df <- as.data.frame (rasterToPoints (stack_file)), que mapeia os números de células da grade do mapa de 1kmx1km para o mapa de cobertura de 300mx300m
  5. Use dplyr para calcular a parcela de células da classe de cobertura do solo nas células de 1kmx1km.
  6. Crie uma nova varredura com base na etapa 5, vinculando-a à grade original de 1kmx1km.

Este procedimento é executado muito rápido e sem problemas de memória no meu PC, quando tentei em um mapa de cobertura do solo com> 15 células de grade de moinho a 300mx300m.

Suponho que a abordagem acima também funcionará se alguém quiser combinar um arquivo de polígono com formas irregulares e dados rasterizados. Talvez você possa combinar as etapas 1 e 2 rasterizando diretamente o arquivo de polígono para uma grade de 300mx300 usando rasterize (raster provavelmente lento) ou gdal_rasterize. No meu caso, eu também precisava reprojetar, então usei o gdalwarp para reprojetar e desagregar ao mesmo tempo.

Michiel van Dijk
fonte
0

Eu tenho que enfrentar esse mesmo problema para extrair valores dentro do polígono de um grande mosaico (50k x 50k). Meu polígono tem apenas 4 pontos. O método mais rápido que encontrei é o cropmosaico no limite do polígono, o polígono triangulado em 2 triângulos e, em seguida, verifique se há pontos no triângulo (o algoritmo mais rápido que encontrei). Compare com a extractfunção, o tempo de execução é reduzido de 20 s para 0,5 s. No entanto, a função cropainda requer cerca de 2 s para cada polígono.

Desculpe, não posso fornecer o exemplo reproduzível completo. Os códigos R abaixo não incluem os arquivos de entrada.

Este método está funcionando apenas para polígonos simples.

par_dsm <- function(i, image_tif_name, field_plys2) {
    library(raster)
    image_tif <- raster(image_tif_name)
    coor <- field_plys2@polygons[[i]]@Polygons[[1]]@coords
    ext <- extent(c(min(coor[,1]), max(coor[,1]), min(coor[,2]), max(coor[,2])))

    extract2 <- function(u, v, us, vs) {
        u1 <- us[2]  - us[1]
        u2 <- us[3]  - us[2]
        u3 <- us[1]  - us[3]
        v1 <- vs[1]  - vs[2]
        v2 <- vs[2]  - vs[3]
        v3 <- vs[3]  - vs[1]
        uv1 <- vs[2] * us[1] - vs[1] * us[2]
        uv2 <- vs[3] * us[2] - vs[2] * us[3]
        uv3 <- vs[1] * us[3] - vs[3] * us[1]

        s1 <- v * u1 + u * v1 + uv1
        s2 <- v * u2 + u * v2 + uv2
        s3 <- v * u3 + u * v3 + uv3
        pos <- s1 * s2 > 0 & s2 * s3 > 0
        pos 
    }

    system.time({
        plot_rect <- crop(image_tif, ext, snap ='out')
        system.time({
        cell_idx <- cellFromXY(plot_rect, coor[seq(1,4),])
        row_idx <- rowFromCell(plot_rect, cell_idx)
        col_idx <- colFromCell(plot_rect, cell_idx)

        rect_idx <- expand.grid(lapply(rev(dim(plot_rect)[1:2]), function(x) seq(length.out = x)))

        pixel_idx1 <- extract2(
            rect_idx[,2], rect_idx[,1], 
            row_idx[c(1,2,3)], col_idx[c(1,2,3)])
        pixel_idx2 <- extract2(
            rect_idx[,2], rect_idx[,1], 
            row_idx[c(1,4,3)], col_idx[c(1,4,3)])
        pixel_idx <- pixel_idx1 | pixel_idx2
        })
    })
    mean(values(plot_rect)[pixel_idx])
}

# field_plys2: An object of polygons
# image_taf_name: file name of mosaic file
library(snowfall)
sfInit(cpus = 14, parallel = TRUE)
system.time(plot_dsm <- sfLapply(
    seq(along = field_plys2), par_dsm, image_tif_name, field_plys2))
sfStop()
Bangyou
fonte