0

我需要计算原始分辨率为 1 x 0.00811 但聚合为 2 度并具有新范围的栅格中 NA 和非 NA 值的数量。

原始栅格（可在此处获得：https ://datadryad.org/stash/dataset/doi:10.5061/dryad.052q5 ，请参见输出 1）与另一个数据集（不幸的是，不是开源的）合并以在输出中生成栅格2：

输出 1

class      : RasterLayer 
dimensions : 19142, 35738, 684096796  (nrow, ncol, ncell)
resolution : 0.01, 0.00811  (x, y)
extent     : -178.6931, 178.6869, -65.29534, 89.94628  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : HumanFootprintWGS84.tif 
names      : HumanFootprintWGS84 
values     : 0, 50  (min, max)

输出 2

dimensions : 19142, 35738, 684096796  (nrow, ncol, ncell)
resolution : 0.01, 0.00811  (x, y)
extent     : -178.6931, 178.6869, -65.29534, 89.94628  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : r_tmp_2022-02-10_145403_42352_01781.grd 
names      : layer 
values     : 0.6730382, 1  (min, max)

我用于重新采样的栅格是一个虚拟栅格，如下所示：

输出 3

class      : RasterLayer 
dimensions : 65, 180, 11700  (nrow, ncol, ncell)
resolution : 2, 2  (x, y)
extent     : -180, 180, -65, 65  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs

不幸的是，合并并因此重新采样的栅格在海岸线和湖泊附近有很多 NA 值；和中的平均方法st_warp将 NA 值转换为 0，这会扭曲某些单元格的平均值的计算。

我决定通过将平均值乘以非 NA 值的数量/（非 NA 的数量 - NA 的数量）来调整平均值产生的值，即：

平均 x {非NA值的数量/（非NA值的数量/NA值的数量）}

为此，我需要知道非 NA 和 NA 值的数量是原始合并栅格（输出 2），但分辨率为 2 度，输出 3 中的范围（-180、180、-65、65）。

我是地图和栅格的新手，如果这是一个基本问题，我深表歉意。

我试图用原始数据集和虚拟网格的坐标进行栅格化，但这不会让我得到 NA 值的数量，只是单元格的数量（加上合并的栅格超过 5.1 Gb）。我试图修剪 NA 值（愚蠢的想法），st_warp但没有平均（无论如何都使用最近的邻居）。

如果有人有任何想法或有更优雅的解决方案，我将不胜感激。

由衷的感谢，

编辑：通过反复试验，我发现无论值是否加载到内存中，重新采样（再次感谢您，@Robert Hijmans）和 st_warp 都会有所不同。

以下是示例输出，加载和不加载值：

terra::resample 进行平均而不在内存中加载数据

terra::resample对内存中加载的数据进行平均

4

1 回答 1

0

示例数据

library(terra)
#terra 1.5.20
library(geodata)
w <- world(path=".")
# input raster
x <- rast(res=3)
x <- rasterize(w, x, field=1)

# output raster
r <- rast(res=10)

解决方案：

rs <- resample(x, r, "sum")

rs
#class       : SpatRaster 
#dimensions  : 18, 36, 1  (nrow, ncol, nlyr)
#resolution  : 10, 10  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 
#source      : memory 
#name        :     layer 
#min value   : 0.1111111 
#max value   :  11.11111 

plot(rs)

在此处输入图像描述

你需要terra 1.5.20这个。那是目前的开发版本。在 Windows 或 OSX 上安装它的最简单方法是使用install.packages('terra', repos='https://rspatial.r-universe.dev')

早期版本terra忽略“sum”选项。

回应您的评论：当我使用文件作为数据源时，sum我得到了相同的结果：average

xx <- writeRaster(x, "test.tif", overwrite=T)
rs <- resample(xx, r, "sum")

我们可以将结果与多边形的精确提取进行比较（这也有效，但会在大数据集上阻塞 R）

p <- as.polygons(r)
e <- extract(x, p, exact=TRUE, fun=sum, na.rm=TRUE)
re <- rast(r)    
re[e[,1]] <- e[,2]

plot(rs, re, xlab="resample", ylab="extract")

在此处输入图像描述

于 2022-02-10T18:24:53.527 回答