From 320f2f18f28e9081b57629942a3a9db4b300798d Mon Sep 17 00:00:00 2001 From: DXC Date: Wed, 10 Jun 2026 17:13:25 +0800 Subject: [PATCH] =?UTF-8?q?feat(step8):=20=E6=96=B0=E5=A2=9E=E6=B0=B4?= =?UTF-8?q?=E8=89=B2=E6=8C=87=E6=95=B0=E5=8F=8D=E6=BC=94=E6=A8=A1=E5=9D=97?= =?UTF-8?q?=20waterindex=5Finversion=20+=20CSV=20=E5=85=AC=E5=BC=8F?= =?UTF-8?q?=E9=A9=B1=E5=8A=A8=E6=9E=B6=E6=9E=84?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/core/algorithms/waterindex_inversion.py | 22 + .../waterindex_inversion/__init__.py | 646 ++++++++++++++++++ src/utils/generate_pure_water_csv.py | 34 + src/utils/pure_water_iops.csv | 402 +++++++++++ src/utils/water_owt_config.py | 114 ++++ 5 files changed, 1218 insertions(+) create mode 100644 src/core/algorithms/waterindex_inversion.py create mode 100644 src/core/algorithms/waterindex_inversion/__init__.py create mode 100644 src/utils/generate_pure_water_csv.py create mode 100644 src/utils/pure_water_iops.csv create mode 100644 src/utils/water_owt_config.py diff --git a/src/core/algorithms/waterindex_inversion.py b/src/core/algorithms/waterindex_inversion.py new file mode 100644 index 0000000..d7ff9ca --- /dev/null +++ b/src/core/algorithms/waterindex_inversion.py @@ -0,0 +1,22 @@ +# -*- coding: utf-8 -*- +""" +水色指数反演模块（包入口） + +从 waterindex.csv 读取公式，对去耀斑 BSQ 高光谱影像进行全图矩阵运算， +输出带完整坐标信息的 GeoTIFF。 + +公式格式（waterindex.csv）： + - 波长占位符：w{nm}，如 w686, w708, w665 + - 支持混合大小写：w686 / W665 均可 + - 示例：NDCI = (w708 - w665) / (w708 + w665) + +输出： + - GeoTIFF (Float32)，LZW 压缩，带 Tile + - 完整克隆原始 BSQ 的 GeoTransform / Projection / NoData + - Step 14 可直接用 rasterio 读取数组和空间范围 +""" + +# 重新导出 WaterIndexProcessor（向后兼容所有已有 import） +from src.core.algorithms.waterindex_inversion import WaterIndexProcessor + +__all__ = ['WaterIndexProcessor'] diff --git a/src/core/algorithms/waterindex_inversion/__init__.py b/src/core/algorithms/waterindex_inversion/__init__.py new file mode 100644 index 0000000..f58eae9 --- /dev/null +++ b/src/core/algorithms/waterindex_inversion/__init__.py @@ -0,0 +1,646 @@ +# -*- coding: utf-8 -*- +""" +水色指数反演模块 + +直接读取去耀斑高光谱 BSQ 影像，应用 waterindex.csv 中的公式， +输出各水质参数指数的 GeoTIFF 栅格图像。 + +公式格式（waterindex.csv）： + - 波长占位符：w{nm}，如 w686, w708, w665 + - 支持混合大小写：w686 / W665 均可 + - 示例：NDCI = (w708 - w665) / (w708 + w665) + BGA_Am09KBBI = (w686 - w658) / (w686 + w658) + +输出： + - GeoTIFF (Float32)，LZW 压缩，带 Tile + - 完整克隆原始 BSQ 的 GeoTransform / Projection / NoData + - Step 14 可直接用 rasterio 读取进行克里金插值 +""" + +from __future__ import annotations + +import csv +import os +import re +import sys +import time +import traceback +from pathlib import Path +from typing import Any, Callable, Dict, List, Optional, Tuple + +import numpy as np +from osgeo import gdal, osr + +# GDAL 驱动注册 +gdal.UseExceptions() + + +# ------------------------------------------------------------------ +# 公共工具 +# ------------------------------------------------------------------ + +def _get_resource_path(relative_path: str) -> str: + """获取 waterindex.csv 等资源的绝对路径，兼容 PyInstaller 打包。""" + if hasattr(sys, '_MEIPASS'): + base = sys._MEIPASS + else: + base = os.path.abspath( + os.path.join(os.path.dirname(os.path.dirname(__file__)), '..', '..', '..') + ) + return os.path.join(base, relative_path) + + +# ------------------------------------------------------------------ +# WaterIndexProcessor +# ------------------------------------------------------------------ + +class WaterIndexProcessor: + """ + 水色指数处理器 + + 读取 waterindex.csv 中的公式，应用于 BSQ 高光谱影像， + 输出带完整坐标信息的 GeoTIFF 指数图。 + + 核心能力： + - 公式解析：w{nm} 占位符 → 实际波段 2D numpy 数组 + - 矩阵运算：全影像批量计算，无需逐点循环 + - 地理信息保持：克隆原始 BSQ 的 GeoTransform / Projection + - NoData 处理：运算中产生的 NaN/Inf 统一标记为 -9999 + """ + + # 内置安全命名空间（公式 eval 白名单） + _SAFE_NS: Dict[str, Any] = { + 'np': np, + 'nan': np.nan, + 'inf': np.inf, + 'pi': np.pi, + 'e': np.e, + } + + def __init__(self, waterindex_csv_path: Optional[str] = None): + """ + Parameters + ---------- + waterindex_csv_path : str, optional + waterindex.csv 路径。 + 若为 None，尝试从默认位置加载： + 1. src/gui/model/waterindex.csv（开发环境） + 2. _MEIPASS/src/gui/model/waterindex.csv（打包环境） + """ + self.csv_path: Optional[str] = None + self.formulas: List[Dict[str, Any]] = [] + + if waterindex_csv_path: + self.csv_path = waterindex_csv_path + else: + candidates = [ + os.path.join(os.path.dirname(__file__), '..', '..', 'gui', 'model', 'waterindex.csv'), + os.path.join(os.path.dirname(__file__), '..', '..', '..', 'gui', 'model', 'waterindex.csv'), + ] + for p in candidates: + if os.path.isfile(p): + self.csv_path = p + break + + if self.csv_path: + self._parse_csv() + else: + self.formulas = [] + + # ------------------------------------------------------------------ + # 公式加载 + # ------------------------------------------------------------------ + + def _parse_csv(self) -> None: + """解析 waterindex.csv，加载所有公式。""" + if not os.path.isfile(self.csv_path): + raise FileNotFoundError(f"公式配置文件不存在: {self.csv_path}") + + # ★★★ 防止多次调用时公式翻倍叠加 ★★★ + self.formulas.clear() + + with open(self.csv_path, 'r', encoding='utf-8-sig') as f: + reader = csv.DictReader(f) + for row in reader: + self.formulas.append(dict(row)) + + print(f"[WaterIndexProcessor] 加载 {len(self.formulas)} 条公式 ← {self.csv_path}") + + def reload(self, waterindex_csv_path: str) -> None: + """重新加载公式配置文件。""" + self.csv_path = waterindex_csv_path + self._parse_csv() + + # ------------------------------------------------------------------ + # 公式查询 + # ------------------------------------------------------------------ + + def list_formulas(self) -> List[Dict[str, Any]]: + """返回所有公式的列表。""" + return list(self.formulas) + + def list_formula_names(self) -> List[str]: + """返回所有公式名称列表。""" + return [f.get('Formula_Name', '') for f in self.formulas] + + def get_formula(self, name: str) -> Optional[Dict[str, Any]]: + """按名称查找单个公式。""" + for f in self.formulas: + if f.get('Formula_Name', '').strip() == name.strip(): + return f + return None + + def list_categories(self) -> List[str]: + """返回所有公式类别（去重排序）。""" + cats = set() + for f in self.formulas: + c = f.get('Category', '').strip() + if c: + cats.add(c) + return sorted(cats) + + def get_formulas_by_category(self, category: str) -> List[Dict[str, Any]]: + """按类别筛选公式。""" + return [f for f in self.formulas + if f.get('Category', '').strip().lower() == category.strip().lower()] + + # ------------------------------------------------------------------ + # 影像元数据 + # ------------------------------------------------------------------ + + def get_image_metadata(self, bsq_path: str, hdr_path: Optional[str] = None) -> Dict[str, Any]: + """获取影像元数据（GDAL + ENVI HDR 双重保障）。 + + Parameters + ---------- + bsq_path : str + BSQ 影像路径 + hdr_path : str, optional + ENVI HDR 路径（None → 自动构造） + + Returns + ------- + dict + 含 keys: width, height, bands, wavelengths, wavelength_range, + geotransform, projection, driver + """ + meta: Dict[str, Any] = {} + + # 1. GDAL 优先（获取空间信息） + try: + ds = gdal.Open(bsq_path, gdal.GA_ReadOnly) + if ds is not None: + meta['width'] = ds.RasterXSize + meta['height'] = ds.RasterYSize + meta['bands'] = ds.RasterCount + meta['driver'] = ds.GetDriver().ShortName + gt = ds.GetGeoTransform() + proj = ds.GetProjection() + if gt and gt != (0, 1, 0, 0, 0, 1): + meta['geotransform'] = gt + if proj: + meta['projection'] = proj + ds = None + except Exception: + pass + + # 2. HDR 补充波长信息 + if hdr_path is None: + hdr_path = os.path.splitext(bsq_path)[0] + '.hdr' + if not os.path.isfile(hdr_path): + hdr_path_alt = os.path.splitext(bsq_path)[0] + '.HDR' + if os.path.isfile(hdr_path_alt): + hdr_path = hdr_path_alt + + if os.path.isfile(hdr_path): + wl = self._parse_wavelengths_from_hdr(hdr_path) + if wl: + meta['wavelengths'] = wl + if len(wl) >= 2: + meta['wavelength_range'] = f"{wl[0]:.1f}–{wl[-1]:.1f} nm ({len(wl)} 波段)" + elif meta.get('bands', 0) > 0: + meta['wavelength_range'] = f"{meta['bands']} 波段（波长信息缺失）" + + return meta + + @staticmethod + def _parse_wavelengths_from_hdr(hdr_path: str) -> Optional[List[float]]: + """从 ENVI .hdr 文件中解析波长列表。""" + try: + with open(hdr_path, 'r', encoding='utf-8', errors='ignore') as f: + content = f.read() + + # 格式1：wavelength = { 400, 401, ... } + m = re.search(r'wavelength\s*=\s*\{([^}]+)\}', content, re.DOTALL) + if m: + vals = [float(v) for v in re.findall(r'[\d.]+', m.group(1)) if v.strip()] + if vals: + return vals + + # 格式2：逐行罗列 + wavelengths: List[float] = [] + in_wl = False + for line in content.split('\n'): + line = line.strip() + if line.startswith('wavelength'): + in_wl = True + continue + if in_wl: + if line.startswith('{'): + continue + try: + wavelengths.append(float(line)) + except ValueError: + if '}' in line: + in_wl = False + return wavelengths if wavelengths else None + except Exception: + return None + + # ------------------------------------------------------------------ + # 公式解析：w{nm} 占位符 → 实际波段数据 + # ------------------------------------------------------------------ + + def _find_nearest_band_index(self, target_wv: float, + wavelengths: List[float]) -> int: + """找到最接近目标波长的 GDAL 波段索引（1-based）。""" + if not wavelengths: + raise ValueError("波长列表为空，无法匹配波段") + nearest = min(range(len(wavelengths)), + key=lambda i: abs(wavelengths[i] - target_wv)) + return nearest + 1 # GDAL 波段从 1 开始 + + def _parse_formula_wavelengths(self, formula: str) -> List[int]: + """从公式字符串中提取所有波长值（去重，int）。""" + raw = re.findall(r'[wW](\d+)', formula) + seen = set() + result: List[int] = [] + for r in raw: + v = int(r) + if v not in seen: + seen.add(v) + result.append(v) + return result + + def _eval_formula_fast(self, formula: str, + band_data: Dict[int, np.ndarray]) -> Optional[np.ndarray]: + """快速公式求值（预处理后直接 eval）。 + + band_data: {波长int: 2D 数组} + formula 示例: "(w708 - w665) / (w708 + w665)" + """ + # 预处理：w708 → _B708（避免与 Python 关键字冲突） + processed = re.sub(r'[wW](\d+)', r'_B\1', formula) + + # 构建局部变量表：_B708 = band_data[708] + local_vars = {f"_B{wv}": arr for wv, arr in band_data.items()} + local_vars.update(self._SAFE_NS) + + try: + result = eval(processed, {"__builtins__": {}}, local_vars) + return result + except Exception as e: + print(f" ⚠ 公式求值失败 [{formula}]: {e}") + return None + + # ------------------------------------------------------------------ + # 单波段读取（带 NoData 处理） + # ------------------------------------------------------------------ + + @staticmethod + def _read_band_as_float(bsq_path: str, band_idx: int) -> np.ndarray: + """读取 BSQ 指定波段（1-based），返回 float64，NaN 替换 NoData。""" + ds = gdal.Open(bsq_path, gdal.GA_ReadOnly) + if ds is None: + raise RuntimeError(f"无法用 GDAL 打开影像: {bsq_path}") + + band = ds.GetRasterBand(band_idx) + arr = band.ReadAsArray() + nodata = band.GetNoDataValue() + ds = None + + arr = arr.astype(np.float64) + if nodata is not None: + arr = np.where(arr == nodata, np.nan, arr) + + return arr + + # ------------------------------------------------------------------ + # 核心处理：逐公式矩阵运算 + GeoTIFF 输出 + # ------------------------------------------------------------------ + + def process_bsq( + self, + bsq_path: str, + hdr_path: Optional[str] = None, + output_dir: Optional[str] = None, + formula_names: Optional[List[str]] = None, + water_mask: Optional[np.ndarray] = None, + nodata_value: float = -9999.0, + progress_callback: Optional[Callable[[str, float], None]] = None, + ) -> Dict[str, str]: + """逐公式处理 BSQ 影像，输出 GeoTIFF。 + + Parameters + ---------- + bsq_path : str + 去耀斑 BSQ 影像路径 + hdr_path : str, optional + ENVI HDR 文件路径（None → 自动构造） + output_dir : str, optional + 输出目录（None → 与 bsq_path 同目录下的 8_WaterIndex_Images/） + formula_names : list, optional + 要处理的公式名列表（None → 处理全部） + water_mask : np.ndarray, optional + 水域掩膜数组（与 BSQ 同形状），掩膜值为 0 表示陆地， + 将被强制赋值为 nodata_value + nodata_value : float + NoData 标记值 + progress_callback : callable, optional + 回调 (msg: str, pct: float) + + Returns + ------- + dict + {公式名: 输出 GeoTIFF 路径} + """ + # ── 自动构造 HDR 路径 ──────────────────────────────────────────── + if hdr_path is None: + hdr_path = os.path.splitext(bsq_path)[0] + '.hdr' + if not os.path.isfile(hdr_path): + hdr_path_alt = os.path.splitext(bsq_path)[0] + '.HDR' + if os.path.isfile(hdr_path_alt): + hdr_path = hdr_path_alt + + # ── 自动构造输出目录 ──────────────────────────────────────────── + if output_dir is None: + output_dir = os.path.join(os.path.dirname(bsq_path), '8_WaterIndex_Images') + os.makedirs(output_dir, exist_ok=True) + + def progress(msg: str, pct: float): + if progress_callback: + progress_callback(msg, pct) + + # ── 获取影像元数据 ─────────────────────────────────────────────── + progress("正在打开影像并读取元数据…", 2) + meta = self.get_image_metadata(bsq_path, hdr_path) + + width = meta.get('width', 0) + height = meta.get('height', 0) + n_bands = meta.get('bands', 0) + wavelengths = meta.get('wavelengths', []) + geotransform = meta.get('geotransform') + projection = meta.get('projection') + + if n_bands == 0 or width == 0 or height == 0: + raise ValueError(f"影像元数据无效，无法处理: {bsq_path}") + + if not wavelengths: + raise ValueError(f"无法从 {hdr_path} 读取波长信息，公式无法解析") + + progress( + f"影像: {width}×{height}像素, {n_bands}波段, " + f"波长 {wavelengths[0]:.1f}–{wavelengths[-1]:.1f}nm", + 5 + ) + + # ── 过滤要处理的公式 ────────────────────────────────────────────── + if formula_names: + formulas_to_run = [ + f for f in self.formulas + if f.get('Formula_Name', '').strip() in formula_names + ] + else: + formulas_to_run = list(self.formulas) + + results: Dict[str, str] = {} + total = len(formulas_to_run) + + # ── 逐公式处理 ─────────────────────────────────────────────────── + for i, formula_row in enumerate(formulas_to_run): + fname = formula_row.get('Formula_Name', '').strip() + fstr = formula_row.get('Formula', '').strip() + category = formula_row.get('Category', '').strip() + ftype = formula_row.get('Formula_Type', '').strip() + + if not fname or not fstr: + continue + + progress( + f"[{i + 1}/{total}] {fname} ({category})", + 5 + 90 * i / total + ) + + try: + # 1) 提取公式所需的波长列表 + required_wvs = self._parse_formula_wavelengths(fstr) + + # 2) 按需读取波段数据（相同波长只读一次） + band_data: Dict[int, np.ndarray] = {} + for wv in required_wvs: + if wv not in band_data: + band_idx = self._find_nearest_band_index(wv, wavelengths) + if not (0 < band_idx <= n_bands): + print(f" ⚠ 公式 '{fname}' 引用波段 {band_idx}，超出范围 ({n_bands})，跳过") + raise ValueError(f"波段 {band_idx} 超出影像范围") + band_data[wv] = self._read_band_as_float(bsq_path, band_idx) + + # 3) 矩阵运算 + index_arr = self._eval_formula_fast(fstr, band_data) + if index_arr is None: + print(f" ⚠ 公式 '{fname}' 计算失败，跳过") + continue + + # 4) NoData 处理：NaN / Inf → nodata_value + index_arr = np.where(np.isfinite(index_arr), index_arr, nodata_value) + + # 4b) 水域掩膜拦截：陆地像素（mask==0）强制赋 NoData + if water_mask is not None: + land_pixels = (water_mask == 0) + land_count = int(land_pixels.sum()) + if land_count > 0: + index_arr = np.where(land_pixels, nodata_value, index_arr) + print(f" 🗺 掩膜处理：陆地像素 {land_count:,} 个已设为 NoData") + + # 5) 输出 GeoTIFF + safe_fname = re.sub(r'[^\w\u4e00-\u9fff-]', '_', fname) + out_tif = os.path.join(output_dir, f"{safe_fname}.tif") + + self._write_geotiff( + out_path=out_tif, + data=index_arr, + reference_bsq=bsq_path, + nodata_value=nodata_value, + description=f"{fname}|{category}|{ftype}|{fstr}", + ) + + results[fname] = out_tif + valid = index_arr[index_arr != nodata_value] + mean_val = float(np.mean(valid)) if valid.size else np.nan + print(f" ✅ {fname} → {out_tif} (mean={mean_val:.4f})") + + except ValueError as ve: + print(f" ⏭ 跳过 '{fname}': {ve}") + continue + except Exception as e: + print(f" ❌ 公式 '{fname}' 失败: {e}\n{traceback.format_exc()}") + continue + + progress(f"完成！共输出 {len(results)} / {total} 个指数图", 100) + return results + + def _write_geotiff( + self, + out_path: str, + data: np.ndarray, + reference_bsq: str, + nodata_value: float = -9999.0, + description: str = "", + ) -> None: + """将数组写入 GeoTIFF，克隆原始 BSQ 的地理信息。 + + Parameters + ---------- + out_path : str + 输出 GeoTIFF 路径 + data : np.ndarray + 2D 数据数组（height, width） + reference_bsq : str + 参考 BSQ 影像路径（用于克隆 GeoTransform / Projection） + nodata_value : float + NoData 标记值 + description : str + GDAL 数据集描述 + """ + height, width = data.shape + + driver = gdal.GetDriverByName('GTiff') + if driver is None: + raise RuntimeError("GDAL GTiff 驱动不可用") + + out_ds = driver.Create( + out_path, + width, height, + 1, + gdal.GDT_Float32, + options=['COMPRESS=LZW', 'TILED=YES', 'BIGTIFF=IF_SAFER'], + ) + if out_ds is None: + raise RuntimeError(f"无法创建 GeoTIFF: {out_path}") + + # 写入数据 + out_band = out_ds.GetRasterBand(1) + out_band.SetNoDataValue(nodata_value) + out_band.WriteArray(data) + out_band.FlushCache() + + # 写入描述 + if description: + out_band.SetDescription(description) + + # ★★★ 克隆原始 BSQ 的 GeoTransform 和 Projection ★★★ + ref_ds = gdal.Open(reference_bsq, gdal.GA_ReadOnly) + if ref_ds is not None: + gt = ref_ds.GetGeoTransform() + proj = ref_ds.GetProjection() + if gt and gt != (0, 1, 0, 0, 0, 1): + out_ds.SetGeoTransform(gt) + if proj: + out_ds.SetProjection(proj) + ref_ds = None + + out_ds = None + + # ------------------------------------------------------------------ + # Pipeline 入口（供 PipelineRunner 调用） + # ------------------------------------------------------------------ + + def run_inversion( + self, + deglint_img_path: str, + work_dir: str, + formula_csv_path: Optional[str] = None, + selected_formulas: Optional[List[str]] = None, + water_mask_path: Optional[str] = None, + nodata_value: float = -9999.0, + callback: Optional[Callable] = None, + **kwargs, + ) -> Dict[str, str]: + """Pipeline 入口方法。 + + Parameters + ---------- + deglint_img_path : str + 去耀斑影像 BSQ 路径 + work_dir : str + 工作目录 + formula_csv_path : str, optional + waterindex.csv 路径（None → 使用初始化时的路径） + selected_formulas : list, optional + 要处理的公式列表 + water_mask_path : str, optional + 水域掩膜路径（如 1_water_mask/water_mask.dat）， + 掩膜中为 0 的像素视为陆地区域，其指数值将被强制设为 NoData。 + nodata_value : float + NoData 标记值，默认 -9999.0 + callback : callable, optional + 进度回调 + + Returns + ------- + dict + {公式名: 输出 GeoTIFF 路径} + """ + # 重新加载公式（如指定了新路径） + if formula_csv_path: + self.reload(formula_csv_path) + elif not self.formulas: + raise RuntimeError("WaterIndexProcessor 未加载公式，请指定 formula_csv_path") + + def notify(msg: str, pct: float): + if callback: + callback(msg, pct) + + notify("开始水色指数反演", 0) + + bsq_path = deglint_img_path + hdr_path = os.path.splitext(bsq_path)[0] + '.hdr' + if not os.path.isfile(hdr_path): + hdr_path_alt = os.path.splitext(bsq_path)[0] + '.HDR' + if os.path.isfile(hdr_path_alt): + hdr_path = hdr_path_alt + + output_dir = os.path.join(work_dir, "8_WaterIndex_Images") + + # ── 加载水域掩膜（可选）─────────────────────────────────────── + water_mask: Optional[np.ndarray] = None + if water_mask_path: + if os.path.isfile(water_mask_path): + try: + import rasterio + with rasterio.open(water_mask_path) as msrc: + water_mask = msrc.read(1) + print(f"[run_inversion] 水域掩膜已加载: {water_mask_path}，" + f"形状={water_mask.shape}，" + f"陆地区域(0)={int((water_mask == 0).sum())}，" + f"水区域(>0)={int((water_mask > 0).sum())}") + except Exception as mask_err: + print(f"[run_inversion] ⚠ 掩膜加载失败，跳过掩膜处理: {mask_err}") + water_mask = None + else: + print(f"[run_inversion] ⚠ 水域掩膜文件不存在: {water_mask_path}，跳过掩膜处理") + + notify("水色指数处理中…", 20) + + results = self.process_bsq( + bsq_path=bsq_path, + hdr_path=hdr_path, + output_dir=output_dir, + formula_names=selected_formulas, + water_mask=water_mask, + nodata_value=nodata_value, + progress_callback=lambda m, p: notify(m, 20 + 70 * p / 100), + ) + + notify("水色指数反演完成", 100) + return results diff --git a/src/utils/generate_pure_water_csv.py b/src/utils/generate_pure_water_csv.py new file mode 100644 index 0000000..bde7868 --- /dev/null +++ b/src/utils/generate_pure_water_csv.py @@ -0,0 +1,34 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +""" +生成纯水 IOPs（吸收与散射系数）插值表 +参考 Pope & Fry (1997) 吸收系数 + Morel & Luukko (1974) 后向散射系数 +""" + +import numpy as np +import pandas as pd +import os + + +def generate_csv(): + anchor_points = { + 400: [0.017, 0.0076], 440: [0.015, 0.0049], 500: [0.026, 0.0029], + 555: [0.060, 0.0019], 600: [0.245, 0.0014], 650: [0.340, 0.0010], + 709: [0.735, 0.0007], 750: [2.470, 0.0005], 800: [2.070, 0.0004] + } + wl_anchors = list(anchor_points.keys()) + aw_anchors = [v[0] for v in anchor_points.values()] + bbw_anchors = [v[1] for v in anchor_points.values()] + + wavelengths = np.arange(400, 801, 1) + aw_interp = np.interp(wavelengths, wl_anchors, aw_anchors) + bbw_interp = np.interp(wavelengths, wl_anchors, bbw_anchors) + + df = pd.DataFrame({'Wavelength': wavelengths, 'aw': aw_interp, 'bbw': bbw_interp}) + save_path = os.path.join(os.path.dirname(__file__), 'pure_water_iops.csv') + df.to_csv(save_path, index=False) + print(f"纯水 IOPs 表已生成: {save_path}") + + +if __name__ == '__main__': + generate_csv() \ No newline at end of file diff --git a/src/utils/pure_water_iops.csv b/src/utils/pure_water_iops.csv new file mode 100644 index 0000000..979ca47 --- /dev/null +++ b/src/utils/pure_water_iops.csv @@ -0,0 +1,402 @@ +Wavelength,aw,bbw +400,0.017,0.0076 +401,0.01695,0.0075325 +402,0.016900000000000002,0.007465 +403,0.01685,0.0073974999999999996 +404,0.016800000000000002,0.00733 +405,0.01675,0.0072625 +406,0.0167,0.007195 +407,0.01665,0.0071275 +408,0.0166,0.00706 +409,0.016550000000000002,0.0069925000000000005 +410,0.0165,0.006925 +411,0.01645,0.0068575 +412,0.0164,0.00679 +413,0.01635,0.0067225 +414,0.016300000000000002,0.006655 +415,0.01625,0.0065875 +416,0.0162,0.00652 +417,0.01615,0.0064525 +418,0.0161,0.006385 +419,0.016050000000000002,0.0063175 +420,0.016,0.00625 +421,0.01595,0.0061825000000000005 +422,0.0159,0.006115 +423,0.01585,0.0060475 +424,0.0158,0.00598 +425,0.01575,0.0059125 +426,0.0157,0.0058449999999999995 +427,0.01565,0.0057775 +428,0.0156,0.00571 +429,0.01555,0.0056425 +430,0.0155,0.005575 +431,0.01545,0.0055075 +432,0.0154,0.00544 +433,0.015349999999999999,0.0053725000000000005 +434,0.0153,0.005305 +435,0.01525,0.0052375 +436,0.0152,0.00517 +437,0.01515,0.005102499999999999 +438,0.015099999999999999,0.0050349999999999995 +439,0.01505,0.0049675 +440,0.015,0.0049 +441,0.015183333333333333,0.004866666666666667 +442,0.015366666666666666,0.004833333333333334 +443,0.01555,0.0048 +444,0.015733333333333332,0.0047666666666666664 +445,0.015916666666666666,0.004733333333333333 +446,0.0161,0.0047 +447,0.016283333333333334,0.004666666666666666 +448,0.016466666666666668,0.004633333333333333 +449,0.016649999999999998,0.0046 +450,0.016833333333333332,0.004566666666666667 +451,0.017016666666666666,0.004533333333333333 +452,0.0172,0.0045 +453,0.017383333333333334,0.0044666666666666665 +454,0.017566666666666668,0.004433333333333333 +455,0.01775,0.004399999999999999 +456,0.017933333333333332,0.004366666666666666 +457,0.018116666666666666,0.004333333333333333 +458,0.0183,0.0043 +459,0.018483333333333334,0.004266666666666667 +460,0.018666666666666665,0.004233333333333333 +461,0.01885,0.0042 +462,0.019033333333333333,0.004166666666666667 +463,0.019216666666666667,0.0041333333333333335 +464,0.0194,0.0040999999999999995 +465,0.019583333333333335,0.004066666666666666 +466,0.019766666666666665,0.004033333333333333 +467,0.01995,0.004 +468,0.020133333333333333,0.003966666666666667 +469,0.020316666666666667,0.003933333333333333 +470,0.020499999999999997,0.0039 +471,0.02068333333333333,0.0038666666666666667 +472,0.020866666666666665,0.003833333333333333 +473,0.02105,0.0037999999999999996 +474,0.021233333333333333,0.0037666666666666664 +475,0.021416666666666667,0.0037333333333333333 +476,0.0216,0.0036999999999999997 +477,0.02178333333333333,0.003666666666666666 +478,0.021966666666666666,0.003633333333333333 +479,0.02215,0.0036 +480,0.022333333333333334,0.0035666666666666663 +481,0.022516666666666667,0.003533333333333333 +482,0.022699999999999998,0.0034999999999999996 +483,0.022883333333333332,0.0034666666666666665 +484,0.023066666666666666,0.0034333333333333334 +485,0.02325,0.0034 +486,0.023433333333333334,0.0033666666666666662 +487,0.023616666666666668,0.003333333333333333 +488,0.0238,0.0033 +489,0.023983333333333332,0.0032666666666666664 +490,0.024166666666666666,0.003233333333333333 +491,0.02435,0.0031999999999999997 +492,0.02453333333333333,0.0031666666666666666 +493,0.024716666666666665,0.003133333333333333 +494,0.0249,0.0030999999999999995 +495,0.025083333333333332,0.0030666666666666663 +496,0.025266666666666666,0.003033333333333333 +497,0.02545,0.003 +498,0.025633333333333334,0.0029666666666666665 +499,0.02581666666666667,0.002933333333333333 +500,0.026,0.0029 +501,0.02661818181818182,0.0028818181818181816 +502,0.027236363636363635,0.0028636363636363633 +503,0.027854545454545455,0.002845454545454545 +504,0.02847272727272727,0.002827272727272727 +505,0.02909090909090909,0.002809090909090909 +506,0.029709090909090907,0.002790909090909091 +507,0.030327272727272727,0.0027727272727272726 +508,0.030945454545454543,0.0027545454545454544 +509,0.03156363636363636,0.002736363636363636 +510,0.03218181818181818,0.002718181818181818 +511,0.032799999999999996,0.0026999999999999997 +512,0.03341818181818182,0.0026818181818181815 +513,0.034036363636363635,0.0026636363636363637 +514,0.03465454545454545,0.0026454545454545455 +515,0.035272727272727275,0.0026272727272727272 +516,0.03589090909090909,0.002609090909090909 +517,0.03650909090909091,0.0025909090909090908 +518,0.037127272727272724,0.0025727272727272725 +519,0.03774545454545454,0.0025545454545454543 +520,0.038363636363636364,0.002536363636363636 +521,0.03898181818181818,0.002518181818181818 +522,0.039599999999999996,0.0025 +523,0.04021818181818182,0.002481818181818182 +524,0.040836363636363636,0.0024636363636363636 +525,0.04145454545454545,0.0024454545454545454 +526,0.042072727272727276,0.002427272727272727 +527,0.042690909090909085,0.002409090909090909 +528,0.04330909090909091,0.0023909090909090907 +529,0.043927272727272725,0.002372727272727273 +530,0.04454545454545454,0.0023545454545454546 +531,0.045163636363636364,0.0023363636363636364 +532,0.04578181818181818,0.002318181818181818 +533,0.0464,0.0023 +534,0.04701818181818182,0.0022818181818181817 +535,0.047636363636363636,0.0022636363636363635 +536,0.04825454545454545,0.0022454545454545453 +537,0.048872727272727276,0.002227272727272727 +538,0.049490909090909085,0.002209090909090909 +539,0.05010909090909091,0.0021909090909090906 +540,0.050727272727272725,0.002172727272727273 +541,0.05134545454545454,0.0021545454545454546 +542,0.051963636363636365,0.0021363636363636363 +543,0.05258181818181818,0.002118181818181818 +544,0.0532,0.0021 +545,0.05381818181818182,0.0020818181818181816 +546,0.05443636363636363,0.002063636363636364 +547,0.05505454545454545,0.0020454545454545456 +548,0.055672727272727277,0.0020272727272727274 +549,0.056290909090909086,0.002009090909090909 +550,0.05690909090909091,0.001990909090909091 +551,0.05752727272727273,0.0019727272727272727 +552,0.05814545454545454,0.0019545454545454545 +553,0.058763636363636365,0.0019363636363636362 +554,0.059381818181818175,0.0019181818181818182 +555,0.06,0.0019 +556,0.0641111111111111,0.001888888888888889 +557,0.06822222222222223,0.0018777777777777777 +558,0.07233333333333333,0.0018666666666666666 +559,0.07644444444444444,0.0018555555555555556 +560,0.08055555555555555,0.0018444444444444443 +561,0.08466666666666667,0.0018333333333333333 +562,0.08877777777777777,0.0018222222222222223 +563,0.09288888888888888,0.001811111111111111 +564,0.097,0.0018 +565,0.10111111111111111,0.001788888888888889 +566,0.10522222222222222,0.0017777777777777779 +567,0.10933333333333334,0.0017666666666666666 +568,0.11344444444444444,0.0017555555555555556 +569,0.11755555555555555,0.0017444444444444445 +570,0.12166666666666667,0.0017333333333333333 +571,0.12577777777777777,0.0017222222222222222 +572,0.1298888888888889,0.0017111111111111112 +573,0.134,0.0017 +574,0.13811111111111113,0.0016888888888888889 +575,0.14222222222222222,0.0016777777777777778 +576,0.14633333333333334,0.0016666666666666666 +577,0.15044444444444444,0.0016555555555555555 +578,0.15455555555555556,0.0016444444444444445 +579,0.15866666666666668,0.0016333333333333334 +580,0.1627777777777778,0.0016222222222222222 +581,0.1668888888888889,0.0016111111111111111 +582,0.17099999999999999,0.0015999999999999999 +583,0.1751111111111111,0.0015888888888888888 +584,0.17922222222222223,0.0015777777777777778 +585,0.18333333333333335,0.0015666666666666667 +586,0.18744444444444444,0.0015555555555555555 +587,0.19155555555555556,0.0015444444444444444 +588,0.19566666666666668,0.0015333333333333334 +589,0.19977777777777778,0.0015222222222222221 +590,0.2038888888888889,0.001511111111111111 +591,0.20800000000000002,0.0015 +592,0.2121111111111111,0.001488888888888889 +593,0.21622222222222223,0.0014777777777777777 +594,0.22033333333333335,0.0014666666666666667 +595,0.22444444444444445,0.0014555555555555554 +596,0.22855555555555557,0.0014444444444444444 +597,0.2326666666666667,0.0014333333333333333 +598,0.23677777777777778,0.0014222222222222223 +599,0.2408888888888889,0.001411111111111111 +600,0.245,0.0014 +601,0.2469,0.001392 +602,0.2488,0.001384 +603,0.2507,0.001376 +604,0.2526,0.001368 +605,0.2545,0.00136 +606,0.2564,0.001352 +607,0.2583,0.001344 +608,0.2602,0.001336 +609,0.2621,0.001328 +610,0.264,0.00132 +611,0.2659,0.001312 +612,0.2678,0.001304 +613,0.2697,0.001296 +614,0.2716,0.001288 +615,0.2735,0.0012799999999999999 +616,0.2754,0.001272 +617,0.2773,0.001264 +618,0.2792,0.001256 +619,0.2811,0.001248 +620,0.28300000000000003,0.00124 +621,0.2849,0.001232 +622,0.2868,0.001224 +623,0.2887,0.001216 +624,0.2906,0.0012079999999999999 +625,0.2925,0.0012000000000000001 +626,0.2944,0.001192 +627,0.2963,0.001184 +628,0.2982,0.001176 +629,0.30010000000000003,0.001168 +630,0.30200000000000005,0.00116 +631,0.3039,0.001152 +632,0.3058,0.001144 +633,0.30770000000000003,0.001136 +634,0.3096,0.001128 +635,0.3115,0.00112 +636,0.3134,0.001112 +637,0.3153,0.001104 +638,0.31720000000000004,0.001096 +639,0.31910000000000005,0.001088 +640,0.321,0.00108 +641,0.3229,0.001072 +642,0.32480000000000003,0.001064 +643,0.3267,0.001056 +644,0.3286,0.0010479999999999999 +645,0.3305,0.0010400000000000001 +646,0.33240000000000003,0.001032 +647,0.33430000000000004,0.001024 +648,0.33620000000000005,0.001016 +649,0.3381,0.001008 +650,0.34,0.001 +651,0.3466949152542373,0.0009949152542372882 +652,0.3533898305084746,0.0009898305084745764 +653,0.3600847457627119,0.0009847457627118643 +654,0.3667796610169492,0.0009796610169491525 +655,0.37347457627118646,0.0009745762711864407 +656,0.38016949152542373,0.0009694915254237288 +657,0.38686440677966105,0.000964406779661017 +658,0.3935593220338983,0.0009593220338983051 +659,0.4002542372881356,0.0009542372881355932 +660,0.4069491525423729,0.0009491525423728814 +661,0.41364406779661017,0.0009440677966101695 +662,0.4203389830508475,0.0009389830508474577 +663,0.42703389830508476,0.0009338983050847457 +664,0.433728813559322,0.0009288135593220339 +665,0.44042372881355935,0.0009237288135593221 +666,0.4471186440677966,0.0009186440677966102 +667,0.45381355932203393,0.0009135593220338983 +668,0.4605084745762712,0.0009084745762711864 +669,0.46720338983050846,0.0009033898305084746 +670,0.4738983050847458,0.0008983050847457628 +671,0.48059322033898305,0.0008932203389830508 +672,0.4872881355932204,0.000888135593220339 +673,0.49398305084745764,0.0008830508474576271 +674,0.5006779661016949,0.0008779661016949153 +675,0.5073728813559322,0.0008728813559322035 +676,0.5140677966101694,0.0008677966101694915 +677,0.5207627118644068,0.0008627118644067796 +678,0.5274576271186441,0.0008576271186440678 +679,0.5341525423728813,0.000852542372881356 +680,0.5408474576271186,0.000847457627118644 +681,0.547542372881356,0.0008423728813559322 +682,0.5542372881355933,0.0008372881355932203 +683,0.5609322033898305,0.0008322033898305085 +684,0.5676271186440678,0.0008271186440677967 +685,0.574322033898305,0.0008220338983050847 +686,0.5810169491525423,0.0008169491525423729 +687,0.5877118644067797,0.000811864406779661 +688,0.594406779661017,0.0008067796610169492 +689,0.6011016949152542,0.0008016949152542373 +690,0.6077966101694916,0.0007966101694915254 +691,0.6144915254237289,0.0007915254237288136 +692,0.6211864406779661,0.0007864406779661017 +693,0.6278813559322034,0.0007813559322033899 +694,0.6345762711864407,0.000776271186440678 +695,0.6412711864406779,0.0007711864406779661 +696,0.6479661016949152,0.0007661016949152542 +697,0.6546610169491525,0.0007610169491525424 +698,0.6613559322033898,0.0007559322033898305 +699,0.6680508474576271,0.0007508474576271187 +700,0.6747457627118645,0.0007457627118644068 +701,0.6814406779661017,0.0007406779661016949 +702,0.688135593220339,0.000735593220338983 +703,0.6948305084745763,0.0007305084745762712 +704,0.7015254237288135,0.0007254237288135593 +705,0.7082203389830508,0.0007203389830508475 +706,0.7149152542372881,0.0007152542372881356 +707,0.7216101694915253,0.0007101694915254237 +708,0.7283050847457627,0.0007050847457627119 +709,0.735,0.0007 +710,0.7773170731707317,0.0006951219512195122 +711,0.8196341463414634,0.0006902439024390243 +712,0.8619512195121951,0.0006853658536585366 +713,0.9042682926829269,0.0006804878048780488 +714,0.9465853658536586,0.000675609756097561 +715,0.9889024390243903,0.0006707317073170731 +716,1.031219512195122,0.0006658536585365854 +717,1.0735365853658538,0.0006609756097560976 +718,1.1158536585365855,0.0006560975609756097 +719,1.158170731707317,0.0006512195121951219 +720,1.200487804878049,0.0006463414634146342 +721,1.2428048780487806,0.0006414634146341463 +722,1.2851219512195122,0.0006365853658536585 +723,1.327439024390244,0.0006317073170731707 +724,1.3697560975609757,0.000626829268292683 +725,1.4120731707317074,0.0006219512195121951 +726,1.4543902439024392,0.0006170731707317073 +727,1.496707317073171,0.0006121951219512195 +728,1.5390243902439025,0.0006073170731707317 +729,1.5813414634146343,0.0006024390243902439 +730,1.6236585365853662,0.0005975609756097561 +731,1.6659756097560976,0.0005926829268292683 +732,1.7082926829268295,0.0005878048780487805 +733,1.7506097560975613,0.0005829268292682927 +734,1.7929268292682927,0.0005780487804878049 +735,1.8352439024390246,0.000573170731707317 +736,1.8775609756097564,0.0005682926829268292 +737,1.9198780487804878,0.0005634146341463415 +738,1.9621951219512197,0.0005585365853658537 +739,2.0045121951219516,0.0005536585365853658 +740,2.0468292682926834,0.0005487804878048781 +741,2.089146341463415,0.0005439024390243903 +742,2.1314634146341467,0.0005390243902439024 +743,2.1737804878048785,0.0005341463414634146 +744,2.21609756097561,0.0005292682926829268 +745,2.258414634146342,0.0005243902439024391 +746,2.3007317073170737,0.0005195121951219512 +747,2.343048780487805,0.0005146341463414634 +748,2.385365853658537,0.0005097560975609757 +749,2.427682926829269,0.0005048780487804878 +750,2.47,0.0005 +751,2.462,0.000498 +752,2.454,0.000496 +753,2.446,0.000494 +754,2.438,0.000492 +755,2.43,0.00049 +756,2.422,0.000488 +757,2.414,0.000486 +758,2.406,0.000484 +759,2.398,0.000482 +760,2.39,0.00048 +761,2.382,0.000478 +762,2.374,0.000476 +763,2.366,0.00047400000000000003 +764,2.358,0.00047200000000000003 +765,2.35,0.00047000000000000004 +766,2.342,0.000468 +767,2.334,0.000466 +768,2.326,0.000464 +769,2.318,0.000462 +770,2.31,0.00046 +771,2.302,0.000458 +772,2.294,0.000456 +773,2.286,0.00045400000000000003 +774,2.278,0.00045200000000000004 +775,2.27,0.00045 +776,2.262,0.000448 +777,2.254,0.000446 +778,2.246,0.000444 +779,2.238,0.000442 +780,2.23,0.00044 +781,2.222,0.000438 +782,2.214,0.000436 +783,2.206,0.00043400000000000003 +784,2.198,0.000432 +785,2.19,0.00043000000000000004 +786,2.182,0.000428 +787,2.174,0.000426 +788,2.166,0.000424 +789,2.158,0.000422 +790,2.15,0.00042 +791,2.142,0.000418 +792,2.134,0.00041600000000000003 +793,2.126,0.00041400000000000003 +794,2.118,0.00041200000000000004 +795,2.11,0.00041 +796,2.102,0.000408 +797,2.094,0.000406 +798,2.086,0.000404 +799,2.078,0.000402 +800,2.07,0.0004 diff --git a/src/utils/water_owt_config.py b/src/utils/water_owt_config.py new file mode 100644 index 0000000..573a226 --- /dev/null +++ b/src/utils/water_owt_config.py @@ -0,0 +1,114 @@ +# -*- coding: utf-8 -*- +""" +水域光学特性配置库 (Water Optical Type Configuration) + +基于 QAA-v5 / QAA-v6 通用参数，为中国典型内陆水体配置参考波段和算法参数。 +参考： + - Lee et al. (2002) QAA-v4 + - Lee et al. (2010) QAA-v5 + - Lee et al. (2014) QAA-v6 +""" + +from typing import Dict, List, Optional + + +LAKE_CONFIGS: Dict[str, Dict] = { + # 太湖：藻型富营养化水体，参考波长 709nm + "太湖 (Lake Taihu)": { + "description": "太湖 — 藻型富营养化水体，参考波长 709nm", + "lambda_0": 709, + "reference_band": 709, + "qaa_version": "QAA-v6", + "S": 0.015, + "gamma": 1.0, + "case": "bloom_dominant", + "notes": "藻类主导水体，b_bp 偏高，叶绿素峰在 700-710nm 区间", + }, + # 滇池：高浊度藻型水体，参考波长 710nm + "滇池 (Lake Dianchi)": { + "description": "滇池 — 高浊度藻型水体，参考波长 710nm", + "lambda_0": 710, + "reference_band": 710, + "qaa_version": "QAA-v5", + "S": 0.020, + "gamma": 1.0, + "case": "turbid_algal", + "notes": "高浊度 + 藻华共存，散射贡献显著", + }, + # 青海湖：寡营养深水湖泊，参考波长 555nm（蓝绿波段） + "青海湖 (Lake Qinghai)": { + "description": "青海湖 — 寡营养深水湖泊，参考波长 555nm", + "lambda_0": 555, + "reference_band": 555, + "qaa_version": "QAA-v6", + "S": 0.006, + "gamma": 0.65, + "case": "oligotrophic_clear", + "notes": "清澈寡营养水体，CDOM 极低，参考 555nm（蓝绿峰）", + }, + # 巢湖：富营养化藻型水体，参考波长 709nm + "巢湖 (Lake Chaohu)": { + "description": "巢湖 — 富营养化藻型水体，参考波长 709nm", + "lambda_0": 709, + "reference_band": 709, + "qaa_version": "QAA-v5", + "S": 0.018, + "gamma": 1.0, + "case": "bloom_dominant", + "notes": "太湖类型水体，蓝藻水华频发", + }, + # 长江中下游典型水体：中等浊度，参考波长 665nm + "长江中下游 (Mid-Lower Yangtze)": { + "description": "长江中下游 — 中等浊度混合型水体，参考波长 665nm", + "lambda_0": 665, + "reference_band": 665, + "qaa_version": "QAA-v5", + "S": 0.012, + "gamma": 0.80, + "case": "turbid_mixed", + "notes": "悬浮物浓度中等，CDOM 与藻类混合贡献", + }, + # 黄河：浑浊高沙河流，参考波长 650nm + "黄河 (Yellow River)": { + "description": "黄河 — 浑浊高沙河流型水体，参考波长 650nm", + "lambda_0": 650, + "reference_band": 650, + "qaa_version": "QAA-v5", + "S": 0.025, + "gamma": 1.0, + "case": "highly_turbid", + "notes": "悬浮泥沙主导，散射系数极高，参考 650nm 避免强吸收干扰", + }, + # 通用水体（默认参数，用于无特定配置时） + "通用水体 (Generic)": { + "description": "通用水体 — 默认 QAA 参数，参考波长 555nm", + "lambda_0": 555, + "reference_band": 555, + "qaa_version": "QAA-v6", + "S": 0.013, + "gamma": 0.75, + "case": "generic", + "notes": "通用参数，适用于大多数内陆水体初步分析", + }, +} + + +def get_lake_config(lake_name: str) -> Optional[Dict]: + """根据水域名称获取配置字典，未找到返回 None""" + return LAKE_CONFIGS.get(lake_name) + + +def get_all_lake_names() -> List[str]: + """返回所有已配置的水域名称列表（按哈希表键序）""" + return list(LAKE_CONFIGS.keys()) + + +def get_default_lake() -> str: + """返回默认水域名称""" + return "通用水体 (Generic)" + + +def get_lambda_0(lake_name: str) -> int: + """快捷方法：获取指定水域的参考波长""" + cfg = get_lake_config(lake_name) + return cfg["lambda_0"] if cfg else 555 \ No newline at end of file