重新采样栅格

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

1. I tried raster::resample()，但由于数据是分类的，我丢失了很多信息：

   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 
   

   正如您所看到的，新栅格具有所需的分辨率，但也有很多零。我认为这是正常的，因为我在中使用了“ngb”选项resample.

2. 第二个策略是使用raster::aggregate()但我发现很难定义一个因子整数，因为分辨率的变化并不简单（如分辨率的双倍或类似）。

   我的高分辨率栅格具有以下分辨率，我希望它将其聚合为0.008333333, 0.008333333 (x, y)分辨率达到相同程度。

   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)
   

   我按照描述尝试了 ~33.33 的系数aggregatehelp: "单元格数是 x 的单元格数除以fact*fact（当事实是单个数字时）。”尽管如此，生成的栅格数据似乎没有与我的其他低分辨率栅格相同的行数和列数。

我从未使用过这种高分辨率数据，而且我在计算上也受到限制（其中一些命令可以使用并行化）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? 

感谢您的帮助！

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 也已经显示 https://stackoverflow.com/a/37956798/2761575.

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 and y是相同的。但它们都可能存在问题，因为您的单元格可能有 5 个类，每个类约占 20%，因此选择一个获胜者是相当不合理的。