首页 文章

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

提问于
浏览 次
2

我正在尝试计算R中空间对象的气候变量的平均值 . 挑战在于我正在尝试为世界上每个二级行政区域（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)
r geospatial shapefile r-raster

1 回答

  • 0

    以下是我下载和镶嵌数据的方法：

    首先，我将使用一个自动下载数据的循环，并为每次下载导出.tif栅格 .

    之后，我将构建一个包含所有导出的.tifs的文件列表，并使用 gdalbuildvrt() 函数创建虚拟马赛克 . 这将为您节省一些时间和硬盘空间 .

    最后，您可以使用 extract() 函数来提取值 . 请注意，提取函数的速度很慢，并且对于像您这样的较大数据集而言需要永久 .

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

    Last word of advice: 尝试将镶嵌和shp拆分为较小的tile / 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)
    回复于

相关问题