如何在Python中实现R terra包的栅格面加权平均功能？

先纠正你R代码的小笔误：Terra中实现该功能的函数是terra::extract()而非terra::rast()，正确代码应为：

library(terra)
terra::extract(my_raster, my_shp, fun = 'mean', na.rm = T, weights = T, touches = T)

针对你的需求，以下是两种Python实现方案：

方法1：使用rasterstats库（最简方案） #

rasterstats是专门处理栅格-矢量统计的工具库，原生支持基于交集面积的加权均值计算，完全对应Terra中weights=T+touches=T的逻辑。

安装依赖 #

pip install rasterstats

实现代码 #

import geopandas as gpd
from rasterstats import zonal_stats

# 读取矢量文件与栅格路径
gdf = gpd.read_file("你的矢量文件路径.shp")
raster_path = "你的栅格文件路径.tif"

# 计算加权均值：weighted=True对应weights=T，all_touched对应touches=T
stats_results = zonal_stats(
    gdf,
    raster_path,
    stats="mean",
    all_touched=True,
    weighted=True,
    nodata=-9999  # 替换为你栅格实际的nodata值
)

# 提取结果（若有多边形，需遍历stats_results）
weighted_mean = stats_results[0]["mean"]

方法2：手动计算（无额外依赖） #

若不想新增库依赖，可通过rasterio+shapely手动计算每个栅格单元与多边形的交集面积，再加权求和得到均值：

import rasterio
import geopandas as gpd
import numpy as np
from shapely.geometry import box
from rasterio.features import geometry_mask

def calculate_weighted_mean(raster_path, shp_path):
    # 读取栅格数据
    with rasterio.open(raster_path) as src:
        raster_data = src.read(1)  # 读取第一波段
        nodata_val = src.nodata
        transform = src.transform
        # 计算单个栅格单元的面积
        cell_width = transform[0]
        cell_height = abs(transform[4])
        cell_area = cell_width * cell_height

    # 读取矢量多边形
    gdf = gpd.read_file(shp_path)
    target_polygon = gdf.geometry.iloc[0]

    # 生成与多边形接触的栅格单元掩码（all_touched=True）
    mask = geometry_mask(
        [target_polygon],
        out_shape=raster_data.shape,
        transform=transform,
        all_touched=True,
        invert=True
    )
    # 获取所有符合条件的栅格单元索引
    y_idx, x_idx = np.where(mask)

    total_weighted_sum = 0.0
    total_intersect_area = 0.0

    for y, x in zip(y_idx, x_idx):
        # 构建当前栅格单元的边界框
        x_min, y_max = transform * (x, y)
        x_max, y_min = transform * (x + 1, y + 1)
        cell_bbox = box(x_min, y_min, x_max, y_max)
        # 计算单元与多边形的交集面积
        intersect_region = cell_bbox.intersection(target_polygon)
        intersect_area = intersect_region.area
        # 跳过nodata值
        current_val = raster_data[y, x]
        if current_val == nodata_val:
            continue
        # 累加加权值与面积
        total_weighted_sum += current_val * intersect_area
        total_intersect_area += intersect_area

    # 返回加权均值，无有效数据时返回NaN
    return total_weighted_sum / total_intersect_area if total_intersect_area > 0 else np.nan

# 调用函数计算结果
weighted_mean = calculate_weighted_mean("你的栅格文件路径.tif", "你的矢量文件路径.shp")

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