Shapely多边形与多多边形交集二进制掩码生成异常排查

目标 #

创建一个二维二进制数组，用于表示正方形Polygon与MultiPolygon的交点对应的点集。

现有代码 #

from typing import List, Tuple
import numpy.typing as npt
from shapely.geometry import Polygon, Point
from shapely.geometry.multipolygon import MultiPolygon
import numpy as np  # 原代码遗漏numpy导入

def binary_mask(tile_poly: Polygon, wkt_mp: MultiPolygon, tile_shape: Tuple[int, int] = (256, 256)) -> npt.NDArray[np.int_]:
    mask = np.zeros(tile_shape)
    for x in range(tile_shape[0]):
        for y in range(tile_shape[1]):
            if tile_poly.contains(Point(x,y)) or wkt_mp.contains(Point(x,y)):
                mask[x, y] = 1

    return mask

问题现象 #

所有生成的数组全为0：

mask = binary_mask(tile_poly, wkt_mp, tile_shape[:2])
print(np.unique(mask))
>>> [0.]

按逻辑，多数tile_poly完全处于wkt_mp内部，应该返回全1数组；边界处的tile_poly，超出wkt_mp的点对应掩码值应为0，但实际结果不符合预期。

数据信息 #

crs='EPSG:4326'

qupath_poly_flip_y[0]
>>> [<shapely.geometry.polygon.Polygon at 0x7fe4fe9f6af0>,
     <shapely.geometry.polygon.Polygon at 0x7fe4fe9f6bb0>,
     ...]

qupath_poly_flip_y[0][0].area  # 对应256x256的瓦片
>>> 65536.0

qupath_poly_flip_y[0][0].bounds
>>> (19200.0, 11336.0, 19456.0, 11592.0)

wkt_multipoly[0].area
>>> 17431020.0

wkt_multipoly[0].bounds
>>> (8474.0, 10026.0, 21732.0, 17267.0)

问题原因与解决方案 #

核心问题 #

代码里的Point(x,y)用的是像素坐标（0-255范围），但你的多边形是基于**地理坐标（如19200、11336这类大数值）**创建的，两者坐标系完全不匹配，因此永远不会产生包含关系，导致掩码全为0。

修复步骤 #

1. 像素坐标转地理坐标：根据tile_poly的地理边界，将每个像素位置映射到对应的地理坐标点。
2. 性能优化：原代码双重循环效率极低，优先采用向量化操作提升速度。

修复后的基础版本代码 #

from typing import List, Tuple
import numpy as np
import numpy.typing as npt
from shapely.geometry import Polygon, Point
from shapely.geometry.multipolygon import MultiPolygon

def binary_mask(tile_poly: Polygon, wkt_mp: MultiPolygon, tile_shape: Tuple[int, int] = (256, 256)) -> npt.NDArray[np.int_]:
    mask = np.zeros(tile_shape, dtype=np.int_)
    # 获取tile_poly的地理边界范围
    min_x, min_y, max_x, max_y = tile_poly.bounds
    # 计算像素到地理坐标的缩放系数
    x_step = (max_x - min_x) / tile_shape[0]
    y_step = (max_y - min_y) / tile_shape[1]
    
    # 遍历每个像素，转换为地理坐标后判断归属
    for x_idx in range(tile_shape[0]):
        for y_idx in range(tile_shape[1]):
            # 取像素中心对应的地理坐标
            geo_x = min_x + x_idx * x_step + x_step/2
            geo_y = min_y + y_idx * y_step + y_step/2
            point = Point(geo_x, geo_y)
            if tile_poly.contains(point) or wkt_mp.contains(point):
                mask[x_idx, y_idx] = 1
    return mask

高性能优化版本 #

如果需要处理大量瓦片，使用shapely.vectorized.contains实现向量化判断，避免循环：

import numpy as np
import numpy.typing as npt
from shapely.geometry import Polygon, MultiPolygon
from shapely.vectorized import contains

def binary_mask_fast(tile_poly: Polygon, wkt_mp: MultiPolygon, tile_shape: Tuple[int, int] = (256, 256)) -> npt.NDArray[np.int_]:
    min_x, min_y, max_x, max_y = tile_poly.bounds
    x_step = (max_x - min_x) / tile_shape[0]
    y_step = (max_y - min_y) / tile_shape[1]
    
    # 生成所有像素中心的地理坐标网格
    xs = np.linspace(min_x + x_step/2, max_x - x_step/2, tile_shape[0])
    ys = np.linspace(min_y + y_step/2, max_y - y_step/2, tile_shape[1])
    xx, yy = np.meshgrid(xs, ys, indexing='ij')
    
    # 向量化判断点是否在多边形内
    in_tile = contains(tile_poly, xx.ravel(), yy.ravel()).reshape(tile_shape)
    in_mp = contains(wkt_mp, xx.ravel(), yy.ravel()).reshape(tile_shape)
    
    return np.where(in_tile | in_mp, 1, 0).astype(np.int_)

验证方法 #

先验证单个点的坐标转换是否正确：

# 取tile_poly中心的地理坐标
min_x, min_y, max_x, max_y = tile_poly.bounds
center_point = Point((min_x+max_x)/2, (min_y+max_y)/2)
print(wkt_mp.contains(center_point))  # 若返回True，说明坐标转换逻辑正确

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