From 757741c6a9a54f3d4e56987f181d6a15ce30771b Mon Sep 17 00:00:00 2001
From: duxin <duxin@example.com>
Date: Wed, 8 Jul 2026 15:02:30 +0800
Subject: [PATCH] =?UTF-8?q?fix:=20=E6=8E=A9=E8=86=9C=E6=94=B9=E7=94=A8=20r?=
 =?UTF-8?q?asterio.features.rasterize=20=E5=85=A8=E5=88=86=E8=BE=A8?=
 =?UTF-8?q?=E7=8E=87=20C=20=E7=BA=A7=E5=85=89=E6=A0=85=E5=8C=96?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

问题: 降采样+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 掩膜, 无降采样无锯齿
---
 src/postprocessing/map.py | 56 +++++++++++++++++++--------------------
 1 file changed, 28 insertions(+), 28 deletions(-)

diff --git a/src/postprocessing/map.py b/src/postprocessing/map.py
index b93c236..dc9558c 100644
--- a/src/postprocessing/map.py
+++ b/src/postprocessing/map.py
@@ -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)