重采样栅格。

Question

我正在尝试将高分辨率（25米）和分类数据（1至13）的森林覆盖栅格重新采样为具有较低分辨率（约1公里）的新RasterLayer。我的想法是将森林覆盖数据与其他低分辨率的栅格数据结合起来：

1. I tried raster::resample(), but since the data is categorical I lost a lot of information:

   summary(as.factor(df$loss_year_mosaic_30m))
     0       1   2   3  4   5   6   7  8   9   10  11   12  13
   3777691  65  101 50 151 145 159 295 291 134 102 126 104  91 
   

   As you can see, the new raster has the desired resolution but have lots of zeros as well. I suppose that is normal since I used the ´ngb´ option in resample.

2. The second strategy was using raster::aggregate() but I find difficult to define a factor integer since the change of resolution is not straightforward (like the double of the resolution or alike).

   My high-resolution raster has the following resolution, and I want it to aggregate it to a 0.008333333, 0.008333333 (x, y) resolution to the same extent.

   loss_year
   class       : RasterLayer 
   dimensions  : 70503, 59566, 4199581698  (nrow, ncol, ncell)
   resolution  : 0.00025, 0.00025  (x, y)
   extent      : -81.73875, -66.84725, -4.2285, 13.39725  (xmin, xmax, ymin, ymax)
   coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
   data source : /Volumes/LaCie/Deforestacion/Hansen/loss_year_mosaic_30m.tif 
   names       : loss_year_mosaic_30m 
   values      : 0, 13  (min, max)
   

   I have tried a factor of ~33.33 following the description of the aggregate help: "The number of cells is the number of cells of x divided by fact*fact (when fact is a single number)." Nonetheless, the resulting raster data do not seem to have the same number of rows and columns as my other low-resolution rasters.

我从未使用过这些高分辨率数据，而且我的计算能力也有限（可以使用clusterR并行化一些命令，但有时它们需要的时间与非并行化的命令相同，特别是对于最近邻计算不起作用的情况）。

我缺乏想法；也许我可以尝试layerize来获得计数栅格，但我必须进行“聚合”，然后出现了factor问题。由于这些过程需要花费我几天的时间，因此我想知道创建低分辨率栅格的最有效方法，同时又不会失去太多信息。

以下是一个可重复的示例：

r_hr <- raster(nrow=70, ncol=70) #High resolution raster with categorical data
set.seed(0)
r_hr[] <- round(runif(1:ncell(r_hr), 1, 5))
r_lr <- raster(nrow=6, ncol=6) #Low resolution raster

第一策略：信息丢失

r <- resample(r_hr, r_lr, method = "ngb") #The raster data is categorical

第二种策略：难以定义聚合因素。

r <- aggregate(r_hr, factor) #How to define a factor to get exactly the same number of cells of h_lr?

另一个选择：layerize

r_brick <- layerize(r_hr)
aggregate(r_brick, factor) #How to define factor to coincide with the r_lr dimensions? 

感谢您的帮助！

- topcat

---

首先要解决这类问题的方法是使用代码创建一个小的内存示例。然后，您可以轻松尝试不同的方法，并在此展示。 - Robert Hijmans

我编辑了问题并尝试创建可重现的示例。谢谢Robert。 - topcat

Answer 1

r_hr <- raster(nrow=70, ncol=70) #High resolution raster with categorical data
set.seed(0)
r_hr[] <- round(runif(1:ncell(r_hr), 1, 5))
r_lr <- raster(nrow=6, ncol=6)

r_hr
#class       : RasterLayer 
#dimensions  : 70, 70, 4900  (nrow, ncol, ncell)
#resolution  : 5.142857, 2.571429  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
#data source : in memory
#names       : layer 
#values      : 1, 5  (min, max)

r_lr
#class       : RasterLayer 
#dimensions  : 6, 6, 36  (nrow, ncol, ncell)
#resolution  : 60, 30  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 

无法进行直接聚合，因为70/6不是整数。

dim(r_hr)[1:2] / dim(r_lr)[1:2]
#[1] 11.66667 11.66667

最近邻重采样也不是一个好主意，因为结果会是任意的。

这里有一个你提出的按层方法，dww 已经展示过了。

b <- layerize(r_hr)
fact <- round(dim(r_hr)[1:2] / dim(r_lr)[1:2])
a <- aggregate(b, fact)
x <- resample(a, r_lr)

现在您有了比例。如果您想要一个单一的类，可以这样做：

y <- which.max(x)

在这种情况下，另一种方法是聚合类。

ag <- aggregate(r_hr, fact, modal) 
agx <- resample(ag, r_lr, method='ngb')

请注意，agx和y是相同的。但它们都可能存在问题，因为您可能会有包含每个约20％的5个类的单元格，这样选择一个获胜者就变得不太合理。

- Robert Hijmans

Answer 2

通常将土地覆盖图聚合成覆盖率层是相当标准的做法。也就是说，您应该致力于生成13个图层，每个图层都类似于网格单元中的覆盖率。这样做可以在保留尽可能多信息的同时降低分辨率。请注意，如果您需要不同的摘要统计信息而不是百分比符号“%”，只需通过更改“aggregate”中的“fun =”函数轻松适应所需的统计信息。

以下方法非常快速（在我的笔记本电脑上处理具有1亿个单元格的栅格仅需几秒钟）：

首先，让我们创建一些虚拟栅格以供使用。

Nhr <- 1e4 # resolution of high-res raster
Nlr <- 333 # resolution of low-res raster
r.hr <- raster(ncols=Nhr, nrows=Nhr)
r.lr <- raster(ncols=Nlr, nrows=Nlr)
r.hr[] <- sample(1:13, Nhr^2, replace=T)

现在，我们首先将高分辨率栅格聚合到与低分辨率栅格的分辨率几乎相同（到最近的整数单元格）。每个生成的图层包含单元格内原始栅格值为N的面积比例。

Nratio <- as.integer(Nhr/Nlr) # ratio of high to low resolutions, to nearest integer value for aggregation
layer1 <- aggregate(r.hr, Nratio, fun=function(x, na.rm=T) {mean(x==1, na.rm=na.rm)})
layer2 <- aggregate(r.hr, Nratio, fun=function(x, na.rm=T) {mean(x==2, na.rm=na.rm)})

最后，将低分辨率栅格重新采样到所需的分辨率。

layer1 <- resample(layer1, r.lr, method = "ngb") 
layer2 <- resample(layer2, r.lr, method = "ngb")

对于每一层，重复此步骤，并将您的图层构建为堆栈或多波段栅格

- dww