如何按像素覆盖比例对含生物量的矢量多边形进行栅格化？

我有一组代表树冠的矢量多边形，每个多边形都关联有生物量值。我希望将该数据集栅格化为10m×10m的栅格，每个像素存储落在其中的生物量总和。

核心需求是根据树冠多边形与像素的重叠比例分配生物量：若某树冠90%与一个像素重叠，该像素应获得其90%的生物量，剩余10%分配给其他重叠像素。部分树冠大于像素，部分小于像素。

目前我通过栅格化多边形质心实现，但这种方式无法反映树冠实际空间范围，导致生物量分配不准确。当前使用R语言的代码如下：

library(terra)
library(sf)

v <- st_read("...tree_vector.geojson")

resolution <- 10
original_ext <- ext(v)
# 生成栅格模板
template <- rast(extent = original_ext, resolution = resolution, crs = st_crs(v)$wkt)

# 提取多边形质心
point <- st_centroid(v)

# 栅格化质心并求和
r_points <- rasterize(point, template, field = "total", fun = "sum")
# 将NA值替换为0
r_final <- classify(r_points, cbind(NA, 0))

要实现按重叠比例分配生物量的栅格化，推荐两种实用方法：

方法一：使用exactextractr包（高效简洁） #

exactextractr专门针对矢量-栅格的精确面积加权计算，能直接输出每个像素的生物量总和：

library(terra)
library(sf)
library(exactextractr)

# 读取矢量数据
v <- st_read("...tree_vector.geojson")

# 生成10m分辨率的栅格模板
resolution <- 10
template <- rast(ext(v), res = resolution, crs = st_crs(v)$wkt)

# 按面积加权计算每个像素的生物量总和
r_final <- exact_extract(v, template, "weighted_sum", weights = "area", column = "total")

# 将NA值替换为0（可选）
r_final <- classify(r_final, cbind(NA, 0))

代码说明 #

• exact_extract的weighted_sum参数会自动计算每个树冠多边形与像素的重叠面积比例，用该比例作为权重乘以生物量值，再对每个像素的所有贡献值求和。
• weights = "area"指定用重叠面积比例作为权重，确保生物量严格按空间占比分配。

方法二：基于terra和sf的手动计算（无需额外包） #

如果不想安装新包，可以通过矢量交集计算实现：

library(terra)
library(sf)

v <- st_read("...tree_vector.geojson")
resolution <- 10

# 生成栅格模板并转换为像素多边形
template <- rast(ext(v), res = resolution, crs = st_crs(v)$wkt)
pixel_polys <- as.polygons(template, dissolve = FALSE) %>% st_as_sf()
pixel_polys$pixel_id <- seq(nrow(pixel_polys)) # 给每个像素分配唯一ID

# 计算树冠与像素的交集
intersections <- st_intersection(v, pixel_polys)

# 计算每个交集的面积占对应树冠总面积的比例
intersections$area_ratio <- st_area(intersections) / st_area(v[intersections$ID, ])

# 计算每个交集对应的生物量分配值
intersections$allocated_biomass <- intersections$total * intersections$area_ratio

# 按像素ID汇总生物量
pixel_biomass <- aggregate(allocated_biomass ~ pixel_id, data = intersections, sum)

# 将汇总结果关联回像素多边形，再转换为栅格
pixel_polys <- merge(pixel_polys, pixel_biomass, by = "pixel_id", all.x = TRUE)
pixel_polys$allocated_biomass[is.na(pixel_polys$allocated_biomass)] <- 0

r_final <- rasterize(pixel_polys, template, field = "allocated_biomass")

代码说明 #

1. 将栅格模板转换为像素多边形，方便后续空间交集计算。
2. 通过st_intersection得到每个树冠与像素的重叠区域，计算该区域占树冠总面积的比例。
3. 用比例乘以树冠生物量得到该像素应分配的生物量，最后按像素汇总并转换回栅格。

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