fix: 补齐缺失的 handler 文件 + IDW 插值退化检测 + V1 代码归档
This commit is contained in:
duxin
99
src/core/handlers/step10_watercolor_handler.py
Normal file
99
src/core/handlers/step10_watercolor_handler.py
Normal file
@ -0,0 +1,99 @@
|
||||
#!/usr/bin/env python
|
||||
# -*- coding: utf-8 -*-
|
||||
"""
|
||||
Step10 处理器：水色指数反演（CSV 散点模式）
|
||||
|
||||
2026-06-30 新建：替代原先错误路由的 step10_qaa（QAA 物理反演）。
|
||||
使用 WaterIndexCsvProcessor 对采样点 CSV 逐行计算选定的水质指数公式，
|
||||
每个公式输出一个独立 CSV 文件（longitude, latitude, <公式值>）。
|
||||
"""
|
||||
|
||||
import os
|
||||
import time
|
||||
from typing import Any, Dict, Optional
|
||||
|
||||
from src.core.handlers.base import BaseStepHandler, PipelineContext
|
||||
|
||||
|
||||
class Step10WatercolorHandler(BaseStepHandler):
|
||||
"""步骤10：水色指数反演（CSV 散点模式）。
|
||||
|
||||
对应 config key: 'step10_watercolor'
|
||||
委托类: WaterIndexCsvProcessor
|
||||
"""
|
||||
|
||||
step_key = 'step10_watercolor'
|
||||
|
||||
def execute(self, context: PipelineContext, config: dict) -> Dict[str, Any]:
|
||||
from src.core.algorithms.waterindex_inversion.csv_processor import (
|
||||
WaterIndexCsvProcessor,
|
||||
)
|
||||
|
||||
step_start_time = time.time()
|
||||
|
||||
# ── 解析参数 ──
|
||||
sampling_csv_path = config.get('sampling_csv_path', '')
|
||||
if not sampling_csv_path or not os.path.isfile(sampling_csv_path):
|
||||
msg = f'采样点 CSV 不存在或路径为空: {sampling_csv_path}'
|
||||
context.notify('step10_watercolor', 'error', msg)
|
||||
return {'error': msg}
|
||||
|
||||
output_dir = config.get('output_dir', '')
|
||||
if not output_dir:
|
||||
output_dir = os.path.join(str(context.work_dir), '10_WaterIndex_CSV')
|
||||
|
||||
# 显式兜底：只有非空 list 才启用过滤，杜绝 falsy 空列表透传
|
||||
raw_selected = config.get('selected_formulas')
|
||||
selected_formulas: Optional[list] = (
|
||||
raw_selected if isinstance(raw_selected, list) and len(raw_selected) > 0 else None
|
||||
)
|
||||
wavelength_offset = float(config.get('wavelength_offset', 0))
|
||||
waterindex_csv = config.get('waterindex_csv', '')
|
||||
|
||||
if not waterindex_csv or not os.path.isfile(waterindex_csv):
|
||||
msg = f'waterindex.csv 不存在（前端未传入或路径无效）: {waterindex_csv!r}'
|
||||
context.notify('step10_watercolor', 'error', msg)
|
||||
return {'error': msg}
|
||||
|
||||
# ── 执行 ──
|
||||
try:
|
||||
processor = WaterIndexCsvProcessor(waterindex_csv)
|
||||
|
||||
from src.gui.core.event_bus import global_event_bus
|
||||
|
||||
def _progress_callback(msg: str, pct: float):
|
||||
context.notify('step10_watercolor', 'info', f'[{pct:.0f}%] {msg}')
|
||||
global_event_bus.publish('ProgressUpdate', {
|
||||
'percentage': int(pct),
|
||||
'message': f'Step10 水色指数: {msg}',
|
||||
})
|
||||
|
||||
out_files = processor.compute_indices_from_csv(
|
||||
sampling_csv_path=sampling_csv_path,
|
||||
output_dir=output_dir,
|
||||
selected_formulas=selected_formulas,
|
||||
progress_callback=_progress_callback,
|
||||
wavelength_offset=wavelength_offset,
|
||||
)
|
||||
|
||||
step_end_time = time.time()
|
||||
context.record_step_time(
|
||||
'步骤10: 水色指数反演', step_start_time, step_end_time
|
||||
)
|
||||
context.notify(
|
||||
'step10_watercolor', 'completed',
|
||||
f'水色指数反演完成，共生成 {len(out_files)} 个 CSV'
|
||||
)
|
||||
|
||||
return {
|
||||
'output_dir': output_dir,
|
||||
'output_files': out_files,
|
||||
}
|
||||
|
||||
except Exception as e:
|
||||
step_end_time = time.time()
|
||||
context.record_step_time(
|
||||
'步骤10: 水色指数反演', step_start_time, step_end_time,
|
||||
status='failed', error=str(e)
|
||||
)
|
||||
raise
|
||||
268
src/core/handlers/step11_map_handler.py
Normal file
268
src/core/handlers/step11_map_handler.py
Normal file
@ -0,0 +1,268 @@
|
||||
#!/usr/bin/env python
|
||||
# -*- coding: utf-8 -*-
|
||||
"""
|
||||
Step11 处理器：专题图生成（CSV 插值 → 分布图 PNG）
|
||||
|
||||
2026-06-30：统一到 Pipeline 体系。
|
||||
- 预矢量化：将 .dat 栅格掩膜转一次为 .shp，避免 63 个 CSV 各自矢量化 1765 万像元。
|
||||
- 串行处理：GDAL 不具有线程安全性，多线程并发 Open/Polygonize 导致 0xC0000005。
|
||||
预矢量化后每个 CSV 处理时间大幅缩短，串行速度已可接受。
|
||||
"""
|
||||
|
||||
import os
|
||||
import shutil
|
||||
import time
|
||||
from pathlib import Path
|
||||
from typing import Any, Dict, List, Tuple
|
||||
|
||||
from src.core.handlers.base import BaseStepHandler, PipelineContext
|
||||
|
||||
|
||||
def _process_one_map(csv_path: str, base_kwargs: dict, output_dir: str
|
||||
) -> Tuple[str, str]:
|
||||
"""处理单个 CSV → 分布图 TIF（子进程内执行）。
|
||||
|
||||
2026-07-01：支持 shared_context 快速通道。
|
||||
当 base_kwargs 包含 'shared_context' 时，绕过 MappingStep，
|
||||
直接用 ContentMapper.process_data(shared_context=ctx) 跳过重复的边界/网格/掩膜计算。
|
||||
"""
|
||||
import matplotlib
|
||||
try:
|
||||
matplotlib.use('Agg', force=True)
|
||||
except Exception:
|
||||
pass
|
||||
|
||||
stem = Path(csv_path).stem
|
||||
output_image_path = str(Path(output_dir) / f'{stem}_distribution.png')
|
||||
|
||||
shared_context = base_kwargs.get('shared_context')
|
||||
|
||||
if shared_context is not None:
|
||||
# ★ 快速通道：直接调 ContentMapper，复用预计算的网格/掩膜
|
||||
from src.postprocessing.map import ContentMapper
|
||||
mapper = ContentMapper(
|
||||
input_crs=base_kwargs.get('input_crs', 'EPSG:32651'),
|
||||
output_crs=base_kwargs.get('output_crs', base_kwargs.get('input_crs', 'EPSG:32651')),
|
||||
)
|
||||
result_path = mapper.process_data(
|
||||
csv_file=csv_path,
|
||||
shp_file=base_kwargs.get('boundary_shp_path'),
|
||||
output_file=output_image_path,
|
||||
resolution=float(base_kwargs.get('resolution', 10.0)),
|
||||
output_format='tif',
|
||||
shared_context=shared_context,
|
||||
)
|
||||
return result_path, csv_path
|
||||
|
||||
# 原有路径（无共享上下文时）
|
||||
kw = dict(base_kwargs)
|
||||
kw.pop('shared_context', None) # 防止透传给 MappingStep
|
||||
kw['prediction_csv_path'] = csv_path
|
||||
kw['output_image_path'] = output_image_path
|
||||
|
||||
from src.core.steps.mapping_step import MappingStep
|
||||
result_path = MappingStep.generate_distribution_map(**kw)
|
||||
return result_path, csv_path
|
||||
|
||||
|
||||
def _pre_vectorize_boundary(boundary_path: str, work_dir: str) -> str:
|
||||
"""将栅格掩膜预转换为矢量 SHP（只做一次），返回 SHP 路径。"""
|
||||
src = Path(boundary_path)
|
||||
if not src.exists():
|
||||
return boundary_path
|
||||
if src.suffix.lower() in ('.shp', '.geojson'):
|
||||
return boundary_path
|
||||
|
||||
cache_dir = Path(work_dir) / '.step11_cache'
|
||||
cache_dir.mkdir(parents=True, exist_ok=True)
|
||||
shp_out = cache_dir / f'{src.stem}_vectorized.shp'
|
||||
if shp_out.exists():
|
||||
return str(shp_out)
|
||||
|
||||
from osgeo import gdal, ogr, osr
|
||||
gdal.UseExceptions()
|
||||
|
||||
ds = gdal.Open(str(src))
|
||||
if ds is None:
|
||||
return boundary_path
|
||||
|
||||
band = ds.GetRasterBand(1)
|
||||
drv = ogr.GetDriverByName('ESRI Shapefile')
|
||||
if drv is None:
|
||||
ds = None
|
||||
return boundary_path
|
||||
|
||||
for ext in ('.shp', '.shx', '.dbf', '.prj'):
|
||||
p = cache_dir / f'{src.stem}_vectorized{ext}'
|
||||
if p.exists():
|
||||
p.unlink()
|
||||
|
||||
out_ds = drv.CreateDataSource(str(shp_out))
|
||||
srs = osr.SpatialReference()
|
||||
proj = ds.GetProjection()
|
||||
if proj:
|
||||
srs.ImportFromWkt(proj)
|
||||
else:
|
||||
srs.ImportFromEPSG(32651)
|
||||
|
||||
layer = out_ds.CreateLayer('boundary', srs, ogr.wkbMultiPolygon)
|
||||
gdal.Polygonize(band, band, layer, 0, [], callback=None)
|
||||
|
||||
out_ds = None
|
||||
ds = None
|
||||
|
||||
if shp_out.exists():
|
||||
return str(shp_out)
|
||||
return boundary_path
|
||||
|
||||
|
||||
class Step11MapHandler(BaseStepHandler):
|
||||
"""步骤11：专题图生成（预矢量化 + 串行）。
|
||||
|
||||
对应 config key: 'step11_map'
|
||||
"""
|
||||
|
||||
step_key = 'step11_map'
|
||||
|
||||
def execute(self, context: PipelineContext, config: dict) -> Dict[str, Any]:
|
||||
step_start_time = time.time()
|
||||
|
||||
prediction_csv_path = config.get('prediction_csv_path', '')
|
||||
prediction_csv_dir = config.get('prediction_csv_dir', '')
|
||||
batch_mode = config.get('batch_mode', False)
|
||||
boundary_shp_path = config.get('boundary_shp_path', '') or ''
|
||||
output_dir = config.get('output_dir', '') or str(context.visualization_dir)
|
||||
output_image_path = config.get('output_image_path', '')
|
||||
|
||||
# ── 收集 CSV ──
|
||||
csv_paths: List[str] = []
|
||||
if prediction_csv_path and os.path.isfile(prediction_csv_path):
|
||||
csv_paths = [prediction_csv_path]
|
||||
elif batch_mode and prediction_csv_dir and os.path.isdir(prediction_csv_dir):
|
||||
import glob
|
||||
csv_paths = sorted(glob.glob(os.path.join(prediction_csv_dir, '*.csv')))
|
||||
|
||||
if not csv_paths:
|
||||
context.notify('step11_map', 'error', '没有找到预测 CSV 文件')
|
||||
return {'error': '没有找到预测 CSV 文件'}
|
||||
|
||||
# ── ★ 预矢量化边界（仅一次）──
|
||||
resolved_boundary = boundary_shp_path
|
||||
if boundary_shp_path and os.path.isfile(boundary_shp_path):
|
||||
try:
|
||||
context.notify('step11_map', 'info', '正在预矢量化水域边界（仅一次）…')
|
||||
resolved_boundary = _pre_vectorize_boundary(
|
||||
boundary_shp_path, str(context.work_dir)
|
||||
)
|
||||
context.notify('step11_map', 'info',
|
||||
f'边界预处理完成，开始生成专题图…')
|
||||
except Exception as e:
|
||||
context.notify('step11_map', 'warning',
|
||||
f'边界预处理失败，使用原始路径: {e}')
|
||||
|
||||
base_kwargs = {
|
||||
'boundary_shp_path': resolved_boundary or None,
|
||||
'resolution': float(config.get('resolution', 10.0)),
|
||||
'input_crs': config.get('input_crs', 'EPSG:32651'),
|
||||
'output_crs': config.get('output_crs', config.get('input_crs', 'EPSG:32651')),
|
||||
'show_sample_points': bool(config.get('show_sample_points', False)),
|
||||
'use_distance_diffusion': bool(config.get('use_distance_diffusion', False)),
|
||||
}
|
||||
if output_image_path and len(csv_paths) == 1:
|
||||
base_kwargs['output_image_path'] = output_image_path
|
||||
|
||||
os.makedirs(output_dir, exist_ok=True)
|
||||
|
||||
total = len(csv_paths)
|
||||
from src.gui.core.event_bus import global_event_bus
|
||||
|
||||
context.notify('step11_map', 'info', f'专题图开始（{total} 个 CSV）')
|
||||
global_event_bus.publish('ProgressUpdate', {
|
||||
'percentage': 0, 'message': f'Step11 专题图: 0/{total}',
|
||||
})
|
||||
|
||||
# ── ★ 2026-07-01：主进程预计算共享空间上下文 ──
|
||||
# 63 个 CSV 坐标一致，边界/网格/掩膜只算一次，塞入 base_kwargs
|
||||
# 子进程通过 shared_context 复用，跳过重复的 ②④⑥ (~1.2s/CSV)
|
||||
if total > 1 and resolved_boundary:
|
||||
try:
|
||||
context.notify('step11_map', 'info',
|
||||
'正在预计算共享空间上下文（边界+网格+掩膜，仅一次）…')
|
||||
from src.postprocessing.map import ContentMapper
|
||||
pre_mapper = ContentMapper(
|
||||
input_crs=base_kwargs['input_crs'],
|
||||
output_crs=base_kwargs['output_crs'],
|
||||
)
|
||||
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,
|
||||
)
|
||||
base_kwargs['shared_context'] = shared_ctx
|
||||
context.notify('step11_map', 'info',
|
||||
f'✓ 共享上下文已就绪，后续 {total} 个 CSV 将直接复用')
|
||||
except Exception as e:
|
||||
context.notify('step11_map', 'warning',
|
||||
f'共享上下文预计算失败: {e}，回退逐个处理')
|
||||
|
||||
# ── 多进程并行（GDAL 线程不安全，但进程隔离下安全）──
|
||||
generated: List[str] = []
|
||||
errors: Dict[str, str] = {}
|
||||
|
||||
import multiprocessing
|
||||
from concurrent.futures import ProcessPoolExecutor, as_completed
|
||||
|
||||
# 留出 1-2 个核心保证电脑不卡死
|
||||
max_workers = max(1, multiprocessing.cpu_count() - 2)
|
||||
context.notify('step11_map', 'info', f'使用 {max_workers} 个进程并行生成')
|
||||
|
||||
with ProcessPoolExecutor(max_workers=max_workers) as executor:
|
||||
future_to_csv = {
|
||||
executor.submit(_process_one_map, csv_p, base_kwargs, output_dir): csv_p
|
||||
for csv_p in csv_paths
|
||||
}
|
||||
done_count = 0
|
||||
for future in as_completed(future_to_csv):
|
||||
csv_p = future_to_csv[future]
|
||||
done_count += 1
|
||||
try:
|
||||
result_path, _ = future.result()
|
||||
generated.append(result_path)
|
||||
except Exception as e:
|
||||
errors[csv_p] = str(e)
|
||||
context.notify('step11_map', 'warning',
|
||||
f'专题图 FAIL: {Path(csv_p).name} — {e}')
|
||||
|
||||
pct = int(done_count / total * 100)
|
||||
global_event_bus.publish('ProgressUpdate', {
|
||||
'percentage': pct,
|
||||
'message': f'Step11 专题图: {done_count}/{total}',
|
||||
})
|
||||
if done_count % max(1, total // 10) == 0 or done_count == total:
|
||||
context.notify('step11_map', 'info',
|
||||
f'专题图 [{done_count}/{total}]')
|
||||
|
||||
step_end_time = time.time()
|
||||
elapsed = step_end_time - step_start_time
|
||||
context.record_step_time('步骤11: 专题图生成', step_start_time, step_end_time)
|
||||
|
||||
result_msg = f'专题图完成：{len(generated)}/{total} 张（{elapsed:.1f}s）'
|
||||
if errors:
|
||||
result_msg += f'，{len(errors)} 个失败'
|
||||
context.notify('step11_map', 'completed', result_msg)
|
||||
|
||||
# 清理临时缓存
|
||||
try:
|
||||
cache_dir = Path(str(context.work_dir)) / '.step11_cache'
|
||||
if cache_dir.exists():
|
||||
shutil.rmtree(cache_dir, ignore_errors=True)
|
||||
except Exception:
|
||||
pass
|
||||
|
||||
return {
|
||||
'generated': generated,
|
||||
'errors': errors,
|
||||
'output_dir': output_dir,
|
||||
'elapsed_seconds': elapsed,
|
||||
}
|
||||
4
src/gui/_legacy/__init__.py
Normal file
4
src/gui/_legacy/__init__.py
Normal file
@ -0,0 +1,4 @@
|
||||
# _legacy — 旧版 V1 代码归档
|
||||
# 2026-07-01：water_quality_gui.py 是原始单体架构（~2000 行上帝类），
|
||||
# 已被 water_quality_gui_v2.py（纯壳 + Manager 模式）取代。
|
||||
# 保留此文件仅供历史参考，不应被任何活跃代码引用。
|
||||
1598
src/gui/_legacy/water_quality_gui.py
Normal file
1598
src/gui/_legacy/water_quality_gui.py
Normal file
File diff suppressed because it is too large
Load Diff
@ -556,7 +556,14 @@ class ContentMapper:
|
||||
return grid_content
|
||||
|
||||
def _perform_interpolation(self, points, values, grid_xx, grid_yy):
|
||||
"""执行空间插值"""
|
||||
"""三级降级插值策略：Kriging → IDW → 最近邻
|
||||
|
||||
2026-07-01 重构：
|
||||
- Kriging 优先（自动拟合球形变异函数，不强制 nugget）
|
||||
- Kriging 退化检测：若结果标准差接近 0（纯色图），自动回退 IDW
|
||||
- IDW 作为首选回退：无需拟合变异函数，不会产生纯色图
|
||||
- scipy linear/nearest 作为最后兜底
|
||||
"""
|
||||
print(f"插值输入检查：")
|
||||
print(f" - 数据点数量: {len(points)}")
|
||||
print(f" - 数据值范围: {values.min():.4f} - {values.max():.4f}")
|
||||
@ -573,149 +580,145 @@ class ContentMapper:
|
||||
if len(points) < 3:
|
||||
raise ValueError(f"有效数据点不足3个（当前：{len(points)}个）")
|
||||
|
||||
# 优先使用Kriging插值
|
||||
# ── 策略 0：值域极窄时直接跳过 Kriging ──
|
||||
value_range = float(values.max() - values.min())
|
||||
value_std = float(np.std(values))
|
||||
|
||||
# ═══════════════════════════════════════════════════════════
|
||||
# 策略 1：Kriging（自动拟合球形变异函数）
|
||||
# ═══════════════════════════════════════════════════════════
|
||||
kriging_degraded = False
|
||||
if PYKRIGE_AVAILABLE:
|
||||
try:
|
||||
print("正在使用Kriging插值（半变异函数模型，块金值=100%）...")
|
||||
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
|
||||
enable_plotting=False,
|
||||
)
|
||||
# ★ 局部邻域 Kriging：只参考最近的 15 个点，避免万阶矩阵求逆
|
||||
z, _ = ok.execute('grid', grid_x, grid_y, backend='loop', n_closest_points=15)
|
||||
z, ss = ok.execute('grid', grid_x, grid_y, backend='loop', n_closest_points=15)
|
||||
grid_content = np.array(z)
|
||||
valid_count = np.sum(~np.isnan(grid_content))
|
||||
print(f"Kriging插值成功，有效点数: {valid_count} / {grid_content.size}")
|
||||
if valid_count > 0:
|
||||
return grid_content
|
||||
else:
|
||||
print("警告：Kriging插值结果为空，将回退到其他插值方法")
|
||||
except Exception as e:
|
||||
print(f"Kriging插值失败: {e}，将回退到其他插值方法")
|
||||
else:
|
||||
print("警告：pykrige未安装，无法使用Kriging插值，将使用其他插值方法")
|
||||
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}, "
|
||||
f"输出std={kriging_std:.6f}, 退化比={degradation_ratio:.3f}")
|
||||
|
||||
if degradation_ratio < 0.05 and value_range > 1e-8:
|
||||
print(f"⚠ Kriging 严重退化（输出 std/输入 std={degradation_ratio:.3f}<5%），"
|
||||
f"判定为纯色图，回退 IDW")
|
||||
kriging_degraded = True
|
||||
else:
|
||||
print(f"Kriging 通过退化检测，直接使用")
|
||||
return grid_content
|
||||
else:
|
||||
print("Kriging 结果全为 NaN，回退")
|
||||
kriging_degraded = True
|
||||
except Exception as e:
|
||||
print(f"Kriging 失败: {e}")
|
||||
kriging_degraded = True
|
||||
else:
|
||||
print("pykrige 未安装，跳过 Kriging")
|
||||
kriging_degraded = True
|
||||
|
||||
# ═══════════════════════════════════════════════════════════
|
||||
# 策略 2：IDW（反距离权重）— 不需要拟合变异函数，绝不纯色
|
||||
# ═══════════════════════════════════════════════════════════
|
||||
if kriging_degraded:
|
||||
try:
|
||||
print("正在使用 IDW 插值（反距离权重, power=2, neighbors=15）...")
|
||||
grid_content = self._idw_interpolation(
|
||||
points, values, grid_xx, grid_yy,
|
||||
power=2, n_neighbors=min(15, len(points)),
|
||||
)
|
||||
valid_count = int(np.sum(~np.isnan(grid_content)))
|
||||
if valid_count > 0:
|
||||
idw_std = float(np.nanstd(grid_content))
|
||||
print(f"IDW 完成: 有效点={valid_count}/{grid_content.size}, 输出std={idw_std:.6f}")
|
||||
if idw_std > 0:
|
||||
return grid_content
|
||||
else:
|
||||
print("IDW std=0（所有输入值完全相同），结果可用")
|
||||
return grid_content
|
||||
else:
|
||||
print("IDW 结果全为 NaN，回退 scipy 插值")
|
||||
except Exception as e:
|
||||
print(f"IDW 失败: {e}，回退 scipy 插值")
|
||||
|
||||
# ═══════════════════════════════════════════════════════════
|
||||
# 策略 3：scipy 线性插值 + 最近邻填充（最终兜底）
|
||||
# ═══════════════════════════════════════════════════════════
|
||||
try:
|
||||
# 首先尝试使用线性插值
|
||||
print("正在尝试线性插值...")
|
||||
print("正在尝试 scipy 线性插值...")
|
||||
grid_content = griddata(
|
||||
points, values, (grid_xx, grid_yy),
|
||||
method='linear', fill_value=np.nan
|
||||
)
|
||||
|
||||
# 检查线性插值结果
|
||||
valid_linear = ~np.isnan(grid_content)
|
||||
valid_count = np.sum(valid_linear)
|
||||
print(f"线性插值结果：有效点数 {valid_count} / {grid_content.size}")
|
||||
valid_count = int(np.sum(~np.isnan(grid_content)))
|
||||
print(f"线性插值: 有效点={valid_count}/{grid_content.size}")
|
||||
|
||||
if valid_count > 0:
|
||||
print(f"线性插值成功，有效区域覆盖率: {valid_count / grid_content.size * 100:.1f}%")
|
||||
|
||||
# 如果有NaN值，用最近邻插值填充
|
||||
nan_count = np.sum(np.isnan(grid_content))
|
||||
nan_count = int(np.sum(np.isnan(grid_content)))
|
||||
if nan_count > 0:
|
||||
print(f"正在用最近邻插值填充 {nan_count} 个缺失值...")
|
||||
print(f"用最近邻填充 {nan_count} 个 NaN...")
|
||||
grid_nearest = griddata(
|
||||
points, values, (grid_xx, grid_yy),
|
||||
method='nearest'
|
||||
points, values, (grid_xx, grid_yy), method='nearest'
|
||||
)
|
||||
# 只填充线性插值的NaN区域
|
||||
nan_mask = np.isnan(grid_content)
|
||||
grid_content[nan_mask] = grid_nearest[nan_mask]
|
||||
print("缺失值填充完成")
|
||||
|
||||
# 最终检查
|
||||
final_valid = ~np.isnan(grid_content)
|
||||
print(f"最终有效点数: {np.sum(final_valid)} / {grid_content.size}")
|
||||
|
||||
grid_content[np.isnan(grid_content)] = grid_nearest[np.isnan(grid_content)]
|
||||
return grid_content
|
||||
else:
|
||||
print("线性插值失败，尝试最近邻插值...")
|
||||
|
||||
except Exception as e:
|
||||
print(f"线性插值失败: {e}")
|
||||
print("尝试最近邻插值...")
|
||||
|
||||
try:
|
||||
# 使用最近邻插值作为备选方案
|
||||
print("执行最近邻插值...")
|
||||
grid_content = griddata(
|
||||
points, values, (grid_xx, grid_yy),
|
||||
method='nearest'
|
||||
)
|
||||
# ═══════════════════════════════════════════════════════════
|
||||
# 策略 4：最近邻（绝对兜底）
|
||||
# ═══════════════════════════════════════════════════════════
|
||||
print("执行最近邻插值（最终兜底）...")
|
||||
grid_content = griddata(
|
||||
points, values, (grid_xx, grid_yy), method='nearest'
|
||||
)
|
||||
if np.sum(~np.isnan(grid_content)) == 0:
|
||||
raise ValueError("所有插值方法均失败")
|
||||
return grid_content
|
||||
|
||||
valid_count = np.sum(~np.isnan(grid_content))
|
||||
print(f"最近邻插值成功，有效点数: {valid_count}")
|
||||
@staticmethod
|
||||
def _idw_interpolation(points, values, grid_xx, grid_yy,
|
||||
power=2, n_neighbors=15):
|
||||
"""IDW（反距离权重）插值 — 不需拟合模型，绝不产生纯色图。
|
||||
|
||||
if valid_count == 0:
|
||||
raise ValueError("最近邻插值也没有产生有效结果")
|
||||
Parameters:
|
||||
points: (N, 2) 采样点坐标
|
||||
values: (N,) 采样点值
|
||||
grid_xx, grid_yy: meshgrid 网格
|
||||
power: 距离衰减幂参数（默认 2）
|
||||
n_neighbors: 每个网格点参考的最近邻数量
|
||||
"""
|
||||
from scipy.spatial import cKDTree
|
||||
grid_shape = grid_xx.shape
|
||||
grid_flat = np.column_stack((grid_xx.ravel(), grid_yy.ravel()))
|
||||
values_flat = values.ravel()
|
||||
|
||||
return grid_content
|
||||
|
||||
except Exception as e:
|
||||
print(f"最近邻插值也失败: {e}")
|
||||
|
||||
# 对于地理坐标系，尝试更简单的方法
|
||||
if self.output_crs == 'EPSG:4326':
|
||||
print("地理坐标系检测到，尝试简化插值...")
|
||||
try:
|
||||
# 创建一个基于距离的简单插值
|
||||
grid_content = np.full(grid_xx.shape, np.nan)
|
||||
|
||||
# 为每个网格点找到最近的数据点
|
||||
for i in range(grid_xx.shape[0]):
|
||||
for j in range(grid_xx.shape[1]):
|
||||
grid_x, grid_y = grid_xx[i, j], grid_yy[i, j]
|
||||
|
||||
# 计算到所有数据点的距离
|
||||
distances = np.sqrt((points[:, 0] - grid_x) ** 2 + (points[:, 1] - grid_y) ** 2)
|
||||
nearest_idx = np.argmin(distances)
|
||||
|
||||
# 如果距离不是太远，就使用该值
|
||||
if distances[nearest_idx] < (grid_xx.max() - grid_xx.min()) * 0.1: # 10%的范围内
|
||||
grid_content[i, j] = values[nearest_idx]
|
||||
|
||||
valid_count = np.sum(~np.isnan(grid_content))
|
||||
print(f"简化插值完成，有效点数: {valid_count}")
|
||||
|
||||
if valid_count > 0:
|
||||
return grid_content
|
||||
else:
|
||||
raise ValueError("简化插值也没有产生有效结果")
|
||||
|
||||
except Exception as e3:
|
||||
print(f"简化插值失败: {e3}")
|
||||
|
||||
print("尝试立方插值作为最后手段...")
|
||||
try:
|
||||
# 最后尝试立方插值
|
||||
grid_content = griddata(
|
||||
points, values, (grid_xx, grid_yy),
|
||||
method='cubic', fill_value=np.nan
|
||||
)
|
||||
|
||||
# 如果立方插值有NaN，用最近邻填充
|
||||
if np.any(np.isnan(grid_content)):
|
||||
print("用最近邻插值填充立方插值的NaN值...")
|
||||
grid_nearest = griddata(
|
||||
points, values, (grid_xx, grid_yy),
|
||||
method='nearest'
|
||||
)
|
||||
nan_mask = np.isnan(grid_content)
|
||||
grid_content[nan_mask] = grid_nearest[nan_mask]
|
||||
|
||||
valid_count = np.sum(~np.isnan(grid_content))
|
||||
print(f"立方插值成功，有效点数: {valid_count}")
|
||||
return grid_content
|
||||
|
||||
except Exception as e4:
|
||||
print(f"立方插值也失败: {e4}")
|
||||
print(f"所有插值方法都失败")
|
||||
raise ValueError("无法完成空间插值，请检查数据点的分布和数值")
|
||||
tree = cKDTree(points)
|
||||
k = min(n_neighbors, len(points))
|
||||
distances, indices = tree.query(grid_flat, k=k)
|
||||
# 防止距离为 0 的除零
|
||||
distances = np.maximum(distances, 1e-12)
|
||||
weights = 1.0 / (distances ** power)
|
||||
# 归一化权重
|
||||
weights /= weights.sum(axis=1, keepdims=True)
|
||||
# 加权求和
|
||||
neighbor_vals = values_flat[indices] if k == 1 else values_flat[indices]
|
||||
if k == 1:
|
||||
result = neighbor_vals
|
||||
else:
|
||||
result = np.sum(weights * neighbor_vals, axis=1)
|
||||
return result.reshape(grid_shape)
|
||||
|
||||
def read_csv_data(self, csv_file, uncertainty_col=None):
|
||||
"""
|
||||
@ -2514,28 +2517,30 @@ class ContentMapper:
|
||||
safe_h = min(float(figsize[1]), _max_inch)
|
||||
fig, ax = plt.subplots(figsize=(safe_w, safe_h))
|
||||
|
||||
# 计算有效值统计（2σ 标准差拉伸，排除长尾异常值干扰）
|
||||
valid = array[~np.isnan(array)]
|
||||
# 1. 明确排除 NaN 以及 nodata_value（与函数参数保持一致）
|
||||
nodata_val = nodata_value
|
||||
valid = array[(~np.isnan(array)) & (array != nodata_val)]
|
||||
|
||||
if valid.size == 0:
|
||||
raise ValueError("GeoTIFF 中没有有效数据（全部为 NoData）")
|
||||
|
||||
mean_val = float(np.nanmean(array))
|
||||
std_val = float(np.nanstd(array))
|
||||
vmin = max(float(np.nanmin(array)), mean_val - 2 * std_val)
|
||||
vmax = min(float(np.nanmax(array)), mean_val + 2 * std_val)
|
||||
# 2. 改用更鲁棒的 2%-98% 百分位拉伸（抗偏态分布）
|
||||
vmin = float(np.percentile(valid, 2))
|
||||
vmax = float(np.percentile(valid, 98))
|
||||
|
||||
if (vmax - vmin) < 1e-9:
|
||||
center = mean_val
|
||||
if (vmax - vmin) < 1e-9: # 防退化：区间过窄 → 取中心 ±1%
|
||||
center = vmin
|
||||
exp = max(abs(center) * 0.01, 1e-9)
|
||||
vmin = center - exp
|
||||
vmax = center + exp
|
||||
|
||||
print(f"[visualize_raster] 2σ 拉伸: vmin={vmin:.4f}, vmax={vmax:.4f}，"
|
||||
f"mean={mean_val:.4f}, std={std_val:.4f}，有效像元: {valid.size}/{array.size}")
|
||||
print(f"[visualize_raster] P2-P98 拉伸: vmin={vmin:.4f}, vmax={vmax:.4f}，"
|
||||
f"有效像元: {valid.size}/{array.size}")
|
||||
|
||||
# ── 栅格绘图 ─────────────────────────────────────────────────
|
||||
# 使用 masked array：NaN 区域自动不显示
|
||||
masked_data = np.ma.masked_invalid(array)
|
||||
# 不仅要屏蔽 NaN，如果 array 中还有 nodata_value，也必须 mask 掉，否则出图时背景会被渲染
|
||||
masked_data = np.ma.masked_where((np.isnan(array)) | (array == nodata_val), array)
|
||||
|
||||
# 【核心修复2】废弃错误的坐标映射逻辑。
|
||||
# 直接使用原生的 imshow，明确告知 matplotlib 第0行在最上方(origin='upper')
|
||||
|
||||
Reference in New Issue
Block a user