fix: 掩膜改用 rasterio.features.rasterize 全分辨率 C 级光栅化

问题: 降采样+np.kron 暴力放大的旧逻辑导致出图边缘出现
  14×14 像素马赛克锯齿。

修复: 彻底删除 _MASK_TARGET_POINTS/_mask_step/shapely.prepared/
  np.kron 升采样链,替换为 rasterio C 底层光栅化:
  - from_bounds 精确计算仿射矩阵
  - rasterize(shapes, out_shape, transform) C 级瞬间盖章
  - flipud 处理 Y 轴对齐
  - 输出 uint8 bool 掩膜, 无降采样无锯齿
This commit is contained in:
duxin
2026-07-08 15:02:30 +08:00
parent e0c6b63446
commit 757741c6a9

View File

@ -2978,39 +2978,39 @@ class ContentMapper:
print(f"[共享上下文] 网格: {nx}×{ny} = {nx*ny}")
# ⑥ 水域掩膜布尔矩阵(降采样到 ~100m 分辨率计算,避免 1m 下千万级 Point 对象 OOM
# ⑥ 水域掩膜布尔矩阵(全分辨率高精度光栅化,消灭锯齿
mask = None
if boundary_gdf is not None:
_MASK_TARGET_POINTS = 50000 # 掩膜降采样目标点数
_mask_step = max(1, int(np.sqrt(grid_xx.size / _MASK_TARGET_POINTS)))
if _mask_step > 1:
mask_xx = grid_xx[::_mask_step, ::_mask_step]
mask_yy = grid_yy[::_mask_step, ::_mask_step]
print(f"[共享上下文] 掩膜降采样 {_mask_step}×"
f"{mask_xx.shape[1]}×{mask_xx.shape[0]} = {mask_xx.size:,}")
else:
mask_xx, mask_yy = grid_xx, grid_yy
import rasterio.features
from rasterio.transform import from_bounds
mask_pts = np.column_stack((mask_xx.ravel(), mask_yy.ravel()))
print(" [共享上下文] 启动 Rasterio全分辨率精确光栅化水体边界...")
# 使用 shapely.prepared 预编译几何体,空间查询加速上百倍
from shapely.prepared import prep
from shapely.geometry import Point
print(" [共享上下文] 预编译水域边界 (C 级空间索引加速)...")
union_poly = boundary_gdf.unary_union
prepared_poly = prep(union_poly)
min_x, max_x = float(grid_xx.min()), float(grid_xx.max())
min_y, max_y = float(grid_yy.min()), float(grid_yy.max())
ny, nx = grid_xx.shape
print(" [共享上下文] 执行快速空间相交测试...")
is_within = [prepared_poly.contains(Point(x, y))
for x, y in mask_pts]
mask_lowres = np.array(is_within).reshape(mask_xx.shape)
# 升采样回原始分辨率(最近邻,掩膜是布尔值)
if _mask_step > 1:
mask = np.kron(mask_lowres, np.ones((_mask_step, _mask_step), dtype=bool))
# 裁剪到精确原始尺寸
mask = mask[:grid_xx.shape[0], :grid_xx.shape[1]]
else:
mask = mask_lowres
# 计算像元大小并生成仿射变换矩阵
dx = (max_x - min_x) / max(1, nx - 1)
dy = (max_y - min_y) / max(1, ny - 1)
transform = from_bounds(min_x - dx / 2, min_y - dy / 2,
max_x + dx / 2, max_y + dy / 2, nx, ny)
# 调用 C 语言底层瞬间完成千万级像素盖章
mask_raster = rasterio.features.rasterize(
shapes=boundary_gdf.geometry.tolist(),
out_shape=(ny, nx),
transform=transform,
fill=0,
default_value=1,
dtype='uint8',
)
# 坐标系 Y 轴方向对齐处理
if grid_yy[0, 0] < grid_yy[-1, 0]:
mask_raster = np.flipud(mask_raster)
mask = mask_raster.astype(bool)
print(f"[共享上下文] 水域掩膜: {int(mask.sum())}/{mask.size} 点在水域内")
return (grid_xx, grid_yy, mask, bounds, boundary_gdf)