大きな集合のラスタと空間オブジェクトの地理空間統計の計算R

Question

Rの空間オブジェクトの気象変数の平均値を計算しようとしています。レベル2の管理エリアごとにこれらの平均を計算しようとしています（www.gadm.org）、統計を効率的に計算する方法が必要です。気候ゾーン/タイルの数が少ない、より小さなエリア定義については問題なしでこれらの統計値を計算しましたが、このタスクを世界全体に拡大しようとすると、物流上の問題が障壁になりました。

gadm.orgのレベル2境界シェイプファイルを全世界に使用し、worldclim.orgの完全な生化学ラスタセット（利用可能な空間解像度の最高値）とゾーン/タイルをインポートしてマージしましたが、リソースをあまりにも多く必要とする。具体的には、ラスタゾーン/タイルの完全なセットを1つのグローバルラスタオブジェクトにマージする動作は、決して終了しない。それは世界全体のラスターゾーンを統合するために、エラーを最小限に抑える最も効率的なアプローチであると思われました。

これらの統計情報を国ごとに計算すると、非常に退屈で非効率的なように、私はここから問題に近づく方法を確信していません。さらに、個々のワールドクリムゾーン/タイルに重なる行政境界層には多数の形状が存在し、単一のゾーン/タイル内に完全には存在しない形状の計算から関連する気候オブジェクトが欠落している場合にはエラーが発生する。

私は、操作の規模がどのように効率的なソリューションを考え出すことができるのだろうかと思います。

レベル2のグローバル管理境界データをダウンロードした後、私は以下のコードを試してみました：ここ

library(raster) 
library(rgdal) 
library(maptools) 
library(foreign) 

#SET WORKING DIRECTORY 
setwd("C:/gadm28") 

#IMPORT GLOBAL ADMINISTRATIVE BOUNDARIES (LEVEL 2) DATA FROM HARD DRIVE 
gadm <- readOGR(dsn="C:/gadm28", layer="gadm28") 

#IMPORT GLOBAL (ALL TILES) BIOCLIMACTIC DATA DIRECTLY FROM WORLDCLIM.ORG 
climatezone00 <- getData('worldclim', var='bio', res=0.5, lon=-180, lat=90) 
climatezone01 <- getData('worldclim', var='bio', res=0.5, lon=-150, lat=90) 
climatezone02 <- getData('worldclim', var='bio', res=0.5, lon=-120, lat=90) 
climatezone03 <- getData('worldclim', var='bio', res=0.5, lon=-90, lat=90) 
climatezone04 <- getData('worldclim', var='bio', res=0.5, lon=-60, lat=90) 
climatezone05 <- getData('worldclim', var='bio', res=0.5, lon=-30, lat=90) 
climatezone06 <- getData('worldclim', var='bio', res=0.5, lon=0, lat=90) 
climatezone07 <- getData('worldclim', var='bio', res=0.5, lon=30, lat=90) 
climatezone08 <- getData('worldclim', var='bio', res=0.5, lon=60, lat=90) 
climatezone09 <- getData('worldclim', var='bio', res=0.5, lon=90, lat=90) 
climatezone010 <- getData('worldclim', var='bio', res=0.5, lon=120, lat=90) 
climatezone011 <- getData('worldclim', var='bio', res=0.5, lon=150, lat=90) 

climatezone10 <- getData('worldclim', var='bio', res=0.5, lon=-180, lat=60) 
climatezone11 <- getData('worldclim', var='bio', res=0.5, lon=-150, lat=60) 
climatezone12 <- getData('worldclim', var='bio', res=0.5, lon=-120, lat=60) 
climatezone13 <- getData('worldclim', var='bio', res=0.5, lon=-90, lat=60) 
climatezone14 <- getData('worldclim', var='bio', res=0.5, lon=-60, lat=60) 
climatezone15 <- getData('worldclim', var='bio', res=0.5, lon=-30, lat=60) 
climatezone16 <- getData('worldclim', var='bio', res=0.5, lon=0, lat=60) 
climatezone17 <- getData('worldclim', var='bio', res=0.5, lon=30, lat=60) 
climatezone18 <- getData('worldclim', var='bio', res=0.5, lon=60, lat=60) 
climatezone19 <- getData('worldclim', var='bio', res=0.5, lon=90, lat=60) 
climatezone110 <- getData('worldclim', var='bio', res=0.5, lon=120, lat=60) 
climatezone111 <- getData('worldclim', var='bio', res=0.5, lon=150, lat=60) 

