在 R 中计算大量栅格和空间对象的地理空间统计数据

Question

我正在尝试计算 R 中空间对象的气候变量的平均值。挑战在于我正在尝试为世界上每个 2 级行政区域 (www.gadm.org) 计算这些平均值，而我需要一种有效的方法来计算统计数据。我计算了这些统计数据，对于跨越较少气候区/区块的较小区域定义没有任何问题，但在尝试将这项任务扩展到整个世界时，后勤问题已成为障碍。

我尝试使用 gadm.org 的 2 级边界 shapefile 为整个世界，然后导入和合并 worldclim.org 的全套生物气候栅格（以最高可用空间分辨率）和区域/图块，但它似乎是对资源要求太高。具体来说，将整套栅格区域/切片合并到一个全局栅格对象中的操作永远不会完成。合并整个世界的栅格区域似乎是最有效的方法，也是最有可能最小化错误的方法。

我不确定如何从这里解决问题，因为逐个国家/地区计算这些统计数据似乎非常繁琐且效率低下。此外，行政边界层中有大量形状与各个 Worldclim 区域/图块重叠，如果在计算不完全位于单个区域/图块内的形状时缺少相关气候对象，则会导致错误.

考虑到手术的规模，我想知道如何提出一个有效的解决方案。

下载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() 函数来提取您的值。请注意，提取功能非常很慢，并且对于像您这样的大型数据集需要很长时间。

我会亲自使用外部软件或其他语言（如 Python、ArcGIS 或 OTB ToolBox）执行此操作。目前，我正在使用 OTB 工具箱中的 otbcli_LSMSVectorization 函数进行大量工作，该函数使您能够基于区域输入栅格和值栅格提取区域统计数据（均值/方差）。

最后的忠告： 尝试将马赛克和 shp 拆分为更小的图块/AOI（尽可能好），然后运行 ​​extract() 函数，可能与 foreach 结合使用循环和%dopar%。这应该会大大减少处理时间。如需更多信息，请查看以下链接。

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