ContentMapper 边界读取支持栅格水掩膜（.dat/.bsq/.tif/.tiff/.img）

src/postprocessing/map.py

@@ -15,7 +15,7 @@ from scipy.spatial.distance import cdist
 from scipy.spatial import ConvexHull
 from shapely.geometry import Point, Polygon
 import rasterio
-from rasterio.features import geometry_mask
+from rasterio.features import geometry_mask, shapes
 from rasterio import windows
 from rasterio.warp import calculate_default_transform, reproject, Resampling
 try:
@@ -785,18 +785,82 @@ class ContentMapper:
         return gdf
 
     def read_boundary_shapefile(self, shp_file):
-        """读取边界shapefile"""
-        print("正在读取边界文件...")
-        boundary = gpd.read_file(shp_file)
+        """读取边界/掩膜文件（同时支持矢量 .shp 与栅格 .dat/.bsq/.tif/.tiff）。
+
+        - .shp → gpd.read_file 读取矢量边界（保持原行为）
+        - .dat/.bsq/.tif/.tiff/.img → rasterio 读取栅格水掩膜 → rasterio.features.shapes
+          矢量化成水体多边形 → gpd.GeoDataFrame 返回
+
+        下游 create_interpolation_grid / create_content_map / visualize_raster
+        始终接收 GeoDataFrame，无需任何改动。
+        """
+        print("正在读取边界/掩膜文件...")
+        suffix = Path(shp_file).suffix.lower()
+        if suffix in (".shp",):
+            boundary = gpd.read_file(shp_file)
+        elif suffix in (".dat", ".bsq", ".tif", ".tiff", ".img"):
+            boundary = self._raster_to_boundary_gdf(shp_file)
+        else:
+            raise ValueError(
+                f"不支持的边界/掩膜文件格式: {suffix}（仅支持 .shp / .dat / .bsq / .tif / .img）"
+            )
+
+        if len(boundary) == 0:
+            raise ValueError(
+                f"边界/掩膜 {shp_file} 矢量化后为空（栅格格式请确认 .dat 包含水体像元 > 0）"
+            )
 
         # 确保边界文件使用目标投影坐标系
-        if boundary.crs != self.output_crs:
-            print(f"正在转换边界文件坐标系到 {self.output_crs}...")
+        if boundary.crs is not None and boundary.crs != self.output_crs:
+            print(f"正在转换边界/掩膜坐标系到 {self.output_crs}...")
             boundary = boundary.to_crs(self.output_crs)
 
-        print(f"边界文件包含 {len(boundary)} 个要素")
+        print(f"边界/掩膜文件包含 {len(boundary)} 个要素")
         return boundary
 
+    def _raster_to_boundary_gdf(self, raster_path):
+        """把栅格二值水掩膜（.dat/.bsq/.tif/.tiff）矢量化成水体多边形 GeoDataFrame。
+
+        修复 Step 11 接收 step1 产出 .dat 水掩膜的兼容性：
+          - rasterio.open 读 band 1（0=非水, 1/任意正数=水）
+          - rasterio.features.shapes 矢量化成多边形
+          - 收集所有 val=1 的多边形 → gpd.GeoDataFrame
+        """
+        try:
+            from shapely.geometry import shape as _shapely_shape
+        except ImportError as e:
+            raise ImportError(
+                "栅格掩膜矢量化需要 shapely（geopandas 自带）。原始错误: " + str(e)
+            )
+
+        with rasterio.open(raster_path) as src:
+            data = src.read(1)
+            transform = src.transform
+            crs = src.crs
+
+        # 二值化：>0 视为水
+        mask_uint8 = (data > 0).astype(np.uint8)
+        if int(mask_uint8.sum()) == 0:
+            raise ValueError(
+                f"栅格掩膜 {raster_path} 中无水体像元（>0），无法矢量化"
+            )
+
+        # 矢量化：shapes 返回 (geom_dict, value) 迭代器
+        geoms = []
+        for geom_dict, val in shapes(mask_uint8, mask=mask_uint8.astype(bool), transform=transform):
+            if int(val) == 1:
+                geoms.append(_shapely_shape(geom_dict))
+
+        if not geoms:
+            raise ValueError(f"栅格掩膜 {raster_path} 矢量化后无有效水体多边形")
+
+        gdf = gpd.GeoDataFrame(geometry=geoms, crs=crs)
+        print(
+            f"栅格掩膜 {Path(raster_path).name} 矢量化完成: "
+            f"{int(mask_uint8.sum())} 个水体像元 → {len(gdf)} 个多边形 (CRS={crs})"
+        )
+        return gdf
+
     def _identify_edge_points(self, points_gdf):
         """
         识别边缘采样点（使用凸包方法）