5

今日は、種の存在/不在記録のために、日付付きの 125,000 の経度/緯度の場所を持つ大規模なデータセットを使用しています。各場所について、日付と日付の 3 か月前の各場所の天気を調べたいと思います。これを行うために、データが取得された 5 年間の特定の気象変数 (最高気温など) の毎日の気象データをダウンロードしました。合計 1,826 個のラスター ファイルがあり、すべて 2 ～ 3 MB です。

すべてのラスター ファイルをスタックし、各ポイントのすべてのラスター (1,826) から値を抽出することを計画していました。これにより、必要な日付を検索するために使用できる巨大なファイルが生成されます。ただし、それほど多くのラスターをスタックできないため、これは不可能です。ラスターを 500 のスタックに分割しようとしましたが、これは機能しますが、生成されるファイルは約 1Gb で、非常に低速です (行数 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)
4

2 に答える 2

5

一度に 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)
}
于 2012-04-07T07:01:56.173 に答える
0

並列処理と、セル数に基づくトリミングの形式を使用しています。この関数は、任意の空間ポイントまたはポリゴンを取得し、大きなラスター スタックから値を返します。これは、大きなポリゴンのコードの良い例です。

私のデータでは、extract を使用すると約 350 秒かかり、16 コアの Linux サーバーでは 32 秒かかります。それが誰かを助けることを願っています!

 # 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)
  }
于 2016-04-11T23:44:12.197 に答える