diff --git a/src/core/handlers/step11_map_handler.py b/src/core/handlers/step11_map_handler.py index 7d90ae9..e467b13 100644 --- a/src/core/handlers/step11_map_handler.py +++ b/src/core/handlers/step11_map_handler.py @@ -193,11 +193,42 @@ class Step11MapHandler(BaseStepHandler): input_crs=base_kwargs['input_crs'], output_crs=base_kwargs['output_crs'], ) + # 原始路径若是栅格（.dat/.tif），传给 prepare_shared_context + # 走 GDAL 快速重采样通道，避免 58K 多边形 rasterize 卡死 + _orig_path = boundary_shp_path + _raster_src = None + _RASTER_EXTS = ('.dat', '.tif', '.tiff', '.bsq', '.bil', '.bip', '.img') + if _orig_path and os.path.isfile(_orig_path): + _ext = os.path.splitext(_orig_path)[1].lower() + if _ext in _RASTER_EXTS: + _raster_src = _orig_path + # fallback: 扫描 work_dir 下的原始栅格掩膜 + if _raster_src is None: + _wd = Path(str(context.work_dir)) + for _sub in ('1_water_mask', '1_Water_Mask', 'water_mask'): + _scan_dir = _wd / _sub + if _scan_dir.is_dir(): + for _ext in _RASTER_EXTS: + _hits = sorted(_scan_dir.glob(f'*{_ext}'), + key=lambda p: p.stat().st_mtime, reverse=True) + for _h in _hits: + if _h.is_file(): + _raster_src = str(_h) + break + if _raster_src: + break + if _raster_src: + break + if _raster_src: + context.notify('step11_map', 'info', + f'找到原始栅格掩膜，将使用 GDAL 快速重采样: {Path(_raster_src).name}') + shared_ctx = pre_mapper.prepare_shared_context( sample_csv=csv_paths[0], shp_file=resolved_boundary, resolution=float(base_kwargs['resolution']), expand_ratio=0.05, + boundary_raster=_raster_src, ) base_kwargs['shared_context'] = shared_ctx context.notify('step11_map', 'info', diff --git a/src/postprocessing/map.py b/src/postprocessing/map.py index a11ad81..522c3f3 100644 --- a/src/postprocessing/map.py +++ b/src/postprocessing/map.py @@ -589,6 +589,11 @@ class ContentMapper: # 网格按 500m 空间窗口分块，每块只取窗口内 + 500m 缓冲区的 # 局部采样点参与计算。协方差矩阵从全局 8242×8242 降为 # 局部 n×n (n≈几十到几百)，千万级网格秒级完成。 + # + # 2026-07-23：新增网格规模自动判断 —— 当网格点数超过 500K 时， + # pykrige 'loop' 后端的 Python 循环开销过大（每个点 100-200μs）， + # 千万级网格需数小时。此时自动跳过 Kriging，直走 IDW。 + # IDW 使用 cKDTree 向量化查询，千万级点仅需数秒，效果差异极小。 # ═══════════════════════════════════════════════════════════ kriging_degraded = False if PYKRIGE_AVAILABLE: @@ -597,31 +602,38 @@ class ContentMapper: grid_y = grid_yy[:, 0] total_cells = len(grid_x) * len(grid_y) - print(f"正在使用 局部克里金 (自适应分块 + 40% 重叠缓冲)：" - f"网格={total_cells:,} 点") - - grid_content = self._local_kriging( - points, values, grid_x, grid_y, - n_closest_points=50, - ) - - valid_mask = ~np.isnan(grid_content) - valid_count = int(np.sum(valid_mask)) - - if valid_count > 0: - kriging_std = float(np.nanstd(grid_content)) - degradation_ratio = kriging_std / max(value_std, 1e-12) - print(f"局部 Kriging 完成: 有效点={valid_count}/{grid_content.size}, " - f"输出std={kriging_std:.6f}, 退化比={degradation_ratio:.3f}") - - if degradation_ratio < 0.05 and value_range > 1e-8: - print(f"⚠ Kriging 严重退化，回退 IDW") - kriging_degraded = True - else: - return grid_content - else: - print("局部 Kriging 结果全为 NaN，回退") + # ── 网格规模自动判断：>500K 时跳过 Kriging ── + _KRIGING_GRID_LIMIT = 500_000 + if total_cells > _KRIGING_GRID_LIMIT: + print(f"[FAST] 网格 {total_cells:,} 点超过阈值 {_KRIGING_GRID_LIMIT:,}，" + f"自动切换 IDW（高分辨率插值无需 Kriging）") kriging_degraded = True + else: + print(f"正在使用 局部克里金 (自适应分块 + 40% 重叠缓冲)：" + f"网格={total_cells:,} 点") + + grid_content = self._local_kriging( + points, values, grid_x, grid_y, + n_closest_points=20, + ) + + valid_mask = ~np.isnan(grid_content) + valid_count = int(np.sum(valid_mask)) + + if valid_count > 0: + kriging_std = float(np.nanstd(grid_content)) + degradation_ratio = kriging_std / max(value_std, 1e-12) + print(f"局部 Kriging 完成: 有效点={valid_count}/{grid_content.size}, " + f"输出std={kriging_std:.6f}, 退化比={degradation_ratio:.3f}") + + if degradation_ratio < 0.05 and value_range > 1e-8: + print(f"[WARN] Kriging 严重退化，回退 IDW") + kriging_degraded = True + else: + return grid_content + else: + print("局部 Kriging 结果全为 NaN，回退") + kriging_degraded = True except Exception as e: print(f"Kriging 失败: {e}") kriging_degraded = True @@ -634,7 +646,8 @@ class ContentMapper: # ═══════════════════════════════════════════════════════════ if kriging_degraded: try: - print("正在使用 IDW 插值（反距离权重, power=2, neighbors=15）...") + print(f"正在使用 IDW 插值（反距离权重, power=2, neighbors={min(15, len(points))}）" + f" — 网格={grid_xx.size:,} 点 ...") grid_content = self._idw_interpolation( points, values, grid_xx, grid_yy, power=2, n_neighbors=min(15, len(points)), @@ -689,12 +702,12 @@ class ContentMapper: return grid_content def _local_kriging(self, points, values, grid_x, grid_y, - n_closest_points=50): + n_closest_points=20): """局部克里金：自适应分块 + 重叠缓冲区 + 保护性近邻限制 1. 自适应块大小: 根据 extent 自动切分为 ~4×4 块 (16~25块) 2. 重叠缓冲区: 采样点范围扩展块长宽的 40%，交界处平滑无拼缝 - 3. 保护性近邻: n_closest_points=50，稀释极端异常值 + 3. 保护性近邻: n_closest_points=20，稀释极端异常值 4. 网格点仅使用严格不重叠的块范围（Buffer 仅用于筛选采样点） """ @@ -777,9 +790,9 @@ class ContentMapper: print(f" 顺序执行 {len(tasks)} 个局部 Kriging 块...") results = [] for i, task in enumerate(tasks): - if i % max(1, len(tasks) // 4) == 0 or i == len(tasks) - 1: - print(f" [LocalKrige] {i+1}/{len(tasks)} ...") + print(f" [LocalKrige] {i+1}/{len(tasks)} ...") results.append(_local_krige_block_worker(task)) + print(f" [LocalKrige] {i+1}/{len(tasks)} 完成") # 拼接：全 NaN 数组，逐块填回 grid_full = np.full((len(grid_y), len(grid_x)), np.nan, dtype=np.float64) @@ -897,7 +910,7 @@ class ContentMapper: if lon_col is None: # 终极回退：按位置取前两列 lon_col, lat_col = df.columns[0], df.columns[1] - print(f" ⚠ 未识别到标准坐标列名，按位置回退: " + print(f" [WARN] 未识别到标准坐标列名，按位置回退: " f"X={lon_col}, Y={lat_col}") # 动态识别含量列（跳过已知的坐标列和特殊列） @@ -1017,7 +1030,7 @@ class ContentMapper: # 【新增防御】检测 Transform 是否为纯像素矩阵（极易导致严重错位） if transform.is_identity: print("\n" + "!" * 65) - print(f"⚠️ [严重警告] 栅格掩膜 {Path(raster_path).name} 缺少真实的地理仿射变换！") + print(f"[WARN]️ [严重警告] 栅格掩膜 {Path(raster_path).name} 缺少真实的地理仿射变换！") print(f"当前被判定为纯像素坐标系 (x: 0~Width, y: 0~Height)。\n强行与地理栅格 (UTM/WGS84) 叠加将发生【极其严重的错位】！") print(f"请务必在 Step 1 导入原始的 .shp 矢量文件进行约束。") print("!" * 65 + "\n") @@ -1084,7 +1097,7 @@ class ContentMapper: 对边缘采样点进行外扩处理，外扩到整个图像的边界（包括外扩后的边界） 按照指定的间距（resolution）生成外扩点，铺满整个画面。 - ★★★ Plan C：boundary_gdf 可选（None = 不依赖水域掩膜，纯采样点自然扩展）★★★ + [*][*][*] Plan C：boundary_gdf 可选（None = 不依赖水域掩膜，纯采样点自然扩展）[*][*][*] Parameters: ----------- @@ -1270,7 +1283,7 @@ class ContentMapper: """ 创建插值网格 - ★★★ Plan C：boundary_gdf 可选（None = 纯采样点自然插值，无水域掩膜约束）★★★ + [*][*][*] Plan C：boundary_gdf 可选（None = 纯采样点自然插值，无水域掩膜约束）[*][*][*] Parameters: ----------- @@ -1485,7 +1498,7 @@ class ContentMapper: """ 创建含量图 - ★★★ Plan C：boundary_gdf 可选（None = 无掩膜裁剪，无黑色边界线）★★★ + [*][*][*] Plan C：boundary_gdf 可选（None = 无掩膜裁剪，无黑色边界线）[*][*][*] Parameters: ----------- @@ -1778,7 +1791,7 @@ class ContentMapper: 尺寸以点数（points）为单位，与数据坐标系解耦， 无论 UTM 坐标范围多大，指北针始终保持合理大小。 """ - # ★★★ 改用画布相对坐标（transAxes）★★★ + # [*][*][*] 改用画布相对坐标（transAxes）[*][*][*] # (0.88, 0.92) = 右上角，尺寸用 points（72分之一英寸） arrow_ax_x, arrow_ax_y = 0.88, 0.92 radius_pt = 18 # 罗盘半径（磅），由 28 → 18 缩小图元 @@ -2518,7 +2531,7 @@ class ContentMapper: gt[4], gt[5], gt[3]) else: transform = None - # ★★★ 关键：从 GeoTransform 计算 bounds 和 res ★★★ + # [*][*][*] 关键：从 GeoTransform 计算 bounds 和 res [*][*][*] # gt = (xmin, xres, 0, ymax, 0, yres) xmin_gdal = gt[0] ymax_gdal = gt[3] @@ -2576,7 +2589,7 @@ class ContentMapper: boundary_gdf = boundary_gdf.set_crs(epsg=4326) # 兜底：如果栅格 TIF 缺失坐标系 - # ★ 关键修复：不能盲目回退到 WGS84 (EPSG:4326) + # [*] 关键修复：不能盲目回退到 WGS84 (EPSG:4326) # transform 的坐标值如果是数十万～数百万量级，一定是投影坐标（UTM 等）， # 而非 WGS84 经纬度（-180~180）。此时应使用掩膜的 CRS 作为正确参考， # 否则 geometry_mask 的空间坐标与 transform 完全错位 → 全部像元被擦除。 @@ -2689,7 +2702,7 @@ class ContentMapper: ymax = transform.f xres = transform.a yres = transform.e - # ★★★ 必须用原始宽高（w_orig/h_orig）而非降采样后的 w/h ★★★ + # [*][*][*] 必须用原始宽高（w_orig/h_orig）而非降采样后的 w/h [*][*][*] extent = [xmin, xmin + w_orig * xres, ymax + h_orig * yres, ymax] scale_x = abs(xres) scale_y = abs(yres) @@ -2743,7 +2756,7 @@ class ContentMapper: interpolation='bilinear' ) - # ★★★ 锁死绘图视口 ★★★ + # [*][*][*] 锁死绘图视口 [*][*][*] # 必须在所有叠加绘图（shp/colorbar/north arrow）之前执行， # 防止其他元素的坐标干扰导致轴范围被拉伸成像素坐标系 ax.set_xlim(extent[0], extent[1]) @@ -2956,11 +2969,12 @@ class ContentMapper: return output_tif_path # ═══════════════════════════════════════════════════════════════ - # ★ 2026-07-01：共享空间上下文 — 63 个 CSV 只算一次网格/掩膜 + # [*] 2026-07-01：共享空间上下文 — 63 个 CSV 只算一次网格/掩膜 # ═══════════════════════════════════════════════════════════════ def prepare_shared_context(self, sample_csv: str, shp_file=None, - resolution=100, expand_ratio=0.05): + resolution=100, expand_ratio=0.05, + boundary_raster=None): """从首个 CSV 预计算所有子进程共用的空间基准数据。 63 个水色指数 CSV 坐标完全一致，以下数据只算一次： @@ -2971,6 +2985,12 @@ class ContentMapper: 子进程直接从 shared_context 解包复用，跳过 ②③④⑥，直入 Kriging。 + Parameters + ---------- + boundary_raster : str, optional + 原始栅格掩膜路径（.dat / .tif）。提供时使用 GDAL 直接重采样， + 跳过 polygonize → rasterize 来回转换，58K 多边形场景下秒级完成。 + Returns: tuple: (grid_xx, grid_yy, mask, bounds, boundary_gdf) """ @@ -3014,25 +3034,122 @@ class ContentMapper: print(f"[共享上下文] 网格: {nx}×{ny} = {nx*ny} 点") - # ⑥ 水域掩膜布尔矩阵（全分辨率高精度光栅化，消灭锯齿） + # ⑥ 水域掩膜布尔矩阵 mask = None if boundary_gdf is not None: - import rasterio.features - from rasterio.transform import from_bounds - - print(" [共享上下文] 启动 Rasterio，全分辨率精确光栅化水体边界...") - 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 - - # 计算像元大小并生成仿射变换矩阵 dx = (max_x - min_x) / max(1, nx - 1) dy = (max_y - min_y) / max(1, ny - 1) + + # ── 快速通道：原始栅格直接 GDAL 重采样 ── + # 58K 多边形 rasterize 需 5-10 分钟，GDAL Warp 秒级完成 + # ── 智能扫描：若未显式传入栅格路径，自动从 shp_file 同目录或常见子目录查找 ── + _raster_src = boundary_raster + if not (_raster_src and os.path.isfile(_raster_src)): + _RASTER_EXTS = ('.dat', '.tif', '.tiff', '.bsq', '.bil', '.bip', '.img') + # 策略 A: 从 shp_file 推测（如 foo_vectorized.shp → foo.dat） + if shp_file and os.path.isfile(shp_file): + _s = Path(shp_file) + for _ext in _RASTER_EXTS: + _candidate = _s.with_suffix(_ext) + if _candidate.is_file(): + _raster_src = str(_candidate) + break + # 策略 B: 扫描 .step11_cache 同级的常见子目录 + if not _raster_src: + for _parent in (_s.parent, _s.parent.parent): + if _parent is None: + continue + for _sub in ('1_water_mask', '1_Water_Mask', 'water_mask', ''): + _scan_dir = _parent / _sub if _sub else _parent + try: + if _scan_dir.is_dir(): + for _ext in _RASTER_EXTS: + _hits = sorted( + _scan_dir.glob(f'*{_ext}'), + key=lambda p: p.stat().st_mtime, reverse=True, + ) + for _h in _hits: + if _h.is_file() and 'vectorized' not in _h.stem.lower(): + _raster_src = str(_h) + break + if _raster_src: + break + except (OSError, PermissionError): + continue + if _raster_src: + break + if _raster_src: + break + # 策略 C: 扫描当前工作目录 + if not _raster_src: + try: + for _ext in _RASTER_EXTS: + _hits = sorted(Path.cwd().glob(f'**/*{_ext}'), + key=lambda p: p.stat().st_mtime, reverse=True) + for _h in _hits: + if _h.is_file() and 'vectorized' not in _h.stem.lower(): + _raster_src = str(_h) + break + if _raster_src: + break + except Exception: + pass + if _raster_src: + print(f" [共享上下文] 自动发现原始栅格: {Path(_raster_src).name}") + + if _raster_src and os.path.isfile(_raster_src): + print(f" [共享上下文] GDAL 快速重采样栅格掩膜 → {nx}×{ny} 网格 ...") + try: + from osgeo import gdal + gdal.UseExceptions() + src_ds = gdal.Open(_raster_src) + if src_ds is not None: + # 目标 SRS：使用 self.output_crs，与网格坐标一致 + from pyproj import CRS as _CRS + _dst_crs = _CRS.from_string(self.output_crs) + _dst_wkt = _dst_crs.to_wkt() + + mem_drv = gdal.GetDriverByName('MEM') + dst_ds = mem_drv.Create('', nx, ny, 1, gdal.GDT_Byte) + dst_ds.SetGeoTransform((min_x, dx, 0, max_y, 0, -dy)) + dst_ds.SetProjection(_dst_wkt) + gdal.ReprojectImage( + src_ds, dst_ds, + src_ds.GetProjection(), _dst_wkt, + gdal.GRA_NearestNeighbour, + ) + mask_raster = dst_ds.ReadAsArray() + src_ds = None + dst_ds = None + + # Y 轴对齐 + if grid_yy[0, 0] < grid_yy[-1, 0]: + mask_raster = np.flipud(mask_raster) + + mask = mask_raster.astype(bool) + print(f" [共享上下文] GDAL 重采样完成: " + f"{int(mask.sum())}/{mask.size} 点在水域内") + return (grid_xx, grid_yy, mask, bounds, boundary_gdf) + except Exception as e: + print(f" [共享上下文] GDAL 重采样失败 ({e})，回退 rasterize ...") + + # ── 常规通道：矢量多边形 rasterize ── + import rasterio.features + from rasterio.transform import from_bounds + + n_polys = len(boundary_gdf) + if n_polys > 1000: + print(f" [共享上下文] {n_polys} 个多边形，先 dissolve 合并再光栅化 ...") + boundary_gdf = boundary_gdf.dissolve() + print(f" [共享上下文] dissolve 完成 → {len(boundary_gdf)} 个要素") + + print(f" [共享上下文] Rasterio 光栅化水体边界 → {nx}×{ny} 网格 ...") 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), @@ -3042,7 +3159,6 @@ class ContentMapper: dtype='uint8', ) - # 坐标系 Y 轴方向对齐处理 if grid_yy[0, 0] < grid_yy[-1, 0]: mask_raster = np.flipud(mask_raster) @@ -3062,7 +3178,7 @@ class ContentMapper: """ 主处理函数 - ★★★ Plan C：shp_file 现在是可选参数（None = 纯采样点插值，不依赖水域掩膜）★★★ + [*][*][*] Plan C：shp_file 现在是可选参数（None = 纯采样点插值，不依赖水域掩膜）[*][*][*] Parameters: ----------- @@ -3084,7 +3200,7 @@ class ContentMapper: # 读取采样点数据 points_gdf = self.read_csv_data(csv_file) - # ── ★ 快速通道：复用预计算的共享上下文 ── + # ── [*] 快速通道：复用预计算的共享上下文 ── if shared_context is not None: grid_xx, grid_yy, mask, bounds, boundary_gdf = shared_context # ③ 仍需边缘扩展（值相关），但跳过 ②④⑥ @@ -3169,7 +3285,7 @@ class ContentMapper: """ 批量处理文件夹中的CSV文件 - ★★★ Plan C：shp_file 可选（None = 不依赖水域掩膜，纯采样点插值）★★★ + [*][*][*] Plan C：shp_file 可选（None = 不依赖水域掩膜，纯采样点插值）[*][*][*] Parameters: ----------- @@ -3251,7 +3367,7 @@ class ContentMapper: ) success_count += 1 - print(f"✓ 成功处理: {csv_basename}.png") + print(f"[OK] 成功处理: {csv_basename}.png") except Exception as e: fail_count += 1