在R语言中统计栅格单元内的多边形数量

问题背景 #

需要统计每个栅格单元内的多边形数量，尝试两种方案均失败，期望用Terra包实现，也接受其他包方案。

问题复现及报错 #

方案1：terra::rasterize #

运行以下代码后出现警告，生成空栅格：

pacman::p_load(terra, raster, rnaturalearth, rts, sf)

# 加载测试多边形
cntrs <- ne_countries(country = c("united kingdom", "argentina", "singapore"), type = "countries")
# 投影转换
cntrs <- project(vect(cntrs), "+proj=moll +lon_0=0 +x_0=0 +y_0=0")
# 创建空栅格
r <- terra::rast(xmin = -180, xmax = 180, ymin = -90, ymax = 90,
          res = 10,
          crs = "+proj=moll +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")
r[] <- 0

# 执行rasterize
terra::rasterize(cntrs, r, fun = "count")

警告信息：

Warning message:
Failed to compute min/max, no valid pixels found in sampling. (GDAL error 1)

方案2：cellFromXY #

转换多边形格式后运行cellFromXY报错：

# 转换多边形格式
cntrs <- as(st_geometry(st_as_sf(cntrs)), "Spatial")
table(cellFromXY(r, cntrs))

错误信息：

Error in (function (classes, fdef, mtable) :
unable to find an inherited method for function ‘cellFromXY’ for signature ‘"SpatRaster", "SpatialPolygons"’

将SpatRaster转为RasterLayer后仍报错：

r <- raster::raster(r)
table(cellFromXY(r, cntrs))

错误信息：

Error: Not compatible with requested type: [type=S4; target=double].

核心原因与解决方法 #

问题出在栅格的范围参数与投影不匹配：Mollweide投影的单位是米，但最初创建栅格时用了经纬度范围（-180到180、-90到90），导致多边形与栅格无重叠。

修正后的栅格创建代码如下：

r <- rast(xmin = -18000000, xmax = 18000000, ymin = -9000000, ymax = 9000000, res = 10000, crs = "+proj=moll +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")

使用该栅格再执行terra::rasterize(cntrs, r, fun = "count")即可得到正确的统计结果。

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