From b59371f441ade06ecae777141e4bb4369b6acc20 Mon Sep 17 00:00:00 2001 From: duxin Date: Tue, 28 Jul 2026 14:59:16 +0800 Subject: [PATCH] =?UTF-8?q?fix:=20=E4=B8=93=E9=A2=98=E5=9B=BE=20CRS=20?= =?UTF-8?q?=E6=8E=A2=E6=B5=8B=E5=A2=9E=E5=BC=BA=20+=20=E5=9C=B0=E7=90=86?= =?UTF-8?q?=E5=9D=90=E6=A0=87=E7=B3=BB=20UTM=20=E4=B8=B4=E6=97=B6=E6=8A=95?= =?UTF-8?q?=E5=BD=B1?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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: 格式清理 --- src/gui/water_quality_gui_v2.py | 1 + src/postprocessing/map.py | 123 +++++++++++++++++++++++++++++++- 2 files changed, 121 insertions(+), 3 deletions(-) diff --git a/src/gui/water_quality_gui_v2.py b/src/gui/water_quality_gui_v2.py index bb746e9..76af3f1 100644 --- a/src/gui/water_quality_gui_v2.py +++ b/src/gui/water_quality_gui_v2.py @@ -362,6 +362,7 @@ class WaterQualityGUI(QMainWindow): self._disable_wheel_for_all_spinboxes() # 第五步：默认选中第一个步骤（延迟执行，确保导航列表和 Tab 均已就位） + QTimer.singleShot(120, self._select_first_nav_item) # 第六步：延迟启动工作目录选择 diff --git a/src/postprocessing/map.py b/src/postprocessing/map.py index 43e8d27..6ac6d73 100644 --- a/src/postprocessing/map.py +++ b/src/postprocessing/map.py @@ -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):