¿Extracción realmente lenta de la trama incluso después de usar el recorte?

Question

Tengo un gran archivo raster (245295396) celdas y pilas de raster con 4 capas cada una que se encuentran en la extensión de este gran raster. Para empezar, estoy tratando de obtener el valor de una pila (3 canales) y para la misma zona de la gran trama. Todo funciona bien, pero la extracción de la trama grande tarda 5 minutos. Si repito este proceso 4000 veces más, tardaré 13 días.

cld<- raster("cdl_30m_r_il_2014_albers.tif") #this is the large raster
r<- stack(paste(path,"/data_robin/", fl,sep="")) #1 stack,I have 4000 similar
mat<-as.data.frame(getValues(r)) # getting values from the stack
xy<-xyFromCell(r,c(1:ncell(r)),spatial = TRUE)
clip1 <- crop(cld, extent(r)) # Tried to crop it to a smaller size
cells<-cellFromXY(clip1,xy)
mat$landuse<- NA
# mat$landuse<-cld[cells]
mat$landuse<- extract(clip1,cells) #this line takes 5 mins based on profiling

> cld
class       : RasterLayer 
dimensions  : 20862, 11758, 245295396  (nrow, ncol, ncell)
resolution  : 30, 30  (x, y)
extent      : 378585, 731325, 1569045, 2194905  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0 
data source : /Users/kaswani/R/Image/cdl_30m_r_il_2014_albers.tif 
names       : cdl_30m_r_il_2014_albers 
values      : 0, 255  (min, max) 

> r
class       : RasterStack 
dimensions  : 9230, 7502, 69243460, 4  (nrow, ncol, ncell, nlayers)
resolution  : 0.7995722, 0.7995722  (x, y)
extent      : 589084.4, 595082.8, 1564504, 1571884  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
names       : m_3608906_ne_16.1, m_3608906_ne_16.2, m_3608906_ne_16.3, m_3608906_ne_16.4 
min values  :                 0,                 0,                 0,                 0 
max values  :               255,               255,               255,               255 

Mis datos están en formato .tiff y soy nuevo en la codificación geoespacial.

También he probado el enfoque en ¿Aumentar la velocidad de recorte, máscara y extracción de raster por muchos polígonos en R? pero durante la parte de enmascaramiento da un error Error in compareRaster(x, mask) : different extent.

Answer 1

1) resample resultados en un 50% de mejora

Pude obtener alrededor de un 50% de mejora al remuestrear directamente desde el cld raster a un nuevo raster con la misma extensión/resolución que r y un método de muestreo del vecino más cercano:

system.time({
  mat<-as.data.frame(getValues(r))
  mat$landuse<- NA
  mat$landuse<-getValues(resample(cld,r,method='ngb'))
})
   user  system elapsed 
   2.60    0.00    2.61

contra.

system.time({
  mat<-as.data.frame(getValues(r)) # getting values from the stack
  xy<-xyFromCell(r,c(1:ncell(r)),spatial = TRUE)
  cells<-cellFromXY(clip1,xy)
  mat$landuse<- NA
  mat$landuse<- extract(clip1,cells) #this line takes 5 mins based on profiling
})
   user  system elapsed 
   4.98    0.00    5.02 

En un conjunto de datos más pequeño, y con una huella de memoria mucho menor

2) Paralelización pourrait conseguir mucho más

Eso mejorará las cosas significativamente, pero se puede obtener una mejora masiva si se puede paralelizar esto. R viene con un par de backends de paralelización y todos ellos se ejecutan a través de foreach . Supongo que va a procesar mat en su lugar o guardarlo para más tarde. Dado que se necesita mucho trabajo para obtener esos datos remuestreados, vamos a suponer que los guardaremos para más adelante. La forma más conveniente es probablemente un raster junto a la data_robin archivos.

Lamentablemente, las opciones de paralelización de Windows y Unix difieren. En linux, utilice doMC En Windows, utilice doSNOW . Suponiendo que empleamos a 4 trabajadores:

inicialización de linux:

library(doMC)
registerDoMC(4) # number of workers should be less than number of CPU cores

inicialización de las ventanas:

library(doSNOW)
cluster<-makeCluster(4, type = "SOCK") # num workers should be < num CPU cores
registerDoSNOW(cluster)

siguiente:

library(foreach)
library(tools)

# Assume you have an array of filenames called 'files'
foreach (i=1:length(filenames), .packages=c('raster')) %dopar% {
    r <- stack(paste0(path, "/data_robin/", files[i]))
    outFilename=paste0(path, "/data_robin/", file_path_sans_ext(files[i]), "_cld.tif")
    cldResampled <- resample(cld,r,method='ngb')
    writeRaster(cldResampled, filename=outFilename, format="GTiff")
}

Uno de los inconvenientes de la foreach es que es difícil saber cuando algo va mal. Sería bueno hacer esto en serie primero reemplazando el %dopar% con %do% hasta que sepas que está funcionando, y luego deja que se ejecute por completo.

Advertencias En mi sencillo ejemplo anterior (cada trama tenía 1/100 de los píxeles en cld y r respectivamente) sólo mejoré un 30% adicional al involucrar a 5 trabajadores respecto a hacerlo en serie con un solo proceso. No pude paralelizar el ejemplo con rasters grandes sin obtener Error in { : task 1 failed - "cannot allocate vector of size xx.x Mb" . Creo que conceptualmente esto debería funcionar, pero no fui capaz de hacerlo funcionar a la escala en la que tú estás trabajando.

Answer 2

@aaryno Gracias por la detallada explicación. Intenté usar tu código pero produce una salida diferente. Tuve que hacer un pequeño cambio sustituyendo en la siguiente línea "cld" mat$landuse<-getValues(resample(cld,r,method='ngb')) por "clip1".Después de eso, funciona bien. He probado a desplegar el código de abajo en una máquina windows y funciona.

path<-"D:/TC/il/il"
out_path<-"D:/TC/out/"
setwd("D:/TC/test/")
fl1 <- list.files(path,pattern="*6.tif|*5.tif")
cld<- raster("cdl_30m_r_il_2014_albers.tif")

worker<- function(i,path,out_path,fl1,cld){ 
  r<- stack(paste0(path, fl1[i]))
  outFilename=paste0(out_path,"cld_",strsplit(fl1[i], "\\.")[[1]][1],".csv")
  mat<- data.frame()
  mat<-as.data.frame(getValues(r))
  clip1 <- crop(cld, extent(r))
  mat$landuse<- NA
  mat$landuse<-getValues(resample(clip1,r,method='ngb'))
  mat$cell<-1:length(mat$landuse)
  rowsToKeep<-which(mat$landuse > 120 & mat$landuse < 125)
  mat<- mat[rowsToKeep,]
  write.csv(mat, file=outFilename, row.names=F)
}

i<- 1:length(fl1)
cl<-makeCluster(8,type="SOCK")
clusterEvalQ(cl, { library(raster); library(stringr);library(rgdal) })
system.time(clusterApply(cl,i,worker,path,out_path,fl1,cld))
stopCluster(cl)

Realmente no tengo que guardarlo en formato csv ni crear un dataframe pero es lo que más cómodo me resulta dentro de R. Agradecería si algún cambio en este código puede suponer una gran mejora.