PAYATHON 2026

如何在 R 中高效地将 25m 分辨率的 raster 聚合为 100m？

技术老齐

可以直接用 terra::aggregate()，通过自定义函数在每个 4 × 4 像元块内统计所有类别的比例。对于 21 类数据，结果是一个包含 21 个图层的 SpatRaster，每层对应一个类别。整个过程无需创建 polygon，也不用先生成 21 个 25m 二值栅格。

这样能减少中间栅格，但运行时间未必更短。自定义 R 函数有调用开销，原方案使用的内置 mean 则经过专门优化。处理大数据时，要根据内存占用和实际耗时选择。

为什么按 4 × 4 聚合就够了

在同一投影坐标系下，目标网格与原网格对齐时：

100 / 25 = 4

每个完整的 100m 单元格包含 16 个 25m 像元。这些像元的面积相同，因此可以直接按数量计算类别比例：

[ p_k=\frac{n_k}{N} ]

其中，(n_k) 是类别 (k) 的像元数量，(N) 是参与统计的像元总数。没有缺失值时，(N=16)。例如，某类别占 6 个像元，其比例就是 6 / 16 = 0.375。

一次聚合，返回所有类别的比例

示例使用两个类别，与问题中的数据对应：

library(terra)

set.seed(42)

# 模拟以米为单位的投影栅格
Land_cover <- rast(
  nrows = 100,
  ncols = 100,
  xmin = 0,
  xmax = 2500,
  ymin = 0,
  ymax = 2500,
  crs = "EPSG:27700"
)

classes <- c(1L, 2L)
values(Land_cover) <- sample(
  classes,
  ncell(Land_cover),
  replace = TRUE
)

class_fraction <- function(v) {
  v <- v[!is.na(v)]

  # 全部缺失时，仍返回固定数量的结果
  if (length(v) == 0L) {
    return(rep(NA_real_, length(classes)))
  }

  # 将类别编码映射到固定的输出图层顺序
  idx <- match(v, classes)

  if (anyNA(idx)) {
    stop("发现未包含在 classes 中的类别编码")
  }

  tabulate(idx, nbins = length(classes)) / length(v)
}

cover100 <- aggregate(
  Land_cover,
  fact = 4,
  fun = class_fraction
)

names(cover100) <- paste0("class_", classes)

cover100

cover100 是一个 25 × 25、分辨率为 100m 的两层栅格。class_1 和 class_2 分别保存两个类别的比例，取值范围为 [0, 1]。

处理实际的 21 类数据时，读取原始文件，并按数据字典设置 classes：

Land_cover <- rast("land_cover_25m.tif")

# 仅当真实类别编码确实是 1 到 21 时使用
classes <- 1:21

cover100 <- aggregate(
  Land_cover,
  fact = 4,
  fun = class_fraction,
  filename = "land_cover_fraction_100m.tif",
  wopt = list(datatype = "FLT4S")
)

names(cover100) <- paste0("class_", classes)

class_fraction() 使用当前的 classes。如果实际编码是 10、20、30 等不连续值，就应填入完整的真实编码列表。match() 会将这些编码转换为 tabulate() 所需的位置编号。

某个聚合单元格中没有出现的类别，其比例会返回 0，所以输出图层的数量和顺序保持一致。

缺失值决定了比例的含义

代码先去掉 NA，计算的是各类别占有效像元的比例，与原方案加上 na.rm = TRUE 后的含义一致。

例如，一个目标单元格中有 8 个有效像元，全部属于类别 1，另外 8 个为 NA，类别 1 的比例就是：

8 / 8 = 1

如果要计算各类别占完整 100m 单元格的比例，分母应固定为 16。将函数最后一行改为：

tabulate(idx, nbins = length(classes)) / 16

同一个例子的结果就变为 0.5，所有类别的比例之和也可能小于 1。全为 NA 的单元格仍会在前面的判断中返回 NA。

还要检查土地覆盖数据是否用 0、255 等编码表示无数据。这些数值不会自动被当作 NA，需要按数据说明处理。

大栅格如何选择更快的方案

自定义函数省去了逐类生成高分辨率二值图层的步骤，但每个目标单元格都要调用一次 R 函数。处理英国全域这样规模的数据时，调用开销可能比较明显。

原先逐类使用 mean 的方法也合理。如果实测更快，可以逐类聚合并写入文件，只保留低分辨率结果：

files <- file.path(
  tempdir(),
  paste0("class_", classes, "_100m.tif")
)

for (i in seq_along(classes)) {
  aggregate(
    Land_cover == classes[i],
    fact = 4,
    fun = "mean",
    na.rm = TRUE,
    filename = files[i],
    overwrite = TRUE,
    wopt = list(datatype = "FLT4S")
  )
}

cover100 <- rast(files)
names(cover100) <- paste0("class_", classes)

这种方法仍需逐类计算，但不用同时保留全部 25m 二值图层。代码使用临时目录，需要长期保存时，应换成正式输出目录。比较耗时时，要使用同一片有代表性的数据，并保持缺失值处理规则和文件输出条件一致。

fact = 4 生成的网格继承原栅格的聚合起点。如果要匹配已有的 100m 模板，还需检查其 CRS、范围和网格起点。源数据的行列数若不能被 4 整除，边缘就会出现不完整的聚合块，解释结果时需要结合前面的分母规则。

目标网格与源网格错位，或统计区域有不规则边界时，才需要考虑覆盖面积加权等方法。题目中的网格规则且对齐，无需先转成 polygon。

备注：内容仅供参考。