使用terra::extract计算栅格统计值结果不一致的原因咨询

背景数据 #

市政多边形（shp_muni） #

shp_muni
class       : SpatVector 
geometry    : polygons 
dimensions  : 1, 1  (geometries, attributes)
extent      : 5.08748, 5.467728, 52.43524, 52.69111  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 
names       :   PC4
type        : <int>
values      :  8244

全球覆盖栅格（temp_rast） #

class       : SpatRaster 
dimensions  : 696, 1440, 1  (nrow, ncol, nlyr)
resolution  : 0.25, 0.25  (x, y)
extent      : -180, 180, -90, 84  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 
source      : memory 
name        : year_value 
min value   :        1 
max value   :      100 

执行的代码与结果 #

方法1：带exact=T的提取 #

cells_score <- terra::extract(temp_rast, shp_muni, exact = T, cells = T, na.rm = TRUE, as.raster = F)
cells_score

输出：

ID year_value   cell   fraction
 1       14      180741 0.24591428
 1       14      180742 0.50551934
 1       17      182181 0.03271280
 1       16      182182 0.09146393

由此计算得：最小值14，最大值17，均值可通过fraction加权计算。

方法2：带"max"/"min"和touches=T的提取 #

terra::extract(temp_rast, shp_muni, "max", touches = T, na.rm = TRUE, as.raster = F)

输出：

ID year_value
1       17

terra::extract(temp_rast, shp_muni, "min", touches = T, na.rm = TRUE, as.raster = F)

输出：

ID year_value
1        0

差异原因 #

两种方法的核心区别在于对「纳入统计的像元」的判定规则：

• 方法1（exact=T）：仅提取与多边形**实际存在重叠区域（即重叠比例fraction>0）**的像元，返回的都是真正和市政区域有交集的像元，因此统计的最值仅基于这些有效像元的数值。
• 方法2（touches=T）：会纳入所有与多边形边界或内部有任何接触的像元——包括那些仅和多边形边界「擦边」但实际重叠面积为0的像元。你得到的最小值0，就是来自这类被错误纳入、但exact=T会过滤掉的像元（这类像元的year_value为0，且在exact=T时因fraction=0不被返回）。

最大值一致是因为两种方法都包含了那个值为17的像元，它在两种判定规则下都属于有效纳入范围。

如果要让touches=T的统计结果和exact=T一致，建议先通过exact=T获取所有重叠比例>0的像元，再在这些像元中计算最值，而不是直接用extract的内置统计参数。

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