Работа с большим количеством данных и множеством растров в R?

Question

G'day, я работаю с большим набором данных с местоположением ~ 125 000 lon/lat с датой, для записей присутствия/отсутствия видов. В каждом месте я хочу выяснить, какая погода была в каждом месте на дату и в течение 3-х минут до даты. Для этого я загрузил ежедневные метеорологические данные для заданной погодной переменной (например, максимальной температуры) в течение 5-го периода, когда данные были сделаны. Я имею в общей сложности 1826 растровых файлов, все между 2-3 МБ.

Я планировал складывать все растровые файлы, а затем извлекать значение из каждого растра (1,826) для каждой точки. Это создаст массивный файл, который я могу использовать для поиска нужных мне дат. Это, однако, невозможно, потому что я не могу собрать столько растров. Я попытался разделить растры на стопки 500, это работает, но файлы, которые он производит, составляют около 1 ГБ и очень медленные (строки, 125 000, столбцы, 500). Кроме того, когда я пытаюсь привести все эти файлы в R, чтобы создать большой фрейм данных, он не работает.

Я хотел бы знать, есть ли способ работать с этим объемом данных в R, или если есть пакет, который я мог бы использовать, чтобы помочь. Могу ли я использовать пакет, например ff? Есть ли у кого-нибудь предложения по менее интенсивному методу делать то, что я хочу делать? Я подумал о чем-то вроде функции лапши, но никогда не использовал его раньше, и я не уверен, с чего начать.

Любая помощь будет действительно замечательной, заблаговременно за ваше время. Код, который я сейчас использую без успеха, приведен ниже.

С наилучшими пожеланиями, Адам

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

Я хотел бы сделать экстракт один растровый файл, в то время, и добавить результаты в файл, как вы идете.

Я обманываю, создавая список матриц, но так как растр может принимать имя файла или матрицу (между прочим), и вы можете индексировать «[[» на символьном векторе, он должен работать примерно так же в вашем случае ,

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

Я использую параллельную обработку и форму обрезки на основе номера ячейки. Эта функция будет принимать любые пространственные точки или многоугольник и возвращать значения из большого растрового стека. Это вариант по коду good example for large polygons.

Для моих данных это занимает около 350 секунд, используя извлечение, или 32 секунды на 16-ядерном сервере linux. Надеюсь, это поможет кому-то!

# 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) 
    }