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
Respostas:
Finalmente consegui melhorar esta função. Descobri que, para meus propósitos, era mais rápido
rasterize()
primeiro no polígono e emgetValues()
vez de usá-loextract()
. 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 acheigetValues()
muito mais rápido que aextract()
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.
Aqui está a função:
Portanto, para usá-lo, ajuste-o
single@data$OWNER
para caber no nome da coluna do seu polígono de identificação (acho que poderia ter sido incorporado à função ...) e insira:fonte
getValues
foi muito mais rápida que aextract
não parece válida, porque se você usarextract
, não precisará fazercrop
erasterize
(oumask
). O código na pergunta original faz as duas coisas, e isso deve dobrar o tempo de processamento.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 ------------------------
fonte
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
raster
pacote contém azonal()
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.fonte
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:
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.
fonte
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
crop
mosaico 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 aextract
função, o tempo de execução é reduzido de 20 s para 0,5 s. No entanto, a funçãocrop
ainda 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.
fonte