fix: 专题图 CRS 探测增强 + 地理坐标系 UTM 临时投影
map.py: - CRS 探测从 Proj4 字符串兜底解析 EPSG（+init=epsg:XXXX / +zone=） - rasterio 支持 .dat/.bsq/.bil/.bip/.img 等 ENVI 格式 - 新增 _is_geographic_crs() / _get_utm_epsg() 方法 - 地理坐标系自动临时投影到 UTM 米制空间进行插值，避免经纬度数值 过小导致 Kriging/IDW 矩阵崩溃、std 退化为 0（纯色图） - 防御：CRS 标记为地理系但坐标值超出经纬度范围时判定为探测错误， 跳过投影直接使用米制坐标 water_quality_gui_v2.py: 格式清理
This commit is contained in:
duxin
@ -362,6 +362,7 @@ class WaterQualityGUI(QMainWindow):
|
||||
self._disable_wheel_for_all_spinboxes()
|
||||
|
||||
# 第五步：默认选中第一个步骤（延迟执行，确保导航列表和 Tab 均已就位）
|
||||
|
||||
QTimer.singleShot(120, self._select_first_nav_item)
|
||||
|
||||
# 第六步：延迟启动工作目录选择
|
||||
|
||||
@ -214,18 +214,35 @@ class ContentMapper:
|
||||
if proj:
|
||||
srs = osr.SpatialReference()
|
||||
srs.ImportFromWkt(proj)
|
||||
# 优先尝试 AuthorityCode
|
||||
epsg = srs.GetAuthorityCode(None)
|
||||
if epsg:
|
||||
return f'EPSG:{epsg}'
|
||||
# ★ 兜底：从 Proj4 字符串解析 EPSG
|
||||
proj4_str = srs.ExportToProj4()
|
||||
if proj4_str:
|
||||
import re
|
||||
# 匹配 +init=epsg:XXXX 或 +init=EPSG:XXXX
|
||||
m = re.search(r'\+init=epsg:(\d+)', proj4_str, re.IGNORECASE)
|
||||
if m:
|
||||
return f'EPSG:{m.group(1)}'
|
||||
# 匹配 UTM zone 信息: +zone=XX +south?
|
||||
zone_m = re.search(r'\+zone=(\d+)', proj4_str)
|
||||
if zone_m and '+proj=utm' in proj4_str:
|
||||
zone = int(zone_m.group(1))
|
||||
south = '+south' in proj4_str
|
||||
return f'EPSG:{32700 + zone if south else 32600 + zone}'
|
||||
except Exception:
|
||||
pass
|
||||
# ── rasterio ──
|
||||
if suffix in ('.tif', '.tiff'):
|
||||
# ── rasterio（扩展支持 .dat / .bsq 等 ENVI 格式）──
|
||||
if suffix in ('.tif', '.tiff', '.dat', '.bsq', '.bil', '.bip', '.img'):
|
||||
try:
|
||||
import rasterio
|
||||
with rasterio.open(file_path) as src:
|
||||
if src.crs is not None:
|
||||
return src.crs.to_string()
|
||||
crs_str = src.crs.to_string()
|
||||
if crs_str:
|
||||
return crs_str
|
||||
except Exception:
|
||||
pass
|
||||
# ── SHP ──
|
||||
@ -239,6 +256,44 @@ class ContentMapper:
|
||||
pass
|
||||
return None
|
||||
|
||||
def _is_geographic_crs(self) -> bool:
|
||||
"""检查 output_crs 是否为地理坐标系（经纬度）。
|
||||
|
||||
用于触发数学空间临时投影：当 CRS 是 EPSG:4326 等地理坐标系时，
|
||||
坐标单位为度，数值极小（如 0.00001°），导致 Kriging/IDW 矩阵
|
||||
计算崩溃、输出 std 退化为 0（纯色图）。
|
||||
"""
|
||||
if self.output_crs is None:
|
||||
return False
|
||||
try:
|
||||
crs_obj = CRS.from_string(self.output_crs)
|
||||
return crs_obj.is_geographic
|
||||
except Exception:
|
||||
return False
|
||||
|
||||
@staticmethod
|
||||
def _get_utm_epsg(lon: float, lat: float) -> int:
|
||||
"""根据经纬度计算所在 UTM 投影带 EPSG 代码。
|
||||
|
||||
Parameters
|
||||
----------
|
||||
lon : float
|
||||
经度（度，WGS84）
|
||||
lat : float
|
||||
纬度（度，WGS84）
|
||||
|
||||
Returns
|
||||
-------
|
||||
int
|
||||
UTM EPSG 代码，如北半球 115°E → 32650，南半球 → 32750
|
||||
"""
|
||||
zone = int((lon + 180) // 6) + 1
|
||||
# 北半球: EPSG:32601–32660; 南半球: EPSG:32701–32760
|
||||
if lat >= 0:
|
||||
return 32600 + zone
|
||||
else:
|
||||
return 32700 + zone
|
||||
|
||||
# ── 内部工具 ─────────────────────────────────────────────────────
|
||||
@staticmethod
|
||||
def _get_chinese_title(stem: str) -> str:
|
||||
@ -646,6 +701,68 @@ class ContentMapper:
|
||||
print(f" - 网格大小: {grid_xx.shape}")
|
||||
print(f" - 坐标系: {self.output_crs}")
|
||||
|
||||
# ═══════════════════════════════════════════════════════════
|
||||
# ★ 数学空间临时投影：地理坐标系 → UTM 米制空间
|
||||
# 经纬度数值极小（如 0.00001°），直接传入 scipy/pykrige
|
||||
# 会导致底层矩阵计算崩溃、输出 std 退化为 0.000000（纯色图）。
|
||||
# 此处将采样点和网格坐标临时转为 UTM 米，插值结果按原网格
|
||||
# 形状填回 — 坐标一一映射，无需空间重采样。
|
||||
#
|
||||
# ★ 防御：若 CRS 标记为地理坐标系，但坐标数值本身超出
|
||||
# 经纬度合理范围（如 CSV 中 longitude 列实际存储 UTM 米），
|
||||
# 则判定 CRS 探测错误，跳过投影，直接用原始米制坐标插值。
|
||||
# ═══════════════════════════════════════════════════════════
|
||||
_orig_grid_shape = grid_xx.shape
|
||||
_is_geo = self._is_geographic_crs()
|
||||
if _is_geo:
|
||||
_mean_x = float(np.mean(points[:, 0]))
|
||||
_mean_y = float(np.mean(points[:, 1]))
|
||||
_x_min, _x_max = float(points[:, 0].min()), float(points[:, 0].max())
|
||||
_y_min, _y_max = float(points[:, 1].min()), float(points[:, 1].max())
|
||||
|
||||
# ★ 防御：坐标值是否在经纬度合理范围内？
|
||||
_x_is_geo = (-180.0 <= _x_min <= 180.0) and (-180.0 <= _x_max <= 180.0)
|
||||
_y_is_geo = (-90.0 <= _y_min <= 90.0) and (-90.0 <= _y_max <= 90.0)
|
||||
|
||||
if not (_x_is_geo and _y_is_geo):
|
||||
# CRS 标记为地理坐标系，但实际坐标值是 UTM 米（数十万～数百万），
|
||||
# 说明 _ensure_crs 未能从边界文件探测到正确投影，CSV 中已存
|
||||
# 投影坐标。跳过投影，直接用米制坐标插值。
|
||||
print(f"[数学投影] ⚠️ CRS 标记为 {self.output_crs}，但坐标值 "
|
||||
f"(X∈[{_x_min:.1f}, {_x_max:.1f}], Y∈[{_y_min:.1f}, {_y_max:.1f}]) "
|
||||
f"超出经纬度范围，判定为已投影的米制坐标，跳过投影直接插值")
|
||||
else:
|
||||
# 坐标值确实是经纬度 → 执行 UTM 临时投影
|
||||
_utm_epsg = self._get_utm_epsg(_mean_x, _mean_y)
|
||||
print(f"[数学投影] 检测到地理坐标系 {self.output_crs}，"
|
||||
f"均值位置 (lon={_mean_x:.4f}°, lat={_mean_y:.4f}°) → "
|
||||
f"临时投影到 EPSG:{_utm_epsg} (UTM 米) 进行插值")
|
||||
|
||||
# 2) 创建临时 transformer: 原始地理 CRS → UTM
|
||||
_trans_to_utm = Transformer.from_crs(
|
||||
CRS.from_string(self.output_crs),
|
||||
CRS.from_string(f'EPSG:{_utm_epsg}'),
|
||||
always_xy=True,
|
||||
)
|
||||
|
||||
# 3) 变换采样点 → UTM 米
|
||||
_pts_x_m, _pts_y_m = _trans_to_utm.transform(
|
||||
points[:, 0], points[:, 1],
|
||||
)
|
||||
points = np.column_stack((_pts_x_m, _pts_y_m))
|
||||
|
||||
# 4) 变换网格坐标 → UTM 米（展平 → 转换 → 重塑）
|
||||
_gx_flat = grid_xx.ravel()
|
||||
_gy_flat = grid_yy.ravel()
|
||||
_gx_m_flat, _gy_m_flat = _trans_to_utm.transform(_gx_flat, _gy_flat)
|
||||
grid_xx = _gx_m_flat.reshape(_orig_grid_shape)
|
||||
grid_yy = _gy_m_flat.reshape(_orig_grid_shape)
|
||||
|
||||
print(f"[数学投影] 变换完成: "
|
||||
f"采样点范围 X[{points[:, 0].min():.1f}, {points[:, 0].max():.1f}]m, "
|
||||
f"Y[{points[:, 1].min():.1f}, {points[:, 1].max():.1f}]m")
|
||||
# ═══════════════════════════════════════════════════════════
|
||||
|
||||
# 检查数据的有效性
|
||||
finite_mask = np.isfinite(values)
|
||||
if not np.all(finite_mask):
|
||||
|
||||
Reference in New Issue
Block a user