climatezone20 <- getData('worldclim', var='bio', res=0.5, lon=-180, lat=30) 
climatezone21 <- getData('worldclim', var='bio', res=0.5, lon=-150, lat=30) 
climatezone22 <- getData('worldclim', var='bio', res=0.5, lon=-120, lat=30) 
climatezone23 <- getData('worldclim', var='bio', res=0.5, lon=-90, lat=30) 
climatezone24 <- getData('worldclim', var='bio', res=0.5, lon=-60, lat=30) 
climatezone25 <- getData('worldclim', var='bio', res=0.5, lon=-30, lat=30) 
climatezone26 <- getData('worldclim', var='bio', res=0.5, lon=0, lat=30) 
climatezone27 <- getData('worldclim', var='bio', res=0.5, lon=30, lat=30) 
climatezone28 <- getData('worldclim', var='bio', res=0.5, lon=60, lat=30) 
climatezone29 <- getData('worldclim', var='bio', res=0.5, lon=90, lat=30) 
climatezone210 <- getData('worldclim', var='bio', res=0.5, lon=120, lat=30) 
climatezone211 <- getData('worldclim', var='bio', res=0.5, lon=150, lat=30) 

climatezone30 <- getData('worldclim', var='bio', res=0.5, lon=-180, lat=0) 
climatezone31 <- getData('worldclim', var='bio', res=0.5, lon=-150, lat=0) 
climatezone32 <- getData('worldclim', var='bio', res=0.5, lon=-120, lat=0) 
climatezone33 <- getData('worldclim', var='bio', res=0.5, lon=-90, lat=0) 
climatezone34 <- getData('worldclim', var='bio', res=0.5, lon=-60, lat=0) 
climatezone35 <- getData('worldclim', var='bio', res=0.5, lon=-30, lat=0) 
climatezone36 <- getData('worldclim', var='bio', res=0.5, lon=0, lat=0) 
climatezone37 <- getData('worldclim', var='bio', res=0.5, lon=30, lat=0) 
climatezone38 <- getData('worldclim', var='bio', res=0.5, lon=60, lat=0) 
climatezone39 <- getData('worldclim', var='bio', res=0.5, lon=90, lat=0) 
climatezone310 <- getData('worldclim', var='bio', res=0.5, lon=120, lat=0) 
climatezone311 <- getData('worldclim', var='bio', res=0.5, lon=150, lat=0) 

climatezone40 <- getData('worldclim', var='bio', res=0.5, lon=-180, lat=-30) 
climatezone41 <- getData('worldclim', var='bio', res=0.5, lon=-150, lat=-30) 
climatezone42 <- getData('worldclim', var='bio', res=0.5, lon=-120, lat=-30) 
climatezone43 <- getData('worldclim', var='bio', res=0.5, lon=-90, lat=-30) 
climatezone44 <- getData('worldclim', var='bio', res=0.5, lon=-60, lat=-30) 
climatezone45 <- getData('worldclim', var='bio', res=0.5, lon=-30, lat=-30) 
climatezone46 <- getData('worldclim', var='bio', res=0.5, lon=0, lat=-30) 
climatezone47 <- getData('worldclim', var='bio', res=0.5, lon=30, lat=-30) 
climatezone48 <- getData('worldclim', var='bio', res=0.5, lon=60, lat=-30) 
climatezone49 <- getData('worldclim', var='bio', res=0.5, lon=90, lat=-30) 
climatezone410 <- getData('worldclim', var='bio', res=0.5, lon=120, lat=-30) 
climatezone411 <- getData('worldclim', var='bio', res=0.5, lon=150, lat=-30) 

