duxingchen/WQ_GUI
2
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

perf(step11): IDW auto-switch for large grids + GDAL fast rasterize

Browse Source 
- 网格 >500K 点自动跳过 Kriging 切 IDW（千万级点从数小时降至数秒）
- n_closest_points 50→20，Block 每块打印进度
- GDAL ReprojectImage 直接重采样栅格掩膜，避免 58K 多边形 rasterize 卡死
- 自动扫描工作目录 .dat 栅格，不依赖 handler 传参
- 特殊 Unicode 字符替换为 ASCII 兼容
This commit is contained in:
duxin
2026-07-24 14:36:51 +08:00
parent d0344b8fa0
commit beed27b1f0
2 changed files with 202 additions and 55 deletions
Download Patch File Download Diff File Expand all files Collapse all files

31
src/core/handlers/step11_map_handler.py
View File

@ -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',

226
src/postprocessing/map.py
View File

@ -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