R中rasterize()生成栅格与多边形值不符，提取结果异常

我尝试栅格化多边形图层，用于基于点缓冲区提取变量平均值的分区统计，期望分辨率约100m。已知分辨率变化会影响栅格值，但希望栅格单元值尽可能贴近原多边形数据。目前遇到两个问题：

• 栅格值远高于多边形实际值：minmax()显示栅格最大值为518，但多边形最大值仅约48
• 裁剪至缓冲区范围后，用extract()或exact_extract()提取的均值与可视化结果完全不符（例：R计算的缓冲区均值为203，和可视化呈现的内容不匹配）

推测问题出在栅格创建环节，尝试创建可复现示例未成功。

栅格创建代码 #

template <- rast(vect(sf_data), resolution=100, crs="EPSG:32610") # sf_data也采用该投影
raster <- terra::rasterize(vect(sf_data), template, field = "Variable")

sf_data为多边形图层，栅格可视化看似正常，但数值存在偏差

均值提取代码 #

site_index <- 10
pt <- sites[site_index, ]
buffer_radius <- 500
buff_terra <- vect(pt) %>% terra::buffer(width = buffer_radius) 

# 裁剪并掩膜栅格
rast_cropped <- crop(raster, buff_terra, snap = "out")

# 提取值并计算均值
extracted_value <- exact_extract(rast_cropped, st_as_sf(buff_terra), fun = "mean", 
                                 weights = "area")
# 也曾尝试：
extracted_value <- terra::extract(rast_cropped, vect(buff_terra), fun = "mean")

1. 栅格化数值异常核心排查 #

• 重叠多边形赋值逻辑问题：terra::rasterize()默认用last规则（后处理的多边形覆盖先处理的），若多边形大量重叠，可能导致异常赋值。先检查重叠情况：

# 检测非自交的多边形重叠
overlaps <- terra::intersect(vect(sf_data), vect(sf_data))
nrow(overlaps[overlaps$ID != overlaps$ID.1, ])

• 投影单位验证：确认sf_data的投影EPSG:32610单位为米，若单位不符，resolution=100会导致实际栅格尺寸偏差：

# 查看投影单位
crs(vect(sf_data), describe=TRUE)$units

• 字段类型检查：确保Variable为数值型，字符型字段会被强制转整数导致数值失真：

class(sf_data$Variable)

2. 栅格化参数优化 #

针对“贴近原多边形数据”的需求，调整栅格化规则：

• 若存在多边形重叠，用mean计算重叠区域的平均值替代默认last：

raster <- terra::rasterize(vect(sf_data), template, field = "Variable", fun = "mean")

• 若多边形面积远小于栅格单元，启用cover=TRUE，仅对完全被多边形覆盖的栅格单元赋值（或适当降低分辨率适配多边形尺寸）：

raster <- terra::rasterize(vect(sf_data), template, field = "Variable", cover=TRUE)

3. 分区统计修正 #

• 裁剪后必须添加掩膜，排除缓冲区外的栅格值干扰：

rast_cropped <- crop(raster, buff_terra, snap = "out") %>% mask(buff_terra)

• 使用terra::extract()时，启用weights=TRUE计算面积加权均值，和exact_extract()逻辑对齐：

extracted_value <- terra::extract(rast_cropped, vect(buff_terra), fun = "mean", weights=TRUE)

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