重新采样栅格

Question

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

1. 我尝试了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 
   

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

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)
   

   我按照aggregate帮助的描述尝试了~33.33的因子:“单元格数量是x的单元格数量除以fact*fact (当事实是一个数字时)。”尽管如此，生成的栅格数据似乎没有与我的其他低分辨率栅格相同的行数和列数。

我从未使用过这种高分辨率数据，而且我的计算能力也有限(其中一些命令可以使用clusterR进行并行化，但有时它们与非并行化命令花费的时间相同，特别是因为它们不适用于最近邻计算)。

我缺乏想法；也许我可以尝试分层来获取计数栅格，但我必须“聚合”，并且出现了因子问题。由于这个过程需要我几天的时间来处理，我确实想知道创建较低分辨率栅格而不丢失太多信息的最有效方法

可重现的示例如下:

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?

另一个选项:分层

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

感谢您的帮助!

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 also showed already 。

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)

现在你有了比例。如果你想要一个单独的类(class)，你可以这样做

y <- which.max(x)

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

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

请注意，agx 和 y 是相同的。但它们都可能存在问题，因为您的单元格可能包含 5 个类，每个类约占 20%，因此选择一个获胜者是相当不合理的。