#COMBINE ZONES TO CREATE ONE COMPLETE CLIMATE OBJECT 
climatemosaic <- mosaic(climatezone01, climatezone02, climatezone03, climatezone04, climatezone05, climatezone06, climatezone07, climatezone08, climatezone09, climatezone010, climatezone011, climatezone10, climatezone11, climatezone12, climatezone13, climatezone14, climatezone15, climatezone16, climatezone17, climatezone18, climatezone19, climatezone110, climatezone111, climatezone20, climatezone21, climatezone22, climatezone23, climatezone24, climatezone25, climatezone26, climatezone27, climatezone28, climatezone29, climatezone210, climatezone211, climatezone30, climatezone31, climatezone32, climatezone33, climatezone34, climatezone35, climatezone36, climatezone37, climatezone38, climatezone39, climatezone310, climatezone311, climatezone40, climatezone41, climatezone42, climatezone43, climatezone44, climatezone45, climatezone46, climatezone47, climatezone48, climatezone49, climatezone410, climatezone411, fun=mean) 

#EXTRACT MEAN VALUES FOR BOUNDARY POLYGONS & ATTACH TO SPDF (WEIGHT AND BUFFER OPTIONS NOT USED HERE) 
gadmMEANS <- extract(climatemosaic, gadm, fun=mean, na.rm=TRUE, small=TRUE, layer=1, nl=19, sp=TRUE) 

Answer 1

は、私がデータをダウンロードし、モザイクする方法を次のとおりです。

まず、私は、ループを使用しますデータは自動的にダウンロードされ、ダウンロードごとに.tifラスタをエクスポートします。

その後、エクスポートされたすべての.tifファイルリストを作成し、gdalbuildvrt()機能を使用して仮想モザイクを作成します。これは、ある程度の時間とハードディスクスペースを安全にします。

最後に、extract()関数を使用して値を抽出することができます。注目すべきは、抽出機能がVERYであり、あなたのような大きなデータセットでは時間がかかります。

私はこのソフトウェアをPython、ArcGIS、OTB ToolBoxのような外部のソフトウェアまたは他の言語で個人的に行います。現在は、OTB Toolboxのotbcli_LSMSVectorization関数を使用して、Zonal Input RasterとValue Rasterに基づいてゾーン統計値（Mean/Var）を抽出することができます。

アドバイスの最後の言葉：多分foreachループと%dopar%と組み合わせて、モザイクあなたを分割するために試してみて、小さなタイル/ののAOI（できるだけ良好）でのごSHP、その後extract()機能を実行します。これは処理時間を大幅に短縮するはずです。詳細は以下のリンクをご覧ください。

library(raster) 
library(gdalUtils) 

lon_vec <- rep(seq(-180,150,30),5) 
lat_vec <- sort(rep(seq(90,-30,-30),12), decreasing=T) 

#Download Worldclim Data and export as Tif 
for(i in 1:length(lon_vec)) { 
    ras <- getData('worldclim', var='bio', res=0.5, lon=lon_vec[i], lat=lat_vec[i]) 
    writeRaster(ras, filename=paste0("YourSubfolder/worldclim_lon_",lon_vec[i],"lat_",lat_vec[i],".tif")) 
} 

#Create list with all exported .tif iles 
ras_lst <- list.files("YourSubfolder/",full.names=T, pattern=".tif$") 

#Build virtual raster mosaic 
gdalbuildvrt(ras_lst, "YourSubfolder/worldclimMosaic.vrt") 

#Load virtual mosaic into R 
climatemosaic <- stack("YourSubfolder/worldclimMosaic.vrt") 

#Extract Mean Values 
gadmMEANS <- extract(climatemosaic, gadm, fun=mean, na.rm=TRUE, small=TRUE, layer=1, nl=19, sp=TRUE) 

• Gdal BuildVRT
• otbcli_LSMSVectorization
• foreach %dopar% with extract