¿Trabajando con muchos datos y muchos rásteres en R?

Question

G'day, estoy trabajando con un gran conjunto de datos con ~ 125,000 ubicaciones lon/lat con fecha, para registros de presencia/ausencia de especies. En cada lugar, quiero saber cómo era el clima en cada lugar en la fecha y durante los 3 meses anteriores a la fecha. Para hacer esto, he descargado datos meteorológicos diarios para una variable climática determinada (por ejemplo, temperatura máxima) durante el período de 5 años que se tomaron los datos. Tengo un total de 1.826 archivos de trama, todos entre 2-3mb.

Había planeado apilar todos los archivos de trama, luego extraer un valor de cada trama (1,826) para cada punto. Esto produciría un archivo masivo que podría usar para buscar las fechas que necesito. Esto, sin embargo, no es posible porque no puedo acumular tantos rásteres. Traté de dividir los rásteres en pilas de 500, esto funciona, pero los archivos que produce son de aproximadamente 1 GB y muy lento (filas, 125 000, columnas, 500). Además, cuando trato de traer todos estos archivos a R para crear un gran marco de datos, no funciona.

Me gustaría saber si hay una manera de trabajar con esta cantidad de datos en R, o si hay un paquete que podría usar para ayudar. ¿Puedo usar un paquete como ff? ¿Alguien tiene alguna sugerencia para un método menos intensivo de energía para hacer lo que quiero hacer? He pensado en algo así como una función de aplicación, pero nunca he usado una antes y no estoy realmente seguro de por dónde empezar.

Cualquier ayuda sería realmente genial, gracias de antemano por su tiempo. El código que estoy usando actualmente sin éxito está debajo.

Saludos cordiales, Adam

library(raster) 
library(rgdal) 
library (maptools) 
library(shapefiles) 

# To create weather data files, first set the working directory to the appropriate location (i.e., maxt) 
# list of raster weather files 
files<- list.files(getwd(), pattern='asc') 
length(files) 

memory.size(4000) 
memory.limit(4000) 

# read in lon/lat data 
X<-read.table(file.choose(), header=TRUE, sep=',') 
SP<- SpatialPoints(cbind(X$lon, X$lat)) 

#separate stacks into mannageable sizes 
s1<- stack(files[1:500]) 
i1 <- extract(s1,SP, cellnumbers = True, layer = 1, nl = 500) 
write.table(i1, file="maxt_vals_all_points_all_dates_1.csv", sep=",", row.names= FALSE, col.names= TRUE) 
rm(s1,i1) 
s2<- stack(files[501:1000]) 
i2 <- extract(s2,SP, cellnumbers = True, layer = 1, nl = 500) 
write.table(i2, file="maxt_vals_all_points_all_dates_2.csv", sep=",", row.names= FALSE, col.names= TRUE) 
rm(s2,i2) 
s3<- stack(files[1001:1500]) 
i3 <- extract(s3,SP, cellnumbers = True, layer = 1, nl = 500) 
write.table(i3, file="maxt_vals_all_points_all_dates_3.csv", sep=",", row.names= FALSE, col.names= TRUE) 
rm(s3,i3) 
s4<- stack(files[1501:1826]) 
i4 <- extract(s4,SP, cellnumbers = True, layer = 1, nl =325) 
write.table(i4, file="maxt_vals_all_points_all_dates_4.csv", sep=",", row.names= FALSE, col.names= TRUE) 
rm(s4,i4) 

# read files back in to bind into final file !!! NOT WORKING FILES ARE TOO BIG!! 
i1<-read.table(file.choose(),header=TRUE,sep=',') 
i2<-read.table(file.choose(),header=TRUE,sep=',') 
i3<-read.table(file.choose(),header=TRUE,sep=',') 
i4<-read.table(file.choose(),header=TRUE,sep=',') 

vals<-data.frame(X, i1, i2, i3 ,i4) 
write.table(vals, file="maxt_master_lookup.csv", sep=",", row.names= FALSE, col.names= TRUE) 

Answer 1

me harían el extracto de un archivo de trama a la vez, y anexar los resultados para presentar a medida que avanza.

Hago trampa haciendo una lista de matrices, pero como el ráster puede tomar un nombre de archivo o una matriz (entre otras cosas) y puede indexar con "[[" en un vector de caracteres, debería funcionar igual en su caso .

files <- list(volcano, volcano * 2, volcano * 3) 
library(sp) 
SP <- SpatialPoints(structure(c(0.455921585146703, 0.237608166502031, 0.397704673508124, 0.678393354622703, 0.342820219769366, 0.554888036966903, 0.777351335399613, 0.654684656824567), .Dim = c(4L, 2L))) 

library(raster) 
for (i in seq_len(length(files))) { 

    r <- raster(files[[i]]) 
    e <- extract(r, SP) 
    ## print(e) ## print for debugging 
    write.table(data.frame(file = i, extract = e),"cellSummary.csv", col.names = i == 1, append = i > 1, sep = ",", row.names = FALSE) 
} 

Answer 2

Estoy utilizando el procesamiento en paralelo y una forma de recortar en función del número de celda. Esta función tomará cualquier punto espacial o polígono y devolverá los valores de una gran pila de ráster. Esta es una variante del código good example for large polygons.

Para mis datos esto toma alrededor de 350 segundos usando extracción, o 32 segundos en un servidor Linux de 16 núcleos. Espero que ayude a alguien!

# Define Functions 
    extract_value_point_polygon = function(point_or_polygon, raster_stack, num_workers){ 
      # Returns list containing values from locations of spatial points or polygons 
      lapply(c('raster','foreach','doParallel'), require, character.only = T) 
      registerDoParallel(num_workers) 
      ply_result = foreach(j = 1:length(point_or_polygon),.inorder=T) %do%{ 
       print(paste('Working on feature: ',j,' out of ',length(point_or_polygon))) 
       get_class= class(point_or_polygon)[1] 
       if(get_class=='SpatialPolygons'|get_class=='SpatialPolygonsDataFrame'){ 
        cell = as.numeric(na.omit(cellFromPolygon(raster_stack, point_or_polygon[j], weights=F)[[1]]))} 
       if(get_class=='SpatialPointsDataFrame'|get_class=='SpatialPoints'){ 
        cell = as.numeric(na.omit(cellFromXY(raster_stack, point_or_polygon[j,])))} 
       if(length(cell)==0)return(NA) 
       r = rasterFromCells(raster_stack, cell,values=F) 
       result = foreach(i = 1:dim(raster_stack)[3],.packages='raster',.inorder=T) %dopar% { 
        crop(raster_stack[[i]],r) 
       } 
       result=as.data.frame(getValues(stack(result))) 
       return(result) 
      } 
      endCluster() 
      return(ply_result) 
    }