diff --git a/src/postprocessing/map.py b/src/postprocessing/map.py index b1ce13f..722530e 100644 --- a/src/postprocessing/map.py +++ b/src/postprocessing/map.py @@ -586,26 +586,54 @@ class ContentMapper: # ═══════════════════════════════════════════════════════════ # 策略 1：Kriging（自动拟合球形变异函数） + # 大网格 (>1M) 自动启用多进程分块 + C 后端加速 # ═══════════════════════════════════════════════════════════ kriging_degraded = False if PYKRIGE_AVAILABLE: try: - print("正在使用 Kriging 插值（球形模型，自动拟合 nugget）...") grid_x = grid_xx[0, :] grid_y = grid_yy[:, 0] - ok = OrdinaryKriging( - points[:, 0], points[:, 1], values, - variogram_model='spherical', - verbose=False, - enable_plotting=False, - ) - z, ss = ok.execute('grid', grid_x, grid_y, backend='loop', n_closest_points=15) - grid_content = np.array(z) + total_cells = len(grid_x) * len(grid_y) + + # 自动检测 pykrige C 后端（比 loop 后端快 50-100×） + _krige_backend = self._detect_kriging_backend() + + # 多进程分块阈值：>1M 网格点时拆分为 N 个独立子网格并行计算 + _CHUNK_THRESHOLD = 1_000_000 + _n_workers = min(os.cpu_count() or 4, 8) + + print(f"正在使用 Kriging 插值（球形模型，backend={_krige_backend}，" + f"网格={total_cells:,} 点，n_closest=15）...") + + if total_cells > _CHUNK_THRESHOLD and _n_workers > 1: + # ── 多进程分块模式 ── + print(f" 网格 > {_CHUNK_THRESHOLD:,}，启用多进程分块 " + f"({_n_workers} 个 worker)...") + grid_content = self._krige_chunked( + points, values, grid_x, grid_y, + n_workers=_n_workers, + backend=_krige_backend, + n_closest_points=15, + ) + else: + # ── 单进程模式 ── + ok = OrdinaryKriging( + points[:, 0], points[:, 1], values, + variogram_model='spherical', + verbose=False, + enable_plotting=False, + ) + z, ss = ok.execute( + 'grid', grid_x, grid_y, + backend=_krige_backend, + n_closest_points=15, + ) + grid_content = np.array(z) + valid_mask = ~np.isnan(grid_content) valid_count = int(np.sum(valid_mask)) if valid_count > 0: - # ★ 退化检测：若插值结果标准差 < 原始数据标准差的 5%，判定为纯色图 kriging_std = float(np.nanstd(grid_content)) degradation_ratio = kriging_std / max(value_std, 1e-12) print(f"Kriging 完成: 有效点={valid_count}/{grid_content.size}, " @@ -687,6 +715,87 @@ class ContentMapper: raise ValueError("所有插值方法均失败") return grid_content + @staticmethod + def _detect_kriging_backend(): + """检测 pykrige 可用的最快后端: 'C' > 'vectorized' > 'loop'""" + try: + from pykrige import __version__ + from pykrige.ok import OrdinaryKriging + # 快速探测：如果 import 正常，大概率 C 后端已编译 + # C 后端的标志：pykrige 安装自 conda-forge 或带编译的 pip wheel + import pykrige.ok as _ok_mod + if hasattr(_ok_mod, 'Okrige') or hasattr(_ok_mod, '_ok'): + print("[Kriging] 检测到 C 编译后端，性能最佳") + return 'C' + except Exception: + pass + # vectorized 比 loop 快，但内存开销大（大网格时可能 OOM） + # 保守起见大网格用 loop，小网格用 vectorized + print("[Kriging] C 后端不可用，使用 loop 后端（多进程补偿）") + return 'loop' + + @staticmethod + def _krige_one_chunk(args): + """单个 Kriging 分块任务（独立进程入口）""" + points, values, grid_x_chunk, grid_y_chunk, chunk_idx, total_chunks = args + + import numpy as np + from pykrige.ok import OrdinaryKriging + + print(f" [Krige Worker] 分块 {chunk_idx + 1}/{total_chunks} " + f"({len(grid_x_chunk)}×{len(grid_y_chunk)} = " + f"{len(grid_x_chunk) * len(grid_y_chunk):,} 点)...") + + ok = OrdinaryKriging( + points[:, 0], points[:, 1], values, + variogram_model='spherical', + verbose=False, + enable_plotting=False, + ) + z, ss = ok.execute( + 'grid', grid_x_chunk, grid_y_chunk, + backend='loop', # 子进程用 loop（C 后端可能跨进程不稳定） + n_closest_points=15, + ) + print(f" [Krige Worker] 分块 {chunk_idx + 1}/{total_chunks} 完成") + return np.array(z), chunk_idx + + def _krige_chunked(self, points, values, grid_x, grid_y, + n_workers=4, backend='loop', n_closest_points=15): + """多进程分块 Kriging：将 Y 轴方向拆分为 N 个 stripe 并行计算 + + 例如 2285×4600 网格 → 4 个 worker， + 每个处理 2285×1150 子网格 → 4× 加速。 + """ + import multiprocessing + + n_y = len(grid_y) + chunk_size = (n_y + n_workers - 1) // n_workers + tasks = [] + for i in range(n_workers): + y_start = i * chunk_size + y_end = min((i + 1) * chunk_size, n_y) + if y_start >= n_y: + break + grid_y_chunk = grid_y[y_start:y_end] + tasks.append(( + points, values, grid_x, grid_y_chunk, i, n_workers, + )) + + if len(tasks) == 1: + # 单 worker，不需要多进程开销 + return self._krige_one_chunk(tasks[0])[0] + + print(f" 启动 {len(tasks)} 个 Kriging worker 进程...") + with multiprocessing.Pool(processes=len(tasks)) as pool: + results = pool.map(self._krige_one_chunk, tasks) + + # 按 chunk_idx 排序后沿 Y 轴拼接 + results.sort(key=lambda r: r[1]) + grid_content = np.vstack([r[0] for r in results]) + print(f" 多进程 Kriging 完成，拼接结果: {grid_content.shape}") + return grid_content + @staticmethod def _idw_interpolation(points, values, grid_xx, grid_yy, power=2, n_neighbors=15): @@ -1175,31 +1284,6 @@ class ContentMapper: total_grid_cells = grid_xx.size print(f"网格大小: {grid_xx.shape[1]} x {grid_xx.shape[0]} (宽 x 高) = {total_grid_cells:,} 个网格点") - # ★ 大网格自动拦截：超过阈值时自动提升分辨率以保证可接受的计算时间 - _MAX_GRID = 500000 # Kriging 在此规模约需 30-120s - _URGENT_MAX_GRID = 2000000 # 绝对上限，超过直接拒绝 - if total_grid_cells > _URGENT_MAX_GRID: - raise ValueError( - f"网格点数 {total_grid_cells:,} 超过绝对上限 {_URGENT_MAX_GRID:,}。" - f"当前分辨率下 Kriging 预计需要数小时。" - f"请将「空间分布图」面板的「分辨率」增大到至少 " - f"{int(resolution * (total_grid_cells / _MAX_GRID) ** 0.5)}m 后重试。" - ) - elif total_grid_cells > _MAX_GRID: - # 自动降分辨率 - scale = (total_grid_cells / _MAX_GRID) ** 0.5 - new_res = int(resolution * scale) - new_nx = max(100, int(grid_points_x / scale)) - new_ny = max(100, int(grid_points_y / scale)) - print(f"⚠ 网格点数 {total_grid_cells:,} 超过推荐上限 {_MAX_GRID:,}，" - f"自动将分辨率从 {resolution}m 提升到 ~{new_res}m " - f"(网格缩减到 ≈{new_nx}×{new_ny} = {new_nx*new_ny:,} 点)") - grid_x = np.linspace(minx, maxx, new_nx) - grid_y = np.linspace(miny, maxy, new_ny) - grid_xx, grid_yy = np.meshgrid(grid_x, grid_y) - total_grid_cells = grid_xx.size - print(f"调整后网格: {grid_xx.shape[1]} x {grid_xx.shape[0]} = {total_grid_cells:,} 个网格点") - if grid_xx.shape[0] < 2 or grid_xx.shape[1] < 2: raise ValueError(f"网格尺寸太小 {grid_xx.shape}，无法进行插值。")