使用gdalUtils::gdalwarp将30m分辨率sf栅格聚合至150m时单元格缺失

• 重采样方法选型错误：通过fasterize生成的30m栅格属于分类/掩膜类栅格（一般为0-1二值值，标识多边形覆盖范围），选用的r="bilinear"双线性内插适用于连续值栅格，用在分类栅格上会产生数值模糊、边界跳变，甚至出现聚合后有效像素缺失的问题。
• 栅格网格未对齐：仅传入目标分辨率tr=c(150,150)时，gdalwarp会自动生成新的像元网格，若新网格原点和30m源栅格的像元角点未形成整数倍偏移，就会出现聚合像元跨边界错位，视觉上表现为单元格跳变缺失。
• 冗余重投影参数：若30m源栅格本身已经是UTM 11N NAD83坐标系，手动传入t_srs参数会触发无意义的重投影流程，过程中的像素对齐误差会进一步放大网格错位问题。
• 未固定输出有效范围：默认状态下gdalwarp会自动按输入范围拉伸裁剪，容易丢失边缘位置的有效像素。

操作时按以下逻辑调整参数即可避免跳变问题，整体运算速度远快于raster::aggregate：

• 去掉冗余的重投影参数，避免不必要的重投影流程
• 强制对齐输出栅格的像元网格，保证每个150m像元严格对应5个连续的30m像元范围，无跨边界错位
• 匹配栅格类型选择合适的重采样方法：
  • 若需要统计每个150m像元内的30m有效像元占比，选择average均值聚合方法
  • 若需要保留分类属性（即每个150m像元取窗口内出现次数最多的分类值），选择mode众数聚合方法
  • 若需要严格保持分类的二值掩膜属性，选择near最邻近重采样方法
• 开启多线程参数进一步提升聚合运算速度

修改后的可直接运行的代码如下：

library(gdalUtils)

input <- ".../globalIP/alp/footprint/pl30m.tif"
output <- ".../globalIP/alp/footprint/pl150m.tif"

gdalwarp(
  srcfile = input,
  dstfile = output,
  tr = c(150, 150),
  t_srs = NULL, # 移除和源栅格重复的投影参数，避免无效重投影
  tap = TRUE, # 自动将输出像元网格对齐到源栅格网格，省去手动计算对齐边界的步骤
  r = "average", # 可根据需求替换为mode/near
  overwrite = TRUE,
  multi = TRUE, # 开启多线程加速
  wo = c("NUM_THREADS=ALL_CPUS")
)

补充提示：如果调整后仍存在边缘像素缺失问题，可以读取源栅格的NoData值，通过dstnodata参数手动指定输出栅格的无值范围，避免gdalwarp自动将边缘有效像素识别为NoData值。

内容的提问来源于stack exchange，提问作者dwiz