Berechnung der Geodatenstatistik für große Mengen von Rastern und Geo-Objekten in R

Question

Ich versuche, Mittelwerte für klimatische Variablen für räumliche Objekte in R zu berechnen. Die Herausforderung besteht darin, diese Mittel für jeden administrativen Bereich der Ebene 2 zu berechnen in der Welt (www.gadm.org), und ich brauche eine effiziente Methode zur Berechnung der Statistiken. Ich habe diese Statistiken ohne Probleme für kleinere Gebietsdefinitionen berechnet, die weniger Klimazonen/Kacheln umfassen, aber die logistischen Probleme sind zu einer Barriere geworden, wenn man versucht hat, diese Aufgabe auf die ganze Welt auszuweiten.

Ich habe versucht, mit gadm.org Level 2 Grenze Shapefile für die ganze Welt, dann importieren und fusionieren Worldclim.org die komplette Reihe von bioklimatischen Rastern (mit der höchsten verfügbaren räumlichen Auflösung) und Zonen/Fliesen, aber es scheint sei zu anspruchsvoll an Ressourcen. Insbesondere endet der Vorgang des Zusammenführens des vollständigen Satzes von Rasterzonen/Kacheln in ein globales Rasterobjekt niemals. Es schien der effizienteste Ansatz zu sein sowie die Wahrscheinlichkeit, Fehler zu minimieren, um die Rasterzonen für die ganze Welt zusammenzuführen.

Ich bin mir nicht sicher, wie ich das Problem von hier aus angehen sollte, da die Berechnung dieser Statistiken von Land zu Land äußerst mühsam und ineffizient ist. Darüber hinaus gibt es eine große Anzahl von Formen in der administrativen Grenzschicht, die einzelne Worldclim-Zonen/Kacheln überlappen, was zu Fehlern führen würde, wenn relevante Klimagegenstände bei den Berechnungen für die Formen fehlen, die nicht vollständig innerhalb einer einzigen Zone/Kachel liegen .

Ich frage mich, wie ich angesichts der Größe der Operation eine effiziente Lösung finden könnte.

Nach der Stufe 2 globale Verwaltungsgrenze des Herunterladen von Daten, habe ich den Code unten versucht:

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

Hier ist, wie ich würde herunterladen und Mosaik die Daten:

Zuerst würde ich eine Schleife verwenden Dadurch werden die Daten automatisch heruntergeladen und ein TIF-Raster für jeden Download exportiert.

Danach würde ich eine Dateiliste mit allen exportierten .tifs erstellen und ein virtuelles Mosaik mit der gdalbuildvrt() Funktion erstellen. Dadurch sparen Sie Zeit und Platz auf der Festplatte.

Schließlich können Sie die extract() Funktion verwenden, um Ihre Werte zu extrahieren. Beachten Sie, dass die Extrakt-Funktion VERY langsam ist und für immer größere Datensätze wie Ihre dauert.

Ich würde diese Operation in einer externen Software oder einer anderen Sprache wie Python, ArcGIS oder OTB ToolBox tun. Momentan arbeite ich viel mit der otbcli_LSMSVectorization Funktion aus der OTB Toolbox, die es ermöglicht, zonale Statistiken (Mean/Var) basierend auf einem Zonal Input Raster und einem Value Raster zu extrahieren.

Letztes Wort der Beratung: Versuchen Sie Mosaik und Ihre shp in kleineren Kacheln/AOI (so gut wie möglich) und führen Sie dann die Funktion extract() aufzuzuspalten, vielleicht mit einer foreach Schleife und %dopar% kombiniert. Dies sollte die Verarbeitungszeit enorm reduzieren. Weitere Informationen finden Sie in den folgenden Links.

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