From 43f50ec07b4cd6ffcd4a818d6bb85bb13618ab9e Mon Sep 17 00:00:00 2001 From: DXC Date: Thu, 25 Jun 2026 15:50:02 +0800 Subject: [PATCH] =?UTF-8?q?=E6=B5=8B=E8=AF=95=E4=BF=AE=E6=94=B9?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../waterindex_inversion/csv_processor.py | 20 + src/core/handlers/step12_kriging.py | 5 +- src/core/steps/mapping_step.py | 11 +- src/gui/core/panel_registry.py | 2 +- src/gui/panels/step11_map_panel.py | 29 +- src/gui/panels/step12_viz_panel.py | 630 ++++++++-------- src/new/main_view.py | 2 +- src/new/services/step11_service.py | 3 +- src/new/services/step12_service.py | 105 ++- src/new/views/step11_view.py | 13 +- src/new/views/step12_view.py | 11 + src/postprocessing/map.py | 712 ++++++++++++------ src/postprocessing/point_map.py | 564 +++----------- src/postprocessing/visualization_reports.py | 10 +- src/utils/band_math.py | 6 +- src/utils/water_index.py | 4 +- 16 files changed, 1098 insertions(+), 1029 deletions(-) diff --git a/src/core/algorithms/waterindex_inversion/csv_processor.py b/src/core/algorithms/waterindex_inversion/csv_processor.py index 58374be..3ea3976 100644 --- a/src/core/algorithms/waterindex_inversion/csv_processor.py +++ b/src/core/algorithms/waterindex_inversion/csv_processor.py @@ -41,6 +41,8 @@ import os import re from typing import Callable, Dict, List, Optional +import numpy as np # P0 修复: 写盘前 inf/-inf 替换与极值截断需要 + class WaterIndexCsvProcessor: """ @@ -200,6 +202,24 @@ class WaterIndexCsvProcessor: for i, name in enumerate(targets): try: per_idx = results_df[name] + # ===== P0 防御: 写盘前清洗（防 Step 11 Kriging / TIN 碎玻璃）===== + # 1) inf / -inf → NaN：pandas 默认会把 inf 写成字面 "Infinity"， + # 下游 ContentMapper 严格按位置读第 3 列时会原样拿到 inf， + # Kriging 变差函数被单一 inf 点拉飞、TIN 出现极长退化三角形。 + # 2) 物理范围截断 ±10：保守大窗，覆盖几乎所有经验公式的合法值域 + # （NDCI 类 ∈ [-1,1]、比值类 ∈ [0,2]、浓度反演类 ∈ [0, ∞)）。 + # 注：保留所有行（包括 NaN），让下游 Step 11 把它当 missing marker； + # dropna 会破坏行号对齐，故不在此调用；如确需剔除 NaN 行， + # 在面板 / service 层读取后自行 .dropna(subset=[col])。 + n_inf = int(np.isinf(per_idx.values).sum()) + per_idx = ( + per_idx + .replace([np.inf, -np.inf], np.nan) + .clip(lower=-10.0, upper=10.0) + ) + if n_inf > 0: + print(f"[WaterIndexCsvProcessor] {name}: 替换 {n_inf} 个 inf/-inf → NaN") + # ===== P0 防御结束 ===== out_df = pd.DataFrame({ "longitude": df["longitude"].values, "latitude": df["latitude"].values, diff --git a/src/core/handlers/step12_kriging.py b/src/core/handlers/step12_kriging.py index d3b8807..8328a10 100644 --- a/src/core/handlers/step12_kriging.py +++ b/src/core/handlers/step12_kriging.py @@ -51,7 +51,10 @@ class Step12KrigingHandler(BaseStepHandler): output_image_path=output_image_path, resolution=config.get('resolution', 30), input_crs=config.get('input_crs', 'EPSG:32651'), - output_crs=config.get('output_crs', 'EPSG:4326'), + # ★★★ 强制 output_crs = input_crs，避免 ContentMapper 把栅格重投影到 EPSG:4326 ★★★ + # 旧实现：output_crs=config.get('output_crs', 'EPSG:4326') + # 重投影会让栅格和基于投影坐标的掩膜在 visualize_raster 叠加时发生仿射变换撕裂 + output_crs=config.get('input_crs', 'EPSG:32651'), show_sample_points=config.get('show_sample_points', False), base_map_tif=config.get('base_map_tif'), use_distance_diffusion=config.get('use_distance_diffusion', True), diff --git a/src/core/steps/mapping_step.py b/src/core/steps/mapping_step.py index fe707df..e48998a 100644 --- a/src/core/steps/mapping_step.py +++ b/src/core/steps/mapping_step.py @@ -16,11 +16,14 @@ class MappingStep: @staticmethod def generate_distribution_map( prediction_csv_path: str, - boundary_shp_path: str, + boundary_shp_path: Optional[str] = None, # ★★★ Plan C: None = 不依赖水域掩膜 ★★★ output_image_path: Optional[str] = None, resolution: float = 30, input_crs: str = "EPSG:32651", - output_crs: str = "EPSG:4326", + # ★★★ 强制默认 output_crs = input_crs，禁止重投影到 EPSG:4326 ★★★ + # 历史默认值 'EPSG:4326' 会让 ContentMapper 将插值栅格从投影坐标系转到经纬度坐标系， + # 与基于 EPSG:32651 的水域掩膜叠加时发生仿射变换撕裂（栅格错位、坐标轴扭曲）。 + output_crs: str = "EPSG:32651", show_sample_points: bool = False, base_map_tif: Optional[str] = None, use_distance_diffusion: bool = True, @@ -31,13 +34,14 @@ class MappingStep: expand_ratio: float = 0.05, output_dir: Union[str, Path] = "./14_visualization", callback: Optional[Callable] = None, + output_format: str = 'tif', # ⭐ 新增：'tif' (默认, 写 GeoTIFF) / 'png' ) -> str: """ 根据采样点的坐标和反演的实测参数，通过插值方法得到水质参数可视化分布图 Args: prediction_csv_path: 预测结果CSV文件路径（前两列为经纬度，第三列为预测值） - boundary_shp_path: 边界shapefile文件路径 + boundary_shp_path: 边界/掩膜文件路径（.shp/.dat/.bsq/.tif 等）。None 时跳过水域掩膜约束。 output_image_path: 输出图片路径（如果为None，自动生成） resolution: 插值网格分辨率（米） input_crs: 输入坐标系 @@ -89,6 +93,7 @@ class MappingStep: "diffusion_power": diffusion_power, "diffusion_n_neighbors": diffusion_n_neighbors, "expand_ratio": expand_ratio, + "output_format": output_format, # ⭐ 透传给 ContentMapper.process_data } optional_kwargs = { diff --git a/src/gui/core/panel_registry.py b/src/gui/core/panel_registry.py index 66c01cb..86b20b4 100644 --- a/src/gui/core/panel_registry.py +++ b/src/gui/core/panel_registry.py @@ -201,7 +201,7 @@ PANEL_REGISTRY = [ 'title': '专题图生成', 'icon': '10.png', 'stage': '模块四 制图与成果汇编', - 'display_name': '11. 专题图生成', + 'display_name': '11. 分布图生成', # 架构解耦（2026-06-22）：dict 键名 = 下游目标控件真实属性名 'dependencies': { # 目标框: self.prediction_csv_dir_edit ← 上游 step9_ml_predict.output_file diff --git a/src/gui/panels/step11_map_panel.py b/src/gui/panels/step11_map_panel.py index fd6629d..132743f 100644 --- a/src/gui/panels/step11_map_panel.py +++ b/src/gui/panels/step11_map_panel.py @@ -290,7 +290,10 @@ class Step11MapPanel(QWidget): params_layout.addRow("输入坐标系:", self.input_crs) self.output_crs = QLineEdit() - self.output_crs.setText("EPSG:4326") + # ★★★ 强制默认输出坐标系与输入一致，禁止从 GUI 误改为 EPSG:4326 ★★★ + # 历史默认值 'EPSG:4326' 会让 ContentMapper 把栅格重投影到经纬度， + # 与基于 EPSG:32651 的水域掩膜叠加时发生仿射变换撕裂（栅格错位、坐标轴扭曲）。 + self.output_crs.setText("EPSG:32651") params_layout.addRow("输出坐标系:", self.output_crs) self.show_points = QCheckBox("显示采样点") @@ -416,7 +419,9 @@ class Step11MapPanel(QWidget): 'boundary_shp_path': self.boundary_file.get_path(), 'resolution': self.resolution.value(), 'input_crs': self.input_crs.text(), - 'output_crs': self.output_crs.text(), + # ★★★ 强制 output_crs = input_crs，不再信任 GUI 的 output_crs 输入框 ★★★ + # 否则 ContentMapper 会把栅格重投影到 EPSG:4326，与水域掩膜叠加时撕裂。 + 'output_crs': self.input_crs.text(), 'show_sample_points': self.show_points.isChecked(), 'use_distance_diffusion': self.use_diffusion.isChecked(), } @@ -437,7 +442,8 @@ class Step11MapPanel(QWidget): 'boundary_shp_path': self.boundary_file.get_path(), 'resolution': self.resolution.value(), 'input_crs': self.input_crs.text(), - 'output_crs': self.output_crs.text(), + # ★★★ 强制 output_crs = input_crs，不再信任 GUI 的 output_crs 输入框 ★★★ + 'output_crs': self.input_crs.text(), 'show_sample_points': self.show_points.isChecked(), 'use_distance_diffusion': self.use_diffusion.isChecked(), } @@ -475,8 +481,9 @@ class Step11MapPanel(QWidget): self.resolution.setValue(config['resolution']) if 'input_crs' in config: self.input_crs.setText(config['input_crs']) - if 'output_crs' in config: - self.output_crs.setText(config['output_crs']) + # ★★★ 反灌入时强制 output_crs = input_crs，避免旧 config 中的 EPSG:4326 回填 ★★★ + if 'output_crs' in config or 'input_crs' in config: + self.output_crs.setText(config.get('input_crs') or config.get('output_crs') or 'EPSG:32651') if 'show_sample_points' in config: self.show_points.setChecked(config['show_sample_points']) if 'use_distance_diffusion' in config: @@ -573,10 +580,8 @@ class Step11MapPanel(QWidget): return boundary_shp_path = self.boundary_file.get_path() - if not boundary_shp_path: - QMessageBox.warning(self, "输入验证失败", "请选择边界文件") - return - if not os.path.exists(boundary_shp_path): + # ── Plan C: 允许 boundary_shp_path 为空（跳过水域掩膜，纯采样点插值）──── + if boundary_shp_path and not os.path.exists(boundary_shp_path): QMessageBox.warning(self, "输入验证失败", "边界文件不存在") return @@ -631,7 +636,8 @@ class Step11MapPanel(QWidget): output_dir=out_dir, boundary_shp_path=boundary_shp_path, input_crs=self.input_crs.text(), - output_crs=self.output_crs.text(), + # ★★★ 强制 output_crs = input_crs，不再信任 GUI 的 output_crs 输入框 ★★★ + output_crs=self.input_crs.text(), ) main_win = parent @@ -671,7 +677,8 @@ class Step11MapPanel(QWidget): boundary_shp_path = self.boundary_file.get_path() input_crs = self.input_crs.text() - output_crs = self.output_crs.text() + # ★★★ 强制 output_crs = input_crs，不再信任 GUI 的 output_crs 输入框 ★★★ + output_crs = self.input_crs.text() # 构造输出路径 out_dir = (self.output_dir.get_path() or "").strip() diff --git a/src/gui/panels/step12_viz_panel.py b/src/gui/panels/step12_viz_panel.py index cbe45ff..5f54234 100644 --- a/src/gui/panels/step12_viz_panel.py +++ b/src/gui/panels/step12_viz_panel.py @@ -440,6 +440,62 @@ class VisualizationWorkerThread(QThread): parts.append(f"浓度统计图: 失败({e})") else: parts.append("浓度统计图: 跳过(无浓度CSV)") + + if self.extra.get("gen_distribution_map"): + dist_dir = wp / "11_Thematic_Map" + out_sub = Path(viz.output_dir) / "distribution_maps" + out_sub.mkdir(parents=True, exist_ok=True) + n_rendered = 0 + if dist_dir.exists(): + import shutil + # 1. 拷贝现成的 PNG + for png in list(dist_dir.glob("*_distribution.png")) + list(dist_dir.glob("*_专题图.png")): + try: + shutil.copy2(png, out_sub / png.name) + n_rendered += 1 + except Exception: + pass + + # 2. 渲染生成的 TIF + tif_files = list(dist_dir.glob("*_distribution.tif")) + list(dist_dir.glob("*_kriging.tif")) + if tif_files: + from src.postprocessing.map import ContentMapper + mapper = ContentMapper() + # 优先级1：直接使用 Step 1 面板中缓存的外部原始 .shp 绝对路径！ + boundary_path = self.extra.get("boundary_shp_path") + + # 优先级2：如果没拿到，全局搜索整个工作目录下的 .shp 文件（放宽限制） + if not boundary_path: + shp_candidates = list(wp.rglob("**/*.shp")) + if shp_candidates: + boundary_path = str(shp_candidates[0]) + + # 优先级3：兜底使用 1_water_mask 下的栅格掩膜 + if not boundary_path: + mask_files = list(wp.rglob("1_water_mask/*")) + other_candidates = [f for f in mask_files if f.suffix.lower() in ('.dat', '.bsq', '.tif', '.tiff')] + if other_candidates: + boundary_path = str(other_candidates[0]) + if not boundary_path: + print(f"[distribution_maps] 未找到水域边界文件，跳过裁剪") + for tif in tif_files: + dst_png = out_sub / f"{tif.stem}_rendered.png" + try: + mapper.visualize_raster( + raster_tif_path=str(tif), + output_file=str(dst_png), + boundary_shp_path=boundary_path, + nodata_value=-9999.0, + figsize=(14, 10), + alpha=0.9 + ) + n_rendered += 1 + except Exception as e: + print(f"渲染 TIF 失败 {tif.name}: {e}") + parts.append(f"空间分布图: {n_rendered} 个") + else: + parts.append("空间分布图: 跳过(无 11_Thematic_Map 目录)") + self.finished_ok.emit({"task": "generate_all_selected", "parts": parts}) else: self.failed.emit(f"未知可视化任务: {self.task}") @@ -574,272 +630,104 @@ class ChartViewerDialog(QDialog): class ImageCategoryTree(QTreeWidget): - """图像分类目录树 - 按真实物理文件夹结构组织图像文件""" - - # 文件名中文翻译映射（key: 文件名前缀 → 中文显示名） - NAME_MAPPING = { - "hsi_preview": "高光谱影像预览", - "hsi_original": "原始高光谱影像", - "hsi_deglint": "去耀斑高光谱影像", - "water_mask_overlay": "水域掩膜叠加图", - "water_mask": "水域掩膜图", - "glint_mask": "耀斑掩膜预览", - "glint_overlay": "耀斑叠加对比图", - "deglint_comparison": "去耀斑前后对比", - "training_spectra": "训练光谱特征", - "spectrum_by_param": "参数光谱图", - "model_evaluation": "模型评估散点图", - "model_scatter": "模型散点图", - "regression": "回归分析图", - "validation": "验证结果图", - "spatial_distribution": "参数空间分布图", - "distribution_map": "分布图", - "thematic_map": "水质专题图", - "water_quality_map": "水质分布图", - "prediction_map": "预测结果图", - "inversion_map": "反演结果图", - "correlation_matrix": "特征相关性矩阵", - "feature_correlation": "特征相关性", - "sampling_point_map": "采样点分布图", - "sampling_points": "采样点图", - "point_locations": "采样位置图", - "boxplot": "箱线图", - "histogram": "直方图", - "statistics": "统计图表", - "statistical_chart": "统计图", - "error_analysis": "误差分析图", - "rmse": "RMSE评估图", - "r2_score": "R²得分图", - "flight": "飞行轨迹图", - "path": "轨迹图", - "trajectory": "轨迹图", - "glint_deglint": "耀斑去耀斑影像", - "enhanced": "增强分布图", - "content": "含量分布图", - "distribution": "分布图", - "prediction": "预测图", - "inversion": "反演图", - "scatter_true_vs_pred": "真值-预测散点图", - "true_vs_pred": "真值-预测散点图", - "correlation_heatmap": "相关性热力图", - "parameter_boxplot": "水质参数箱线图", - "spectrum_comparison": "光谱曲线对比图", - "scatter": "散点图", - } - - # 目录层级中文翻译 - DIR_MAPPING = { - "14_visualization": "统计与分析报表", - "1_water_mask": "水域掩膜识别", - "2_Glint_Detection": "耀斑区域检测", - "3_deglint": "去耀斑影像结果", - "5_training_spectra": "训练光谱特征", - "8_Regression_Modeling": "回归建模分析", - "9_water_quality_prediction": "水质预测结果", - "10_feature_construction": "特征构建散点", - "11_12_13_predictions": "空间分布专题图", - "glint_deglint_previews": "耀斑处理预览", - "sampling_maps": "采样点空间分布", - "flight_maps": "无人机飞行轨迹", - "9_ML_Prediction": "机器学习预测", - "Non_Empirical_Prediction": "非经验模型预测", - "Custom_Regression_Prediction": "自定义回归预测", - "boxplot_dir": "水质参数箱线图", - "boxplot": "水质参数箱线图", - "output_dir": "输出目录", - "8_spatial_inversion": "空间反演", - "4_processed_data": "处理数据", - "9_Concentration": "物理反演浓度分布", - } + """现代化的图像分类目录树 - 支持智能归类、筛选和高颜值样式""" def __init__(self, parent=None): super().__init__(parent) - self._dir_node_map: dict = {} # 目录路径字符串 → QTreeWidgetItem - self._work_path: Optional[Path] = None - self.setHeaderLabel("图像目录") - self.setMaximumWidth(300) - self.setMinimumWidth(250) + self._work_path = None + self._all_image_files = [] # 缓存所有扫描到的图片路径 + + self.setHeaderHidden(True) # 隐藏表头，显得更清爽 + self.setAlternatingRowColors(True) # 斑马纹交替背景 + self.setMaximumWidth(340) + self.setMinimumWidth(280) + + # 核心：高颜值现代化 CSS 样式 self.setStyleSheet(""" QTreeWidget { - border: 1px solid #ddd; - border-radius: 5px; - background-color: #f8f9fa; + border: 1px solid #E2E8F0; + border-radius: 8px; + background-color: #FFFFFF; + alternate-background-color: #F8FAFC; + font-family: "Microsoft YaHei", "Segoe UI"; + font-size: 13px; + padding: 4px; } QTreeWidget::item { - padding: 5px; - border-radius: 3px; - } - QTreeWidget::item:selected { - background-color: #0078D4; - color: white; + height: 32px; + border-radius: 6px; + margin: 2px 4px; } QTreeWidget::item:hover { - background-color: #e3f2fd; + background-color: #F1F5F9; + } + QTreeWidget::item:selected { + background-color: #E0F2FE; + color: #0369A1; + font-weight: bold; + } + QTreeWidget::branch:has-children:!has-siblings:closed, + QTreeWidget::branch:closed:has-children:has-siblings { + border-image: none; + image: none; + } + QTreeWidget::branch:open:has-children:!has-siblings, + QTreeWidget::branch:open:has-children:has-siblings { + border-image: none; + image: none; } """) - def clear_all_images(self): - """清除所有图像项""" - try: - self.invisibleRootItem().takeChildren() - if hasattr(self, '_dir_node_map'): - self._dir_node_map.clear() - except Exception as e: - print(f"清空树状图出错: {e}") - import traceback - traceback.print_exc() + def _parse_file_info(self, file_path: Path): + """智能解析文件名，提取【水质参数】和【图表类型】""" + name_upper = file_path.name.upper() - def _translate_dir_name(self, dir_name: str) -> str: - """翻译目录名为中文""" - return self.DIR_MAPPING.get(dir_name, dir_name) + # 1. 提取图表类型 + chart_type = "其他图表" + if "DISTRIBUTION" in name_upper or "专题图" in name_upper or "RENDERED" in name_upper: + chart_type = "空间分布图" + elif "SCATTER" in name_upper or "散点" in name_upper: + chart_type = "模型散点图" + elif "SPECTRUM" in name_upper or "光谱" in name_upper: + chart_type = "光谱曲线图" + elif "HEATMAP" in name_upper or "热力图" in name_upper: + chart_type = "相关性热力图" + elif "BOXPLOT" in name_upper or "HISTOGRAM" in name_upper or "箱线" in name_upper or "直方" in name_upper: + chart_type = "统计箱线图" + elif "SAMPLING" in name_upper or "采样" in name_upper: + chart_type = "采样点地图" + elif "GLINT" in name_upper or "MASK" in name_upper or "PREVIEW" in name_upper: + chart_type = "掩膜与预览" - def _translate_filename(self, filename: str) -> str: - # 1. 后缀替换 (图表类型) - type_mapping = { - '_scatter_true_vs_pred': ' 真值预测散点图', - '_spectrum_comparison': ' 光谱曲线对比图', - '_spectrum': ' 光谱特征图', - '_histogram': ' 分布直方图', - '_boxplot_seaborn': ' Seaborn箱线图', - '_boxplot': ' 箱线图', - '_distribution_enhanced': ' 增强空间分布图', - '_distribution': ' 空间分布图', - '_sampling_map': ' 采样点地图', - '_flight_paths': ' 飞行轨迹图', - '_preview': ' 效果预览图', - 'water_mask_overlay': '水域掩膜叠加图', - 'hsi_preview': '原始影像预览', - 'correlation_heatmap': '特征相关性热力图', - 'parameter_boxplot': '水质参数汇总箱线图', - 'all_parameters_boxplot': '全参数汇总箱线图', - 'content_map': '含量分布专题图', - '_scatter_with_confidence': ' 置信区间散点图' + # 2. 提取参数名 + param_name = "综合/未分类" + # 常见水质参数字典映射 + params_map = { + 'CHLOROPHYLL': 'Chlorophyll (叶绿素)', 'CHL_A': 'Chlorophyll (叶绿素)', 'CHLA': 'Chlorophyll (叶绿素)', + 'COD': 'COD (化学需氧量)', 'DO': 'DO (溶解氧)', 'PH': 'pH', + 'TEMPERATURE': 'Temperature (温度)', 'SPCOND': 'spCond (电导率)', + 'TURBIDITY': 'Turbidity (浊度)', 'TDS': 'TDS (总溶解固体)', + 'CL-': 'Cl- (氯离子)', 'NO3-N': 'NO3-N (硝态氮)', 'NH3-N': 'NH3-N (氨氮)', + 'BGA': 'BGA (蓝绿藻)', 'TT': 'TT (透明度)' } + for key, display_name in params_map.items(): + if key in name_upper or key.replace('-', '') in name_upper: + param_name = display_name + break - name = filename - for eng, chn in type_mapping.items(): - if eng in name: - name = name.replace(eng, chn) - - # 2. 常见水质参数前缀替换 - param_mapping = { - 'Chlorophyll': '叶绿素', 'Chl_a': '叶绿素a', 'Chla': '叶绿素a', - 'Turbidity': '浊度', 'Temperature': '温度', 'spCond': '电导率', - 'COD': '化学需氧量', 'DO': '溶解氧', 'PH': 'pH值', 'TDS': '总溶解固体', - 'BGA': '蓝绿藻', 'TT': '透明度', 'NH3-N': '氨氮', 'NO3-N': '硝酸盐氮', - 'glint_severe_glint_area': '重度耀斑区域', - 'severe_glint_area': '重度耀斑区域', - 'deglint_goodman': 'Goodman算法去耀斑', - 'deglint_Goodman': 'Goodman算法去耀斑', - 'glint_': '耀斑检测_', - 'deglint_': '耀斑去除_', - } - - for eng, chn in param_mapping.items(): - if name.startswith(eng + ' ') or name.startswith(eng + '_'): - name = name.replace(eng, chn, 1) - elif eng in name: - name = name.replace(eng, chn) - - return name.strip('_') - - def add_image_by_dir(self, file_path: Path, work_path: Path): - """按真实物理目录层级挂载图片节点 - - Args: - file_path: 图片文件的完整路径 - work_path: 工作目录根路径 - """ - # 计算相对路径 - try: - rel_path = file_path.relative_to(work_path) - except ValueError: - rel_path = Path(file_path.name) - - # 分离父目录链和文件名 - parts = rel_path.parts - if len(parts) <= 1: - parent_key = "__root__" - parent_display = "根目录" - else: - # 父目录路径（相对于work_path） - parent_key = str(Path(*parts[:-1])) - # 取最后一层目录名作为显示名 - parent_display = self._translate_dir_name(parts[-2]) - - # 根目录节点特殊处理 - root_display = self._translate_dir_name(parts[0]) if parts else "根目录" - - # 获取或创建根目录节点 - if root_display not in self._dir_node_map: - root_item = QTreeWidgetItem(self) - root_item.setText(0, f"📁 {root_display}") - root_item.setData(0, Qt.UserRole, {"type": "root_dir", "path": str(work_path / parts[0])}) - root_item.setExpanded(True) - self._dir_node_map[root_display] = root_item - self._dir_node_map[f"__root__{root_display}"] = root_item - root_item = self._dir_node_map.get(f"__root__{root_display}") - - if len(parts) > 1: - # 获取或创建子目录节点 - if parent_key not in self._dir_node_map: - dir_item = QTreeWidgetItem(root_item) - dir_item.setText(0, f" 📂 {parent_display}") - dir_item.setData(0, Qt.UserRole, {"type": "sub_dir", "path": str(work_path / parent_key)}) - dir_item.setExpanded(True) - self._dir_node_map[parent_key] = dir_item - parent_item = self._dir_node_map[parent_key] - else: - parent_item = root_item - - # 创建图片节点（根据翻译后的名称分配图标） - display_name = self._translate_filename(file_path.stem) + file_path.suffix - icon = "🖼️" # 默认 - if "散点" in display_name: - icon = "📊" - elif "光谱" in display_name or "曲线" in display_name: - icon = "📈" - elif "箱线" in display_name or "直方" in display_name: - icon = "📉" - elif "分布" in display_name or "地图" in display_name or "轨迹" in display_name: - icon = "🗺️" - image_item = QTreeWidgetItem(parent_item) - image_item.setText(0, f" {icon} {display_name}") - image_item.setData(0, Qt.UserRole, {"type": "image", "path": str(file_path), "display_name": display_name}) - image_item.setToolTip(0, str(file_path)) - - return image_item + return param_name, chart_type def scan_directory(self, work_dir: str): - """扫描目录中的所有图像文件（深度递归扫描）—— 按真实物理目录结构挂载""" + """全量扫描文件并缓存（不直接构建树，而是交给 rebuild_tree 渲染）""" try: - if not work_dir: - print("可视化面板：工作目录为空，跳过扫描") - return - + if not work_dir: return self._work_path = Path(work_dir) + if not self._work_path.exists(): return - # 阻塞信号，防止在清空树状图时触发 selected 槽函数导致崩溃 - # 因为当前类继承自 QTreeWidget，所以 self 本身就是树 - self.blockSignals(True) - self.clear_all_images() - self.blockSignals(False) - - if not self._work_path.exists(): - return - except Exception as e: - import traceback - print(f"可视化面板初始化扫描出错: {e}") - traceback.print_exc() - # 确保信号锁被解开 - self.blockSignals(False) - return - - try: - image_extensions = ['*.png', '*.jpg', '*.jpeg', '*.tif', '*.tiff', '*.bmp'] - - # 拓宽扫描根目录列表（新增多个遗漏目录） - scan_roots: List[Path] = [ + # 仅扫描用于视觉展示的常规图片格式，屏蔽科学栅格 TIF 以免无法渲染报错 + image_extensions = ['*.png', '*.jpg', '*.jpeg', '*.bmp'] + # 扩展扫描路径 + scan_roots = [ Path(resolve_subdir(str(self._work_path), 'visualization')), Path(resolve_subdir(str(self._work_path), 'prediction_dir')), Path(resolve_subdir(str(self._work_path), 'regression_modeling')), @@ -850,51 +738,98 @@ class ImageCategoryTree(QTreeWidget): Path(resolve_subdir(str(self._work_path), 'water_mask')), self._work_path / "9_water_quality_prediction", self._work_path / "9_Concentration", + self._work_path / "11_Thematic_Map" ] - - # 只保留存在的目录，并补充工作根目录作为兜底 scan_roots = [p for p in scan_roots if p.is_dir()] - if not scan_roots: - scan_roots.append(self._work_path) + if not scan_roots: scan_roots.append(self._work_path) + + seen_norm = set() + self._all_image_files = [] - seen_norm: set = set() - image_files: List[Path] = [] for root in scan_roots: for ext in image_extensions: for p in root.rglob(ext): key = os.path.normcase(os.path.normpath(str(p.resolve()))) - if key in seen_norm: - continue + if key in seen_norm: continue seen_norm.add(key) - image_files.append(p) + if p.name.startswith('.') or 'thumb' in p.name.lower(): continue + self._all_image_files.append(p) - for img_file in sorted(image_files): - if img_file.name.startswith('.') or 'thumb' in img_file.name.lower(): - continue - self.add_image_by_dir(img_file, self._work_path) - - # 更新目录节点计数 - for key, item in self._dir_node_map.items(): - if key.startswith("__root__"): - continue - if item.data(0, Qt.UserRole).get("type") == "sub_dir": - count = item.childCount() - name = item.text(0) - if count > 0 and f"({count})" not in name: - # 从目录名中提取显示名并附加计数 - display = name.strip() - item.setText(0, f" 📂 {display} ({count})") + self._all_image_files.sort(key=lambda x: x.name) + # 默认构建模式 + self.rebuild_tree(group_mode='parameter', filter_type='all') except Exception as e: - import traceback - print(f"可视化面板图片挂载出错: {e}") - traceback.print_exc() + print(f"目录扫描出错: {e}") + + def rebuild_tree(self, group_mode='parameter', filter_type='all'): + """根据下拉框的【模式】和【筛选条件】实时重新构建 UI 树""" + self.blockSignals(True) + self.clear() + root_nodes = {} + + for img_file in self._all_image_files: + param, chart_type = self._parse_file_info(img_file) + + # 1. 应用筛选器逻辑 + if filter_type != 'all' and filter_type != chart_type: + continue + + # 2. 决定分组基准 + if group_mode == 'parameter': + group_key = param + group_icon = "💧" if param != "综合/未分类" else "📁" + display_name = f"[{chart_type}] {img_file.name}" + elif group_mode == 'type': + group_key = chart_type + group_icon = "📊" + display_name = f"[{param.split(' ')[0]}] {img_file.name}" + else: + # 物理文件夹原样模式 + try: + rel_path = img_file.relative_to(self._work_path) + group_key = str(rel_path.parent) if len(rel_path.parts) > 1 else "根目录" + except: + group_key = "其他" + group_icon = "📂" + display_name = img_file.name + + # 3. 创建父节点 + if group_key not in root_nodes: + root_item = QTreeWidgetItem(self) + root_item.setText(0, f"{group_icon} {group_key}") + root_item.setExpanded(True) + font = root_item.font(0) + font.setBold(True) + root_item.setFont(0, font) + root_item.setData(0, Qt.UserRole, {"type": "root"}) + root_nodes[group_key] = root_item + + parent_item = root_nodes[group_key] + + # 4. 挂载子节点及专属图标 + icon = "🖼️" + if "散点" in chart_type: icon = "📌" + elif "光谱" in chart_type: icon = "📈" + elif "分布" in chart_type: icon = "🗺️" + elif "箱线" in chart_type: icon = "📉" + + image_item = QTreeWidgetItem(parent_item) + image_item.setText(0, f" {icon} {display_name}") + image_item.setData(0, Qt.UserRole, {"type": "image", "path": str(img_file)}) + image_item.setToolTip(0, str(img_file)) + + # 统计数量 + for i in range(self.topLevelItemCount()): + root_item = self.topLevelItem(i) + count = root_item.childCount() + old_text = root_item.text(0) + root_item.setText(0, f"{old_text} ({count})") + + self.blockSignals(False) def get_selected_image_path(self) -> Optional[str]: - """获取当前选中的图像路径""" selected_item = self.currentItem() - if not selected_item: - return None - + if not selected_item: return None data = selected_item.data(0, Qt.UserRole) if data and data.get("type") == "image": return data.get("path") @@ -1401,17 +1336,19 @@ class Step12VizPanel(QWidget): QMessageBox.critical(self, "错误", f"可视化任务失败:\n{err[:1200]}") def init_ui(self): - """初始化UI - 使用左右分栏布局""" + """初始化UI - 使用全新的三列布局（控制参数 | 独立满高目录树 | 图像查看器）""" main_layout = QHBoxLayout() - main_layout.setSpacing(10) + main_layout.setSpacing(12) main_layout.setContentsMargins(10, 10, 10, 10) - # ===== 左侧面板 ===== - left_panel = QWidget() - left_layout = QVBoxLayout() - left_layout.setContentsMargins(0, 0, 0, 0) + # ========================================== + # 第一列：控制面板（目录选择 + 生成配置） + # ========================================== + control_panel = QWidget() + control_layout = QVBoxLayout() + control_layout.setContentsMargins(0, 0, 0, 0) - # 工作目录选择 + # 1. 工作目录选择 dir_group = QGroupBox("工作目录") dir_layout = QHBoxLayout() self.work_dir_edit = QLineEdit() @@ -1422,9 +1359,9 @@ class Step12VizPanel(QWidget): dir_layout.addWidget(self.work_dir_edit, 1) dir_layout.addWidget(dir_browse_btn) dir_group.setLayout(dir_layout) - left_layout.addWidget(dir_group) + control_layout.addWidget(dir_group) - # 图像目录选择（优先指向预测结果目录） + # 2. 图像目录选择 img_dir_group = QGroupBox("图像目录") img_dir_layout = QHBoxLayout() self.img_dir_edit = QLineEdit() @@ -1435,18 +1372,9 @@ class Step12VizPanel(QWidget): img_dir_layout.addWidget(self.img_dir_edit, 1) img_dir_layout.addWidget(img_dir_browse_btn) img_dir_group.setLayout(img_dir_layout) - left_layout.addWidget(img_dir_group) + control_layout.addWidget(img_dir_group) - # 图像目录树 - tree_group = QGroupBox("图像目录") - tree_layout = QVBoxLayout() - self.image_tree = ImageCategoryTree() - self.image_tree.itemClicked.connect(self.on_tree_item_clicked) - tree_layout.addWidget(self.image_tree) - tree_group.setLayout(tree_layout) - left_layout.addWidget(tree_group, 1) - - # 可视化配置 + # 3. 可视化配置 config_group = QGroupBox("可视化配置") config_layout = QVBoxLayout() @@ -1470,6 +1398,11 @@ class Step12VizPanel(QWidget): self.gen_sampling_map.setChecked(True) config_layout.addWidget(self.gen_sampling_map) + self.gen_distribution_map = QCheckBox("空间分布图 (Step 11 产物)") + self.gen_distribution_map.setChecked(True) + self.gen_distribution_map.setToolTip("渲染并汇总 Step 11 生成的 TIF 分布图") + config_layout.addWidget(self.gen_distribution_map) + config_layout.addSpacing(10) line = QFrame() line.setFrameShape(QFrame.HLine) @@ -1479,23 +1412,73 @@ class Step12VizPanel(QWidget): self.gen_all_btn = QPushButton("🚀 生成全部") self.gen_all_btn.setToolTip("生成所有类型的可视化图表") - self.gen_all_btn.setStyleSheet("background-color: #4CAF50; color: white; font-weight: bold;") + self.gen_all_btn.setStyleSheet("background-color: #4CAF50; color: white; font-weight: bold; padding: 8px; border-radius: 4px;") self.gen_all_btn.clicked.connect(self.generate_all_visualizations) config_layout.addWidget(self.gen_all_btn) self.scan_btn = QPushButton("📁 扫描目录") self.scan_btn.setToolTip("扫描工作目录中的图像文件") + self.scan_btn.setStyleSheet("padding: 6px; border-radius: 4px;") self.scan_btn.clicked.connect(self.scan_work_directory) config_layout.addWidget(self.scan_btn) config_group.setLayout(config_layout) - left_layout.addWidget(config_group) + control_layout.addWidget(config_group) - left_panel.setLayout(left_layout) - left_panel.setMaximumWidth(350) - main_layout.addWidget(left_panel, 0) + control_layout.addStretch() # 把控制面板的内容往上顶 + control_panel.setLayout(control_layout) + control_panel.setMaximumWidth(280) # 稍微收窄第一列 + control_panel.setMinimumWidth(230) + main_layout.addWidget(control_panel, 0) # stretch=0，不横向拉伸 - # ===== 右侧面板 ===== + # ========================================== + # 第二列：满高独立的目录树（选图与筛选面板） + # ========================================== + from PyQt5.QtWidgets import QComboBox + tree_panel = QWidget() + tree_layout = QVBoxLayout() + tree_layout.setContentsMargins(0, 0, 0, 0) + + tree_group = QGroupBox("图像浏览与筛选") + group_layout = QVBoxLayout() + group_layout.setSpacing(8) + + # 添加过滤控制栏 + filter_layout = QFormLayout() + filter_layout.setContentsMargins(0, 0, 0, 0) + + # 分组模式下拉框 + self.view_mode_cb = QComboBox() + self.view_mode_cb.addItems(["按水质参数归类", "按图表类型归类", "按物理文件夹"]) + self.view_mode_cb.currentIndexChanged.connect(self.update_image_tree_view) + self.view_mode_cb.setStyleSheet("QComboBox { padding: 4px; border-radius: 4px; border: 1px solid #ccc; }") + + # 图表类型筛选下拉框 + self.chart_filter_cb = QComboBox() + self.chart_filter_cb.addItems(["全部图表", "空间分布图", "模型散点图", "光谱曲线图", "统计箱线图", "相关性热力图", "掩膜与预览", "采样点地图"]) + self.chart_filter_cb.currentIndexChanged.connect(self.update_image_tree_view) + self.chart_filter_cb.setStyleSheet("QComboBox { padding: 4px; border-radius: 4px; border: 1px solid #ccc; }") + + filter_layout.addRow("视图模式:", self.view_mode_cb) + filter_layout.addRow("类型筛选:", self.chart_filter_cb) + group_layout.addLayout(filter_layout) + + # 挂载新版大树 + self.image_tree = ImageCategoryTree() + self.image_tree.itemClicked.connect(self.on_tree_item_clicked) + group_layout.addWidget(self.image_tree, 1) # stretch=1 让树垂直填满 + + tree_group.setLayout(group_layout) + tree_layout.addWidget(tree_group, 1) # stretch=1 让GroupBox垂直填满 + + tree_panel.setLayout(tree_layout) + tree_panel.setMaximumWidth(320) + tree_panel.setMinimumWidth(260) + main_layout.addWidget(tree_panel, 0) # stretch=0，不抢占右侧图片的宽度 + + # ========================================== + # 第三列：图像查看器 + # ========================================== right_panel = QWidget() right_layout = QVBoxLayout() right_layout.setContentsMargins(0, 0, 0, 0) @@ -1503,10 +1486,33 @@ class Step12VizPanel(QWidget): self.image_viewer.refresh_btn.clicked.connect(self.scan_work_directory) right_layout.addWidget(self.image_viewer, 1) right_panel.setLayout(right_layout) - main_layout.addWidget(right_panel, 1) + main_layout.addWidget(right_panel, 1) # stretch=1，右侧画板填满所有剩余宽度 self.setLayout(main_layout) + def update_image_tree_view(self): + """响应下拉框改变，重新渲染树状图""" + if not hasattr(self, 'image_tree') or not self.image_tree._all_image_files: + return + + # 1. 提取当前选中的分组模式 + mode_idx = self.view_mode_cb.currentIndex() + if mode_idx == 0: + group_mode = 'parameter' + elif mode_idx == 1: + group_mode = 'type' + else: + group_mode = 'folder' + + # 2. 提取当前的图表筛选条件 + filter_type = self.chart_filter_cb.currentText() + if filter_type == "全部图表": + filter_type = 'all' + + # 3. 触发重绘 + self.image_tree.rebuild_tree(group_mode, filter_type) + self._load_first_image_from_tree() + def set_work_dir(self, work_dir): """设置工作目录""" self.work_dir = work_dir @@ -1678,7 +1684,7 @@ class Step12VizPanel(QWidget): return if not (self.gen_scatter.isChecked() or self.gen_spectrum.isChecked() or self.gen_boxplots.isChecked() or self.gen_mask_glint.isChecked() or - self.gen_sampling_map.isChecked()): + self.gen_sampling_map.isChecked() or self.gen_distribution_map.isChecked()): QMessageBox.information(self, "提示", "请至少勾选一项可视化配置选项以生成图表。") return reply = QMessageBox.question( @@ -1694,9 +1700,20 @@ class Step12VizPanel(QWidget): "gen_boxplots": self.gen_boxplots.isChecked(), "gen_mask_glint": self.gen_mask_glint.isChecked(), "gen_sampling_map": self.gen_sampling_map.isChecked(), + "gen_distribution_map": self.gen_distribution_map.isChecked(), } main_window = self.window() factory = getattr(main_window, '_panel_factory', None) if main_window else None + + # [新增] 直接从 Step 1 面板读取原始 .shp 的绝对路径，突破工作目录限制 + step1_panel = factory.get_panel('step1_mask') if factory else None + if step1_panel: + s1_conf = step1_panel.get_config() + s1_mask = s1_conf.get('mask_path') + # 确保文件存在且是shp格式，存入extra透传给后台线程 + if s1_mask and Path(s1_mask).is_file() and str(s1_mask).lower().endswith('.shp'): + extra["boundary_shp_path"] = str(s1_mask) + step6_panel = factory.get_panel('step6_feature') if factory else None if step6_panel and getattr(step6_panel, 'output_file', None): _resolved_csv = step6_panel.output_file.get_path() @@ -1886,6 +1903,7 @@ class Step12VizPanel(QWidget): 'generate_spectrum': self.gen_spectrum.isChecked(), 'generate_glint_previews': self.gen_mask_glint.isChecked(), 'generate_sampling_maps': self.gen_sampling_map.isChecked(), + 'generate_distribution_maps': self.gen_distribution_map.isChecked(), 'scatter_config': { 'metric': 'test_r2', 'feature_start_column': 13, 'test_size': 0.2, 'random_state': 42 @@ -1913,3 +1931,5 @@ class Step12VizPanel(QWidget): self.gen_mask_glint.setChecked(config['generate_glint_previews']) if 'generate_sampling_maps' in config: self.gen_sampling_map.setChecked(config.get('generate_sampling_maps', True)) + if 'generate_distribution_maps' in config: + self.gen_distribution_map.setChecked(config.get('generate_distribution_maps', True)) diff --git a/src/new/main_view.py b/src/new/main_view.py index 6d5766e..35e38ba 100644 --- a/src/new/main_view.py +++ b/src/new/main_view.py @@ -166,7 +166,7 @@ ROUTES = [ { # data/icons/ 没有 12.png/13.png；Step11/12/13 暂时复用 9.png（11 个 png 对 13 个 step 必然有共用） "id": "step11", - "name": "11. 专题图生成", + "name": "11. 分布图生成", "icon": "9.png", "view_module": "src.new.views.step11_view", "view_class": "Step11View", diff --git a/src/new/services/step11_service.py b/src/new/services/step11_service.py index 1c8af7b..e906033 100644 --- a/src/new/services/step11_service.py +++ b/src/new/services/step11_service.py @@ -25,7 +25,8 @@ Step11 后端计算服务（专题图生成 / 克里金插值） "boundary_shp_path": "D:/boundary.shp", # 边界 shp（可选） "resolution": 30.0, # 空间分辨率（米） "input_crs": "EPSG:32651", - "output_crs": "EPSG:4326", + # ★★★ 强制默认 output_crs = input_crs，禁止从 service 配置误改为 EPSG:4326 ★★★ + "output_crs": "EPSG:32651", "output_dir": "D:/11_Thematic_Map", # 输出目录 "enabled": True, "work_dir": "D:/workspace", # 工作目录 diff --git a/src/new/services/step12_service.py b/src/new/services/step12_service.py index fb7f188..98ed3a8 100644 --- a/src/new/services/step12_service.py +++ b/src/new/services/step12_service.py @@ -46,7 +46,7 @@ Step12 后端计算服务（数据可视化——散点/光谱/箱线/掩膜缩 from __future__ import annotations from pathlib import Path -from typing import Any, Dict +from typing import Any, Dict, List, Optional from src.new.services._output_resolver import get_user_output_path, is_user_specified, resolve_output_dir @@ -168,18 +168,115 @@ def _try_sampling_maps(work_dir: str, output_dir: Path) -> Dict[str, Any]: return {"status": "completed" if p else "error", "path": p, "output_dir": str(output_dir / "sampling_maps")} +def _resolve_thematic_map_dir(work_dir: str) -> "Optional[Path]": + """自动推断 Step 11 分布图输出目录 + + 与 step11_service._resolve_output_dir 的 default 分支完全镜像—— + 用户未指定 output_dir 时默认为 ``work_dir/11_Thematic_Map``。 + """ + if not work_dir: + return None + cand = Path(work_dir) / "11_Thematic_Map" + return cand if cand.is_dir() else None + + +def _try_distribution_maps(work_dir: str, output_dir: Path) -> Dict[str, Any]: + """渲染/汇集 Step 11 生成的分布图到 14_visualization/distribution_maps + + 处理两类产物（独立子步骤，单个失败不影响另一个）： + + 1. **PNG 直拷** —— ``11_Thematic_Map/*_专题图.png`` + （Step 11 GeoTIFF 栅格模式的产物，已带外框/图例/采样点，直接复制即可） + 2. **GeoTIFF 渲染** —— ``11_Thematic_Map/*_kriging.tif`` + （Step 11 CSV 插值模式的纯 GeoTIFF 产物，调用 ``ContentMapper.visualize_raster`` + 渲染为同风格的 PNG） + + Returns: + {"status": "completed"|"skipped", "count": int, + "png_copied": int, "tif_rendered": int, + "output_dir": str, "details": [{"src", "dst", "kind"}, ...]} + """ + import shutil + + src_dir = _resolve_thematic_map_dir(work_dir) + if src_dir is None: + raise FileNotFoundError( + f"Step 11 输出目录不存在: {work_dir}/11_Thematic_Map" + f"（请先运行 Step 11 生成分布图）" + ) + + out_sub = output_dir / "distribution_maps" + out_sub.mkdir(parents=True, exist_ok=True) + + details: List[Dict[str, str]] = [] + n_png_copy = 0 + n_tif_render = 0 + + # --- 1. PNG 直拷（GeoTIFF 模式已有产物） --- + for src_png in sorted(src_dir.glob("*_专题图.png")): + dst_png = out_sub / src_png.name + try: + shutil.copy2(src_png, dst_png) + details.append({"src": str(src_png), "dst": str(dst_png), "kind": "png_copy"}) + n_png_copy += 1 + except Exception as copy_err: # noqa: BLE001 + print(f"[distribution_maps] ⚠ 复制失败 {src_png.name}: {copy_err}") + + # --- 2. GeoTIFF 渲染（CSV 模式纯栅格产物） --- + tif_paths = sorted(src_dir.glob("*_kriging.tif")) + if tif_paths: + import matplotlib + matplotlib.use("Agg") + from src.postprocessing.map import ContentMapper + + mapper = ContentMapper() + for tif_path in tif_paths: + stem = tif_path.stem + chinese_title = mapper._get_chinese_title(stem) + dst_png = out_sub / f"{chinese_title}_分布图.png" + try: + mapper.visualize_raster( + raster_tif_path=str(tif_path), + output_file=str(dst_png), + boundary_shp_path=None, + nodata_value=-9999.0, + figsize=(14, 10), + alpha=0.9, + ) + details.append({"src": str(tif_path), "dst": str(dst_png), "kind": "tif_render"}) + n_tif_render += 1 + except Exception as render_err: # noqa: BLE001 + print(f"[distribution_maps] ⚠ 渲染失败 {tif_path.name}: {render_err}") + + total = n_png_copy + n_tif_render + if total == 0: + raise FileNotFoundError( + f"Step 11 输出目录 {src_dir} 中无 *_专题图.png 也无 *_kriging.tif" + ) + + return { + "status": "completed", + "count": total, + "png_copied": n_png_copy, + "tif_rendered": n_tif_render, + "output_dir": str(out_sub), + "details": details, + } + + def execute_step12(config: Dict[str, Any]) -> Dict[str, Any]: """Step 12 后端计算入口——纯函数""" work_dir: str = config.get("work_dir") or "" img_dir: str = config.get("img_dir") or "" enabled: bool = bool(config.get("enabled", True)) output_dir: str = config.get("output_dir") or "" - # 5 个开关：缺省默认 True（与旧 panel 行为一致） + # 6 个开关：缺省默认 True（与旧 panel 行为一致；新增 distribution_maps 用于闭环 Step 11） gen_scatter = bool(config.get("generate_scatter", True)) gen_spectrum = bool(config.get("generate_spectrum", True)) gen_boxplots = bool(config.get("generate_boxplots", True)) gen_glint = bool(config.get("generate_glint_previews", True)) gen_sampling = bool(config.get("generate_sampling_maps", True)) + gen_distribution = bool(config.get("generate_distribution_maps", True)) output_path, _source = _resolve_output_dir(config, work_dir) mode = "viz_generate" @@ -228,12 +325,14 @@ def execute_step12(config: Dict[str, Any]) -> Dict[str, Any]: tasks.append(("glint_previews", _try_glint_previews)) if gen_sampling: tasks.append(("sampling_maps", _try_sampling_maps)) + if gen_distribution: + tasks.append(("distribution_maps", _try_distribution_maps)) if not tasks: return { "status": "completed", "output_path": str(output_path).replace("\\", "/"), - "message": "无可视化任务（5 个开关全部 False）", + "message": "无可视化任务（6 个开关全部 False）", "mode": mode, } diff --git a/src/new/views/step11_view.py b/src/new/views/step11_view.py index 3299cfc..b035fcc 100644 --- a/src/new/views/step11_view.py +++ b/src/new/views/step11_view.py @@ -181,7 +181,10 @@ class Step11View(BaseView): params_layout.addRow("输入坐标系:", self.input_crs) self.output_crs = QLineEdit() - self.output_crs.setText("EPSG:4326") + # ★★★ 强制默认输出坐标系与输入一致，禁止从 GUI 误改为 EPSG:4326 ★★★ + # 历史默认值 'EPSG:4326' 会让 ContentMapper 把栅格重投影到经纬度， + # 与基于 EPSG:32651 的水域掩膜叠加时发生仿射变换撕裂（栅格错位、坐标轴扭曲）。 + self.output_crs.setText("EPSG:32651") params_layout.addRow("输出坐标系:", self.output_crs) self.show_points = QCheckBox("显示采样点") @@ -302,7 +305,8 @@ class Step11View(BaseView): "boundary_shp_path": self.boundary_file.get_path(), "resolution": self.resolution.value(), "input_crs": self.input_crs.text(), - "output_crs": self.output_crs.text(), + # ★★★ 强制 output_crs = input_crs，不再信任 GUI 的 output_crs 输入框 ★★★ + "output_crs": self.input_crs.text(), "show_sample_points": self.show_points.isChecked(), "use_distance_diffusion": self.use_diffusion.isChecked(), "enabled": self.enable_checkbox.isChecked(), @@ -342,8 +346,9 @@ class Step11View(BaseView): self.resolution.setValue(config["resolution"]) if "input_crs" in config: self.input_crs.setText(config["input_crs"]) - if "output_crs" in config: - self.output_crs.setText(config["output_crs"]) + # ★★★ 反灌入时强制 output_crs = input_crs，避免旧 config 中的 EPSG:4326 回填 ★★★ + if "output_crs" in config or "input_crs" in config: + self.output_crs.setText(config.get("input_crs") or config.get("output_crs") or "EPSG:32651") if "show_sample_points" in config: self.show_points.setChecked(config["show_sample_points"]) if "use_distance_diffusion" in config: diff --git a/src/new/views/step12_view.py b/src/new/views/step12_view.py index 2e96f52..a653509 100644 --- a/src/new/views/step12_view.py +++ b/src/new/views/step12_view.py @@ -100,6 +100,14 @@ class Step12View(BaseView): self.gen_sampling_map.setChecked(True) config_layout.addWidget(self.gen_sampling_map) + self.gen_distribution_map = QCheckBox("空间分布图（Step 11 产物）") + self.gen_distribution_map.setChecked(True) + self.gen_distribution_map.setToolTip( + "渲染/汇集 Step 11 生成的分布图（PNG 直拷 + GeoTIFF 渲染）" + "到 14_visualization/distribution_maps" + ) + config_layout.addWidget(self.gen_distribution_map) + config_layout.addSpacing(10) line = QFrame() line.setFrameShape(QFrame.HLine) @@ -186,6 +194,7 @@ class Step12View(BaseView): "generate_boxplots": self.gen_boxplots.isChecked(), "generate_glint_previews": self.gen_mask_glint.isChecked(), "generate_sampling_maps": self.gen_sampling_map.isChecked(), + "generate_distribution_maps": self.gen_distribution_map.isChecked(), "enabled": True, } @@ -207,6 +216,8 @@ class Step12View(BaseView): self.gen_mask_glint.setChecked(config["generate_glint_previews"]) if "generate_sampling_maps" in config: self.gen_sampling_map.setChecked(config.get("generate_sampling_maps", True)) + if "generate_distribution_maps" in config: + self.gen_distribution_map.setChecked(config.get("generate_distribution_maps", True)) def update_work_directory(self, work_dir: str): super().update_work_directory(work_dir) diff --git a/src/postprocessing/map.py b/src/postprocessing/map.py index 458b12c..9c8576e 100644 --- a/src/postprocessing/map.py +++ b/src/postprocessing/map.py @@ -838,6 +838,14 @@ class ContentMapper: transform = src.transform crs = src.crs + # 【新增防御】检测 Transform 是否为纯像素矩阵（极易导致严重错位） + if transform.is_identity: + print("\n" + "!" * 65) + print(f"⚠️ [严重警告] 栅格掩膜 {Path(raster_path).name} 缺少真实的地理仿射变换！") + print(f"当前被判定为纯像素坐标系 (x: 0~Width, y: 0~Height)。\n强行与地理栅格 (UTM/WGS84) 叠加将发生【极其严重的错位】！") + print(f"请务必在 Step 1 导入原始的 .shp 矢量文件进行约束。") + print("!" * 65 + "\n") + # 二值化：>0 视为水 mask_uint8 = (data > 0).astype(np.uint8) if int(mask_uint8.sum()) == 0: @@ -895,56 +903,72 @@ class ContentMapper: print(f"识别边缘点时出错: {e}，将使用所有点作为边缘点") return np.arange(len(points)) - def _expand_edge_points(self, points_gdf, boundary_gdf, resolution=100, expand_ratio=0.05): + def _expand_edge_points(self, points_gdf, boundary_gdf=None, resolution=100, expand_ratio=0.05): """ 对边缘采样点进行外扩处理，外扩到整个图像的边界（包括外扩后的边界） - 按照指定的间距（resolution）生成外扩点，铺满整个画面 - + 按照指定的间距（resolution）生成外扩点，铺满整个画面。 + + ★★★ Plan C：boundary_gdf 可选（None = 不依赖水域掩膜，纯采样点自然扩展）★★★ + Parameters: ----------- points_gdf : gpd.GeoDataFrame 原始采样点GeoDataFrame - boundary_gdf : gpd.GeoDataFrame - 水域掩膜边界GeoDataFrame + boundary_gdf : gpd.GeoDataFrame, optional + 水域掩膜边界GeoDataFrame。None 时基于采样点自身范围做外扩 resolution : float, default=100 外扩点的间距（单位与坐标相同），与插值网格分辨率一致 expand_ratio : float, default=0.05 边界外扩比例（与create_interpolation_grid中的expand_ratio一致） - + Returns: -------- expanded_gdf : gpd.GeoDataFrame 外扩后的采样点GeoDataFrame """ + # ── Plan C: 无水域掩膜时，基于采样点自身范围外扩 ────────────── + if boundary_gdf is None: + print(f"[Plan C] 无水域掩膜，基于采样点范围做自然外扩（expand_ratio={expand_ratio}）...") + points = np.column_stack( + (points_gdf['proj_x'].values, points_gdf['proj_y'].values) + ) + p_minx, p_miny = points.min(axis=0) + p_maxx, p_maxy = points.max(axis=0) + width = p_maxx - p_minx + height = p_maxy - p_miny + expand_x = width * expand_ratio + expand_y = height * expand_ratio + image_minx = p_minx - expand_x + image_maxx = p_maxx + expand_x + image_miny = p_miny - expand_y + image_maxy = p_maxy + expand_y + print(f"[Plan C] 采样点范围: X[{p_minx:.2f}, {p_maxx:.2f}], Y[{p_miny:.2f}, {p_maxy:.2f}]") + print(f"[Plan C] 外扩后范围: X[{image_minx:.2f}, {image_maxx:.2f}], Y[{image_miny:.2f}, {image_maxy:.2f}]") + else: + boundary_bounds = boundary_gdf.total_bounds + mask_minx, mask_miny, mask_maxx, mask_maxy = boundary_bounds + width = mask_maxx - mask_minx + height = mask_maxy - mask_miny + expand_x = width * expand_ratio + expand_y = height * expand_ratio + image_minx = mask_minx - expand_x + image_maxx = mask_maxx + expand_x + image_miny = mask_miny - expand_y + image_maxy = mask_maxy + expand_y + print(f"正在对边缘采样点进行外扩处理（按照 {resolution} 的间距外扩到整个图像边界）...") - + # 识别边缘点 edge_indices = self._identify_edge_points(points_gdf) - + if len(edge_indices) == 0: print("未识别到边缘点，跳过外扩处理") return points_gdf.copy() - # 获取水域掩膜的边界范围 - boundary_bounds = boundary_gdf.total_bounds # [minx, miny, maxx, maxy] - mask_minx, mask_miny, mask_maxx, mask_maxy = boundary_bounds - - # 计算范围大小 - width = mask_maxx - mask_minx - height = mask_maxy - mask_miny - - # 外扩边界，与create_interpolation_grid中的逻辑一致，确保外扩到整个图像范围 - expand_x = width * expand_ratio - expand_y = height * expand_ratio - image_minx = mask_minx - expand_x - image_maxx = mask_maxx + expand_x - image_miny = mask_miny - expand_y - image_maxy = mask_maxy + expand_y - # 获取所有点的坐标 points = np.column_stack((points_gdf['proj_x'].values, points_gdf['proj_y'].values)) - - # 计算点集的范围和中心 + + # 计算点集的范围和中心（Plan C 两个分支都要用，提前算） x_min, x_max = points[:, 0].min(), points[:, 0].max() y_min, y_max = points[:, 1].min(), points[:, 1].max() center = np.array([(x_min + x_max) / 2, (y_min + y_max) / 2]) @@ -1055,7 +1079,8 @@ class ContentMapper: print(f"外扩完成：原始点 {len(points_gdf)} 个，边缘点 {len(edge_indices)} 个，" f"新增外扩点 {len(new_data_list)} 个（间距 {resolution}），总计 {len(result_gdf)} 个点") - print(f"水域掩膜范围: X[{mask_minx:.2f}, {mask_maxx:.2f}], Y[{mask_miny:.2f}, {mask_maxy:.2f}]") + if boundary_gdf is not None: + print(f"水域掩膜范围: X[{mask_minx:.2f}, {mask_maxx:.2f}], Y[{mask_miny:.2f}, {mask_maxy:.2f}]") print(f"图像范围（含外扩）: X[{image_minx:.2f}, {image_maxx:.2f}], Y[{image_miny:.2f}, {image_maxy:.2f}]") return result_gdf @@ -1063,42 +1088,52 @@ class ContentMapper: print("未生成外扩点，返回原始点集") return points_gdf.copy() - def create_interpolation_grid(self, points_gdf, boundary_gdf, resolution=100, expand_ratio=0.05, + def create_interpolation_grid(self, points_gdf, boundary_gdf=None, resolution=100, expand_ratio=0.05, use_distance_diffusion=True, max_diffusion_distance=None, diffusion_power=2, diffusion_n_neighbors=15): """ 创建插值网格 + ★★★ Plan C：boundary_gdf 可选（None = 纯采样点自然插值，无水域掩膜约束）★★★ + Parameters: ----------- + boundary_gdf : gpd.GeoDataFrame, optional + 水域掩膜边界GeoDataFrame。None 时基于采样点自身范围插值，不做掩膜裁剪和填充。 expand_ratio : float, default=0.05 - 边界外扩比例（5%），确保图像边界不完全挨着地图 + 边界外扩比例（5%），用于从采样点范围外扩出图像边界。 use_distance_diffusion : bool, default=True - 是否使用距离扩散方法填充边界空白区域 + 是否使用距离扩散方法填充边界空白区域（仅在 boundary_gdf 有值时生效）。 max_diffusion_distance : float, optional - 最大扩散距离（单位与坐标相同）。如果为None，自动计算为网格分辨率的5倍 + 最大扩散距离（单位与坐标相同）。如果为None，自动计算为网格分辨率的5倍。 diffusion_power : float, default=2 - 距离扩散的IDW幂参数，值越大，距离衰减越快 + 距离扩散的IDW幂参数，值越大，距离衰减越快。 diffusion_n_neighbors : int, default=15 - 距离扩散使用的最近邻点数 - + 距离扩散使用的最近邻点数。 + Returns: -------- grid_xx, grid_yy, grid_content, bounds : tuple """ print("正在创建插值网格...") - # 获取边界范围 - bounds = boundary_gdf.total_bounds - minx, miny, maxx, maxy = bounds + # ── Plan C: 无水域掩膜时，基于采样点自身范围计算边界 ───────── + if boundary_gdf is None: + print("[Plan C] 无水域掩膜，基于采样点范围创建插值网格...") + points = np.column_stack((points_gdf['proj_x'], points_gdf['proj_y'])) + minx = points[:, 0].min() + maxx = points[:, 0].max() + miny = points[:, 1].min() + maxy = points[:, 1].max() + print(f"采样点范围: X({minx:.6f} - {maxx:.6f}), Y({miny:.6f} - {maxy:.6f})") + else: + bounds = boundary_gdf.total_bounds + minx, miny, maxx, maxy = bounds + print(f"水域掩膜范围: X({minx:.6f} - {maxx:.6f}), Y({miny:.6f} - {maxy:.6f})") - print(f"原始边界范围: X({minx:.6f} - {maxx:.6f}), Y({miny:.6f} - {maxy:.6f})") - - # 计算范围大小 + # 计算范围大小并外扩 width = maxx - minx height = maxy - miny - - # 外扩边界，确保图像不完全挨着地图 expand_x = width * expand_ratio expand_y = height * expand_ratio minx -= expand_x @@ -1111,84 +1146,77 @@ class ContentMapper: if self.output_crs == 'EPSG:4326': print(f"区域尺寸: 宽度={width:.6f}°, 高度={height:.6f}°") - # 对于地理坐标系，需要调整分辨率单位（度） - # 1度约等于111公里，所以100米约等于0.0009度 - resolution_deg = resolution / 111000.0 # 将米转换为度 + resolution_deg = resolution / 111000.0 print(f"网格分辨率: {resolution}m ≈ {resolution_deg:.6f}°") else: print(f"区域尺寸: 宽度={width:.2f}m, 高度={height:.2f}m") resolution_deg = resolution - # 检查分辨率是否合理 - min_grid_points = 50 # 增加最少网格点数以获得更平滑的插值效果 - + # 计算网格点数 + min_grid_points = 50 if self.output_crs == 'EPSG:4326': - # 地理坐标系的网格点计算 grid_points_x = max(int(width / resolution_deg), min_grid_points) grid_points_y = max(int(height / resolution_deg), min_grid_points) else: - # 投影坐标系的网格点计算 grid_points_x = max(int(width / resolution), min_grid_points) grid_points_y = max(int(height / resolution), min_grid_points) - # 确保网格足够密集以获得平滑效果 grid_points_x = max(grid_points_x, 100) grid_points_y = max(grid_points_y, 100) - # 创建网格 grid_x = np.linspace(minx, maxx, grid_points_x) grid_y = np.linspace(miny, maxy, grid_points_y) grid_xx, grid_yy = np.meshgrid(grid_x, grid_y) print(f"网格大小: {grid_xx.shape[1]} x {grid_xx.shape[0]} (宽 x 高)") - # 检查网格大小 if grid_xx.shape[0] < 2 or grid_xx.shape[1] < 2: - raise ValueError(f"网格尺寸太小 {grid_xx.shape}，无法进行插值。请检查数据范围和分辨率设置。") + raise ValueError(f"网格尺寸太小 {grid_xx.shape}，无法进行插值。") - # 准备插值数据（使用原始点+外扩点的合并数据） + # 准备插值数据 points = np.column_stack((points_gdf['proj_x'], points_gdf['proj_y'])) values = points_gdf['content'].values - print(f"插值数据点数量: {len(points)}（包含原始采样点和外扩点）") - print(f"数据点范围: X({points[:, 0].min():.6f} - {points[:, 0].max():.6f}), " - f"Y({points[:, 1].min():.6f} - {points[:, 1].max():.6f})") + print(f"插值数据点数量: {len(points)}") print(f"含量值范围: {values.min():.4f} - {values.max():.4f}") - print(f"含量值统计: 平均={values.mean():.4f}, 标准差={values.std():.4f}") - # 检查数据点数量 if len(points) < 3: raise ValueError("插值需要至少3个数据点") - # 检查数据点的几何分布 self._check_point_distribution(points) - # 执行插值（先对整个网格插值，包括边界外） - print("正在执行空间插值（整个网格，包括边界外）...") + # 执行插值 + print("正在执行空间插值...") grid_content = self._perform_interpolation(points, values, grid_xx, grid_yy) - # 创建边界掩膜（用于识别边界内外） + # ── Plan C: 无水域掩膜时，跳过所有掩膜裁剪和边缘填充逻辑 ────── + if boundary_gdf is None: + print("[Plan C] 无水域掩膜，跳过边缘填充，保留插值空白区域（NaN）") + valid_data = ~np.isnan(grid_content) + valid_count = np.sum(valid_data) + print(f"有效插值点数量: {valid_count} / {grid_content.size}") + if valid_count > 0: + valid_values = grid_content[valid_data] + print(f"插值后数据统计: 最小值={valid_values.min():.4f}, " + f"最大值={valid_values.max():.4f}, 平均值={valid_values.mean():.4f}") + expanded_bounds = np.array([minx, miny, maxx, maxy]) + return grid_xx, grid_yy, grid_content, expanded_bounds + + # ── 以下为原有水域掩膜逻辑（boundary_gdf 有值时执行）──────────── print("正在识别边界区域...") - # 创建掩膜 mask_points = np.column_stack((grid_xx.ravel(), grid_yy.ravel())) mask_geometry = [Point(x, y) for x, y in mask_points] mask_gdf = gpd.GeoDataFrame(geometry=mask_geometry, crs=self.output_crs) - - # 检查哪些点在边界内 within_boundary = mask_gdf.within(boundary_gdf.unary_union) mask = within_boundary.values.reshape(grid_xx.shape) - # 找到边界边缘上的点（在边界内，但靠近边界） print("正在提取边界边缘值并填充边界外区域...") - # 方法：找到边界内有效值的边缘点，然后填充到边界外 - # 1. 先填充边界内的NaN（使用距离扩散方法） nan_mask = np.isnan(grid_content) within_boundary_nan = nan_mask & mask if np.any(within_boundary_nan): if use_distance_diffusion: - # 使用距离扩散方法填充边界内的空白区域 grid_content = self._fill_boundary_blanks_with_distance_diffusion( grid_content, grid_xx, grid_yy, mask, boundary_gdf, max_diffusion_distance=max_diffusion_distance, @@ -1196,72 +1224,49 @@ class ContentMapper: n_neighbors=diffusion_n_neighbors ) else: - # 使用传统的最近邻插值方法 print(f"填充边界内的 {np.sum(within_boundary_nan)} 个NaN点（使用最近邻插值）...") valid_mask = ~nan_mask & mask if np.sum(valid_mask) > 0: valid_points = np.column_stack((grid_xx[valid_mask], grid_yy[valid_mask])) valid_values = grid_content[valid_mask] nan_points = np.column_stack((grid_xx[within_boundary_nan], grid_yy[within_boundary_nan])) - - filled_values = griddata( - valid_points, valid_values, nan_points, - method='nearest' - ) + filled_values = griddata(valid_points, valid_values, nan_points, method='nearest') grid_content[within_boundary_nan] = filled_values - print(f"边界内填充完成") + print("边界内填充完成") - # 2. 找到边界边缘的值（边界内但靠近边界外的点） - # 使用形态学操作找到边界边缘 boundary_mask_binary = mask.astype(int) - # 创建边界外掩膜 outside_mask = ~mask - - # 找到边界边缘（在边界内，但相邻有边界外的点） - # 对边界外区域进行膨胀，找到边界边缘 kernel = np.ones((3, 3), dtype=bool) dilated_outside = ndimage.binary_dilation(outside_mask, structure=kernel) - edge_mask = mask & dilated_outside # 边界内但靠近边界外的点 + edge_mask = mask & dilated_outside - # 3. 提取边缘值，填充到边界外 if np.any(edge_mask): edge_values = grid_content[edge_mask] edge_valid = ~np.isnan(edge_values) if np.any(edge_valid): - # 使用边缘的有效值填充边界外 edge_mean = np.nanmean(edge_values) print(f"边界边缘平均值: {edge_mean:.4f}") - - # 将边缘值填充到边界外的所有NaN点 outside_nan = outside_mask & np.isnan(grid_content) if np.any(outside_nan): - # 使用最近邻插值从边缘值填充 edge_points = np.column_stack((grid_xx[edge_mask & ~np.isnan(grid_content)], grid_yy[edge_mask & ~np.isnan(grid_content)])) if len(edge_points) > 0: edge_vals = grid_content[edge_mask & ~np.isnan(grid_content)] outside_points = np.column_stack((grid_xx[outside_nan], grid_yy[outside_nan])) - - outside_filled = griddata( - edge_points, edge_vals, outside_points, - method='nearest' - ) + outside_filled = griddata(edge_points, edge_vals, outside_points, method='nearest') grid_content[outside_nan] = outside_filled print(f"已填充边界外的 {np.sum(~np.isnan(outside_filled))} 个点") else: - # 如果没有边缘值，使用边缘平均值填充 grid_content[outside_nan] = edge_mean print(f"使用边缘平均值填充边界外的 {np.sum(outside_nan)} 个点") else: print("边界外区域已全部填充") else: - # 如果边缘没有有效值，使用全局平均值填充边界外 global_mean = np.nanmean(grid_content[mask]) if not np.isnan(global_mean): grid_content[outside_mask & np.isnan(grid_content)] = global_mean print(f"使用全局平均值 {global_mean:.4f} 填充边界外") else: - # 如果没有找到边缘，直接使用边界内的平均值填充边界外 mean_in_boundary = np.nanmean(grid_content[mask]) if not np.isnan(mean_in_boundary): grid_content[outside_mask & np.isnan(grid_content)] = mean_in_boundary @@ -1269,26 +1274,16 @@ class ContentMapper: print("整个画面已铺满，边界外区域已用边缘值填充") - # 最终检查：确保边界内所有区域都有值 final_check_nan = np.isnan(grid_content) & mask if np.any(final_check_nan): - print(f"警告: 仍有 {np.sum(final_check_nan)} 个边界内的点未填充，使用平均值填充...") + print(f"警告: 仍有 {np.sum(final_check_nan)} 个边界内的点未填充...") if np.sum(~np.isnan(grid_content) & mask) > 0: mean_value = np.nanmean(grid_content[mask]) grid_content[final_check_nan] = mean_value - print(f" 使用平均值 {mean_value:.4f} 填充剩余 {np.sum(final_check_nan)} 个点") else: - # 如果边界内完全没有有效值，使用全局平均值 global_mean = np.nanmean(grid_content) - if not np.isnan(global_mean): - grid_content[final_check_nan] = global_mean - else: - grid_content[final_check_nan] = 0 - print(" 使用全局平均值填充") - else: - print("边界内所有区域已完全填充") + grid_content[final_check_nan] = global_mean if not np.isnan(global_mean) else 0 - # 检查插值结果 valid_data = ~np.isnan(grid_content) valid_count = np.sum(valid_data) print(f"有效插值点数量: {valid_count} / {grid_content.size}") @@ -1299,23 +1294,22 @@ class ContentMapper: if valid_count < 4: print("警告：有效数据点很少，可能影响绘图效果") - # 输出插值结果的统计信息 valid_values = grid_content[valid_data] - print( - f"插值后数据统计: 最小值={valid_values.min():.4f}, 最大值={valid_values.max():.4f}, 平均值={valid_values.mean():.4f}") + print(f"插值后数据统计: 最小值={valid_values.min():.4f}, " + f"最大值={valid_values.max():.4f}, 平均值={valid_values.mean():.4f}") - # 返回外扩后的bounds expanded_bounds = np.array([minx, miny, maxx, maxy]) - return grid_xx, grid_yy, grid_content, expanded_bounds - def create_content_map(self, points_gdf, boundary_gdf, grid_xx, grid_yy, - grid_content, bounds, output_file='content_map.png', + def create_content_map(self, points_gdf, boundary_gdf=None, grid_xx=None, grid_yy=None, + grid_content=None, bounds=None, output_file='content_map.png', show_sample_points=False, base_map_tif=None, cmap='viridis'): """ 创建含量图 + ★★★ Plan C：boundary_gdf 可选（None = 无掩膜裁剪，无黑色边界线）★★★ + Parameters: ----------- base_map_tif : str, optional @@ -1330,19 +1324,22 @@ class ContentMapper: # 创建边界掩膜（用于绘图时只显示边界内） print("创建边界掩膜用于绘图...") - try: - # 创建网格点的GeoDataFrame - grid_points = gpd.GeoDataFrame( - geometry=[Point(x, y) for x, y in zip(grid_xx.flatten(), grid_yy.flatten())], - crs=points_gdf.crs - ) - # 检查哪些点在边界内 - within_boundary = grid_points.within(boundary_gdf.unary_union) - mask = within_boundary.values.reshape(grid_xx.shape) - print(f"边界内点数: {np.sum(mask)} / {mask.size}") - except Exception as e: - print(f"创建边界掩膜时出现错误: {e}，继续绘图...") - mask = np.ones_like(grid_content, dtype=bool) # 如果失败，显示全部 + # ── Plan C: 无水域掩膜时，显示全部插值区域（NaN 区域本身就不填充）───────── + if boundary_gdf is None: + print("[Plan C] 无水域掩膜，显示全部插值区域，不做边界裁剪") + mask = np.ones_like(grid_content, dtype=bool) + else: + try: + grid_points = gpd.GeoDataFrame( + geometry=[Point(x, y) for x, y in zip(grid_xx.flatten(), grid_yy.flatten())], + crs=points_gdf.crs + ) + within_boundary = grid_points.within(boundary_gdf.unary_union) + mask = within_boundary.values.reshape(grid_xx.shape) + print(f"边界内点数: {np.sum(mask)} / {mask.size}") + except Exception as e: + print(f"创建边界掩膜时出现错误: {e}，继续绘图...") + mask = np.ones_like(grid_content, dtype=bool) # 如果失败，显示全部 valid_data = ~np.isnan(grid_content) if np.sum(valid_data) == 0: @@ -1495,12 +1492,13 @@ class ContentMapper: print(f"所有绘图方法都失败: {e3}") raise ValueError("无法生成颜色图，请检查数据") - # 绘制边界（黑色） - try: - boundary_gdf.boundary.plot(ax=ax, color='black', linewidth=2, alpha=1.0) - print("边界绘制成功（黑色）") - except Exception as e: - print(f"边界绘制失败: {e}") + # 绘制边界（黑色）—— Plan C: 无掩膜时不绘制边界线 + if boundary_gdf is not None: + try: + boundary_gdf.boundary.plot(ax=ax, color='black', linewidth=2, alpha=1.0) + print("边界绘制成功（黑色）") + except Exception as e: + print(f"边界绘制失败: {e}") # 可选择性绘制采样点（默认不绘制，以显示平滑的颜色分布） if show_sample_points: @@ -1606,19 +1604,20 @@ class ContentMapper: # ★★★ 改用画布相对坐标（transAxes）★★★ # (0.88, 0.92) = 右上角，尺寸用 points（72分之一英寸） arrow_ax_x, arrow_ax_y = 0.88, 0.92 - radius_pt = 28 # 罗盘半径（磅），固定大小 + radius_pt = 18 # 罗盘半径（磅），由 28 → 18 缩小图元 # 统一在数据坐标系下绘制（transform=ax.transData） # 但 position 由 axes 坐标决定，radius 用固定点数 # 将 axes 坐标转为数据坐标：取右上角 + 偏移 xlim = ax.get_xlim() ylim = ax.get_ylim() - dx = (xlim[1] - xlim[0]) * 0.08 - dy = (ylim[1] - ylim[0]) * 0.08 + # 偏移系数 0.08 → 0.05 让指北针更靠中心，避免紧贴角落被裁切 + dx = (xlim[1] - xlim[0]) * 0.05 + dy = (ylim[1] - ylim[0]) * 0.05 arrow_x = xlim[1] - dx arrow_y = ylim[1] - dy - # radius 转为数据坐标单位（近似） - radius = min(dx, dy) * 0.6 + # radius 系数 0.6 → 0.42 缩小指北针整体半径 + radius = min(dx, dy) * 0.42 # 绘制圆形背景（外圈） circle_outer = patches.Circle( @@ -1626,7 +1625,7 @@ class ContentMapper: radius=radius, facecolor='white', edgecolor='black', - linewidth=2.5, + linewidth=1.5, zorder=10, transform=ax.transData, ) @@ -1638,7 +1637,7 @@ class ContentMapper: radius=radius * 0.7, facecolor='none', edgecolor='gray', - linewidth=1.5, + linewidth=0.8, linestyle='--', zorder=11, transform=ax.transData, @@ -1646,7 +1645,7 @@ class ContentMapper: ax.add_patch(circle_inner) # 绘制四个方向的刻度线 - tick_width = 1.5 + tick_width = 1.0 # 北方向刻度（主刻度） ax.plot([arrow_x, arrow_x], [arrow_y, arrow_y + radius * 0.85], @@ -1683,7 +1682,7 @@ class ContentMapper: arrow_points, facecolor='black', edgecolor='black', - linewidth=2, + linewidth=1.2, zorder=13, transform=ax.transData, ) @@ -1700,30 +1699,31 @@ class ContentMapper: south_arrow_points, facecolor='white', edgecolor='black', - linewidth=1.5, + linewidth=1.0, zorder=13, transform=ax.transData, ) ax.add_patch(south_arrow_poly) # 添加方向标记（N, S, E, W） - label_offset = radius * 1.15 - font_size = 9 + label_offset = radius * 1.1 + # 字号 9 → 7（与缩小的指北针半径相匹配） + font_size = 7 ax.text(arrow_x, arrow_y + label_offset, 'N', fontsize=font_size, fontweight='bold', ha='center', va='bottom', color='black', zorder=14) ax.text(arrow_x, arrow_y - label_offset, 'S', - fontsize=font_size * 0.8, fontweight='bold', ha='center', va='top', + fontsize=font_size, fontweight='bold', ha='center', va='top', color='black', zorder=14) ax.text(arrow_x + label_offset, arrow_y, 'E', - fontsize=font_size * 0.8, fontweight='bold', ha='left', va='center', + fontsize=font_size, fontweight='bold', ha='left', va='center', color='black', zorder=14) ax.text(arrow_x - label_offset, arrow_y, 'W', - fontsize=font_size * 0.8, fontweight='bold', ha='right', va='center', + fontsize=font_size, fontweight='bold', ha='right', va='center', color='black', zorder=14) def add_scale_bar(self, ax, scale_x=None, scale_y=None): @@ -1748,7 +1748,7 @@ class ContentMapper: location='lower left', box_alpha=0.8, color='black', - font_properties={'size': 10}, + font_properties={'size': 8}, label_loc='bottom', ) ax.add_artist(scalebar) @@ -1760,7 +1760,7 @@ class ContentMapper: location='lower left', box_alpha=0.8, color='black', - font_properties={'size': 10}, + font_properties={'size': 8}, label_loc='bottom' ) ax.add_artist(scalebar) @@ -1768,7 +1768,7 @@ class ContentMapper: else: scalebar = ScaleBar(1, units='m', location='lower left', box_alpha=0.8, color='black', - font_properties={'size': 10}) + font_properties={'size': 8}) ax.add_artist(scalebar) print("投影坐标系比例尺添加成功") except Exception as e: @@ -1801,16 +1801,16 @@ class ContentMapper: scale_length_deg = distance_km / 111.0 # 转换为度数 # 绘制比例尺线 - ax.plot([scale_x, scale_x + scale_length_deg], [scale_y, scale_y], - 'k-', linewidth=3) - ax.plot([scale_x, scale_x], [scale_y - y_range * 0.01, scale_y + y_range * 0.01], + ax.plot([scale_x, scale_x + scale_length_deg], [scale_y, scale_y], 'k-', linewidth=2) - ax.plot([scale_x + scale_length_deg, scale_x + scale_length_deg], - [scale_y - y_range * 0.01, scale_y + y_range * 0.01], 'k-', linewidth=2) - + ax.plot([scale_x, scale_x], [scale_y - y_range * 0.01, scale_y + y_range * 0.01], + 'k-', linewidth=1.5) + ax.plot([scale_x + scale_length_deg, scale_x + scale_length_deg], + [scale_y - y_range * 0.01, scale_y + y_range * 0.01], 'k-', linewidth=1.5) + # 添加文字标注 ax.text(scale_x + scale_length_deg / 2, scale_y + y_range * 0.02, - f'{distance_km} km', ha='center', va='bottom', fontsize=10, + f'{distance_km} km', ha='center', va='bottom', fontsize=8, bbox=dict(boxstyle='round,pad=0.3', facecolor='white', alpha=0.8)) else: # 投影坐标系：使用米为单位 @@ -1824,18 +1824,18 @@ class ContentMapper: else: # 小于2km scale_length = 500 # 500m scale_text = '500 m' - + # 绘制比例尺线 - ax.plot([scale_x, scale_x + scale_length], [scale_y, scale_y], - 'k-', linewidth=3) - ax.plot([scale_x, scale_x], [scale_y - y_range * 0.01, scale_y + y_range * 0.01], + ax.plot([scale_x, scale_x + scale_length], [scale_y, scale_y], 'k-', linewidth=2) - ax.plot([scale_x + scale_length, scale_x + scale_length], - [scale_y - y_range * 0.01, scale_y + y_range * 0.01], 'k-', linewidth=2) - + ax.plot([scale_x, scale_x], [scale_y - y_range * 0.01, scale_y + y_range * 0.01], + 'k-', linewidth=1.5) + ax.plot([scale_x + scale_length, scale_x + scale_length], + [scale_y - y_range * 0.01, scale_y + y_range * 0.01], 'k-', linewidth=1.5) + # 添加文字标注 ax.text(scale_x + scale_length / 2, scale_y + y_range * 0.02, - scale_text, ha='center', va='bottom', fontsize=10, + scale_text, ha='center', va='bottom', fontsize=8, bbox=dict(boxstyle='round,pad=0.3', facecolor='white', alpha=0.8)) def _add_base_map(self, ax, base_map_tif, bounds, mask, grid_xx, grid_yy, boundary_gdf): @@ -2367,7 +2367,86 @@ class ContentMapper: # 保存原始宽高：transform 回退分支需用原始尺寸计算 extent w_orig, h_orig = w, h + # ── 全面 NoData 清洗：-9999.0 / NaN / Inf → 统一转为 np.nan ── + # 这一步确保陆地像素（无论来自掩膜还是原始 NoData）均被清除， + # 使 nanpercentile 分位数拉伸 100% 精准锁定水体内部 + array = np.where( + (array == nodata_value) | np.isnan(array) | np.isinf(array), + np.nan, + array + ) + + # ====== 新增：矢量掩膜物理擦除（必须在降采样之前，否则 array.shape 与 transform 错位）====== + # 把水域多边形外部的陆地像素物理擦除为 NaN，让下游 mean/std 统计 100% 干净（无陆地假数据污染） + # 同时保留 boundary_gdf_plotted 给末尾描边复用，避免重复读取 + 重复矢量化 + boundary_gdf_plotted: Optional[Any] = None + if boundary_shp_path and os.path.isfile(boundary_shp_path) and transform is not None: + try: + boundary_ext = Path(boundary_shp_path).suffix.lower() + if boundary_ext in ('.shp',): + # 矢量：直接读取 + boundary_gdf = gpd.read_file(boundary_shp_path) + elif boundary_ext in ('.dat', '.bsq', '.tif', '.tiff'): + # 栅格：复用 ContentMapper._raster_to_boundary_gdf 矢量化 + boundary_gdf = self._raster_to_boundary_gdf(boundary_shp_path) + else: + raise ValueError(f"不支持的边界文件格式: {boundary_ext}") + + # 兜底：如果 SHP 缺少投影文件（如 .prj 丢失），默认赋予 WGS84 (EPSG:4326) + if boundary_gdf.crs is None: + print(f"[visualize_raster] 警告: 掩膜 SHP 缺失坐标系，" + f"默认按 WGS84 (EPSG:4326) 处理 ({Path(boundary_shp_path).name})。") + boundary_gdf = boundary_gdf.set_crs(epsg=4326) + + # 兜底：如果栅格 TIF 缺失坐标系，默认也赋予 WGS84（CRS 顶部已导入，无需局部 import） + if crs_obj is None: + crs_obj = CRS.from_epsg(4326) + print(f"[visualize_raster] 警告: 栅格 TIF 缺失坐标系，" + f"默认按 WGS84 渲染 ({Path(raster_tif_path).name})。") + + # 坐标系对齐到当前栅格的 CRS（不是 self.output_crs，必须与 transform 保持一致） + if boundary_gdf.crs != crs_obj: + boundary_gdf = boundary_gdf.to_crs(crs_obj) + + # ==================================================================== + # 【核心修复：仅提取实心外部轮廓，忽略内部的耀斑孔洞和数据缺失缺口】 + # 彻底解决密密麻麻的黑点（孔洞描边）以及被异常挖空的方形区域 + from shapely.geometry import Polygon, MultiPolygon + exteriors = [] + for geom in boundary_gdf.geometry: + if geom is None or geom.is_empty: + continue + if geom.geom_type == 'Polygon': + # 只取外壳（exterior），丢弃所有内部孔洞（interiors） + exteriors.append(Polygon(geom.exterior)) + elif geom.geom_type == 'MultiPolygon': + for part in geom.geoms: + exteriors.append(Polygon(part.exterior)) + + if exteriors: + boundary_gdf = gpd.GeoDataFrame(geometry=exteriors, crs=boundary_gdf.crs) + # ==================================================================== + + boundary_gdf_plotted = boundary_gdf # 留给末尾描边代码复用 + + # 物理擦除：多边形内为 True（保留原值），外部为 False（强制 NaN） + geom_mask = geometry_mask( + geometries=boundary_gdf.geometry, + out_shape=array.shape, + transform=transform, + invert=True, + ) + array = np.where(geom_mask, array, np.nan) + kept = int((~np.isnan(array)).sum()) + print(f"[visualize_raster] 矢量掩膜物理擦除完成: 陆地背景 → NaN " + f"(boundary={Path(boundary_shp_path).name}, 擦除后有效像元: " + f"{kept}/{array.size})") + except Exception as e: + print(f"[visualize_raster] 矢量掩膜物理擦除失败 ({boundary_shp_path}): {e}") + # ==================================================================== + # ── 极速降采样：>400 万像元时，将矩阵降维至约 200 万像素 ───────── + # 必须放在物理擦除之后，否则 geometry_mask 的 out_shape 会与降采样后 array 不对齐 # extent 使用原始 bounds（与降采样无关），保证坐标轴 UTM 米精确 # 降采样切片仅影响绘图渲染，可将 1 亿像素图在 1 秒内降至 ~200 万像素 _MAX_VIZ_PIXELS = 4_000_000 @@ -2379,15 +2458,6 @@ class ContentMapper: f"(step={step})，节省内存并加速渲染") w, h = w_downsampled, h_downsampled - # ── 全面 NoData 清洗：-9999.0 / NaN / Inf → 统一转为 np.nan ── - # 这一步确保陆地像素（无论来自掩膜还是原始 NoData）均被清除， - # 使 nanpercentile 分位数拉伸 100% 精准锁定水体内部 - array = np.where( - (array == nodata_value) | np.isnan(array) | np.isinf(array), - np.nan, - array - ) - # ── 从描述推断参数名和 colormap ─────────────────────────────── # 描述格式：Formula_Name|Category|Formula_Type|Formula param_name: Optional[str] = None @@ -2437,13 +2507,11 @@ class ContentMapper: scale_y = 1.0 # ── 准备图形 ───────────────────────────────────────────────── - # 画布大小保护：超大图像（如 40000×40000 px）在 DPI=300 输出时会导致 - # MemoryError；限制每维最大 100 英寸，防止内存爆炸 - _max_inch = 100 - safe_w = min(w / 100, _max_inch) # 像素 / 100 = 英寸，向上封顶 - safe_h = min(h / 100, _max_inch) - safe_figsize = (safe_w, safe_h) - fig, ax = plt.subplots(figsize=safe_figsize) + # 尊重用户的 figsize 参数（同时设置内存安全上限防止 DPI=300 下爆内存） + _max_inch = 60 # 每维最大 60 英寸，远超常规打印需求 + safe_w = min(float(figsize[0]), _max_inch) + safe_h = min(float(figsize[1]), _max_inch) + fig, ax = plt.subplots(figsize=(safe_w, safe_h)) # 计算有效值统计（2σ 标准差拉伸，排除长尾异常值干扰） valid = array[~np.isnan(array)] @@ -2468,27 +2536,17 @@ class ContentMapper: # 使用 masked array：NaN 区域自动不显示 masked_data = np.ma.masked_invalid(array) - try: - # 优先：pcolormesh（矢量输出，平滑颜色过渡） - im = ax.pcolormesh( - extent[0], extent[2], masked_data, - cmap=cmap or 'viridis', - vmin=vmin, vmax=vmax, - alpha=alpha, - shading='gouraud', # 颜色插值，平滑 - ) - except Exception: - # 备选：contourf - x_coords = np.linspace(extent[0], extent[1], w) - y_coords = np.linspace(extent[2], extent[3], h) - xx, yy = np.meshgrid(x_coords, y_coords) - im = ax.contourf( - xx, yy, masked_data, - levels=100, - cmap=cmap or 'viridis', - vmin=vmin, vmax=vmax, - alpha=alpha, - ) + # 【核心修复2】废弃错误的坐标映射逻辑。 + # 直接使用原生的 imshow，明确告知 matplotlib 第0行在最上方(origin='upper') + im = ax.imshow( + masked_data, + extent=[extent[0], extent[1], extent[2], extent[3]], + origin='upper', + cmap=cmap or 'viridis', + vmin=vmin, vmax=vmax, + alpha=alpha, + interpolation='bilinear' + ) # ★★★ 锁死绘图视口 ★★★ # 必须在所有叠加绘图（shp/colorbar/north arrow）之前执行， @@ -2496,34 +2554,37 @@ class ContentMapper: ax.set_xlim(extent[0], extent[1]) ax.set_ylim(extent[2], extent[3]) - # ── 边界 shapefile（叠加水域边界线）────────────────────────── - if boundary_shp_path and os.path.isfile(boundary_shp_path): + # ── 边界描边（直接复用上面的 boundary_gdf_plotted，不再重复读取/矢量化）───────── + if boundary_gdf_plotted is not None: try: - boundary_gdf = gpd.read_file(boundary_shp_path) - # 坐标系转换 - if crs_obj is not None: - target_crs = CRS.from_string(self.output_crs) - if boundary_gdf.crs != target_crs: - boundary_gdf = boundary_gdf.to_crs(target_crs) - boundary_gdf.boundary.plot(ax=ax, color='black', linewidth=1.5) + boundary_gdf_plotted.boundary.plot(ax=ax, color='black', linewidth=1.5) except Exception as e: - print(f"[visualize_raster] 边界 shapefile 叠加失败: {e}") + print(f"[visualize_raster] 边界描边失败: {e}") - # ── 坐标轴标签（固定 UTM 米，无条件覆盖）───────────────────── - ax.set_xlabel('X (UTM Meters)', fontsize=11) - ax.set_ylabel('Y (UTM Meters)', fontsize=11) + # ── 坐标轴标签（动态：根据 bounds 阈值判断经纬度 vs UTM 米）─── + # 经纬度坐标系下 bounds.left 在 [-180, 180]，UTM 投影坐标在百万级 + if _src_bounds is not None and _src_bounds.left < 180 and _src_bounds.bottom < 90: + ax.set_xlabel('Longitude', fontsize=14) + ax.set_ylabel('Latitude', fontsize=14) + else: + ax.set_xlabel('X (UTM Meters)', fontsize=14) + ax.set_ylabel('Y (UTM Meters)', fontsize=14) + # 显式设置刻度字号（避免 matplotlib 默认 10pt 在大画布下偏小） + ax.tick_params(axis='both', which='major', labelsize=12) + ax.tick_params(axis='both', which='minor', labelsize=10) ax.grid(True, linestyle='--', linewidth=0.5, alpha=0.4, color='gray') ax.set_axisbelow(True) # ── 标题（中文）────────────────────────────────────────────── - ax.set_title(chinese_title, fontsize=13, fontweight='bold', pad=10) + ax.set_title(chinese_title, fontsize=18, fontweight='bold', pad=15) # ── 颜色条（工业级样式：extend 三角 + MaxNLocator 刻度防重叠）───────── if show_colorbar and im is not None: try: - cbar = fig.colorbar(im, ax=ax, shrink=0.55, aspect=35, pad=0.02, extend='both') - cbar.set_label('Index Value', fontsize=10) + cbar = fig.colorbar(im, ax=ax, shrink=0.6, aspect=30, pad=0.03, extend='both') + cbar.set_label('Index Value', fontsize=14) + cbar.ax.tick_params(labelsize=12) cbar.locator = MaxNLocator(nbins=6) cbar.update_ticks() print("[visualize_raster] 颜色条添加成功") @@ -2544,7 +2605,8 @@ class ContentMapper: print(f"[visualize_raster] 指北针添加失败: {e}") # ── 紧凑布局并保存 ─────────────────────────────────────────── - plt.tight_layout() + # 给 tight_layout 显式 pad=2.0 防止指北针/比例尺/标题互相重叠 + plt.tight_layout(pad=2.0, h_pad=1.5, w_pad=1.5) try: plt.savefig( @@ -2567,20 +2629,159 @@ class ContentMapper: plt.close(fig) return output_file - def process_data(self, csv_file, shp_file, output_file='content_map.png', + # ------------------------------------------------------------------ + # Step 11 改造：插值网格 → GeoTIFF 物理落盘（替代 PNG 渲染） + # ------------------------------------------------------------------ + @staticmethod + def _redirect_png_to_tif_path(output_file: str) -> str: + """ + 智能路径重定向：把 14_visualization/visualization 改写到 11_Thematic_Map；.png 换 .tif + + 触发条件（与 panel 默认输出目录约定保持一致）： + - 原路径最后一级目录名为 ``14_visualization`` 或 ``visualization`` → 重定向到 ``11_Thematic_Map`` + - 其它情况（如用户自定义路径）→ 仅替换后缀，保留原目录 + """ + from pathlib import Path + p = Path(output_file) + parent = p.parent + + if parent.name in ('14_visualization', 'visualization'): + new_parent = parent.parent / '11_Thematic_Map' + else: + new_parent = parent + + new_filename = p.stem + '.tif' + return str(new_parent / new_filename) + + def _save_as_geotiff(self, grid_content, grid_xx, grid_yy, + output_tif_path, nodata_value=-9999.0): + """ + 将插值网格矩阵落盘为带坐标系的 GeoTIFF 文件（GDAL 实现）。 + + 与 ``src/utils/kriging.py:KrigingInterpolator.save_raster`` 风格保持一致： + GTiff + LZW + TILED + BIGTIFF=IF_SAFER + Float64。 + + Parameters + ---------- + grid_content : np.ndarray, shape (rows, cols) + 插值后的二维网格矩阵（含 NaN） + grid_xx : np.ndarray, shape (rows, cols) + X 坐标 meshgrid（与 grid_content 同 shape） + grid_yy : np.ndarray, shape (rows, cols) + Y 坐标 meshgrid（与 grid_content 同 shape） + output_tif_path : str + 输出 GeoTIFF 完整路径 + nodata_value : float + NoData 值（默认 -9999.0；NaN 将被替换为该值） + """ + try: + from osgeo import gdal, osr + except ImportError as e: + raise ImportError( + "需要 osgeo (GDAL) 库来写 GeoTIFF，请检查 conda 环境: " + str(e) + ) + + # ── 1. 计算 GeoTransform 6 参数 ────────────────────────────── + # GDAL 约定：transform = (x_min, dx, 0, y_max, 0, -abs(dy)) + # y 轴向下方为正（影像行列与地理坐标对应） + x_min = float(grid_xx[0, 0]) + y_max = float(grid_yy[-1, 0]) # meshgrid 末尾 = 北方 + dx = float(grid_xx[0, 1] - grid_xx[0, 0]) if grid_xx.shape[1] > 1 else 0.0 + dy_raw = float(grid_yy[1, 0] - grid_yy[0, 0]) if grid_yy.shape[0] > 1 else 0.0 + # dy 在地理上为正（向北递增），GDAL 用负值表示"行索引向下走 Y 增大" + dy = -abs(dy_raw) if dy_raw != 0 else -abs(dx) + + if abs(dx) < 1e-12 or abs(dy) < 1e-12: + raise ValueError( + f"网格分辨率异常: dx={dx}, dy={dy}，请检查 create_interpolation_grid 输入" + ) + + rows, cols = grid_content.shape + geotransform = (x_min, dx, 0, y_max, 0, dy) + + # ── 2. NaN → nodata；统一 Float64（与 kriging.py 一致）──────── + data_clean = np.where(np.isnan(grid_content), nodata_value, grid_content) + data_clean = data_clean.astype(np.float64) + + # 【核心修复1】上下翻转矩阵 + # Python网格第0行是南方，但GeoTIFF规范第0行是北方。 + # 必须翻转，否则落盘的 TIF 永远是上下颠倒的！ + data_clean = np.flipud(data_clean) + + # ── 3. CRS 字符串 → WKT（output_crs 是 EPSG 字符串或 WKT 均可）──── + try: + srs = osr.SpatialReference() + srs.SetFromUserInput(self.output_crs) + proj_wkt = srs.ExportToWkt() + except Exception as e: + print(f"[_save_as_geotiff] CRS 解析失败 ({e})，兜底使用 EPSG:4326") + srs = osr.SpatialReference() + srs.SetFromUserInput("EPSG:4326") + proj_wkt = srs.ExportToWkt() + + # ── 4. 创建输出目录 ───────────────────────────────────────── + out_dir = os.path.dirname(output_tif_path) + if out_dir: + os.makedirs(out_dir, exist_ok=True) + + # ── 5. GDAL 写盘 ───────────────────────────────────────────── + driver = gdal.GetDriverByName("GTiff") + if driver is None: + raise RuntimeError("GDAL GTiff 驱动不可用，请检查 osgeo 安装") + + dataset = driver.Create( + output_tif_path, cols, rows, 1, gdal.GDT_Float64, + options=[ + "COMPRESS=LZW", + "TILED=YES", + "BIGTIFF=IF_SAFER", # ⭐ 与 Step 8/10 Kriging 写盘保持一致 + ], + ) + if dataset is None: + raise RuntimeError(f"无法创建输出 GeoTIFF: {output_tif_path}") + + dataset.SetGeoTransform(geotransform) + dataset.SetProjection(proj_wkt) + + band = dataset.GetRasterBand(1) + band.WriteArray(data_clean) + band.SetNoDataValue(nodata_value) + band.ComputeStatistics(0) + band.FlushCache() + dataset.FlushCache() + del dataset + + valid_mask = ~np.isnan(grid_content) + print(f"[_save_as_geotiff] ✅ GeoTIFF 已保存: {output_tif_path}") + print(f" 分辨率: dx={dx:.6f}, dy={abs(dy):.6f}") + print(f" 范围: x=[{x_min:.4f}, {x_min + cols * dx:.4f}], " + f"y=[{y_max + rows * dy:.4f}, {y_max:.4f}]") + print(f" 尺寸: {rows} × {cols}, CRS: {self.output_crs}") + print(f" NoData={nodata_value}, 有效像元: {int(valid_mask.sum())}/{grid_content.size}") + return output_tif_path + + def process_data(self, csv_file, shp_file=None, output_file='content_map.png', resolution=100, show_sample_points=False, base_map_tif=None, use_distance_diffusion=True, max_diffusion_distance=None, diffusion_power=2, diffusion_n_neighbors=15, cmap=None, - expand_ratio=0.05): + expand_ratio=0.05, + output_format='tif'): """ 主处理函数 + ★★★ Plan C：shp_file 现在是可选参数（None = 纯采样点插值，不依赖水域掩膜）★★★ + Parameters: ----------- + csv_file : str + CSV文件路径 + shp_file : str, optional + 水域掩膜/边界文件路径（.shp / .dat / .bsq / .tif 等）。 + ★★★ None 时跳过所有边界相关逻辑，插值仅基于采样点自然扩展 ★★★ base_map_tif : str, optional TIF正射底图文件路径。如果提供，将在水域掩膜外显示底图 use_distance_diffusion : bool, default=True - 是否使用距离扩散方法填充边界空白区域 + 是否使用距离扩散方法填充边界空白区域（shp_file=None 时无效） max_diffusion_distance : float, optional 最大扩散距离（单位与坐标相同）。如果为None，自动计算为网格分辨率的5倍 diffusion_power : float, default=2 @@ -2590,7 +2791,9 @@ class ContentMapper: cmap : str, optional 颜色映射。如果为None，将从CSV文件名或内容中自动识别参数并选择对应的colormap expand_ratio : float, default=0.05 - 边界外扩比例（5%），确保图像边界不完全挨着地图 + 边界外扩比例（5%），用于从采样点范围外扩出图像边界 + output_format : str, default='tif' + 输出格式：'tif'（GeoTIFF）或 'png'（渲染图） """ try: # 自动识别参数名称并获取colormap @@ -2599,15 +2802,23 @@ class ContentMapper: cmap = self._get_colormap(param_name) else: print(f"使用指定的颜色映射: {cmap}") - - # 读取数据 - points_gdf = self.read_csv_data(csv_file) - boundary_gdf = self.read_boundary_shapefile(shp_file) - - # 对边缘采样点进行外扩处理（外扩到整个图像边界，按照resolution间距） - points_gdf = self._expand_edge_points(points_gdf, boundary_gdf, resolution=resolution, expand_ratio=expand_ratio) - # 创建插值网格 + # 读取采样点数据 + points_gdf = self.read_csv_data(csv_file) + + # ── Plan C: shp_file=None 时跳过所有水域掩膜逻辑 ─────────── + if shp_file is None: + print("[Plan C] shp_file=None，跳过水域掩膜读取，插值不依赖边界约束") + boundary_gdf = None + else: + boundary_gdf = self.read_boundary_shapefile(shp_file) + + # 对边缘采样点进行外扩处理（boundary_gdf=None 时基于采样点自身范围外扩） + points_gdf = self._expand_edge_points( + points_gdf, boundary_gdf, resolution=resolution, expand_ratio=expand_ratio + ) + + # 创建插值网格（boundary_gdf=None 时纯采样点插值，无掩膜裁剪） grid_xx, grid_yy, grid_content, bounds = self.create_interpolation_grid( points_gdf, boundary_gdf, resolution, expand_ratio=expand_ratio, @@ -2617,39 +2828,46 @@ class ContentMapper: diffusion_n_neighbors=diffusion_n_neighbors ) - # 生成含量图（包含不确定性叠加） - self.create_content_map( - points_gdf, boundary_gdf, grid_xx, grid_yy, - grid_content, bounds, output_file, show_sample_points, base_map_tif, - cmap=cmap - ) + # ── 按 output_format 分发落盘方式 ─────────────────────────── + if output_format == 'tif': + output_tif_path = self._redirect_png_to_tif_path(output_file) + self._save_as_geotiff(grid_content, grid_xx, grid_yy, output_tif_path) + actual_output = output_tif_path + else: + self.create_content_map( + points_gdf, boundary_gdf, grid_xx, grid_yy, + grid_content, bounds, output_file, show_sample_points, + base_map_tif, cmap=cmap + ) + actual_output = output_file print("处理完成！") - - # 输出统计信息 print(f"\n统计信息：") print(f"数据点数量: {len(points_gdf)}") print(f"含量值范围: {points_gdf['content'].min():.2f} - {points_gdf['content'].max():.2f}") print(f"含量值平均: {points_gdf['content'].mean():.2f}") print(f"含量值标准差: {points_gdf['content'].std():.2f}") + print(f"输出文件: {actual_output}") except Exception as e: print(f"处理过程中出现错误: {str(e)}") raise - def process_batch(self, csv_folder, shp_file, output_folder=None, + def process_batch(self, csv_folder, shp_file=None, output_folder=None, resolution=100, show_sample_points=False, base_map_tif=None, use_distance_diffusion=True, max_diffusion_distance=None, diffusion_power=2, diffusion_n_neighbors=15): """ 批量处理文件夹中的CSV文件 - + + ★★★ Plan C：shp_file 可选（None = 不依赖水域掩膜，纯采样点插值）★★★ + Parameters: ----------- csv_folder : str 包含CSV文件的文件夹路径 - shp_file : str - 边界shapefile文件路径 + shp_file : str, optional + 水域掩膜/边界文件路径。None 时跳过掩膜约束。 output_folder : str, optional 输出文件夹路径。如果为None，将在CSV文件所在文件夹创建'map_output'子文件夹 resolution : int, default=100 diff --git a/src/postprocessing/point_map.py b/src/postprocessing/point_map.py index b2edd92..a77c6bc 100644 --- a/src/postprocessing/point_map.py +++ b/src/postprocessing/point_map.py @@ -2,13 +2,6 @@ # -*- coding: utf-8 -*- """ 采样点地图生成模块 - 在高光谱假彩色影像上标注采样点 - -支持功能： -1. 读取高光谱影像并生成假彩色RGB图像 -2. 读取CSV文件中的采样点坐标（前两列为纬度、经度） -3. 在影像上标注红色采样点 -4. 添加指北针、图例和比例尺 -5. 支持地理坐标系转换 """ import numpy as np @@ -21,13 +14,13 @@ from matplotlib.patches import FancyArrowPatch import matplotlib.patheffects as path_effects # 性能优化配置 -plt.rcParams['agg.path.chunksize'] = 10000 # 提高矢量渲染性能 +plt.rcParams['agg.path.chunksize'] = 10000 plt.rcParams['path.simplify'] = True plt.rcParams['path.simplify_threshold'] = 0.1 -# 导入GDAL用于影像读写 try: from osgeo import gdal, osr + GDAL_AVAILABLE = True except ImportError: GDAL_AVAILABLE = False @@ -35,26 +28,15 @@ except ImportError: class SamplingPointMap: - """采样点地图生成类 - 在高光谱假彩色影像上标注采样点""" - def __init__(self, output_dir: str = "./point_maps", fast_mode: bool = False): - """ - 初始化采样点地图生成器 - - Args: - output_dir: 输出目录，用于保存生成的地图 - fast_mode: 是否启用快速模式（降低质量换取速度） - """ self.output_dir = Path(output_dir) self.output_dir.mkdir(parents=True, exist_ok=True) self.fast_mode = fast_mode - # 设置中文字体 plt.rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei', 'DejaVu Sans', 'Arial Unicode MS'] plt.rcParams['axes.unicode_minus'] = False plt.rcParams['font.size'] = 12 - # 性能优化设置 if fast_mode: plt.rcParams['figure.dpi'] = 150 plt.rcParams['savefig.dpi'] = 150 @@ -62,86 +44,45 @@ class SamplingPointMap: else: plt.rcParams['figure.dpi'] = 300 plt.rcParams['savefig.dpi'] = 300 - warnings.filterwarnings('ignore') - def create_sampling_point_map(self, - hyperspectral_path: str, - csv_path: str, - output_filename: Optional[str] = None, - rgb_bands: Optional[List[int]] = None, - point_color: str = 'red', - point_size: int = 80, - point_alpha: float = 0.8, - show_north_arrow: bool = True, - show_scale_bar: bool = True, - show_legend: bool = True, - dpi: int = None, - downsample: bool = False) -> str: - """ - 创建采样点地图：在高光谱假彩色影像上标注采样点 - - Args: - hyperspectral_path: 高光谱影像文件路径 (.dat, .bsq, .tif等) - csv_path: 采样点CSV文件路径（前两列为纬度、经度） - output_filename: 输出文件名（如果为None则自动生成） - rgb_bands: 用于RGB合成的三个波段索引 [R, G, B]，默认为None自动选择 - point_color: 采样点颜色 - point_size: 采样点大小 - point_alpha: 采样点透明度 - show_north_arrow: 是否显示指北针 - show_scale_bar: 是否显示比例尺 - show_legend: 是否显示图例 - dpi: 输出图像分辨率（None时使用fast_mode设置） - downsample: 是否对图像进行下采样以加快速度（大影像推荐启用） - - Returns: - 生成的地图文件路径 - """ + def create_sampling_point_map(self, hyperspectral_path: str, csv_path: str, + output_filename: Optional[str] = None, rgb_bands: Optional[List[int]] = None, + point_color: str = 'red', point_size: int = 80, point_alpha: float = 0.8, + show_north_arrow: bool = True, show_scale_bar: bool = True, + show_legend: bool = True, dpi: int = None, downsample: bool = False) -> str: if not GDAL_AVAILABLE: raise ImportError("GDAL未安装，无法处理地理坐标转换") print(f"正在生成采样点地图...{' (快速模式)' if self.fast_mode else ''}") - # 读取高光谱影像 - 优化：仅读取需要的RGB波段 hyperspectral_img, geotransform, projection, width, height, sample_factor = self._read_hyperspectral( hyperspectral_path, rgb_bands, downsample) - # 读取采样点 sampling_points = self._read_sampling_points(csv_path) - - # 生成假彩色图像 - 应用线性拉伸 rgb_image = self._create_false_color_image(hyperspectral_img) - - # 将地理坐标转换为像素坐标 - 支持投影系转换和下采样 pixel_coords = self._geo_to_pixel(sampling_points, geotransform, width, height, projection, sample_factor) - # 创建地图 if output_filename is None: csv_name = Path(csv_path).stem hs_name = Path(hyperspectral_path).stem output_filename = f"{hs_name}_{csv_name}_sampling_map.png" output_path = self.output_dir / output_filename - - # 使用更优化的绘图设置 if dpi is None: dpi = 150 if self.fast_mode else 200 self._create_map_visualization( - rgb_image, pixel_coords, sampling_points, - str(output_path), point_color, point_size, point_alpha, - show_north_arrow, show_scale_bar, show_legend, dpi, - geotransform, width, height, downsample, projection, sample_factor + rgb_image, pixel_coords, sampling_points, str(output_path), point_color, point_size, point_alpha, + show_north_arrow, show_scale_bar, show_legend, dpi, geotransform, width, height, downsample, projection, + sample_factor ) print(f"采样点地图已保存: {output_path}") return str(output_path) - def _read_hyperspectral(self, hyperspectral_path: str, - rgb_bands: Optional[List[int]] = None, - downsample: bool = False) -> Tuple[np.ndarray, tuple, str, int, int]: - """优化版：读取高光谱影像 - 仅读取需要的RGB波段""" + def _read_hyperspectral(self, hyperspectral_path: str, rgb_bands: Optional[List[int]] = None, + downsample: bool = False) -> Tuple[np.ndarray, tuple, str, int, int]: dataset = gdal.Open(hyperspectral_path) if dataset is None: raise ValueError(f"无法打开高光谱影像: {hyperspectral_path}") @@ -150,487 +91,222 @@ class SamplingPointMap: height = dataset.RasterYSize band_count = dataset.RasterCount - # 确定要读取的波段 - 优先使用指定波长 (650nm, 550nm, 460nm) if rgb_bands is None: if band_count >= 3: try: - # 使用find_band_number根据波长查找RGB波段 from src.utils.util import find_band_number rgb_bands = [ - find_band_number(650.0, hyperspectral_path), # Red ~650nm - find_band_number(550.0, hyperspectral_path), # Green ~550nm - find_band_number(460.0, hyperspectral_path) # Blue ~460nm + find_band_number(650.0, hyperspectral_path), + find_band_number(550.0, hyperspectral_path), + find_band_number(460.0, hyperspectral_path) ] - print(f" 根据波长选择RGB波段: R={rgb_bands[0]}, G={rgb_bands[1]}, B={rgb_bands[2]}") - except Exception as e: - print(f" 波长查找失败 ({e})，使用默认索引") - # 回退到基于索引的选择 - rgb_bands = [min(band_count-1, int(band_count*0.25)), - min(band_count-1, int(band_count*0.15)), - min(band_count-1, int(band_count*0.05))] + except Exception: + rgb_bands = [min(band_count - 1, int(band_count * 0.25)), + min(band_count - 1, int(band_count * 0.15)), + min(band_count - 1, int(band_count * 0.05))] else: rgb_bands = [0, 0, 0] - # 下采样控制 - 用户反馈下采样读取会导致像素值全为0 if downsample and (width > 2000 or height > 2000): - print(f" ⚠ 下采样暂被禁用（会导致像素值全0），使用原始分辨率: {width}x{height}") + print(f" ⚠ 下采样暂被禁用，使用原始分辨率: {width}x{height}") sample_factor = 1 - target_width = width - target_height = height else: sample_factor = 1 - target_width = width - target_height = height - # 只读取需要的RGB波段（性能关键优化） rgb_data = [] for band_idx in rgb_bands: band = dataset.GetRasterBand(band_idx + 1) - # 直接使用完整分辨率读取，避免下采样导致像素值为0的问题 band_data = band.ReadAsArray().astype(np.float32) rgb_data.append(band_data) - # 堆叠为RGB图像 (height, width, 3) if len(rgb_data) == 3: image_array = np.stack(rgb_data, axis=2) else: - # 如果只有1个波段，复制为RGB - image_array = np.stack([rgb_data[0]]*3, axis=2) + image_array = np.stack([rgb_data[0]] * 3, axis=2) geotransform = dataset.GetGeoTransform() projection = dataset.GetProjection() - - # 释放数据集 dataset = None - # 更新尺寸信息 - final_width = target_width if sample_factor > 1 else width - final_height = target_height if sample_factor > 1 else height - - print(f" 读取影像: {final_width}x{final_height}x{image_array.shape[2]} (RGB)") - if projection: - proj_type = "投影坐标系" if "PROJCS" in projection else "地理坐标系" - print(f" 影像投影: {proj_type}") - if sample_factor > 1: - print(f" 下采样因子: {sample_factor}") - - return image_array, geotransform, projection, final_width, final_height, sample_factor + return image_array, geotransform, projection, width, height, sample_factor def _read_sampling_points(self, csv_path: str) -> pd.DataFrame: - """读取采样点CSV文件""" - if not Path(csv_path).exists(): - raise FileNotFoundError(f"CSV文件不存在: {csv_path}") - + """智能读取采样点，自动识别模糊列名，允许UTM坐标，自动修复颠倒坐标""" df = pd.read_csv(csv_path) - - # 检查前两列是否为纬度和经度 if len(df.columns) < 2: - raise ValueError("CSV文件至少需要两列（纬度、经度）") + raise ValueError("CSV文件至少需要两列（经度、纬度 或 X、Y）") - # 假设前两列是纬度和经度 - lat_col = df.columns[0] - lon_col = df.columns[1] + # 智能子串匹配 + lat_aliases = ['lat', 'y', '纬'] + lon_aliases = ['lon', 'lng', 'x', '经'] + + lat_col = None + lon_col = None + cols_lower = {c: str(c).strip().lower() for c in df.columns} + + for c, lc in cols_lower.items(): + if lat_col is None and any(a in lc for a in lat_aliases): + lat_col = c + elif lon_col is None and any(a in lc for a in lon_aliases): + lon_col = c + + # 兜底：取前两列，默认列0=X(lon), 列1=Y(lat) + if lat_col is None or lon_col is None: + c0, c1 = df.columns[0], df.columns[1] + lon_col, lat_col = c0, c1 - # 重命名列 df = df.rename(columns={lat_col: 'latitude', lon_col: 'longitude'}) - - # 确保数值类型 df['latitude'] = pd.to_numeric(df['latitude'], errors='coerce') df['longitude'] = pd.to_numeric(df['longitude'], errors='coerce') + n_nan = int(df[['latitude', 'longitude']].isna().any(axis=1).sum()) + df = df.dropna(subset=['latitude', 'longitude']).reset_index(drop=True) - # 删除无效的坐标 - df = df.dropna(subset=['latitude', 'longitude']) + if len(df) > 0: + lat_max = df['latitude'].abs().max() + lon_max = df['longitude'].abs().max() - print(f"读取到 {len(df)} 个采样点") + # 智能对调：如果纬度 > 90，且经度 <= 90，说明用户把经纬度两列搞反了 + if lat_max > 90 and lon_max <= 90 and lat_max <= 180: + print(" ⚠ 检测到经纬度数值颠倒 (纬度>90, 经度<=90)，系统已自动对调坐标列") + df['latitude'], df['longitude'] = df['longitude'], df['latitude'] + # UTM 投影坐标判定：只要数值远大于180，就是米级别的投影系统 + elif lat_max > 180 or lon_max > 180: + print(f" ℹ 检测到坐标值远超180 (X:{lon_max:.1f}, Y:{lat_max:.1f})，判定为投影坐标(UTM)") + + print(f" CSV 列匹配: lat_col='{lat_col}', lon_col='{lon_col}'") + if n_nan: + print(f" 剔除 {n_nan} 个无效(NaN)行") + print(f" 读取到 {len(df)} 个有效采样点 (不再拦截越界拦截)") return df - def _create_false_color_image(self, image_array: np.ndarray, - rgb_bands: Optional[List[int]] = None) -> np.ndarray: - """创建假彩色RGB图像 - 应用线性拉伸和Gamma校正""" - # 由于_read_hyperspectral已返回RGB图像，这里仅进行最终处理 + def _create_false_color_image(self, image_array: np.ndarray, rgb_bands: Optional[List[int]] = None) -> np.ndarray: if image_array.shape[2] != 3: - # 确保是3通道 if len(image_array.shape) == 2 or image_array.shape[2] == 1: - if len(image_array.shape) == 2: - image_array = np.stack([image_array]*3, axis=2) - else: - image_array = np.repeat(image_array, 3, axis=2) + image_array = np.stack([image_array] * 3, axis=2) if len(image_array.shape) == 2 else np.repeat( + image_array, 3, axis=2) - print(f" 处理前图像范围: R[{image_array[:,:,0].min():.3f}-{image_array[:,:,0].max():.3f}], " - f"G[{image_array[:,:,1].min():.3f}-{image_array[:,:,1].max():.3f}], " - f"B[{image_array[:,:,2].min():.3f}-{image_array[:,:,2].max():.3f}]") - - # 增强型线性拉伸 - 解决图像太暗的问题 def simple_linear_stretch(data, min_percent=1, max_percent=99): - """增强对比度的线性拉伸""" valid_data = data[np.isfinite(data)] - if len(valid_data) == 0: - return np.zeros_like(data, dtype=np.float32) - - # 计算百分位数，使用更激进的拉伸 (1%-99%) + if len(valid_data) == 0: return np.zeros_like(data, dtype=np.float32) p_low = np.percentile(valid_data, min_percent) p_high = np.percentile(valid_data, max_percent) - if p_high - p_low < 1e-8: - # 如果数据范围太小，使用最小最大值归一化 - data_min = valid_data.min() - data_max = valid_data.max() - if data_max > data_min: - stretched = (data - data_min) / (data_max - data_min) - else: - stretched = np.zeros_like(data, dtype=np.float32) - else: - stretched = (data - p_low) / (p_high - p_low) + d_min, d_max = valid_data.min(), valid_data.max() + return (data - d_min) / (d_max - d_min) if d_max > d_min else np.zeros_like(data, dtype=np.float32) + stretched = (data - p_low) / (p_high - p_low) + return np.clip(stretched, 0.0, 1.0) - # 允许轻微过饱和以增加对比度 - stretched = np.clip(stretched, 0.0, 1.05) - stretched = np.clip(stretched, 0.0, 1.0) # 最终确保在[0,1] - return stretched - - # 对每个通道进行拉伸 r_stretched = simple_linear_stretch(image_array[:, :, 0]) g_stretched = simple_linear_stretch(image_array[:, :, 1]) b_stretched = simple_linear_stretch(image_array[:, :, 2]) - - # 合成为RGB图像 - rgb_image = np.stack([r_stretched, g_stretched, b_stretched], axis=2) - rgb_image = np.nan_to_num(rgb_image, nan=0.0) - - # 最终确保范围在[0,1]，并轻微增强对比度 + rgb_image = np.nan_to_num(np.stack([r_stretched, g_stretched, b_stretched], axis=2), nan=0.0) rgb_image = np.clip(rgb_image, 0.0, 1.0) + return (rgb_image * 255).astype(np.uint8) - # 可选：Gamma校正增加亮度（解决太暗问题） - gamma = 1 # <1会增加亮度 - rgb_image = np.power(rgb_image, gamma) - - # 映射到0-255范围（uint8），这样imshow显示效果更好 - rgb_image = (rgb_image * 255).astype(np.uint8) - - print(f" 处理后图像范围: [0-255] (Gamma={gamma})") - - return rgb_image - - def _geo_to_pixel(self, sampling_points: pd.DataFrame, - geotransform: tuple, width: int, height: int, - projection: str = "", sample_factor: int = 1) -> List[Tuple[float, float]]: - """ - 使用GDAL进行地理坐标到像素坐标的投影变换 - 支持下采样 - - 原始点位坐标格式: 41.66054612 124.2208338 (WGS84地理坐标: 纬度,经度) - 高光谱影像通常使用UTM或其他投影坐标系 - 当图像下采样时，sample_factor > 1，需要相应缩放坐标 - """ + def _geo_to_pixel(self, sampling_points: pd.DataFrame, geotransform: tuple, width: int, height: int, + projection: str = "", sample_factor: int = 1) -> List[Tuple[float, float]]: if geotransform is None or len(sampling_points) == 0: - # 如果没有地理变换信息，使用图像中心 - return [(width/2, height/2) for _ in range(len(sampling_points))] + return [(width / 2, height / 2) for _ in range(len(sampling_points))] pixel_coords = [] gt = geotransform + needs_transform = projection and ("PROJCS" in projection or "GEOGCS" in projection) - # 检查是否需要投影转换 - needs_transform = False - if projection and ("PROJCS" in projection or "GEOGCS" in projection): - needs_transform = True - print(f" 检测到影像投影: {projection[:80]}...") + # 智能判定是否为 WGS84 + sample_lon = float(sampling_points['longitude'].iloc[0]) + sample_lat = float(sampling_points['latitude'].iloc[0]) + is_wgs84 = (abs(sample_lon) <= 180) and (abs(sample_lat) <= 90) - # 创建坐标转换对象（WGS84 -> 影像投影） transform = None - if needs_transform and GDAL_AVAILABLE: + if needs_transform and is_wgs84 and GDAL_AVAILABLE: try: - # 源坐标系: WGS84 (EPSG:4326) src_srs = osr.SpatialReference() - src_srs.ImportFromEPSG(4326) # WGS84 - - # 目标坐标系: 影像的投影 + src_srs.ImportFromEPSG(4326) dst_srs = osr.SpatialReference() dst_srs.ImportFromWkt(projection) - - # 创建坐标转换 transform = osr.CoordinateTransformation(src_srs, dst_srs) - print(" ✓ 已创建WGS84到影像投影的坐标转换") except Exception as e: - print(f" ⚠ 坐标转换创建失败: {e}，使用简化变换") transform = None + elif not is_wgs84: + print(" ℹ 采样点为投影坐标(UTM)，跳过WGS84投影转换，直接使用放射变换映射") for _, row in sampling_points.iterrows(): - lon = float(row['longitude']) # 经度 (WGS84) - lat = float(row['latitude']) # 纬度 (WGS84) + lon, lat = float(row['longitude']), float(row['latitude']) if transform is not None: - # 使用GDAL进行投影转换: (经度, 纬度) -> (投影X, 投影Y) try: - proj_x, proj_y, _ = transform.TransformPoint(lat, lon) - # 再转换为像素坐标 - x = (proj_x - gt[0]) / gt[1] - y = (proj_y - gt[3]) / gt[5] - except Exception as e: - # 转换失败时回退到直接计算 - x = (lon - gt[0]) / gt[1] - y = (lat - gt[3]) / gt[5] + proj_x, proj_y, _ = transform.TransformPoint(lon, lat) + x, y = (proj_x - gt[0]) / gt[1], (proj_y - gt[3]) / gt[5] + except Exception: + x, y = width / 2, height / 2 else: - # 直接使用仿射变换（坐标系一致的情况） - x = (lon - gt[0]) / gt[1] - y = (lat - gt[3]) / gt[5] + x, y = (lon - gt[0]) / gt[1], (lat - gt[3]) / gt[5] - # 如果图像进行了下采样，需要相应缩放坐标 if sample_factor > 1: - x = x / sample_factor - y = y / sample_factor + x, y = x / sample_factor, y / sample_factor - # 限制在图像范围内（使用下采样后的尺寸） - x = max(0, min(x, width - 1)) - y = max(0, min(y, height - 1)) - - pixel_coords.append((x, y)) - - if transform is not None: - print(f" ✓ 使用GDAL投影变换处理 {len(pixel_coords)} 个采样点") - else: - print(f" 使用直接仿射变换处理 {len(pixel_coords)} 个采样点") + pixel_coords.append((max(0, min(x, width - 1)), max(0, min(y, height - 1)))) return pixel_coords - def _create_map_visualization(self, rgb_image: np.ndarray, - pixel_coords: List[Tuple[float, float]], - sampling_points: pd.DataFrame, - output_path: str, - point_color: str, - point_size: int, - point_alpha: float, - show_north_arrow: bool, - show_scale_bar: bool, - show_legend: bool, - dpi: int, - geotransform: tuple, - width: int, - height: int, - downsample: bool = False, - projection: str = "", - sample_factor: int = 1): - """创建地图可视化 - 优化版""" - # 使用更小的figure尺寸加快渲染 + def _create_map_visualization(self, rgb_image: np.ndarray, pixel_coords: List[Tuple[float, float]], + sampling_points: pd.DataFrame, output_path: str, point_color: str, point_size: int, + point_alpha: float, show_north_arrow: bool, show_scale_bar: bool, show_legend: bool, + dpi: int, geotransform: tuple, width: int, height: int, downsample: bool = False, + projection: str = "", sample_factor: int = 1): figsize = (10, 8) if self.fast_mode or downsample else (12, 10) fig, ax = plt.subplots(figsize=figsize, dpi=100 if self.fast_mode else 150) - - # 显示假彩色图像 - 现在已经是0-255的uint8格式 - print(f" 最终图像数据范围: [{rgb_image.min()}, {rgb_image.max()}] (uint8)") ax.imshow(rgb_image, interpolation='nearest' if self.fast_mode else 'bilinear') - # 绘制采样点 - 优化：使用scatter代替循环plot if pixel_coords: - x_coords = [p[0] for p in pixel_coords] - y_coords = [p[1] for p in pixel_coords] - ax.scatter(x_coords, y_coords, c=point_color, s=point_size, - alpha=point_alpha, edgecolors='white', linewidth=1.5) + x_coords, y_coords = [p[0] for p in pixel_coords], [p[1] for p in pixel_coords] + ax.scatter(x_coords, y_coords, c=point_color, s=point_size, alpha=point_alpha, edgecolors='white', + linewidth=1.5) - # 添加指北针 - if show_north_arrow: - self._add_north_arrow(ax, width, height, position='bottom-left', direction='down') + if show_north_arrow: self._add_north_arrow(ax, width, height, position='bottom-left', direction='down') + if show_scale_bar and geotransform is not None: self._add_scale_bar(ax, geotransform, width, height) - # 添加比例尺 - if show_scale_bar and geotransform is not None: - self._add_scale_bar(ax, geotransform, width, height) - - # 添加图例 if show_legend: - legend_text = f'采样点 (n={len(sampling_points)})' - ax.plot([], [], 'o', color=point_color, markersize=8, label=legend_text) + ax.plot([], [], 'o', color=point_color, markersize=8, label=f'采样点 (n={len(sampling_points)})') ax.legend(loc='lower right', frameon=True, facecolor='white', edgecolor='gray') - # 设置标题和标签 ax.set_title('高光谱影像采样点分布图', fontsize=16, fontweight='bold', pad=20) - - - # 隐藏坐标轴刻度 ax.set_xticks([]) ax.set_yticks([]) - - # 添加网格 ax.grid(True, alpha=0.2, linestyle='--') - plt.tight_layout() - # 保存参数 - 避免传递不兼容的参数 - save_kwargs = { - 'dpi': dpi, - 'bbox_inches': 'tight', - 'pad_inches': 0.05, - 'facecolor': 'white' - } - - # 仅添加matplotlib支持的参数 - if self.fast_mode: - save_kwargs['dpi'] = min(dpi, 180) # 快速模式降低DPI - + save_kwargs = {'dpi': min(dpi, 180) if self.fast_mode else dpi, 'bbox_inches': 'tight', 'pad_inches': 0.05, + 'facecolor': 'white'} plt.savefig(output_path, **save_kwargs) plt.close(fig) - def _add_north_arrow(self, ax, width: int, height: int, position='top-right', direction='down', - size=0.08, color='white', n_color='white', outline_color='black'): - """ - 添加指北针，可配置位置、方向、大小、颜色。 + def _add_north_arrow(self, ax, width: int, height: int, position='top-right', direction='down', size=0.08, + color='white', n_color='white', outline_color='black'): + pos_map = {'top-left': (0.08, 0.88), 'top-right': (0.92, 0.88), 'bottom-left': (0.08, 0.12), + 'bottom-right': (0.92, 0.12)} + arrow_x, arrow_y = width * pos_map.get(position, (0.92, 0.88))[0], height * pos_map.get(position, (0.92, 0.88))[ + 1] + dx, dy = {'up': (0, size), 'down': (0, -size), 'left': (-size, 0), 'right': (size, 0)}.get(direction, + (0, -size)) - 参数： - ax: matplotlib Axes对象 - width, height: 图像宽高（用于相对定位） - position: 'top-left', 'top-right', 'bottom-left', 'bottom-right' - direction: 'up', 'down', 'left', 'right' 箭头指向 - size: 箭头长度相对于高度的比例（0.05~0.12） - color: 箭头颜色 - n_color: 'N' 文字颜色 - outline_color: 文字描边颜色 - """ - # 位置映射（偏移系数） - pos_map = { - 'top-left': (0.08, 0.88), - 'top-right': (0.92, 0.88), - 'bottom-left': (0.08, 0.12), - 'bottom-right': (0.92, 0.12), - } - arrow_x_ratio, arrow_y_ratio = pos_map.get(position, (0.92, 0.88)) - arrow_x = width * arrow_x_ratio - arrow_y = height * arrow_y_ratio - - # 方向映射（箭头终点偏移） - direction_map = { - 'up': (0, +size), - 'down': (0, -size), - 'left': (-size, 0), - 'right': (+size, 0), - } - dx, dy = direction_map.get(direction, (0, -size)) - end_x = arrow_x + dx * width # 注意：dx是比例，乘以宽度/高度保持比例一致 - end_y = arrow_y + dy * height - - # 箭头绘制 - arrow = FancyArrowPatch((arrow_x, arrow_y), (end_x, end_y), - color=color, linewidth=3, - arrowstyle='->', mutation_scale=20) + arrow = FancyArrowPatch((arrow_x, arrow_y), (arrow_x + dx * width, arrow_y + dy * height), color=color, + linewidth=3, arrowstyle='->', mutation_scale=20) ax.add_patch(arrow) - - # N 文字位置：在箭头尾部或头部？通常放在箭头指向的反方向末端 - # 这里放在箭头尾部向外偏移一点（便于阅读） - # 偏移系数根据方向决定 - offset_scale = 0.02 # 偏移量比例 - if direction == 'up': - text_x = arrow_x - text_y = arrow_y - height * offset_scale # 放在箭头下方 - elif direction == 'down': - text_x = arrow_x - text_y = arrow_y + height * offset_scale # 放在箭头上方 - elif direction == 'left': - text_x = arrow_x + width * offset_scale - text_y = arrow_y - else: # right - text_x = arrow_x - width * offset_scale - text_y = arrow_y - - ax.text(text_x, text_y, 'N', fontsize=14, fontweight='bold', - color=n_color, ha='center', va='center', + text_y = arrow_y - height * 0.02 if direction == 'up' else arrow_y + height * 0.02 + ax.text(arrow_x, text_y, 'N', fontsize=14, fontweight='bold', color=n_color, ha='center', va='center', path_effects=[path_effects.withStroke(linewidth=3, foreground=outline_color)]) def _add_scale_bar(self, ax, geotransform: tuple, width: int, height: int): - """添加比例尺""" - if geotransform is None: - return - - # 计算图像实际宽度（米） + if geotransform is None: return pixel_size_x = abs(geotransform[1]) - image_width_meters = width * pixel_size_x - - # 选择合适的比例尺长度（图像宽度的1/4） - scale_length_m = image_width_meters / 4 - scale_length_pixels = width / 4 - - # 找到合适的刻度 - scale_options = [1000, 500, 200, 100, 50, 20, 10, 5, 2, 1] - scale_meters = next((s for s in scale_options if s <= scale_length_m), 1) - + scale_length_m = (width * pixel_size_x) / 4 + scale_meters = next((s for s in [1000, 500, 200, 100, 50, 20, 10, 5, 2, 1] if s <= scale_length_m), 1) scale_pixels = int(scale_meters / pixel_size_x) + bar_x, bar_y = width * 0.08, height * 0.92 - # 在左下角添加比例尺 - bar_x = width * 0.08 - bar_y = height * 0.92 - - # 绘制比例尺线 ax.plot([bar_x, bar_x + scale_pixels], [bar_y, bar_y], color='white', linewidth=4) - - # 添加刻度线 ax.plot([bar_x, bar_x], [bar_y, bar_y + 8], color='white', linewidth=2) ax.plot([bar_x + scale_pixels, bar_x + scale_pixels], [bar_y, bar_y + 8], color='white', linewidth=2) - - # 添加文字 - ax.text(bar_x + scale_pixels/2, bar_y , f'{scale_meters} m', - fontsize=11, ha='center', va='bottom', fontweight='bold', - bbox=dict(facecolor='white', alpha=0.8, edgecolor='none', pad=1)) - - def batch_create_maps(self, hyperspectral_path: str, - csv_folder: str, - output_subdir: str = "sampling_maps", - fast_mode: bool = True) -> Dict[str, str]: - """ - 批量创建采样点地图 - - Args: - hyperspectral_path: 高光谱影像路径 - csv_folder: 包含多个CSV文件的文件夹 - output_subdir: 输出子目录 - - Returns: - 生成的地图文件路径字典 - """ - csv_folder_path = Path(csv_folder) - if not csv_folder_path.exists(): - raise FileNotFoundError(f"CSV文件夹不存在: {csv_folder}") - - # 创建输出目录 - output_dir = self.output_dir / output_subdir - output_dir.mkdir(parents=True, exist_ok=True) - - map_paths = {} - - # 查找所有CSV文件 - csv_files = list(csv_folder_path.glob("*.csv")) - - print(f"找到 {len(csv_files)} 个CSV文件，开始批量生成采样点地图... (快速模式: {fast_mode})") - - for csv_file in csv_files: - try: - output_filename = f"{Path(hyperspectral_path).stem}_{csv_file.stem}_sampling_map.png" - map_path = self.create_sampling_point_map( - hyperspectral_path=hyperspectral_path, - csv_path=str(csv_file), - output_filename=output_filename, - downsample=True, # 批量模式默认下采样 - dpi=120 if fast_mode else 200 - ) - map_paths[csv_file.name] = map_path - print(f"✓ 生成: {csv_file.name}") - - except Exception as e: - print(f"✗ 处理 {csv_file.name} 失败: {e}") - - print(f"批量生成完成，共生成 {len(map_paths)} 个采样点地图") - return map_paths - - -# 测试代码 -if __name__ == "__main__": - # 示例用法 - map_generator = SamplingPointMap(output_dir="./point_maps") - - # 测试代码已禁用，避免直接运行时出错 - map_generator_fast = SamplingPointMap(output_dir="./point_maps", fast_mode=True) - map_path = map_generator_fast.create_sampling_point_map( - hyperspectral_path=r"D:\BaiduNetdiskDownload\yaobao\result3.bsq", - csv_path=r"E:\code\WQ\pipeline_result\work_dir\4_processed_data\processed_data.csv", - downsample=True, - dpi=150 - ) - print("测试代码已注释，请通过GUI或手动调用使用。") - - print("SamplingPointMap类已创建，可以用于生成带采样点的地图。") - print("性能优化功能：") - print(" - fast_mode=True: 快速模式 (推荐用于预览)") - print(" - downsample=True: 对大影像下采样 (推荐用于>2000x2000影像)") - print(" - 使用: SamplingPointMap(fast_mode=True).create_sampling_point_map(...)") + ax.text(bar_x + scale_pixels / 2, bar_y, f'{scale_meters} m', fontsize=11, ha='center', va='bottom', + fontweight='bold', bbox=dict(facecolor='white', alpha=0.8, edgecolor='none', pad=1)) \ No newline at end of file diff --git a/src/postprocessing/visualization_reports.py b/src/postprocessing/visualization_reports.py index 0011225..932f9f2 100644 --- a/src/postprocessing/visualization_reports.py +++ b/src/postprocessing/visualization_reports.py @@ -385,7 +385,7 @@ class WaterQualityVisualization: return output_paths def plot_distribution_map_enhanced(self, prediction_csv_path: str, - boundary_shp_path: str, + boundary_shp_path: Optional[str] = None, # ★★★ Plan C: None = 不依赖水域掩膜 ★★★ parameter_column: str = 'prediction', output_path: Optional[str] = None, resolution: float = 30, @@ -394,12 +394,12 @@ class WaterQualityVisualization: colormap: str = 'viridis') -> str: """ 生成增强的含量分布图（彩色填充图） - - 这是对step9的增强版本，使用更丰富的颜色映射 - + + ★★★ Plan C：boundary_shp_path 可选（None = 不依赖水域掩膜）★★★ + Args: prediction_csv_path: 预测结果CSV文件路径 - boundary_shp_path: 边界shapefile文件路径 + boundary_shp_path: 边界/掩膜文件路径。None 时跳过水域掩膜约束。 parameter_column: 参数值列名 output_path: 输出图片路径 resolution: 插值网格分辨率 diff --git a/src/utils/band_math.py b/src/utils/band_math.py index ffdb166..331bffc 100644 --- a/src/utils/band_math.py +++ b/src/utils/band_math.py @@ -99,8 +99,10 @@ class BandMathCalculator: # 【新增安全防护】引入 numpy 命名空间，让 eval 引擎安全识别 nan 与 inf import numpy as np try: - # 即使 calc_expression 含有纯字符 nan，也能被 np.nan 安全接管 - result = eval(calc_expression, {"__builtins__": None}, {"nan": np.nan, "inf": np.inf, "np": np}) + # 【P0 修复】包 np.errstate 抑制除零 / 无效操作产生的 RuntimeWarning 洪水 + with np.errstate(divide='ignore', invalid='ignore'): + # 即使 calc_expression 含有纯字符 nan，也能被 np.nan 安全接管 + result = eval(calc_expression, {"__builtins__": None}, {"nan": np.nan, "inf": np.inf, "np": np}) except Exception as e: print(f"⚠️ 警告：公式计算异常 ({e})，该点赋值为 nan") result = np.nan diff --git a/src/utils/water_index.py b/src/utils/water_index.py index 57d3279..8660e28 100644 --- a/src/utils/water_index.py +++ b/src/utils/water_index.py @@ -117,7 +117,9 @@ class WaterQualityIndexCalculator: calc_expr, ) try: - r = eval(calc_expr, {"__builtins__": None}, {"nan": np.nan, "inf": np.inf, "np": np}) + # 【P0 修复】包 np.errstate 抑制除零 / 无效操作产生的 RuntimeWarning 洪水 + with np.errstate(divide='ignore', invalid='ignore'): + r = eval(calc_expr, {"__builtins__": None}, {"nan": np.nan, "inf": np.inf, "np": np}) except Exception: r = np.nan results.append(r)