import pandas as pd import numpy as np import geopandas as gpd from osgeo import gdal from pathlib import Path from typing import Optional, Tuple from pyproj import CRS, Transformer import matplotlib.pyplot as plt import matplotlib.patches as patches from matplotlib.ticker import FuncFormatter, MaxNLocator from matplotlib_scalebar.scalebar import ScaleBar from scipy.interpolate import griddata from scipy import ndimage from scipy.spatial.distance import cdist from scipy.spatial import ConvexHull from shapely.geometry import Point, Polygon import rasterio from rasterio.features import geometry_mask, shapes from rasterio import windows from rasterio.warp import calculate_default_transform, reproject, Resampling try: from affine import Affine except ImportError: try: from rasterio.transform import Affine except ImportError: Affine = None import warnings import math import os import random import glob # 尝试导入pykrige（可选依赖） try: from pykrige.ok import OrdinaryKriging PYKRIGE_AVAILABLE = True except ImportError: PYKRIGE_AVAILABLE = False print("警告: pykrige未安装，Kriging不确定性计算将不可用") warnings.filterwarnings('ignore') # 设置中文字体 plt.rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei'] plt.rcParams['axes.unicode_minus'] = False # 参数到颜色映射的字典 PARAMS_CMAP = { "Chlorophyll": "YlGnBu_r", "COD": "coolwarm", "DO": "RdYlBu", "pH": "Spectral", "Temperature": "turbo", "spCond": "cividis", "Turbidity": "YlOrBr", "TDS": "inferno", "Cl-": "RdYlBu_r", "NO3-N": "YlOrRd", "NH3-N": "magma", "BGA": "viridis", "TT": "RdYlBu_r" } # ── 水色指数英文名 → 中文标题映射（精确整词匹配）──────────────────── # 优先于动态拼接；当文件名恰好命中完整 key 时使用 INDEX_TITLE_MAP = { # Chl_a 叶绿素a "Chl_Conc_NDCI": "叶绿素a浓度估算_NDCI模型", "Chl_MM12NDCI": "叶绿素a相对指数_Matthews12模型", "Chl_Conc_Gao": "叶绿素a浓度估算_Gao模型", "Chl_Conc_QAA": "叶绿素a浓度估算_QAA模型", # BGA 蓝藻/藻蓝蛋白 "BGA_Go04MCI": "蓝藻相对指数_Gower04模型", "BGA_PC_Conc_Mishra": "藻蓝蛋白浓度估算_Mishra模型", "BGA_Conc": "蓝藻浓度估算", "BGA_Am09KBBI": "蓝藻相对指数_Am09模型", # Turb 浊度 "Turb_Dox02NIRoverRed": "水体浊度指数_Doxaran02模型", "Turbidity": "水体浊度", "Turb_Conc": "浊度估算", # TSM 悬浮物 "TSM_Conc_Bowling": "总悬浮物浓度估算_Bowling模型", "TSM_Conc": "总悬浮物浓度估算", # CDOM 有色溶解有机物 "CDOM_Conc": "有色溶解有机物浓度估算", # WI 水色指数综合 "WaterIndex": "水色指数", "NDCI": "归一化叶绿素差值指数_NDCI", "MCI": "最大Chlorophyll指数_MCI", } # ── 关键词 → 中文词根映射（用于动态拼接）──────────────────────────── # 顺序即优先级：长词优先，短词兜底 PART_NAME_MAP = [ # 1. 指数/浓度类（最高优先级，描述输出类型） ("Conc", "浓度估算"), ("_Conc", "浓度估算"), ("Concentration", "浓度"), ("Index", "指数"), ("_Index", "指数"), # 2. 模型/方法标识（次高优先级） ("NDCI", "NDCI模型"), ("MCI", "MCI模型"), ("FLH", "FLH荧光基线"), ("QAA", "QAA模型"), ("Go04", "Gower04"), ("MM12", "Matthews12"), ("Gao", "Gao模型"), ("Dox02", "Doxaran02"), ("Bowling", "Bowling模型"), ("Mishra", "Mishra模型"), ("Am09", "Am09"), ("KBBI", "KBBI"), ("PC", "藻蓝蛋白"), # 3. 参数大类（核心水质指标） ("BGA", "蓝藻相对指数"), ("Chl", "叶绿素a"), ("Chl_a", "叶绿素a"), ("Turb", "浊度"), ("TSM", "总悬浮物"), ("CDOM", "有色溶解有机物"), ("DO", "溶解氧"), ("pH", "pH值"), ("NH3", "氨氮"), ("NO3", "硝态氮"), ("TDS", "溶解性总固体"), ] class ContentMapper: def __init__(self, input_crs='EPSG:32651', output_crs='EPSG:4326'): """ 初始化ContentMapper - 生成平滑的含量分布图 本类专门用于生成平滑、均匀的颜色分布图，而不是显示离散的采样点。 通过高密度网格插值和多级颜色映射，创建连续的颜色过渡效果。 Parameters: ----------- input_crs : str 输入坐标系，默认为'EPSG:32651' (WGS_1984_UTM_Zone_51N) output_crs : str 输出坐标系，默认为'EPSG:4326' (WGS84) """ # 定义坐标转换器 self.input_crs = input_crs self.output_crs = output_crs self.transformer = Transformer.from_crs( CRS.from_string(input_crs), CRS.from_string(output_crs), always_xy=True ) # 参数到颜色映射的字典 self.params_cmap = PARAMS_CMAP.copy() # 所有可用的matplotlib colormap列表（用于随机选择） self.available_cmaps = ['viridis', 'plasma', 'inferno', 'magma', 'cividis', 'coolwarm', 'RdYlBu', 'Spectral', 'YlGnBu_r', 'YlOrBr', 'YlOrRd', 'turbo', 'RdYlBu_r', 'cool', 'hot', 'jet'] print(f"坐标转换设置: {input_crs} -> {output_crs}") # ── 内部工具 ───────────────────────────────────────────────────── @staticmethod def _get_chinese_title(stem: str) -> str: """ 根据 GeoTIFF 文件名 stem 返回中文图表标题（绝对唯一）。 匹配策略： 1. 精确整词命中 INDEX_TITLE_MAP 2. 中文分类前缀 + 原始模型后缀（绝对唯一保证） 3. 未匹配任何关键词 → 返回原英文 stem Parameters ---------- stem : str GeoTIFF 文件名（不含路径和扩展名） Returns ------- str 中文标题；若未匹配则返回英文 stem """ # 策略1：精确整词匹配（优先级最高） if stem in INDEX_TITLE_MAP: return INDEX_TITLE_MAP[stem] # 策略2：中文分类前缀 + 原始模型后缀（确保绝对唯一） category = "" suffix = "" if stem.startswith("BGA_PC_Conc"): category = "藻蓝蛋白浓度估算" suffix = stem.replace("BGA_PC_Conc_", "") elif stem.startswith("BGA"): category = "蓝藻相对指数" suffix = stem.replace("BGA_", "") elif stem.startswith("Chl_Conc"): category = "叶绿素a浓度估算" suffix = stem.replace("Chl_Conc_", "") elif stem.startswith("Chl"): category = "叶绿素a相对指数" suffix = stem.replace("Chl_", "") elif stem.startswith("Turb_Conc"): category = "浊度浓度估算" suffix = stem.replace("Turb_Conc_", "") elif stem.startswith("Turb"): category = "浊度相对指数" suffix = stem.replace("Turb_", "") elif stem.startswith("TSM_Conc"): category = "总悬浮物浓度估算" suffix = stem.replace("TSM_Conc_", "") else: # 兜底机制：未知分类直接使用原名 category = "水质参数" suffix = stem return f"{category}_{suffix}" def _extract_param_name(self, csv_file): """ 从CSV文件名或内容中提取参数名称 Parameters: ----------- csv_file : str CSV文件路径 Returns: -------- param_name : str or None 提取的参数名称，如果未找到则返回None """ print(f"[调试] 开始从文件 {csv_file} 中提取参数名称") print(f"[调试] 字典中的参数键: {list(self.params_cmap.keys())}") # 从文件名中提取（去除路径和扩展名） file_name = os.path.basename(csv_file) file_name_no_ext = os.path.splitext(file_name)[0] print(f"[调试] 文件名（不含扩展名）: {file_name_no_ext}") # 尝试从文件名中匹配参数名称（不区分大小写） file_name_upper = file_name_no_ext.upper() for param in self.params_cmap.keys(): param_upper = param.upper() if param_upper in file_name_upper: print(f"从文件名中识别到参数: {param} (匹配到 '{param_upper}' 在 '{file_name_upper}' 中)") return param # 返回字典中的原始键（保持大小写） # 如果文件名中没有找到，尝试从CSV内容中提取（检查列名） try: df = pd.read_csv(csv_file, encoding='utf-8', nrows=0) # 只读取列名 columns = [col.upper() for col in df.columns] print(f"[调试] CSV列名: {list(df.columns)}") for param in self.params_cmap.keys(): param_upper = param.upper() # 检查列名中是否包含参数名称 for col in columns: if param_upper in col or col in param_upper: print(f"从CSV列名中识别到参数: {param} (匹配到列名 '{col}')") return param # 返回字典中的原始键（保持大小写） except Exception as e: print(f"读取CSV列名时出错: {e}") print(f"未能在文件 {csv_file} 中识别参数名称") print(f"[调试] 可用的参数名: {list(self.params_cmap.keys())}") return None def _get_colormap(self, param_name=None): """ 根据参数名称获取对应的colormap（支持精确匹配、模糊匹配） Parameters ---------- param_name : str, optional 参数名称。如果为None或不在映射中，则随机选择一个colormap Returns ------- cmap : str 颜色映射名称 """ print(f"[调试] _get_colormap 被调用，param_name={param_name}") print(f"[调试] 当前字典中的键: {list(self.params_cmap.keys())}") if param_name: # 精确匹配（区分大小写） if param_name in self.params_cmap: cmap = self.params_cmap[param_name] print(f"使用参数 '{param_name}' 对应的颜色映射: {cmap}") return cmap # 不区分大小写匹配 param_name_upper = param_name.upper() for key in self.params_cmap.keys(): if key.upper() == param_name_upper: cmap = self.params_cmap[key] print(f"使用参数 '{key}' (不区分大小写匹配 '{param_name}') 对应的颜色映射: {cmap}") return cmap # ── 模糊匹配（关键字包含检测）─────────────────────────── pn_upper = param_name.upper() pn_lower = param_name.lower() # 蓝藻 / BGA / Phycocyanin → YlGn（蓝绿色系） if any(k in pn_upper for k in ('BGA', 'PHYCO', 'CYAN', '蓝藻', '藻蓝')): cmap = self.params_cmap.get('BGA', 'YlGn') print(f"模糊匹配 BGA/Phycocyanin → '{cmap}'") return cmap # 叶绿素 / Chlorophyll / Chl → YlGn（绿色系） if any(k in pn_upper for k in ('CHL', '叶绿素', 'CHLORO')): cmap = self.params_cmap.get('Chl_a', 'YlGn') print(f"模糊匹配 Chl/叶绿素 → '{cmap}'") return cmap # CDOM / 有色溶解有机物 if any(k in pn_upper for k in ('CDOM', '色DOM', '有色溶解')): cmap = self.params_cmap.get('CDOM', 'BrBG') print(f"模糊匹配 CDOM → '{cmap}'") return cmap # 悬浮物 / TSM / SS → YlOrBr（黄棕系） if any(k in pn_upper for k in ('TSM', 'SS', '悬浮物', '总悬浮')): cmap = self.params_cmap.get('TSM', 'YlOrBr') print(f"模糊匹配 TSM/悬浮物 → '{cmap}'") return cmap # 透明度 / SD / Secchi → Blues（蓝色系） if any(k in pn_upper for k in ('SD', 'SECCHI', '透明度', '透明')): cmap = self.params_cmap.get('SD', 'Blues') print(f"模糊匹配 SD/透明度 → '{cmap}'") return cmap # 氨氮 / NH4 / NH3 → Oranges if any(k in pn_upper for k in ('NH4', 'NH3', '氨氮', '氨')): cmap = 'Oranges' print(f"模糊匹配 NH4/氨氮 → '{cmap}'") return cmap # 总磷 / TP / 总氮 / TN → RdYlGn if any(k in pn_upper for k in ('TP', '总磷')): cmap = 'RdYlGn_r' print(f"模糊匹配 TP/总磷 → '{cmap}'") return cmap if any(k in pn_upper for k in ('TN', '总氮')): cmap = 'RdYlGn_r' print(f"模糊匹配 TN/总氮 → '{cmap}'") return cmap # 高浊度 / Turbidity → PuBu（紫蓝系） if any(k in pn_upper for k in ('TURBIDITY', '浊度', 'TURB')): cmap = 'PuBu' print(f"模糊匹配 Turbidity/浊度 → '{cmap}'") return cmap # 溶解氧 / DO → cool（蓝白冷色） if any(k in pn_upper for k in ('DO', '溶解氧', 'DISSOLVED')): cmap = 'cool' print(f"模糊匹配 DO/溶解氧 → '{cmap}'") return cmap # 仍不匹配 → 随机 cmap = random.choice(self.available_cmaps) print(f"警告: 参数 '{param_name}' 不在映射中，随机选择颜色映射: {cmap}") print(f"可用的参数名: {list(self.params_cmap.keys())}") return cmap else: cmap = random.choice(self.available_cmaps) print(f"未指定参数名称，随机选择颜色映射: {cmap}") return cmap def _check_point_distribution(self, points): """检查数据点的几何分布""" print("正在检查数据点分布...") # 检查是否有重复点 unique_points = np.unique(points, axis=0) if len(unique_points) < len(points): print(f"警告：发现 {len(points) - len(unique_points)} 个重复数据点") # 检查点是否共线 if len(unique_points) >= 3: # 计算前三个点构成的三角形面积 p1, p2, p3 = unique_points[:3] area = 0.5 * abs((p2[0] - p1[0]) * (p3[1] - p1[1]) - (p3[0] - p1[0]) * (p2[1] - p1[1])) if area < 1e-10: # 面积太小，可能共线 print("警告：前三个数据点可能共线") # 尝试找到不共线的点 for i in range(3, len(unique_points)): p4 = unique_points[i] area = 0.5 * abs((p2[0] - p1[0]) * (p4[1] - p1[1]) - (p4[0] - p1[0]) * (p2[1] - p1[1])) if area > 1e-10: print(f"找到非共线点，使用点 {i}") break else: print("警告：所有数据点可能都共线，这会导致插值失败") # 检查坐标范围 x_range = points[:, 0].max() - points[:, 0].min() y_range = points[:, 1].max() - points[:, 1].min() if x_range < 1e-6 or y_range < 1e-6: print(f"警告：坐标范围很小 (X范围: {x_range:.2e}, Y范围: {y_range:.2e})") print("这可能导致插值数值不稳定") return unique_points def _fill_boundary_blanks_with_distance_diffusion(self, grid_content, grid_xx, grid_yy, mask, boundary_gdf, max_diffusion_distance=None, power=2, n_neighbors=15): """ 使用距离扩散方法填充边界附近的空白区域 Parameters: ----------- grid_content : np.ndarray 插值网格数据 grid_xx : np.ndarray 网格X坐标 grid_yy : np.ndarray 网格Y坐标 mask : np.ndarray 边界掩膜（True表示边界内） boundary_gdf : gpd.GeoDataFrame 边界几何数据 max_diffusion_distance : float, optional 最大扩散距离（单位与坐标相同）。如果为None，自动计算为网格分辨率的5倍 power : float, default=2 IDW距离衰减幂参数 n_neighbors : int, default=15 使用的最近邻点数 Returns: -------- grid_content : np.ndarray 填充后的网格数据 """ print("正在使用距离扩散方法填充边界空白区域...") # 找到边界内的空白区域 nan_mask = np.isnan(grid_content) within_boundary_nan = nan_mask & mask if not np.any(within_boundary_nan): print("边界内没有空白区域需要填充") return grid_content blank_count = np.sum(within_boundary_nan) print(f"发现 {blank_count} 个边界内的空白点，开始距离扩散填充...") # 找到边界内有效值的点 valid_mask = ~nan_mask & mask if np.sum(valid_mask) == 0: print("警告：边界内没有有效值，无法进行距离扩散") return grid_content # 计算网格分辨率（用于确定最大扩散距离） if max_diffusion_distance is None: # 自动计算：使用网格点之间的平均距离 dx = np.abs(grid_xx[0, 1] - grid_xx[0, 0]) if grid_xx.shape[1] > 1 else 1.0 dy = np.abs(grid_yy[1, 0] - grid_yy[0, 0]) if grid_xx.shape[0] > 1 else 1.0 avg_resolution = (dx + dy) / 2.0 max_diffusion_distance = avg_resolution * 5.0 # 5倍网格分辨率 print(f"自动计算最大扩散距离: {max_diffusion_distance:.6f}") # 准备有效数据点 valid_points = np.column_stack((grid_xx[valid_mask], grid_yy[valid_mask])) valid_values = grid_content[valid_mask] # 准备空白点 blank_points = np.column_stack((grid_xx[within_boundary_nan], grid_yy[within_boundary_nan])) print(f"使用 {len(valid_points)} 个有效点填充 {len(blank_points)} 个空白点...") # 使用向量化计算距离矩阵 distances = cdist(blank_points, valid_points) # 对每个空白点，找到最近的n_neighbors个有效点 n_neighbors = min(n_neighbors, len(valid_points)) # 应用最大扩散距离限制 if max_diffusion_distance > 0: # 只考虑在最大扩散距离内的点 # 对于每个空白点，找到在扩散距离内的最近邻 filled_values = np.full(len(blank_points), np.nan) global_mean = np.nanmean(valid_values) for i in range(len(blank_points)): point_distances = distances[i, :] valid_idx = point_distances <= max_diffusion_distance if np.any(valid_idx): # 找到最近的n_neighbors个点（在扩散距离内） valid_dist = point_distances[valid_idx] valid_vals = valid_values[valid_idx] # 如果有效点数量超过n_neighbors，只取最近的n_neighbors个 if len(valid_dist) > n_neighbors: nearest_idx = np.argpartition(valid_dist, n_neighbors-1)[:n_neighbors] valid_dist = valid_dist[nearest_idx] valid_vals = valid_vals[nearest_idx] # 避免除零 valid_dist = np.maximum(valid_dist, 1e-10) # 计算IDW权重 weights = 1.0 / (valid_dist ** power) weight_sum = np.sum(weights) if weight_sum > 0: # 距离加权平均 filled_values[i] = np.sum(weights * valid_vals) / weight_sum else: filled_values[i] = global_mean else: # 如果该点不在任何有效点的扩散距离内，使用全局平均值 filled_values[i] = global_mean else: # 不使用距离限制，对所有点进行IDW插值（批量处理以提高效率） # 对每个空白点，找到最近的n_neighbors个点 nearest_indices = np.argpartition(distances, n_neighbors-1, axis=1)[:, :n_neighbors] # 批量提取距离和值 nearest_dists = np.take_along_axis(distances, nearest_indices, axis=1) nearest_vals = valid_values[nearest_indices] # 避免除零 nearest_dists = np.maximum(nearest_dists, 1e-10) # 批量计算IDW权重 weights = 1.0 / (nearest_dists ** power) weight_sums = np.sum(weights, axis=1) # 批量计算加权平均值 filled_values = np.sum(weights * nearest_vals, axis=1) / weight_sums # 处理可能的NaN（如果weight_sum为0） nan_mask = np.isnan(filled_values) | (weight_sums == 0) if np.any(nan_mask): filled_values[nan_mask] = np.nanmean(valid_values) # 填充空白点 grid_content[within_boundary_nan] = filled_values # 检查填充结果 filled_count = np.sum(~np.isnan(filled_values)) print(f"距离扩散填充完成：成功填充 {filled_count} / {blank_count} 个空白点") return grid_content def _perform_interpolation(self, points, values, grid_xx, grid_yy): """三级降级插值策略：Kriging → IDW → 最近邻 2026-07-01 重构： - Kriging 优先（自动拟合球形变异函数，不强制 nugget） - Kriging 退化检测：若结果标准差接近 0（纯色图），自动回退 IDW - IDW 作为首选回退：无需拟合变异函数，不会产生纯色图 - scipy linear/nearest 作为最后兜底 """ print(f"插值输入检查：") print(f" - 数据点数量: {len(points)}") print(f" - 数据值范围: {values.min():.4f} - {values.max():.4f}") print(f" - 网格大小: {grid_xx.shape}") print(f" - 坐标系: {self.output_crs}") # 检查数据的有效性 finite_mask = np.isfinite(values) if not np.all(finite_mask): print(f"警告：发现 {np.sum(~finite_mask)} 个无效数据值，将被移除") points = points[finite_mask] values = values[finite_mask] if len(points) < 3: raise ValueError(f"有效数据点不足3个（当前：{len(points)}个）") # ── 策略 0：值域极窄时直接跳过 Kriging ── value_range = float(values.max() - values.min()) value_std = float(np.std(values)) # ═══════════════════════════════════════════════════════════ # 策略 1：Kriging（自动拟合球形变异函数） # ═══════════════════════════════════════════════════════════ kriging_degraded = False if PYKRIGE_AVAILABLE: try: print("正在使用 Kriging 插值（球形模型，自动拟合 nugget）...") grid_x = grid_xx[0, :] grid_y = grid_yy[:, 0] ok = OrdinaryKriging( points[:, 0], points[:, 1], values, variogram_model='spherical', verbose=False, enable_plotting=False, ) z, ss = ok.execute('grid', grid_x, grid_y, backend='loop', n_closest_points=15) grid_content = np.array(z) valid_mask = ~np.isnan(grid_content) valid_count = int(np.sum(valid_mask)) if valid_count > 0: # ★ 退化检测：若插值结果标准差 < 原始数据标准差的 5%，判定为纯色图 kriging_std = float(np.nanstd(grid_content)) degradation_ratio = kriging_std / max(value_std, 1e-12) print(f"Kriging 完成: 有效点={valid_count}/{grid_content.size}, " f"输出std={kriging_std:.6f}, 退化比={degradation_ratio:.3f}") if degradation_ratio < 0.05 and value_range > 1e-8: print(f"⚠ Kriging 严重退化（输出 std/输入 std={degradation_ratio:.3f}<5%），" f"判定为纯色图，回退 IDW") kriging_degraded = True else: print(f"Kriging 通过退化检测，直接使用") return grid_content else: print("Kriging 结果全为 NaN，回退") kriging_degraded = True except Exception as e: print(f"Kriging 失败: {e}") kriging_degraded = True else: print("pykrige 未安装，跳过 Kriging") kriging_degraded = True # ═══════════════════════════════════════════════════════════ # 策略 2：IDW（反距离权重）— 不需要拟合变异函数，绝不纯色 # ═══════════════════════════════════════════════════════════ if kriging_degraded: try: print("正在使用 IDW 插值（反距离权重, power=2, neighbors=15）...") grid_content = self._idw_interpolation( points, values, grid_xx, grid_yy, power=2, n_neighbors=min(15, len(points)), ) valid_count = int(np.sum(~np.isnan(grid_content))) if valid_count > 0: idw_std = float(np.nanstd(grid_content)) print(f"IDW 完成: 有效点={valid_count}/{grid_content.size}, 输出std={idw_std:.6f}") if idw_std > 0: return grid_content else: print("IDW std=0（所有输入值完全相同），结果可用") return grid_content else: print("IDW 结果全为 NaN，回退 scipy 插值") except Exception as e: print(f"IDW 失败: {e}，回退 scipy 插值") # ═══════════════════════════════════════════════════════════ # 策略 3：scipy 线性插值 + 最近邻填充（最终兜底） # ═══════════════════════════════════════════════════════════ try: print("正在尝试 scipy 线性插值...") grid_content = griddata( points, values, (grid_xx, grid_yy), method='linear', fill_value=np.nan ) valid_count = int(np.sum(~np.isnan(grid_content))) print(f"线性插值: 有效点={valid_count}/{grid_content.size}") if valid_count > 0: nan_count = int(np.sum(np.isnan(grid_content))) if nan_count > 0: print(f"用最近邻填充 {nan_count} 个 NaN...") grid_nearest = griddata( points, values, (grid_xx, grid_yy), method='nearest' ) grid_content[np.isnan(grid_content)] = grid_nearest[np.isnan(grid_content)] return grid_content except Exception as e: print(f"线性插值失败: {e}") # ═══════════════════════════════════════════════════════════ # 策略 4：最近邻（绝对兜底） # ═══════════════════════════════════════════════════════════ print("执行最近邻插值（最终兜底）...") grid_content = griddata( points, values, (grid_xx, grid_yy), method='nearest' ) if np.sum(~np.isnan(grid_content)) == 0: raise ValueError("所有插值方法均失败") return grid_content @staticmethod def _idw_interpolation(points, values, grid_xx, grid_yy, power=2, n_neighbors=15): """IDW（反距离权重）插值 — 不需拟合模型，绝不产生纯色图。 Parameters: points: (N, 2) 采样点坐标 values: (N,) 采样点值 grid_xx, grid_yy: meshgrid 网格 power: 距离衰减幂参数（默认 2） n_neighbors: 每个网格点参考的最近邻数量 """ from scipy.spatial import cKDTree grid_shape = grid_xx.shape grid_flat = np.column_stack((grid_xx.ravel(), grid_yy.ravel())) values_flat = values.ravel() tree = cKDTree(points) k = min(n_neighbors, len(points)) distances, indices = tree.query(grid_flat, k=k) # 防止距离为 0 的除零 distances = np.maximum(distances, 1e-12) weights = 1.0 / (distances ** power) # 归一化权重 weights /= weights.sum(axis=1, keepdims=True) # 加权求和 neighbor_vals = values_flat[indices] if k == 1 else values_flat[indices] if k == 1: result = neighbor_vals else: result = np.sum(weights * neighbor_vals, axis=1) return result.reshape(grid_shape) def read_csv_data(self, csv_file, uncertainty_col=None): """ 读取CSV文件并进行坐标转换 Parameters: ----------- csv_file : str CSV文件路径 uncertainty_col : str, optional 不确定性数据列名。如果为None，将自动检测包含'variance'、'uncertainty'、'std'、'sigma'、'var'、'mc_dropout'的列 Returns: -------- gdf : gpd.GeoDataFrame 包含坐标和含量数据的GeoDataFrame，如果找到不确定性列，会包含'uncertainty'列 """ print("正在读取CSV文件...") df = pd.read_csv(csv_file, encoding='utf-8') # 假设前三列分别是经度、纬度、含量 if df.shape[1] < 3: raise ValueError("CSV文件必须至少包含3列：经度、纬度、含量") # 获取列名 lon_col = df.columns[0] lat_col = df.columns[1] content_col = df.columns[2] print(f"检测到列名：经度({lon_col})，纬度({lat_col})，含量({content_col})") # 自动检测不确定性列 if uncertainty_col is None: uncertainty_keywords = ['variance', 'uncertainty', 'std', 'sigma', 'var', 'mc_dropout'] for col in df.columns: col_lower = col.lower() if any(keyword in col_lower for keyword in uncertainty_keywords): uncertainty_col = col print(f"自动检测到不确定性列: {uncertainty_col}") break # 坐标转换 print(f"正在进行坐标转换: {self.input_crs} -> {self.output_crs}") transformed_x, transformed_y = self.transformer.transform( df[lon_col].values, df[lat_col].values ) # 创建GeoDataFrame geometry = [Point(x, y) for x, y in zip(transformed_x, transformed_y)] gdf = gpd.GeoDataFrame( df, geometry=geometry, crs=self.output_crs ) gdf['proj_x'] = transformed_x gdf['proj_y'] = transformed_y gdf['content'] = df[content_col] # 如果找到不确定性列，添加到GeoDataFrame if uncertainty_col and uncertainty_col in df.columns: gdf['uncertainty'] = df[uncertainty_col].values print(f"已加载不确定性数据列: {uncertainty_col}") print(f"不确定性值范围: {gdf['uncertainty'].min():.4f} - {gdf['uncertainty'].max():.4f}") print(f"成功读取 {len(gdf)} 个数据点") return gdf def read_boundary_shapefile(self, shp_file): """读取边界/掩膜文件（同时支持矢量 .shp 与栅格 .dat/.bsq/.tif/.tiff）。 - .shp → gpd.read_file 读取矢量边界（保持原行为） - .dat/.bsq/.tif/.tiff/.img → rasterio 读取栅格水掩膜 → rasterio.features.shapes 矢量化成水体多边形 → gpd.GeoDataFrame 返回 下游 create_interpolation_grid / create_content_map / visualize_raster 始终接收 GeoDataFrame，无需任何改动。 """ print("正在读取边界/掩膜文件...") suffix = Path(shp_file).suffix.lower() if suffix in (".shp",): boundary = gpd.read_file(shp_file) elif suffix in (".dat", ".bsq", ".tif", ".tiff", ".img"): boundary = self._raster_to_boundary_gdf(shp_file) else: raise ValueError( f"不支持的边界/掩膜文件格式: {suffix}（仅支持 .shp / .dat / .bsq / .tif / .img）" ) if len(boundary) == 0: raise ValueError( f"边界/掩膜 {shp_file} 矢量化后为空（栅格格式请确认 .dat 包含水体像元 > 0）" ) # 确保边界文件使用目标投影坐标系 if boundary.crs is not None and boundary.crs != self.output_crs: print(f"正在转换边界/掩膜坐标系到 {self.output_crs}...") boundary = boundary.to_crs(self.output_crs) print(f"边界/掩膜文件包含 {len(boundary)} 个要素") return boundary def _raster_to_boundary_gdf(self, raster_path): """把栅格二值水掩膜（.dat/.bsq/.tif/.tiff）矢量化成水体多边形 GeoDataFrame。 修复 Step 11 接收 step1 产出 .dat 水掩膜的兼容性： - rasterio.open 读 band 1（0=非水, 1/任意正数=水） - rasterio.features.shapes 矢量化成多边形 - 收集所有 val=1 的多边形 → gpd.GeoDataFrame """ try: from shapely.geometry import shape as _shapely_shape except ImportError as e: raise ImportError( "栅格掩膜矢量化需要 shapely（geopandas 自带）。原始错误: " + str(e) ) with rasterio.open(raster_path) as src: data = src.read(1) 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: raise ValueError( f"栅格掩膜 {raster_path} 中无水体像元（>0），无法矢量化" ) # 矢量化：shapes 返回 (geom_dict, value) 迭代器 geoms = [] for geom_dict, val in shapes(mask_uint8, mask=mask_uint8.astype(bool), transform=transform): if int(val) == 1: geoms.append(_shapely_shape(geom_dict)) if not geoms: raise ValueError(f"栅格掩膜 {raster_path} 矢量化后无有效水体多边形") gdf = gpd.GeoDataFrame(geometry=geoms, crs=crs) print( f"栅格掩膜 {Path(raster_path).name} 矢量化完成: " f"{int(mask_uint8.sum())} 个水体像元 → {len(gdf)} 个多边形 (CRS={crs})" ) return gdf def _identify_edge_points(self, points_gdf): """ 识别边缘采样点（使用凸包方法） Parameters: ----------- points_gdf : gpd.GeoDataFrame 采样点GeoDataFrame Returns: -------- edge_indices : np.ndarray 边缘点的索引数组 """ print("正在识别边缘采样点...") # 获取所有点的坐标 points = np.column_stack((points_gdf['proj_x'].values, points_gdf['proj_y'].values)) if len(points) < 3: print("警告：采样点数量少于3个，无法识别边缘点") return np.array([]) try: # 使用凸包识别边缘点 hull = ConvexHull(points) edge_indices = hull.vertices print(f"识别到 {len(edge_indices)} 个边缘采样点（共 {len(points)} 个点）") return edge_indices except Exception as e: print(f"识别边缘点时出错: {e}，将使用所有点作为边缘点") return np.arange(len(points)) def _expand_edge_points(self, points_gdf, boundary_gdf=None, resolution=100, expand_ratio=0.05): """ 对边缘采样点进行外扩处理，外扩到整个图像的边界（包括外扩后的边界） 按照指定的间距（resolution）生成外扩点，铺满整个画面。 ★★★ Plan C：boundary_gdf 可选（None = 不依赖水域掩膜，纯采样点自然扩展）★★★ Parameters: ----------- points_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() # 获取所有点的坐标 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]) # 存储新添加的点 new_points_list = [] new_data_list = [] # 对每个边缘点进行外扩 for edge_idx in edge_indices: edge_point = points[edge_idx] # 计算从中心到边缘点的方向向量 direction = edge_point - center distance_to_center = np.linalg.norm(direction) if distance_to_center < 1e-10: # 如果边缘点就是中心点，跳过 continue # 归一化方向向量 direction_unit = direction / distance_to_center # 计算该方向与水域掩膜边界的交点 # 使用射线法：从边缘点沿方向延伸，找到与边界框的交点 max_distance = 0 # 检查与四个边界的交点（使用整个图像的范围，包括外扩后的边界） # 上边界 (y = image_maxy) if direction_unit[1] > 1e-10: # 向上 t = (image_maxy - edge_point[1]) / direction_unit[1] if t > 0: intersect_x = edge_point[0] + direction_unit[0] * t if image_minx <= intersect_x <= image_maxx: max_distance = max(max_distance, t) # 下边界 (y = image_miny) if direction_unit[1] < -1e-10: # 向下 t = (image_miny - edge_point[1]) / direction_unit[1] if t > 0: intersect_x = edge_point[0] + direction_unit[0] * t if image_minx <= intersect_x <= image_maxx: max_distance = max(max_distance, t) # 右边界 (x = image_maxx) if direction_unit[0] > 1e-10: # 向右 t = (image_maxx - edge_point[0]) / direction_unit[0] if t > 0: intersect_y = edge_point[1] + direction_unit[1] * t if image_miny <= intersect_y <= image_maxy: max_distance = max(max_distance, t) # 左边界 (x = image_minx) if direction_unit[0] < -1e-10: # 向左 t = (image_minx - edge_point[0]) / direction_unit[0] if t > 0: intersect_y = edge_point[1] + direction_unit[1] * t if image_miny <= intersect_y <= image_maxy: max_distance = max(max_distance, t) # 如果找到了边界交点，按照resolution间距创建外扩点 if max_distance > 1e-10: # 从边缘点到边界的距离 distance_to_boundary = max_distance # 计算需要生成的外扩点数量（按照resolution间距） # 使用ceil确保能铺满到边界 n_points = int(np.ceil(distance_to_boundary / resolution)) # 从边缘点开始，按照resolution间距生成点，直到边界 for i in range(1, n_points + 1): # 计算外扩距离（从边缘点开始，按照resolution间距） expand_distance = i * resolution # 如果超过边界距离，使用边界距离作为最后一个点 if expand_distance >= distance_to_boundary: expand_distance = distance_to_boundary # 计算新点位置 new_point = edge_point + direction_unit * expand_distance # 确保新点在图像范围内（包括外扩后的边界） new_point[0] = np.clip(new_point[0], image_minx, image_maxx) new_point[1] = np.clip(new_point[1], image_miny, image_maxy) # 创建新点的数据（复制边缘点的所有属性） new_row = points_gdf.iloc[edge_idx].copy() new_row['proj_x'] = new_point[0] new_row['proj_y'] = new_point[1] # 更新geometry new_row['geometry'] = Point(new_point[0], new_point[1]) new_points_list.append(new_point) new_data_list.append(new_row) # 如果已经到达边界，停止生成 if expand_distance >= distance_to_boundary: break # 合并原始点和外扩点 if len(new_data_list) > 0: # 创建新点的GeoDataFrame expanded_gdf = gpd.GeoDataFrame(new_data_list, crs=points_gdf.crs) # 合并原始点和外扩点（使用gpd.concat以确保geometry列正确处理） result_gdf = gpd.GeoDataFrame(pd.concat([points_gdf, expanded_gdf], ignore_index=True), crs=points_gdf.crs) print(f"外扩完成：原始点 {len(points_gdf)} 个，边缘点 {len(edge_indices)} 个，" f"新增外扩点 {len(new_data_list)} 个（间距 {resolution}），总计 {len(result_gdf)} 个点") 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 else: print("未生成外扩点，返回原始点集") return points_gdf.copy() 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%），用于从采样点范围外扩出图像边界。 use_distance_diffusion : bool, default=True 是否使用距离扩散方法填充边界空白区域（仅在 boundary_gdf 有值时生效）。 max_diffusion_distance : float, optional 最大扩散距离（单位与坐标相同）。如果为None，自动计算为网格分辨率的5倍。 diffusion_power : float, default=2 距离扩散的IDW幂参数，值越大，距离衰减越快。 diffusion_n_neighbors : int, default=15 距离扩散使用的最近邻点数。 Returns: -------- grid_xx, grid_yy, grid_content, bounds : tuple """ print("正在创建插值网格...") # ── 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})") # 计算范围大小并外扩 width = maxx - minx height = maxy - miny expand_x = width * expand_ratio expand_y = height * expand_ratio minx -= expand_x maxx += expand_x miny -= expand_y maxy += expand_y print(f"外扩后边界范围: X({minx:.6f} - {maxx:.6f}), Y({miny:.6f} - {maxy:.6f})") print(f"外扩比例: {expand_ratio * 100:.1f}%") if self.output_crs == 'EPSG:4326': print(f"区域尺寸: 宽度={width:.6f}°, 高度={height:.6f}°") 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 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) total_grid_cells = grid_xx.size print(f"网格大小: {grid_xx.shape[1]} x {grid_xx.shape[0]} (宽 x 高) = {total_grid_cells:,} 个网格点") # ★ 大网格自动拦截：超过阈值时自动提升分辨率以保证可接受的计算时间 _MAX_GRID = 500000 # Kriging 在此规模约需 30-120s _URGENT_MAX_GRID = 2000000 # 绝对上限，超过直接拒绝 if total_grid_cells > _URGENT_MAX_GRID: raise ValueError( f"网格点数 {total_grid_cells:,} 超过绝对上限 {_URGENT_MAX_GRID:,}。" f"当前分辨率下 Kriging 预计需要数小时。" f"请将「空间分布图」面板的「分辨率」增大到至少 " f"{int(resolution * (total_grid_cells / _MAX_GRID) ** 0.5)}m 后重试。" ) elif total_grid_cells > _MAX_GRID: # 自动降分辨率 scale = (total_grid_cells / _MAX_GRID) ** 0.5 new_res = int(resolution * scale) new_nx = max(100, int(grid_points_x / scale)) new_ny = max(100, int(grid_points_y / scale)) print(f"⚠ 网格点数 {total_grid_cells:,} 超过推荐上限 {_MAX_GRID:,}，" f"自动将分辨率从 {resolution}m 提升到 ~{new_res}m " f"(网格缩减到 ≈{new_nx}×{new_ny} = {new_nx*new_ny:,} 点)") grid_x = np.linspace(minx, maxx, new_nx) grid_y = np.linspace(miny, maxy, new_ny) grid_xx, grid_yy = np.meshgrid(grid_x, grid_y) total_grid_cells = grid_xx.size print(f"调整后网格: {grid_xx.shape[1]} x {grid_xx.shape[0]} = {total_grid_cells:,} 个网格点") if grid_xx.shape[0] < 2 or grid_xx.shape[1] < 2: 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"含量值范围: {values.min():.4f} - {values.max():.4f}") if len(points) < 3: raise ValueError("插值需要至少3个数据点") self._check_point_distribution(points) # 执行插值 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("正在提取边界边缘值并填充边界外区域...") 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, power=diffusion_power, 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') grid_content[within_boundary_nan] = filled_values print("边界内填充完成") 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 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}") 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') 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 print(f"使用边界内平均值 {mean_in_boundary:.4f} 填充边界外") print("整个画面已铺满，边界外区域已用边缘值填充") final_check_nan = np.isnan(grid_content) & mask if np.any(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 else: global_mean = np.nanmean(grid_content) 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}") if valid_count == 0: raise ValueError("边界掩膜后没有有效数据点，请检查数据点是否在边界范围内") if valid_count < 4: print("警告：有效数据点很少，可能影响绘图效果") 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 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 TIF正射底图文件路径。如果提供，将在水域掩膜外显示底图 cmap : str, default='viridis' 含量数据的颜色映射 """ print("正在生成含量图...") # 检查网格数据 print(f"网格形状: {grid_content.shape}") # 创建边界掩膜（用于绘图时只显示边界内） print("创建边界掩膜用于绘图...") # ── 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: raise ValueError("没有有效的插值数据用于绘图") # 计算数据统计 valid_values = grid_content[valid_data] print( f"插值结果统计: 最小值={valid_values.min():.4f}, 最大值={valid_values.max():.4f}, 平均值={valid_values.mean():.4f}") print(f"有效数据点数量: {np.sum(valid_data)} / {grid_content.size}") # 检查数据范围 data_range = valid_values.max() - valid_values.min() print(f"数据范围: {data_range:.6f}") if data_range == 0: print("警告：所有数据值都相同，将使用单一颜色显示") # 创建图形 fig, ax = plt.subplots(figsize=(12, 10)) # 如果提供了底图，先绘制底图（在水域掩膜外） if base_map_tif is not None: try: print(f"正在加载底图: {base_map_tif}") self._add_base_map(ax, base_map_tif, bounds, mask, grid_xx, grid_yy, boundary_gdf) print("底图加载成功") except Exception as e: print(f"加载底图失败: {e}，将跳过底图显示") # 设置颜色映射参数 im = None try: if data_range > 0: # 设置颜色范围，确保有足够的对比度 vmin = valid_values.min() vmax = valid_values.max() # 如果范围很小，稍微扩展一下以增加对比度 if data_range < 1e-6: center = valid_values.mean() expansion = max(abs(center) * 0.01, 1e-6) # 扩展1%或最小值 vmin = center - expansion vmax = center + expansion print(f"颜色映射范围: {vmin:.6f} - {vmax:.6f}") # 方法1：尝试使用contourf try: print("尝试使用contourf绘制...") # 使用掩膜数组：边界外的数据被掩膜掉，只显示边界内 # mask已经在前面创建好了 masked_data = np.ma.masked_where(~mask, grid_content) # 创建更多等级数以获得更平滑的颜色过渡 levels = np.linspace(vmin, vmax, 100) # 创建100个等级以获得平滑效果 im = ax.contourf(grid_xx, grid_yy, masked_data, levels=levels, cmap=cmap, alpha=0.9, vmin=vmin, vmax=vmax, extend='both') print("contourf绘制成功") # 可选择性添加等值线（默认不添加，以保持平滑效果） # 如果需要等值线，可以取消注释下面的代码 # try: # contour_levels = np.linspace(vmin, vmax, 11) # contours = ax.contour(grid_xx, grid_yy, grid_content, # levels=contour_levels, colors='white', # alpha=0.3, linewidths=0.5) # ax.clabel(contours, inline=True, fontsize=8, fmt='%.3f') # print("等值线添加成功") # except Exception as e: # print(f"等值线绘制失败: {e}") except Exception as e: print(f"contourf失败: {e}") # 方法2：使用pcolormesh try: print("尝试使用pcolormesh绘制...") # 使用掩膜数组：边界外的数据被掩膜掉，只显示边界内 # mask已经在前面创建好了 masked_data = np.ma.masked_where(~mask, grid_content) im = ax.pcolormesh(grid_xx, grid_yy, masked_data, cmap=cmap, alpha=0.9, vmin=vmin, vmax=vmax, shading='gouraud') # 使用gouraud平滑着色 print("pcolormesh绘制成功") except Exception as e2: print(f"pcolormesh也失败: {e2}") raise e2 else: # 所有值相同的情况 print("使用单一颜色填充（所有值相同）") # 创建一个简单的填充 single_value = valid_values[0] im = ax.contourf(grid_xx, grid_yy, grid_content, levels=[single_value - 0.001, single_value + 0.001], cmap=cmap, alpha=0.8) except Exception as e: print(f"主要绘图方法失败，尝试备选方案: {e}") # 备选方案1：imshow try: print("尝试使用imshow...") # 处理NaN值 display_data = grid_content.copy() nan_mask = np.isnan(display_data) if np.any(nan_mask): # 用平均值填充NaN display_data[nan_mask] = valid_values.mean() if data_range > 0: vmin = valid_values.min() vmax = valid_values.max() im = ax.imshow(display_data, extent=[grid_xx.min(), grid_xx.max(), grid_yy.min(), grid_yy.max()], cmap=cmap, alpha=0.8, origin='lower', vmin=vmin, vmax=vmax, aspect='auto') else: im = ax.imshow(display_data, extent=[grid_xx.min(), grid_xx.max(), grid_yy.min(), grid_yy.max()], cmap=cmap, alpha=0.8, origin='lower', aspect='auto') print("imshow绘制成功") except Exception as e2: print(f"imshow也失败: {e2}") # 备选方案2：散点图 try: print("尝试使用散点图...") valid_x = grid_xx[valid_data] valid_y = grid_yy[valid_data] valid_z = grid_content[valid_data] if data_range > 0: im = ax.scatter(valid_x, valid_y, c=valid_z, cmap=cmap, alpha=0.8, s=10, vmin=valid_values.min(), vmax=valid_values.max()) else: im = ax.scatter(valid_x, valid_y, c=valid_z, cmap=cmap, alpha=0.8, s=10) print("散点图绘制成功") except Exception as e3: print(f"所有绘图方法都失败: {e3}") raise ValueError("无法生成颜色图，请检查数据") # 绘制边界（黑色）—— 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: try: points_gdf.plot(ax=ax, color='black', markersize=6, alpha=0.7, marker='+', edgecolors='white', linewidth=1) print("采样点绘制成功") except Exception as e: print(f"采样点绘制失败: {e}") # 设置坐标轴标签和格式 # 由于输入是投影坐标系，输出是地理坐标系，始终显示为地理坐标 ax.set_xlabel('经度 (°)', fontsize=12) ax.set_ylabel('纬度 (°)', fontsize=12) # 格式化坐标轴刻度为经纬度格式（保留3位小数） def lon_formatter(x, p): return f'{x:.3f}°' def lat_formatter(x, p): return f'{x:.3f}°' ax.xaxis.set_major_formatter(FuncFormatter(lon_formatter)) ax.yaxis.set_major_formatter(FuncFormatter(lat_formatter)) # 添加格网线 ax.grid(True, linestyle='--', linewidth=0.5, alpha=0.5, color='gray') ax.set_axisbelow(True) # 将格网线放在图层下方 # ax.set_title('含量分布图', fontsize=16, fontweight='bold', pad=20) # 已去除标题 # 添加颜色条 try: if im is not None: cbar = plt.colorbar(im, ax=ax, shrink=0.5, aspect=40, pad=0.02) cbar.set_label('含量值', fontsize=10) # 设置颜色条的刻度 if data_range > 0: tick_values = np.linspace(valid_values.min(), valid_values.max(), 6) cbar.set_ticks(tick_values) cbar.set_ticklabels([f'{val:.3f}' for val in tick_values]) cbar.ax.tick_params(labelsize=8) # 缩小刻度标签字体 print("颜色条添加成功") else: print("警告：无法添加颜色条，im对象为None") except Exception as e: print(f"颜色条添加失败: {e}") # 添加指北针 try: self.add_north_arrow(ax, bounds) except Exception as e: print(f"指北针添加失败: {e}") # 添加比例尺 try: self.add_scale_bar(ax) except Exception as e: print(f"比例尺添加失败: {e}") # 添加图例 try: self.add_legend(ax) except Exception as e: print(f"图例添加失败: {e}") # 设置图形边界（进一步外扩1%以确保边界不完全挨着地图） try: x_range = bounds[2] - bounds[0] y_range = bounds[3] - bounds[1] display_expand = 0.01 # 显示时再外扩1% ax.set_xlim(bounds[0] - x_range * display_expand, bounds[2] + x_range * display_expand) ax.set_ylim(bounds[1] - y_range * display_expand, bounds[3] + y_range * display_expand) except Exception as e: print(f"设置图形边界失败: {e}") # 调整布局 plt.tight_layout() # 保存图片 try: plt.savefig(output_file, dpi=300, bbox_inches='tight', facecolor='white', edgecolor='none') print(f"含量图已保存为：{output_file}") except Exception as e: print(f"图片保存失败: {e}") # 显示图片 try: plt.show() except Exception as e: print(f"图片显示失败: {e}") def add_north_arrow(self, ax, bounds): """添加指北针（右上角）- 画布相对坐标，不依赖数据坐标系。 使用 ax.transAxes 将指北针固定在右上角， 尺寸以点数（points）为单位，与数据坐标系解耦， 无论 UTM 坐标范围多大，指北针始终保持合理大小。 """ # ★★★ 改用画布相对坐标（transAxes）★★★ # (0.88, 0.92) = 右上角，尺寸用 points（72分之一英寸） arrow_ax_x, arrow_ax_y = 0.88, 0.92 radius_pt = 18 # 罗盘半径（磅），由 28 → 18 缩小图元 # 统一在数据坐标系下绘制（transform=ax.transData） # 但 position 由 axes 坐标决定，radius 用固定点数 # 将 axes 坐标转为数据坐标：取右上角 + 偏移 xlim = ax.get_xlim() ylim = ax.get_ylim() # 偏移系数 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 系数 0.6 → 0.42 缩小指北针整体半径 radius = min(dx, dy) * 0.42 # 绘制圆形背景（外圈） circle_outer = patches.Circle( (arrow_x, arrow_y), radius=radius, facecolor='white', edgecolor='black', linewidth=1.5, zorder=10, transform=ax.transData, ) ax.add_patch(circle_outer) # 绘制内圈（装饰） circle_inner = patches.Circle( (arrow_x, arrow_y), radius=radius * 0.7, facecolor='none', edgecolor='gray', linewidth=0.8, linestyle='--', zorder=11, transform=ax.transData, ) ax.add_patch(circle_inner) # 绘制四个方向的刻度线 tick_width = 1.0 # 北方向刻度（主刻度） ax.plot([arrow_x, arrow_x], [arrow_y, arrow_y + radius * 0.85], 'k-', linewidth=tick_width * 2, zorder=12) # 南方向刻度 ax.plot([arrow_x, arrow_x], [arrow_y, arrow_y - radius * 0.85], 'k-', linewidth=tick_width, zorder=12) # 东方向刻度 ax.plot([arrow_x, arrow_x + radius * 0.85], [arrow_y, arrow_y], 'k-', linewidth=tick_width, zorder=12) # 西方向刻度 ax.plot([arrow_x, arrow_x - radius * 0.85], [arrow_y, arrow_y], 'k-', linewidth=tick_width, zorder=12) # 绘制次要刻度（45度方向） for angle in [45, 135, 225, 315]: angle_rad = math.radians(angle) x_end = arrow_x + radius * 0.7 * math.cos(angle_rad) y_end = arrow_y + radius * 0.7 * math.sin(angle_rad) ax.plot([arrow_x, x_end], [arrow_y, y_end], 'k-', linewidth=tick_width * 0.5, alpha=0.6, zorder=12) # 绘制指北箭头（三角形，填充） arrow_size = radius * 0.6 arrow_points = np.array([ [arrow_x, arrow_y + radius * 0.9], # 顶点（北） [arrow_x - arrow_size * 0.3, arrow_y + radius * 0.3], # 左下 [arrow_x + arrow_size * 0.3, arrow_y + radius * 0.3] # 右下 ]) arrow_poly = patches.Polygon( arrow_points, facecolor='black', edgecolor='black', linewidth=1.2, zorder=13, transform=ax.transData, ) ax.add_patch(arrow_poly) # 绘制指南箭头（三角形，填充，但较小） south_arrow_size = radius * 0.4 south_arrow_points = np.array([ [arrow_x, arrow_y - radius * 0.6], # 顶点（南） [arrow_x - south_arrow_size * 0.2, arrow_y - radius * 0.2], # 左上 [arrow_x + south_arrow_size * 0.2, arrow_y - radius * 0.2] # 右上 ]) south_arrow_poly = patches.Polygon( south_arrow_points, facecolor='white', edgecolor='black', linewidth=1.0, zorder=13, transform=ax.transData, ) ax.add_patch(south_arrow_poly) # 添加方向标记（N, S, E, W） 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, fontweight='bold', ha='center', va='top', color='black', zorder=14) ax.text(arrow_x + label_offset, arrow_y, 'E', 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, fontweight='bold', ha='right', va='center', color='black', zorder=14) def add_scale_bar(self, ax, scale_x=None, scale_y=None): """添加比例尺 Parameters ---------- ax : matplotlib Axes 绘图坐标轴 scale_x : float, optional X 方向像素分辨率（米），由 visualize_raster 从 src.res 传入。 若传入则直接作为 ScaleBar 的 scale 值，忽略 self.output_crs 判断。 scale_y : float, optional Y 方向像素分辨率（米），同 scale_x。 """ try: if scale_x is not None and scale_y is not None: # visualize_raster 传入真实像素分辨率，直接用米为单位 scalebar = ScaleBar( scale_x, units='m', location='lower left', box_alpha=0.8, color='black', font_properties={'size': 8}, label_loc='bottom', ) ax.add_artist(scalebar) print(f"比例尺添加成功（像素分辨率: {scale_x:.4f} m）") elif self.output_crs == 'EPSG:4326': scalebar = ScaleBar( 111000, # 1度 = 111000米 units='m', location='lower left', box_alpha=0.8, color='black', font_properties={'size': 8}, label_loc='bottom' ) ax.add_artist(scalebar) print("地理坐标系比例尺添加成功") else: scalebar = ScaleBar(1, units='m', location='lower left', box_alpha=0.8, color='black', font_properties={'size': 8}) ax.add_artist(scalebar) print("投影坐标系比例尺添加成功") except Exception as e: print(f"比例尺添加失败: {e}") # 如果matplotlib-scalebar失败，尝试手动添加简单的比例尺 try: self._add_manual_scale_bar(ax) print("手动比例尺添加成功") except Exception as e2: print(f"手动比例尺也失败: {e2}") def _add_manual_scale_bar(self, ax): """手动添加简单的比例尺""" # 获取当前坐标轴的范围 xlim = ax.get_xlim() ylim = ax.get_ylim() # 计算比例尺的位置和长度 x_range = xlim[1] - xlim[0] y_range = ylim[1] - ylim[0] # 比例尺位置（左下角） scale_x = xlim[0] + x_range * 0.05 scale_y = ylim[0] + y_range * 0.1 if self.output_crs == 'EPSG:4326': # 地理坐标系：计算合适的比例尺长度（度） # 选择一个合理的距离，比如1公里、5公里或10公里 distance_km = 5 # 5公里 scale_length_deg = distance_km / 111.0 # 转换为度数 # 绘制比例尺线 ax.plot([scale_x, scale_x + scale_length_deg], [scale_y, scale_y], '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=8, bbox=dict(boxstyle='round,pad=0.3', facecolor='white', alpha=0.8)) else: # 投影坐标系：使用米为单位 # 选择合适的比例尺长度 if x_range > 10000: # 大于10km scale_length = 5000 # 5km scale_text = '5 km' elif x_range > 2000: # 大于2km scale_length = 1000 # 1km scale_text = '1 km' else: # 小于2km scale_length = 500 # 500m scale_text = '500 m' # 绘制比例尺线 ax.plot([scale_x, scale_x + scale_length], [scale_y, scale_y], '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=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): """添加正射底图（在水域掩膜外显示） Parameters: ----------- ax : matplotlib.axes.Axes 绘图轴对象 base_map_tif : str TIF底图文件路径 bounds : np.ndarray 显示范围 [minx, miny, maxx, maxy] mask : np.ndarray 水域掩膜（True表示水域内） grid_xx : np.ndarray 网格X坐标 grid_yy : np.ndarray 网格Y坐标 boundary_gdf : gpd.GeoDataFrame 边界几何数据 """ print("正在读取底图文件...") with rasterio.open(base_map_tif) as src: # 获取底图的坐标系 tif_crs = src.crs tif_bounds = src.bounds print(f"底图坐标系: {tif_crs}") print(f"底图范围: {tif_bounds}") print(f"目标范围: {bounds}") # 检查是否需要投影转换 target_crs = CRS.from_string(self.output_crs) need_reproject = tif_crs != target_crs # 读取底图数据 if need_reproject: print(f"底图坐标系({tif_crs})与目标坐标系({target_crs})不同，正在转换...") # 计算转换后的变换参数和尺寸 transform, width, height = calculate_default_transform( tif_crs, target_crs, src.width, src.height, left=bounds[0], bottom=bounds[1], right=bounds[2], top=bounds[3] ) # 创建目标数组 if src.count == 1: # 单波段 base_map_data = np.zeros((height, width), dtype=src.dtypes[0]) reproject( source=src.read(1), destination=base_map_data, src_transform=src.transform, src_crs=tif_crs, dst_transform=transform, dst_crs=target_crs, resampling=Resampling.bilinear ) else: # 多波段（RGB），取前3个波段 num_bands = min(3, src.count) base_map_data = np.zeros((num_bands, height, width), dtype=src.dtypes[0]) for i in range(num_bands): reproject( source=src.read(i + 1), destination=base_map_data[i], src_transform=src.transform, src_crs=tif_crs, dst_transform=transform, dst_crs=target_crs, resampling=Resampling.bilinear ) # 如果是RGB，转换为(height, width, 3)格式 if num_bands == 3: base_map_data = np.transpose(base_map_data, (1, 2, 0)) # 创建extent用于显示 extent = [bounds[0], bounds[2], bounds[1], bounds[3]] else: # 不需要投影转换，直接读取对应范围的数据 print("底图坐标系与目标坐标系一致，直接读取...") # 计算需要读取的窗口 row_min, col_min = src.index(bounds[0], bounds[3]) # 左上角 row_max, col_max = src.index(bounds[2], bounds[1]) # 右下角 # 确保索引在有效范围内 row_min = max(0, row_min) row_max = min(src.height, row_max + 1) col_min = max(0, col_min) col_max = min(src.width, col_max + 1) window = windows.Window.from_slices( (row_min, row_max), (col_min, col_max) ) # 读取数据 if src.count == 1: base_map_data = src.read(1, window=window) else: # 多波段，取前3个波段 num_bands = min(3, src.count) base_map_data = src.read(list(range(1, num_bands + 1)), window=window) if num_bands == 3: # 转换为(height, width, 3)格式 base_map_data = np.transpose(base_map_data, (1, 2, 0)) # 计算extent window_transform = windows.transform(window, src.transform) left = window_transform[2] top = window_transform[5] right = left + window_transform[0] * base_map_data.shape[1] bottom = top + window_transform[4] * base_map_data.shape[0] # 确保extent不超过bounds extent = [ max(bounds[0], left), min(bounds[2], right), max(bounds[1], bottom), min(bounds[3], top) ] # 将底图数据缩放到网格大小以便显示 # 创建底图的显示掩膜：只在边界外显示 print("正在创建底图显示掩膜...") # 创建底图网格（与显示范围对齐） base_map_height, base_map_width = base_map_data.shape[:2] # 性能优化：如果底图分辨率过高，进行降采样以提高处理速度 # 限制最大边长为2000像素（保持足够清晰度的同时提高速度） max_display_size = 2000 scale_factor = 1.0 if max(base_map_height, base_map_width) > max_display_size: scale_factor = max_display_size / max(base_map_height, base_map_width) new_height = int(base_map_height * scale_factor) new_width = int(base_map_width * scale_factor) print( f"底图分辨率较高 ({base_map_width}x{base_map_height})，降采样到 {new_width}x{new_height} 以提高速度") # 使用scipy的zoom进行降采样 if base_map_data.ndim == 2: base_map_data = ndimage.zoom(base_map_data, scale_factor, order=1) else: base_map_data = ndimage.zoom(base_map_data, (scale_factor, scale_factor, 1), order=1) base_map_height, base_map_width = base_map_data.shape[:2] # 更新extent以匹配新的分辨率 extent_width = extent[1] - extent[0] extent_height = extent[3] - extent[2] extent = [ extent[0], extent[0] + extent_width, extent[2], extent[2] + extent_height ] # 使用rasterio的geometry_mask快速生成掩膜（比创建大量Point对象快得多） # 创建底图的变换矩阵 if need_reproject: # 如果进行了投影转换，使用计算得到的transform base_map_transform = transform else: # 如果没有投影转换，使用窗口变换 base_map_transform = window_transform # 如果进行了降采样，需要调整transform if scale_factor < 1.0: # 调整transform以适应新的分辨率 # rasterio的transform是6元素tuple或Affine对象，需要调整像素大小 # 获取transform的6个参数 (a, b, c, d, e, f) # 其中a和e是像素大小，需要除以scale_factor if Affine is not None: # 获取6个参数 if hasattr(base_map_transform, '__iter__') and len(base_map_transform) == 6: a, b, c, d, e, f = base_map_transform else: a, b, c, d, e, f = base_map_transform[0], base_map_transform[1], base_map_transform[2], \ base_map_transform[3], base_map_transform[4], base_map_transform[5] # 创建新的transform，调整像素大小（a和e是像素大小） base_map_transform = Affine(a / scale_factor, b, c, d, e / scale_factor, f) else: # 降级方案：使用tuple if hasattr(base_map_transform, '__iter__') and len(base_map_transform) == 6: a, b, c, d, e, f = base_map_transform else: a, b, c, d, e, f = base_map_transform[0], base_map_transform[1], base_map_transform[2], \ base_map_transform[3], base_map_transform[4], base_map_transform[5] base_map_transform = (a / scale_factor, b, c, d, e / scale_factor, f) # 调试信息：检查边界数据和底图范围 print(f"底图显示范围 (extent): {extent}") print(f"底图分辨率: {base_map_width}x{base_map_height}") print(f"底图transform: {base_map_transform}") if boundary_gdf is not None and len(boundary_gdf) > 0: boundary_bounds = boundary_gdf.total_bounds print(f"边界数据范围: {boundary_bounds}") print(f"边界数据坐标系: {boundary_gdf.crs}") print(f"边界要素数量: {len(boundary_gdf)}") # 检查边界是否与底图范围重叠 overlap_x = not (boundary_bounds[2] < extent[0] or boundary_bounds[0] > extent[1]) overlap_y = not (boundary_bounds[3] < extent[2] or boundary_bounds[1] > extent[3]) if not (overlap_x and overlap_y): print("警告: 边界数据范围与底图显示范围不重叠！") print(" 将不应用掩膜，显示整个底图") # 创建全True的掩膜（显示所有区域） base_map_mask = np.ones((base_map_height, base_map_width), dtype=bool) else: # 使用geometry_mask生成掩膜（True表示在几何体内，即水域内） # 注意：geometry_mask返回True表示需要掩膜的区域（在几何体内）， # 但我们想要的是边界外的区域（不在几何体内），所以需要反转 try: within_boundary_mask = geometry_mask( boundary_gdf.geometry, out_shape=(base_map_height, base_map_width), transform=base_map_transform, invert=False # False表示掩膜几何体内的区域（水域内） ) # 反转掩膜：True表示边界外（需要显示的区域） base_map_mask = ~within_boundary_mask except Exception as e: print(f"生成掩膜时出错: {e}") print(" 将不应用掩膜，显示整个底图") import traceback traceback.print_exc() # 创建全True的掩膜（显示所有区域） base_map_mask = np.ones((base_map_height, base_map_width), dtype=bool) else: print("警告: 边界数据为空，将不应用掩膜，显示整个底图") # 创建全True的掩膜（显示所有区域） base_map_mask = np.ones((base_map_height, base_map_width), dtype=bool) # 调试信息：检查掩膜状态 mask_ratio = np.sum(base_map_mask) / base_map_mask.size print( f"底图掩膜状态: 可显示区域占比 {mask_ratio * 100:.2f}% ({np.sum(base_map_mask)}/{base_map_mask.size} 像素)") # 如果掩膜后没有可显示区域，警告并显示整个底图 if mask_ratio == 0.0: print("警告: 掩膜后没有可显示区域，将显示整个底图（不应用掩膜）") base_map_mask = np.ones((base_map_height, base_map_width), dtype=bool) # 归一化数据以便显示（如果是数值型） # 注意：先归一化整个数据，再应用掩膜，这样可以保证归一化范围正确 if base_map_data.dtype != np.uint8: if base_map_data.ndim == 2: # 单波段：归一化到0-1 # 使用整个数据集的范围进行归一化（不仅仅是掩膜区域） data_min = np.nanmin(base_map_data) data_max = np.nanmax(base_map_data) print(f"底图数据范围: [{data_min}, {data_max}], dtype: {base_map_data.dtype}") if data_max > data_min: # 先归一化整个数组 base_map_normalized = (base_map_data - data_min) / (data_max - data_min) # 然后应用掩膜：只显示边界外的区域 base_map_display = np.ma.masked_where(~base_map_mask, base_map_normalized) else: # 如果数据范围无效，创建全0的掩膜数组 print("警告: 底图数据范围无效，所有值相同") base_map_display = np.ma.masked_where(~base_map_mask, np.zeros_like(base_map_data)) else: # RGB：每个波段单独归一化 base_map_normalized = base_map_data.copy().astype(np.float32) for i in range(base_map_data.shape[2]): band_data = base_map_data[:, :, i] data_min = np.nanmin(band_data) data_max = np.nanmax(band_data) print(f"底图波段 {i} 数据范围: [{data_min}, {data_max}]") if data_max > data_min: # 归一化整个波段 base_map_normalized[:, :, i] = (band_data - data_min) / (data_max - data_min) else: print(f"警告: 底图波段 {i} 数据范围无效，所有值相同") base_map_normalized[:, :, i] = np.zeros_like(band_data) # 应用掩膜：只显示边界外的区域 mask_3d = np.broadcast_to(~base_map_mask[..., np.newaxis], base_map_data.shape) base_map_display = np.ma.masked_where(mask_3d, base_map_normalized) else: # uint8类型：直接使用，但可能需要归一化到0-1用于imshow if base_map_data.ndim == 2: # 单波段：uint8通常已经是0-255范围，归一化到0-1 base_map_normalized = base_map_data.astype(np.float32) / 255.0 base_map_display = np.ma.masked_where(~base_map_mask, base_map_normalized) else: # RGB：uint8归一化到0-1 base_map_normalized = base_map_data.astype(np.float32) / 255.0 mask_3d = np.broadcast_to(~base_map_mask[..., np.newaxis], base_map_data.shape) base_map_display = np.ma.masked_where(mask_3d, base_map_normalized) # 检查归一化后的数据范围 if isinstance(base_map_display, np.ma.MaskedArray): valid_data = base_map_display[~base_map_display.mask] if len(valid_data) > 0: print( f"归一化后有效数据范围: [{np.nanmin(valid_data):.3f}, {np.nanmax(valid_data):.3f}], 有效像素数: {len(valid_data)}") else: print("警告: 归一化后没有有效数据显示区域") # 绘制底图 print("正在绘制底图...") # 注意：extent格式为 [left, right, bottom, top] # 对于地理坐标系，y轴通常向上为正，所以使用origin='lower' try: if base_map_data.ndim == 2: # 单波段：使用灰度图 # 确保数据在0-1范围内 if isinstance(base_map_display, np.ma.MaskedArray): # 对于masked array，确保数据范围正确 if np.ma.max(base_map_display) > 1.0 or np.ma.min(base_map_display) < 0.0: base_map_display = np.ma.clip(base_map_display, 0.0, 1.0) else: base_map_display = np.clip(base_map_display, 0.0, 1.0) im = ax.imshow(base_map_display, extent=extent, origin='lower', cmap='gray', alpha=0.8, zorder=0, interpolation='bilinear', vmin=0.0, vmax=1.0) else: # RGB：直接显示 # 确保数据格式正确（需要在0-1范围内） if isinstance(base_map_display, np.ma.MaskedArray): # 对于masked array，确保数据在0-1范围内 if np.ma.max(base_map_display) > 1.0 or np.ma.min(base_map_display) < 0.0: base_map_display = np.ma.clip(base_map_display, 0.0, 1.0) else: base_map_display = np.clip(base_map_display, 0.0, 1.0) # 确保是float32类型，imshow期望0-1范围的float数组 if base_map_display.dtype != np.float32 and base_map_display.dtype != np.float64: base_map_display = base_map_display.astype(np.float32) im = ax.imshow(base_map_display, extent=extent, origin='lower', alpha=0.8, zorder=0, interpolation='bilinear') print(f"底图绘制成功") except Exception as e: print(f"底图绘制出错: {e}") import traceback traceback.print_exc() # 如果绘制失败，至少尝试绘制一个简单的占位图 print("尝试使用备用方法绘制底图...") try: if base_map_data.ndim == 2: # 使用简单的numpy数组，不应用掩膜 simple_display = np.clip(base_map_data.astype(np.float32) / np.nanmax(base_map_data), 0, 1) ax.imshow(simple_display, extent=extent, origin='lower', cmap='gray', alpha=0.5, zorder=0) else: simple_display = np.clip(base_map_data.astype(np.float32) / 255.0, 0, 1) ax.imshow(simple_display, extent=extent, origin='lower', alpha=0.5, zorder=0) print("备用方法绘制成功") except Exception as e2: print(f"备用方法也失败: {e2}") print(f"底图已绘制，显示范围: {extent}") def add_legend(self, ax): """ 添加图例 Parameters: ----------- """ legend_elements = [ # 移除边界标签 # plt.Line2D([0], [0], color='red', linewidth=2, label='边界'), # 移除采样点和等值线图例项，以突出平滑的颜色分布效果 # plt.Line2D([0], [0], marker='+', color='w', markerfacecolor='black', # markersize=8, label='采样点'), ] # 如果图例为空，则不显示图例 if legend_elements: ax.legend(handles=legend_elements, loc='upper left', framealpha=0.9, fontsize=10) # ------------------------------------------------------------------ # Step 14 适配：水色指数 GeoTIFF 可视化（绕过 CSV 插值） # ------------------------------------------------------------------ def visualize_raster( self, raster_tif_path: str, output_file: Optional[str] = None, boundary_shp_path: Optional[str] = None, cmap: Optional[str] = None, nodata_value: float = -9999.0, show_colorbar: bool = True, figsize: Tuple[int, int] = (12, 10), title: Optional[str] = None, alpha: float = 0.9, ) -> str: """直接读取 GeoTIFF 栅格数据，生成水质指数专题图。 适用场景： - WaterIndexProcessor 输出的水色指数 GeoTIFF - Step 14 接收 GeoTIFF 路径后直接可视化（不通过 CSV 插值） Parameters ---------- raster_tif_path : str 水色指数 GeoTIFF 文件路径（由 WaterIndexProcessor 输出） output_file : str, optional 输出图片路径（None → 自动从 GeoTIFF 文件名派生） boundary_shp_path : str, optional 边界 shapefile 路径（None → 纯栅格显示，无水域掩膜裁切） cmap : str, optional 颜色映射（None → 自动从 GeoTIFF 描述或文件名推断） nodata_value : float NoData 标记值（GeoTIFF 中存储的无效值） show_colorbar : bool 是否显示颜色条 figsize : tuple 图形尺寸（英寸） title : str, optional 图形标题（None → 从 GeoTIFF 描述推断或使用文件名） alpha : float 透明度（0-1） Returns ------- str 输出图片路径 """ # ── 始终从路径提取 stem（供后续中文标题和文件派生使用）────────── stem = Path(raster_tif_path).stem # ── 输出路径自动派生（中文文件名）────────────────────────────── if output_file is None: chinese_title = self._get_chinese_title(stem) out_dir = Path(raster_tif_path).parent / 'visualization' out_dir.mkdir(parents=True, exist_ok=True) output_file = str(out_dir / f"{chinese_title}_专题图.png") # ── 读取 GeoTIFF（优先 rasterio，备选 GDAL）────────────────── tif_path = Path(raster_tif_path) if not tif_path.is_file(): raise FileNotFoundError(f"GeoTIFF 文件不存在: {raster_tif_path}") array: Optional[np.ndarray] = None transform: Optional[Any] = None crs_obj: Optional[Any] = None nodata_read: Optional[float] = None desc: str = "" # 方式1：rasterio（推荐，坐标系信息更完整） _src_bounds = None # rasterio 原生边界（优先用于 extent） _src_res = None # rasterio 像素分辨率 (xres, yres) try: with rasterio.open(raster_tif_path) as src: array = src.read(1).astype(np.float64) transform = src.transform crs_obj = src.crs nodata_read = src.nodata desc = src.descriptions[0] if src.descriptions else "" # 保存原生边界和分辨率，供后续 extent/scale_bar 使用 _src_bounds = src.bounds # left, bottom, right, top _src_res = src.res # (xres, yres) # 替换 NoData 为 NaN（用于绘图） nd = nodata_read if nodata_read is not None else nodata_value if nd is not None: array = np.where(array == nd, np.nan, array) else: array = np.where(np.isnan(array), np.nan, array) print(f"[visualize_raster] rasterio 读取成功: {raster_tif_path}") use_rasterio = True except Exception as rio_err: print(f"[visualize_raster] rasterio 失败 ({rio_err})，回退到 GDAL") use_rasterio = False # 方式2：GDAL（备选） if array is None: try: ds = gdal.Open(raster_tif_path, gdal.GA_ReadOnly) if ds is None: raise RuntimeError("GDAL 无法打开文件") array = ds.GetRasterBand(1).ReadAsArray().astype(np.float64) gt = ds.GetGeoTransform() proj = ds.GetProjection() nodata_read = ds.GetRasterBand(1).GetNoDataValue() desc = ds.GetDescription() or "" if nodata_read is not None: array = np.where(array == nodata_read, np.nan, array) else: array = np.where(np.isnan(array), np.nan, array) # 从 GeoTransform 构造仿射变换（用于计算 extent） if gt and gt != (0, 1, 0, 0, 0, 1): if Affine is not None: transform = Affine(gt[1], gt[2], gt[0], gt[4], gt[5], gt[3]) else: transform = None # ★★★ 关键：从 GeoTransform 计算 bounds 和 res ★★★ # gt = (xmin, xres, 0, ymax, 0, yres) xmin_gdal = gt[0] ymax_gdal = gt[3] xres_gdal = gt[1] yres_gdal = gt[5] width_gdal = ds.RasterXSize height_gdal = ds.RasterYSize xmax_gdal = xmin_gdal + width_gdal * xres_gdal ymin_gdal = ymax_gdal + height_gdal * yres_gdal _src_bounds = rasterio.coords.BoundingBox(xmin_gdal, ymin_gdal, xmax_gdal, ymax_gdal) _src_res = (abs(xres_gdal), abs(yres_gdal)) else: transform = None ds = None except Exception as gdal_err: raise RuntimeError( f"无法读取 GeoTIFF（rasterio 和 GDAL 均失败）: {gdal_err}" ) # ── 宽高变量（供 extent 计算和 figsize 保护使用）───────────── w, h = array.shape[1], array.shape[0] # 保存原始宽高：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 (EPSG:4326) # transform 的坐标值如果是数十万～数百万量级，一定是投影坐标（UTM 等）， # 而非 WGS84 经纬度（-180~180）。此时应使用掩膜的 CRS 作为正确参考， # 否则 geometry_mask 的空间坐标与 transform 完全错位 → 全部像元被擦除。 if crs_obj is None: if (boundary_gdf.crs is not None and boundary_gdf.crs.is_projected and transform is not None): # transform 坐标量级检测：UTM 坐标通常在 10^5~10^7 范围 x_magnitude = max(abs(transform.c), abs(transform.c + transform.a * array.shape[1])) if x_magnitude > 180.0: crs_obj = boundary_gdf.crs print(f"[visualize_raster] 栅格 TIF 缺失坐标系，" f"根据 transform 坐标量级 (x≈{x_magnitude:.0f}) " f"判定为投影坐标系，使用掩膜 CRS ({crs_obj.to_epsg()}) " f"({Path(raster_tif_path).name})") 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 if array.size > _MAX_VIZ_PIXELS: step = int(np.ceil(np.sqrt(array.size / _MAX_VIZ_PIXELS))) array = array[::step, ::step] w_downsampled, h_downsampled = array.shape[1], array.shape[0] print(f"[visualize_raster] 极速降采样: {w}×{h} → {w_downsampled}×{h_downsampled} " f"(step={step})，节省内存并加速渲染") w, h = w_downsampled, h_downsampled # ── 从描述推断参数名和 colormap ─────────────────────────────── # 描述格式：Formula_Name|Category|Formula_Type|Formula param_name: Optional[str] = None if desc and '|' in desc: parts = desc.split('|') param_name = parts[0].strip() if len(parts) >= 2: category = parts[1].strip() if not cmap: cmap = self._get_colormap(category) elif not cmap: # 从文件名推断 stem = tif_path.stem param_name = self._extract_param_name(str(tif_path)) cmap = self._get_colormap(param_name) # ── 中文标题（文件名汉化 + 绘图标题）────────────────────────── # 用户显式传入 title 时直接使用；否则用中文映射 chinese_title = self._get_chinese_title(stem) if not title else title # ── 计算空间范围（extent）────────────────────────────────────── # 优先使用 rasterio 原生 bounds，保证坐标轴为真实 UTM 米 # GDAL 回退使用 GeoTransform 计算 if _src_bounds is not None: extent = [ _src_bounds.left, # xmin _src_bounds.right, # xmax _src_bounds.bottom, # ymin _src_bounds.top, # ymax ] # 从 bounds 推导分辨率（取绝对值，正数用于比例尺） scale_x = abs(_src_res[0]) if _src_res else 1.0 scale_y = abs(_src_res[1]) if _src_res else 1.0 elif transform is not None: xmin = transform.c ymax = transform.f xres = transform.a yres = transform.e # ★★★ 必须用原始宽高（w_orig/h_orig）而非降采样后的 w/h ★★★ extent = [xmin, xmin + w_orig * xres, ymax + h_orig * yres, ymax] scale_x = abs(xres) scale_y = abs(yres) else: # 回退到像素索引范围（使用原始尺寸） extent = [0, w_orig, 0, h_orig] scale_x = 1.0 scale_y = 1.0 # ── 准备图形 ───────────────────────────────────────────────── # 尊重用户的 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)) # 1. 明确排除 NaN 以及 nodata_value（与函数参数保持一致） nodata_val = nodata_value valid = array[(~np.isnan(array)) & (array != nodata_val)] if valid.size == 0: raise ValueError("GeoTIFF 中没有有效数据（全部为 NoData）") # 2. 改用更鲁棒的 2%-98% 百分位拉伸（抗偏态分布） vmin = float(np.percentile(valid, 2)) vmax = float(np.percentile(valid, 98)) if (vmax - vmin) < 1e-9: # 防退化：区间过窄 → 取中心 ±1% center = vmin exp = max(abs(center) * 0.01, 1e-9) vmin = center - exp vmax = center + exp print(f"[visualize_raster] P2-P98 拉伸: vmin={vmin:.4f}, vmax={vmax:.4f}，" f"有效像元: {valid.size}/{array.size}") # ── 栅格绘图 ───────────────────────────────────────────────── # 使用 masked array：NaN 区域自动不显示 # 不仅要屏蔽 NaN，如果 array 中还有 nodata_value，也必须 mask 掉，否则出图时背景会被渲染 masked_data = np.ma.masked_where((np.isnan(array)) | (array == nodata_val), array) # 【核心修复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）之前执行， # 防止其他元素的坐标干扰导致轴范围被拉伸成像素坐标系 ax.set_xlim(extent[0], extent[1]) ax.set_ylim(extent[2], extent[3]) # ── 边界描边（直接复用上面的 boundary_gdf_plotted，不再重复读取/矢量化）───────── if boundary_gdf_plotted is not None: try: boundary_gdf_plotted.boundary.plot(ax=ax, color='black', linewidth=1.5) except Exception as e: print(f"[visualize_raster] 边界描边失败: {e}") # ── 坐标轴标签（动态：根据 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=18, fontweight='bold', pad=15) # ── 颜色条（工业级样式：extend 三角 + MaxNLocator 刻度防重叠）───────── if show_colorbar and im is not None: try: 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] 颜色条添加成功") except Exception as e: print(f"[visualize_raster] 颜色条添加失败: {e}") # ── 比例尺 ─────────────────────────────────────────────────── try: self.add_scale_bar(ax, scale_x=scale_x, scale_y=scale_y) except Exception as e: print(f"[visualize_raster] 比例尺添加失败: {e}") # ── 指北针 ─────────────────────────────────────────────────── try: bounds_arr = np.array(extent) self.add_north_arrow(ax, bounds_arr) except Exception as e: print(f"[visualize_raster] 指北针添加失败: {e}") # ── 紧凑布局并保存 ─────────────────────────────────────────── # 给 tight_layout 显式 pad=2.0 防止指北针/比例尺/标题互相重叠 plt.tight_layout(pad=2.0, h_pad=1.5, w_pad=1.5) try: plt.savefig( output_file, dpi=300, bbox_inches='tight', facecolor='white', edgecolor='none', ) print(f"[visualize_raster] ✅ 专题图已保存: {output_file}") except Exception as e: print(f"[visualize_raster] 保存失败: {e}") raise try: plt.show() except Exception: pass plt.close(fig) return output_file # ------------------------------------------------------------------ # 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 # ═══════════════════════════════════════════════════════════════ # ★ 2026-07-01：共享空间上下文 — 63 个 CSV 只算一次网格/掩膜 # ═══════════════════════════════════════════════════════════════ def prepare_shared_context(self, sample_csv: str, shp_file=None, resolution=100, expand_ratio=0.05): """从首个 CSV 预计算所有子进程共用的空间基准数据。 63 个水色指数 CSV 坐标完全一致，以下数据只算一次： - boundary_gdf (水域边界) - grid_xx, grid_yy (插值网格) - mask (水域掩膜布尔矩阵) - bounds (空间范围) 子进程直接从 shared_context 解包复用，跳过 ②③④⑥，直入 Kriging。 Returns: tuple: (grid_xx, grid_yy, mask, bounds, boundary_gdf) """ print(f"[共享上下文] 从 {Path(sample_csv).name} 预计算空间基准...") # ② 读边界（只此一次） if shp_file is None: boundary_gdf = None else: boundary_gdf = self.read_boundary_shapefile(shp_file) # ③ 边缘外扩（只此一次）—— 需要读第一个CSV获取坐标结构 points_gdf = self.read_csv_data(sample_csv) points_gdf = self._expand_edge_points( points_gdf, boundary_gdf, resolution=resolution, expand_ratio=expand_ratio ) # ④ 计算网格几何 if boundary_gdf is None: pts = np.column_stack((points_gdf['proj_x'], points_gdf['proj_y'])) minx, maxx = pts[:, 0].min(), pts[:, 0].max() miny, maxy = pts[:, 1].min(), pts[:, 1].max() else: bnd = boundary_gdf.total_bounds minx, miny, maxx, maxy = bnd width = maxx - minx height = maxy - miny minx -= width * expand_ratio maxx += width * expand_ratio miny -= height * expand_ratio maxy += height * expand_ratio res = resolution / 111000.0 if self.output_crs == 'EPSG:4326' else resolution nx = max(int(width / res), 100) ny = max(int(height / res), 100) grid_x = np.linspace(minx, maxx, nx) grid_y = np.linspace(miny, maxy, ny) grid_xx, grid_yy = np.meshgrid(grid_x, grid_y) bounds = np.array([minx, miny, maxx, maxy]) print(f"[共享上下文] 网格: {nx}×{ny} = {nx*ny} 点") # ⑥ 水域掩膜布尔矩阵（只此一次） mask = None if boundary_gdf is not None: mask_pts = np.column_stack((grid_xx.ravel(), grid_yy.ravel())) mask_gdf = gpd.GeoDataFrame( geometry=[Point(x, y) for x, y in mask_pts], crs=self.output_crs ) mask = mask_gdf.within(boundary_gdf.unary_union).values.reshape(grid_xx.shape) print(f"[共享上下文] 水域掩膜: {int(mask.sum())}/{mask.size} 点在水域内") return (grid_xx, grid_yy, mask, bounds, boundary_gdf) 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, output_format='tif', shared_context=None): """ 主处理函数 ★★★ Plan C：shp_file 现在是可选参数（None = 纯采样点插值，不依赖水域掩膜）★★★ Parameters: ----------- csv_file : str CSV文件路径 shp_file : str, optional 水域掩膜/边界文件路径（.shp / .dat / .bsq / .tif 等）。 shared_context : tuple, optional (2026-07-01 批量优化) 由 prepare_shared_context() 返回的 (grid_xx, grid_yy, mask, bounds, boundary_gdf)。 提供时跳过 读边界/边缘外扩/建网格/算掩膜，直入 Kriging 插值阶段。 ... (其他参数同上) """ try: # 自动识别参数名称并获取colormap if cmap is None: param_name = self._extract_param_name(csv_file) cmap = self._get_colormap(param_name) # 读取采样点数据 points_gdf = self.read_csv_data(csv_file) # ── ★ 快速通道：复用预计算的共享上下文 ── if shared_context is not None: grid_xx, grid_yy, mask, bounds, boundary_gdf = shared_context # ③ 仍需边缘扩展（值相关），但跳过 ②④⑥ points_gdf = self._expand_edge_points( points_gdf, boundary_gdf, resolution=resolution, expand_ratio=expand_ratio ) # ⑤ 直接用共享网格执行 Kriging pts = np.column_stack((points_gdf['proj_x'], points_gdf['proj_y'])) vals = points_gdf['content'].values grid_content = self._perform_interpolation(pts, vals, grid_xx, grid_yy) # ⑥ 复用共享掩膜裁剪 if mask is not None: grid_content[~mask] = np.nan # 边界内 NaN 填充 nan_mask = np.isnan(grid_content) within_nan = nan_mask & mask if np.any(within_nan): valid_m = ~nan_mask & mask if np.sum(valid_m) > 0: v_pts = np.column_stack((grid_xx[valid_m], grid_yy[valid_m])) v_vals = grid_content[valid_m] n_pts = np.column_stack((grid_xx[within_nan], grid_yy[within_nan])) try: from scipy.interpolate import griddata grid_content[within_nan] = griddata( v_pts, v_vals, n_pts, method='nearest' ) except Exception: grid_content[within_nan] = np.nanmean(grid_content[mask]) else: # ── 原有完整流程（单图模式）────────── if shp_file is None: print("[Plan C] shp_file=None，跳过水域掩膜读取") boundary_gdf = None else: boundary_gdf = self.read_boundary_shapefile(shp_file) points_gdf = self._expand_edge_points( points_gdf, boundary_gdf, resolution=resolution, expand_ratio=expand_ratio ) grid_xx, grid_yy, grid_content, bounds = self.create_interpolation_grid( points_gdf, boundary_gdf, resolution, expand_ratio=expand_ratio, use_distance_diffusion=use_distance_diffusion, max_diffusion_distance=max_diffusion_distance, diffusion_power=diffusion_power, diffusion_n_neighbors=diffusion_n_neighbors ) # ── 按 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=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, optional 水域掩膜/边界文件路径。None 时跳过掩膜约束。 output_folder : str, optional 输出文件夹路径。如果为None，将在CSV文件所在文件夹创建'map_output'子文件夹 resolution : int, default=100 网格分辨率（米） show_sample_points : bool, default=False 是否显示采样点 base_map_tif : str, optional TIF正射底图文件路径 use_distance_diffusion : bool, default=True 是否使用距离扩散方法 max_diffusion_distance : float, optional 最大扩散距离 diffusion_power : float, default=2 距离扩散的IDW幂参数 diffusion_n_neighbors : int, default=15 距离扩散使用的最近邻点数 """ print("=" * 60) print("开始批量处理CSV文件") print("=" * 60) # 检查输入文件夹是否存在 if not os.path.isdir(csv_folder): raise ValueError(f"输入文件夹不存在: {csv_folder}") # 获取所有CSV文件 csv_files = glob.glob(os.path.join(csv_folder, "*.csv")) if len(csv_files) == 0: raise ValueError(f"在文件夹 {csv_folder} 中未找到CSV文件") print(f"找到 {len(csv_files)} 个CSV文件") # 创建输出文件夹 if output_folder is None: output_folder = os.path.join(csv_folder, "map_output") if not os.path.exists(output_folder): os.makedirs(output_folder) print(f"创建输出文件夹: {output_folder}") else: print(f"使用输出文件夹: {output_folder}") # 统计信息 success_count = 0 fail_count = 0 failed_files = [] # 批量处理每个CSV文件 for i, csv_file in enumerate(csv_files, 1): print("\n" + "=" * 60) print(f"处理文件 {i}/{len(csv_files)}: {os.path.basename(csv_file)}") print("=" * 60) try: # 生成输出文件名（使用CSV文件名，但扩展名为.png） csv_basename = os.path.splitext(os.path.basename(csv_file))[0] output_file = os.path.join(output_folder, f"{csv_basename}.png") # 处理单个文件（自动识别参数并选择colormap） self.process_data( csv_file=csv_file, shp_file=shp_file, output_file=output_file, resolution=resolution, show_sample_points=show_sample_points, base_map_tif=base_map_tif, use_distance_diffusion=use_distance_diffusion, max_diffusion_distance=max_diffusion_distance, diffusion_power=diffusion_power, diffusion_n_neighbors=diffusion_n_neighbors, cmap=None # 自动识别 ) success_count += 1 print(f"✓ 成功处理: {csv_basename}.png") except Exception as e: fail_count += 1 failed_files.append((os.path.basename(csv_file), str(e))) print(f"✗ 处理失败: {os.path.basename(csv_file)}") print(f" 错误信息: {e}") import traceback traceback.print_exc() # 输出批量处理结果统计 print("\n" + "=" * 60) print("批量处理完成") print("=" * 60) print(f"总文件数: {len(csv_files)}") print(f"成功: {success_count}") print(f"失败: {fail_count}") print(f"输出文件夹: {output_folder}") if failed_files: print("\n失败的文件列表:") for file_name, error in failed_files: print(f" - {file_name}: {error}") return { 'total': len(csv_files), 'success': success_count, 'failed': fail_count, 'output_folder': output_folder, 'failed_files': failed_files } def main(): """主函数 - 使用示例""" # 创建处理器实例 mapper = ContentMapper() # 示例1：处理单个文件 csv_file = r"E:\code\WQ\pipeline_result\tests1\11_12_13_predictions\BGA.csv" # 采样点的预测值 shp_file = r"D:\BaiduNetdiskDownload\yaobao\roi\roi.shp" # 水体边界shapefile路径 output_file = r"E:\code\WQ\pipeline_result\work_dir\11_12_13_predictions\BGA.png" # 输出图片路径 # mapper.process_data( csv_file=csv_file, shp_file=shp_file, output_file=output_file, resolution=30, # 网格分辨率（米），更小的值产生更平滑的效果 show_sample_points=False, # 设置为False以显示平滑的颜色分布，True则显示采样点位置 base_map_tif=None, # 正射底图路径（可选） cmap=None # 自动从文件名或内容中识别参数并选择对应的colormap ) # # 示例2：批量处理文件夹中的所有CSV文件 # csv_folder = r"E:\code\WQ\xiaogujia\使用腰堡模型\predict\TT.csv" # CSV文件所在文件夹 # shp_file = r"E:\code\WQ\xiaogujia\SHP\shp\watemask.shp" # 水体边界shapefile路径 # output_folder = r"E:\code\WQ\xiaogujia\使用腰堡模型\map\TT.png" # 输出文件夹（可选，如果为None则在CSV文件夹下创建map_output） # 批量处理（会自动识别每个CSV文件的参数名称并选择对应的colormap） # result = mapper.process_batch( # csv_folder=csv_folder, # shp_file=shp_file, # output_folder=output_folder, # 如果为None，将在CSV文件夹下创建map_output子文件夹 # resolution=30, # 网格分辨率（米） # show_sample_points=False, # 是否显示采样点 # base_map_tif=None, # 正射底图路径（可选） # ) # # print(f"\n批量处理结果: {result}") if __name__ == "__main__": # 使用示例 print("含量分布图生成器") print("=" * 50) # 如果要直接运行，请取消下面的注释并修改文件路径 main() # 或者交互式使用 # print("使用方法：") # print("1. 准备CSV文件（前两列为WGS84经纬度，第三列为含量数据）") # print("2. 准备边界Shapefile文件") # print("3. 调用以下代码：") # print(""" # mapper = ContentMapper() # mapper.process_data( # csv_file='your_data.csv', # shp_file='your_boundary.shp', # output_file='output_map.png', # resolution=50 # ) # """)