如何计算numpy数组掩码与shapely多边形的交集面积

共有两种成熟方案可实现需求，你可以根据精度要求、代码复杂度偏好选择： #

方案1：修正现有掩码转多边形的坐标映射问题 #

你当前遇到的尺寸不匹配问题，本质是mask_to_polygons返回的像素坐标范围和输入的thr数组尺寸不匹配，需要先做缩放再映射到实际地理坐标：

import numpy as np
from shapely.affinity import affine_transform, translate

# 1. 计算单像素对应的实际地理尺寸
pix_width = (x_right - x_left) / thr.shape[1]
pix_height = (y_top - y_bottom) / thr.shape[0]

# 2. 转换掩码为多边形后先缩放，把0~511的坐标范围缩放到匹配256x256的0~255范围
raw_mask_pol = mask_to_polygons(np.flip(thr, 0))
scaled_mask_pol = affine_transform(raw_mask_pol, [0.5, 0, 0, 0.5, 0, 0])

# 3. 平移到实际地理坐标范围
mask_pol = translate(scaled_mask_pol, xoff=x_left, yoff=y_bottom)

# 4. 直接求交集面积
target_area = building.intersection(mask_pol).area

该方案的优势是精度更高，适合对面积计算精度要求极高的场景。 #

方案2：将shapely多边形栅格化为掩码计算交集（更推荐，代码更简洁不易出错） #

如果没有极致的精度要求，直接把shapely多边形转成和thr同尺寸的掩码，通过像素统计计算面积即可，完全不需要处理复杂的坐标对齐问题：

import numpy as np
from rasterio.features import rasterize
from affine import Affine

# 1. 计算单像素对应的实际地理尺寸，定义地理坐标到像素坐标的转换矩阵
pix_width = (x_right - x_left) / thr.shape[1]
pix_height = (y_top - y_bottom) / thr.shape[0]
transform = Affine(pix_width, 0, x_left, 0, -pix_height, y_top)

# 2. 将building多边形栅格化为和thr同尺寸的二值掩码
building_mask = rasterize(
    [(building, 1)],
    out_shape=thr.shape,
    transform=transform,
    fill=0,
    dtype=np.uint8
)

# 3. 统计两个掩码同时为真的像素数，乘单像素面积得到最终交集面积
intersection_count = ((thr > 0) & (building_mask > 0)).sum()
target_area = intersection_count * pix_width * pix_height

该方案代码量更小，不会出现多边形转换的坐标偏移问题，适合绝大多数工程应用场景。 #

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