Initial commit

ocbrdf/.idea/.gitignore

@@ -0,0 +1,8 @@
+# Default ignored files
+/shelf/
+/workspace.xml
+# 基于编辑器的 HTTP 客户端请求
+/httpRequests/
+# Datasource local storage ignored files
+/dataSources/
+/dataSources.local.xml

ocbrdf/.idea/inspectionProfiles/Project_Default.xml

@@ -0,0 +1,44 @@
+<component name="InspectionProjectProfileManager">
+  <profile version="1.0">
+    <option name="myName" value="Project Default" />
+    <inspection_tool class="PyPackageRequirementsInspection" enabled="true" level="WARNING" enabled_by_default="true">
+      <option name="ignoredPackages">
+        <value>
+          <list size="31">
+            <item index="0" class="java.lang.String" itemvalue="catboost" />
+            <item index="1" class="java.lang.String" itemvalue="opencv-python" />
+            <item index="2" class="java.lang.String" itemvalue="xgboost" />
+            <item index="3" class="java.lang.String" itemvalue="opencv-python-headless" />
+            <item index="4" class="java.lang.String" itemvalue="zstandard" />
+            <item index="5" class="java.lang.String" itemvalue="scikit-image" />
+            <item index="6" class="java.lang.String" itemvalue="scipy" />
+            <item index="7" class="java.lang.String" itemvalue="tensorflow_gpu" />
+            <item index="8" class="java.lang.String" itemvalue="h5py" />
+            <item index="9" class="java.lang.String" itemvalue="matplotlib" />
+            <item index="10" class="java.lang.String" itemvalue="numpy" />
+            <item index="11" class="java.lang.String" itemvalue="opencv_python" />
+            <item index="12" class="java.lang.String" itemvalue="Pillow" />
+            <item index="13" class="java.lang.String" itemvalue="tensorflow" />
+            <item index="14" class="java.lang.String" itemvalue="jupyter" />
+            <item index="15" class="java.lang.String" itemvalue="ipykernel" />
+            <item index="16" class="java.lang.String" itemvalue="pandas" />
+            <item index="17" class="java.lang.String" itemvalue="Werkzeug" />
+            <item index="18" class="java.lang.String" itemvalue="cellpose" />
+            <item index="19" class="java.lang.String" itemvalue="torchvision" />
+            <item index="20" class="java.lang.String" itemvalue="Flask" />
+            <item index="21" class="java.lang.String" itemvalue="fiona" />
+            <item index="22" class="java.lang.String" itemvalue="joblib" />
+            <item index="23" class="java.lang.String" itemvalue="pytest" />
+            <item index="24" class="java.lang.String" itemvalue="geopandas" />
+            <item index="25" class="java.lang.String" itemvalue="PyYAML" />
+            <item index="26" class="java.lang.String" itemvalue="pybaselines" />
+            <item index="27" class="java.lang.String" itemvalue="molmass" />
+            <item index="28" class="java.lang.String" itemvalue="noise" />
+            <item index="29" class="java.lang.String" itemvalue="scikit-gstat" />
+            <item index="30" class="java.lang.String" itemvalue="plotly" />
+          </list>
+        </value>
+      </option>
+    </inspection_tool>
+  </profile>
+</component>

ocbrdf/.idea/inspectionProfiles/profiles_settings.xml

@@ -0,0 +1,6 @@
+<component name="InspectionProjectProfileManager">
+  <settings>
+    <option name="USE_PROJECT_PROFILE" value="false" />
+    <version value="1.0" />
+  </settings>
+</component>

ocbrdf/.idea/misc.xml

@@ -0,0 +1,7 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<project version="4">
+  <component name="Black">
+    <option name="sdkName" value="WQ" />
+  </component>
+  <component name="ProjectRootManager" version="2" project-jdk-name="hypro" project-jdk-type="Python SDK" />
+</project>

ocbrdf/.idea/modules.xml

@@ -0,0 +1,8 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<project version="4">
+  <component name="ProjectModuleManager">
+    <modules>
+      <module fileurl="file://$PROJECT_DIR$/.idea/ocbrdf.iml" filepath="$PROJECT_DIR$/.idea/ocbrdf.iml" />
+    </modules>
+  </component>
+</project>

ocbrdf/.idea/ocbrdf.iml

@@ -0,0 +1,12 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<module type="PYTHON_MODULE" version="4">
+  <component name="NewModuleRootManager">
+    <content url="file://$MODULE_DIR$" />
+    <orderEntry type="jdk" jdkName="hypro" jdkType="Python SDK" />
+    <orderEntry type="sourceFolder" forTests="false" />
+  </component>
+  <component name="PyDocumentationSettings">
+    <option name="format" value="PLAIN" />
+    <option name="myDocStringFormat" value="Plain" />
+  </component>
+</module>

ocbrdf/Raman.py

@@ -0,0 +1,36 @@
+import numpy as np
+import xarray as xr
+
+''' 来自 Lee 等人 (2013) 的拉曼校正：
+      《全球海洋中紫外-可见太阳辐射的穿透：来自海洋水色遥感的见解》
+      J Geophys Res Oceans 2013;118:4241–55
+
+    光谱插值在最近波段进行
+'''
+
+""" 拉曼校正查找表的类 """ # TODO 存储并在 ADF 中读取
+class Raman():
+    def __init__(self):
+        bands = [412,443,488,531,551,667]
+        coords = {'bands': bands}
+        self.alpha = xr.DataArray([0.003,0.004,0.011,0.015,0.017,0.018], dims='bands', coords=coords)
+        self.beta1 = xr.DataArray([0.014,0.015,0.010,0.010,0.010,0.010], dims='bands', coords=coords)
+        self.beta2 = xr.DataArray([-0.022,-0.023,-0.051,-0.070,-0.080,-0.081], dims='bands', coords=coords)
+
+    def correct(self, Rrs):
+        # Interpolate coefficients at nearest bands
+        interp_opt = {'method':'nearest', 'kwargs':{'fill_value':'extrapolate'}}
+        alpha = self.alpha.interp(bands=Rrs.bands, **interp_opt)
+        beta1 = self.beta1.interp(bands=Rrs.bands, **interp_opt)
+        beta2 = self.beta2.interp(bands=Rrs.bands, **interp_opt)
+
+        # Select reference bands closest to 440 and 550
+        Rrs440 = Rrs.sel(bands=440, method='nearest')
+        Rrs550 = Rrs.sel(bands=550, method='nearest')
+
+        # Compute Raman factor
+        RF = alpha*Rrs440/Rrs550 + beta1*np.power(Rrs550,beta2)
+
+        # Correct
+        return Rrs/(1+RF)
+

ocbrdf/aviris.py

@@ -0,0 +1,978 @@
+
+import gc
+import os
+import numpy as np
+import sys
+import xarray as xr
+import rasterio
+from rasterio.windows import Window
+import spectral
+
+from brdf_model_M02 import M02
+from brdf_model_M02SeaDAS import M02SeaDAS
+from brdf_model_L11 import L11
+from brdf_model_O25 import O25
+from brdf_utils import ADF_OCP, squeeze_trivial_dims
+
+# 尝试导入 tqdm 进度条库，如果不可用则设为 None
+try:
+    from tqdm import tqdm
+except ImportError:
+    tqdm = None
+
+UNC_LUT_CACHE = {}
+
+"""
+主 BRDF 校正模块
+输入为 xarray 数据集
+所需的光谱维度为 "bands"，其他维度自由
+输入数据集中必填的字段：
+Rw: 方向性海洋反射率
+sza: 太阳天顶角
+vza: 观测天顶角
+raa: 相对方位角（当太阳和观测在同一侧时，raa=0）
+输入数据集中可选的字段：
+Rw_unc: Rw 的不确定度（若缺失，设为零）
+输出数据集中的字段：
+nrrs: 完全归一化的遥感反射率
+rho_ex_w: nrrs * PI
+omega_b: bb/(a+bb)
+eta_b: bbw/bb
+C_brdf: BRDF 校正因子
+brdf_unc: C_brdf 的不确定度
+nrrs_unc: nrrs 的不确定度
+
+使用和适配 brdf_hypercp 模块需保留的信息：
+Brdf_hypercp 是 EUMETSAT 研究项目的一部分
+"BRDF correction of S3 OLCI water reflectance products"（S3 OLCI 水体反射率产品的 BRDF 校正），
+合同号：RB_EUM-CO-21-4600002626-JIG。
+研究团队成员：Davide D'Alimonte (davide.dalimonte@aequora.org)，
+Tamito Kajiyama (tamito.kajiyama@aequora.org)，
+Jaime Pitarch (jaime.pitarchportero@artov.ismar.cnr.it)，
+Vittorio Brando (vittorio.brando@cnr.it)，
+Marco Talone (talone@icm.csic.es) 以及
+Constant Mazeran (constant.mazeran@solvo.fr)。
+BRDF 查找表中的相对方位角遵循 OLCI 约定。详见 https://www.eumetsat.int/media/50720，图 6。
+"""
+
+
+# ---------------------------------------------------------------------------
+# File I/O helpers
+# ---------------------------------------------------------------------------
+
+def _resolve_envi_header(filepath):
+    """Return the ENVI header path accepted by spectral."""
+    if filepath.lower().endswith(".hdr"):
+        return filepath
+
+    hdr_path = filepath + ".hdr"
+    if os.path.exists(hdr_path):
+        return hdr_path
+
+    base, _ext = os.path.splitext(filepath)
+    hdr_path = base + ".hdr"
+    if os.path.exists(hdr_path):
+        return hdr_path
+
+    return filepath
+
+
+def read_bsq(filepath, wavelengths=None, dtype=np.float32):
+    """Read an ENVI BSQ hyperspectral image.
+
+    Parameters
+    ----------
+    filepath : str
+        Path to the ENVI header (.hdr) or image file.
+    wavelengths : array-like, optional
+        Band centre wavelengths in nm.  When *None*, wavelengths are read
+        from the ENVI header if available.
+
+    Returns
+    -------
+    data : ndarray, shape (bands, rows, cols), float64
+    wavelengths : ndarray, shape (bands,)
+    img : spectral image object (metadata / transform access)
+    """
+    img = spectral.open_image(_resolve_envi_header(filepath))
+    data = np.asarray(img.open_memmap(interleave="source"), dtype=dtype)
+    # data = np.transpose(data, (2, 0, 1)) # -> (bands, rows, cols)
+
+    if wavelengths is None:
+        if hasattr(img, 'bands') and img.bands.centers is not None:
+            wavelengths = np.array(img.bands.centers, dtype=np.float32)
+        else:
+            n_bands = data.shape[0]
+            wavelengths = np.arange(1, n_bands + 1, dtype=np.float32)
+
+    return data, np.asarray(wavelengths, dtype=dtype), img
+
+
+def relative_azimuth_olci(saa_deg, vaa_deg):
+    """Relative azimuth (RAA / delta_phi) in [0, 180] deg, OLCI convention.
+
+    Same convention as the module docstring: sun and sensor on the same side -> RAA = 0.
+    """
+    saa = np.asarray(saa_deg, dtype=np.float64)
+    vaa = np.asarray(vaa_deg, dtype=np.float64)
+    d = np.abs(saa - vaa)
+    out = np.minimum(d, 360.0 - d)
+    return out.astype(np.float32)
+
+
+def read_bip_angles(filepath, dtype=np.float32):
+    """Read an ENVI BIP geometry file (10 bands) and build BRDF angles.
+
+    Band layout (1-indexed), last dimension = band (BIP):
+        1  path length
+        2  to-sensor azimuth (VAA)
+        3  to-sensor zenith (VZA)
+        4  to-sun azimuth (SAA)
+        5  to-sun zenith (SZA)
+        6  phase angle
+        7  slope
+        8  aspect
+        9  cosine i
+        10 UTC time
+
+    RAA is not stored; it is computed from SAA and VAA via :func:`relative_azimuth_olci`.
+
+    Parameters
+    ----------
+    filepath : str
+        Path to the ENVI BIP file (header or image file).
+
+    Returns
+    -------
+    dict with keys 'sza', 'saa', 'vza', 'vaa', 'raa'.
+    Each value is a 2-D ndarray of shape (rows, cols).
+    """
+    img = spectral.open_image(_resolve_envi_header(filepath))
+    data = np.asarray(img.open_memmap(interleave="source"), dtype=dtype)
+
+    n_bands = data.shape[2]
+    if n_bands != 10:
+        raise ValueError(
+            f"Geometry file must have exactly 10 bands (got {n_bands}): "
+            "path length, VAA, VZA, SAA, SZA, phase, slope, aspect, cos i, UTC."
+        )
+
+    vaa = data[:, :, 1].copy()
+    vza = data[:, :, 2].copy()
+    saa = data[:, :, 3].copy()
+    sza = data[:, :, 4].copy()
+    raa = relative_azimuth_olci(saa, vaa)
+
+    return {
+        'sza': sza,
+        'saa': saa,
+        'vza': vza,
+        'vaa': vaa,
+        'raa': raa,
+    }
+
+
+def read_water_mask(filepath):
+    """Read a single-band GeoTIFF water mask.
+
+    Parameters
+    ----------
+    filepath : str
+        Path to the GeoTIFF file.
+
+    Returns
+    -------
+    mask : ndarray, shape (rows, cols), int8
+        Pixels with value 1 indicate water; all other values are non-water.
+    """
+    with rasterio.open(filepath) as src:
+        mask = src.read(1).astype(np.int8)
+    return mask
+
+
+def _block_windows(src, block_size=512):
+    """Yield raster windows for block-wise processing."""
+    rows, cols = src.height, src.width
+    for row in range(0, rows, block_size):
+        row_end = min(row + block_size, rows)
+        for col in range(0, cols, block_size):
+            col_end = min(col + block_size, cols)
+            window = Window(col, row, col_end - col, row_end - row)
+            yield window
+
+
+def _resolve_output_envi_path(output_file):
+    """Return the ENVI image path used for rasterio writes."""
+    if output_file.lower().endswith(".img"):
+        return output_file
+    if output_file.lower().endswith(".hdr"):
+        return output_file[:-4] + ".img"
+    return output_file + ".img"
+
+
+def _read_wavelengths_from_header(hyperspectral_file):
+    """Read wavelength metadata from ENVI header if available."""
+    img = spectral.open_image(_resolve_envi_header(hyperspectral_file))
+    if hasattr(img, 'bands') and img.bands.centers is not None:
+        wavelengths = img.bands.centers
+    else:
+        wavelengths = img.metadata.get('wavelength')
+    if wavelengths is None:
+        return None
+    return np.asarray(wavelengths, dtype=np.float32)
+
+
+def save_brdf_result(
+    ds_out,
+    output_file,
+    source_file=None,
+    output_var="Rw_brdf",
+    output_format="ENVI",
+    dtype=np.float32,
+):
+    """Save the BRDF result to ENVI or NetCDF."""
+    output_format = output_format.upper()
+
+    if output_var not in ds_out:
+        raise KeyError(f"Output variable '{output_var}' not found in result dataset.")
+
+    if output_format == "NETCDF":
+        if not output_file.lower().endswith(".nc"):
+            output_file = output_file + ".nc"
+        ds_out.to_netcdf(output_file)
+        return output_file
+
+    if output_format != "ENVI":
+        raise ValueError("output_format must be 'ENVI' or 'NETCDF'.")
+
+    cube = np.asarray(ds_out[output_var].values, dtype=dtype)
+    cube = np.transpose(cube, (1, 2, 0))
+
+    if output_file.lower().endswith(".hdr"):
+        hdr_path = output_file
+    else:
+        hdr_path = output_file + ".hdr"
+
+    metadata = {
+        "interleave": "bsq",
+        "description": f"BRDF corrected result: {output_var}",
+        "bands": cube.shape[2],
+        "lines": cube.shape[0],
+        "samples": cube.shape[1],
+    }
+    if "bands" in ds_out.coords:
+        metadata["wavelength"] = [str(v) for v in ds_out.coords["bands"].values.tolist()]
+
+    if source_file is not None:
+        metadata["source_file"] = source_file
+
+    spectral.envi.save_image(
+        hdr_path,
+        cube,
+        dtype=dtype,
+        interleave="bsq",
+        metadata=metadata,
+        force=True,
+    )
+    return hdr_path
+
+
+def _get_output_filename(base_output_file: str, var_name: str) -> str:
+    """Generate output filename with variable suffix.
+    
+    Examples:
+        input.hdr + Rw_brdf     -> input_Rw_brdf.hdr
+        input.hdr + brdf_unc    -> input_brdf_unc.hdr
+    """
+    base, ext = os.path.splitext(base_output_file)
+    if ext.lower() in ('.hdr', '.img'):
+        base = base  # keep base without extension
+    suffix = f"_{var_name}" if var_name.lower() != "rw_brdf" else ""
+    return f"{base}{suffix}{ext}"
+
+
+def save_multiple_brdf_results(
+    ds_out,
+    output_file,
+    source_file=None,
+    output_vars=None,
+    output_format="ENVI",
+    dtype=np.float32,
+):
+    """Save multiple output variables to separate files."""
+    if output_vars is None or len(output_vars) == 0:
+        output_vars = ["Rw_brdf"]
+    
+    saved_files = []
+    for var in output_vars:
+        out_file = _get_output_filename(output_file, var)
+        try:
+            saved_path = save_brdf_result(
+                ds_out=ds_out,
+                output_file=out_file,
+                source_file=source_file,
+                output_var=var,
+                output_format=output_format,
+                dtype=dtype,
+            )
+            saved_files.append((var, saved_path))
+            print(f"  Saved {var} to: {saved_path}")
+        except KeyError as e:
+            print(f"  Warning: {e}. Skipping {var}.")
+    
+    return saved_files
+
+
+def run_brdf_correction_block(
+    hyperspectral_file,
+    angle_file,
+    output_file,
+    mask_file=None,
+    wavelengths=None,
+    brdf_model='L11',
+    output_vars=None,           # 支持列表
+    block_size=512,
+    compute_uncertainty=False,
+    progress_bar=True,
+):
+    """
+    Run BRDF correction block-wise and write multiple variables to separate ENVI files.
+
+    ``angle_file`` must be a 10-band raster (same grid as hyperspectral); RAA is derived
+    from bands 2 (VAA) and 4 (SAA); see :func:`read_bip_angles` for band order.
+
+    If ``mask_file`` is None, every pixel is corrected (full scene); otherwise only pixels
+    with mask value 1 are corrected (other pixels pass through unchanged for Rw_brdf).
+    """
+    if output_vars is None:
+        output_vars = ["Rw_brdf"]
+    if isinstance(output_vars, str):
+        output_vars = [output_vars]
+
+    # 1. Open all inputs
+    src_rw = rasterio.open(hyperspectral_file)
+    n_bands = src_rw.count
+    n_rows = src_rw.height
+    n_cols = src_rw.width
+
+    src_ang = rasterio.open(angle_file)
+    if src_ang.height != n_rows or src_ang.width != n_cols:
+        raise ValueError("Angle file dimensions do not match hyperspectral file.")
+    if src_ang.count != 10:
+        raise ValueError(
+            f"Geometry file must have 10 bands (got {src_ang.count}). "
+            "Expected: path length, VAA, VZA, SAA, SZA, phase, slope, aspect, cos i, UTC."
+        )
+
+    use_mask = mask_file is not None
+    src_mask = rasterio.open(mask_file) if use_mask else None
+    if use_mask and (src_mask.height != n_rows or src_mask.width != n_cols):
+        raise ValueError("Mask file dimensions do not match hyperspectral file.")
+
+    # 2. Wavelengths
+    if wavelengths is None:
+        wavelengths = _read_wavelengths_from_header(hyperspectral_file)
+    if wavelengths is None:
+        wavelengths = np.arange(1, n_bands + 1, dtype=np.float32)
+    bands_xr = xr.DataArray(wavelengths, dims=('bands',), coords={'bands': wavelengths}, name='bands')
+
+    # 3. Prepare multiple output files (one per variable)
+    out_meta = src_rw.meta.copy()
+    out_meta.update({
+        'driver': 'ENVI',
+        'dtype': 'float32',
+        'count': n_bands,
+        'interleave': 'bsq',
+    })
+
+    # Open one output file for each variable
+    output_dsts = {}
+    output_paths = {}
+    for var_name in output_vars:
+        var_file = _get_output_filename(output_file, var_name)
+        out_path = _resolve_output_envi_path(var_file)
+        output_dsts[var_name] = rasterio.open(out_path, 'w', **out_meta)
+        output_paths[var_name] = out_path
+        print(f"  Writing {var_name} to {out_path}")
+
+    # 4. Band indices in geometry file (1-based); RAA computed from SAA and VAA
+    band_vaa = 2
+    band_vza = 3
+    band_saa = 4
+    band_sza = 5
+
+    # 5. Map output_var to ds_corr field (support uncertainty variables)
+    var_map = {
+        'Rw_brdf': 'rho_ex_w',
+        'rho_ex_w': 'rho_ex_w',
+        'nrrs': 'nrrs',
+        'C_brdf': 'C_brdf',
+        'brdf_unc': 'brdf_unc',
+        'nrrs_unc': 'nrrs_unc',
+    }
+
+    # Validate all requested output variables
+    for v in output_vars:
+        if v not in var_map and v not in ['Rw_brdf']:
+            raise ValueError(
+                f"Unknown output_var '{v}'. Allowed: {list(var_map.keys())}"
+            )
+
+    # 6. Progress bar setup
+    iterator = _block_windows(src_rw, block_size=block_size)
+    if progress_bar and tqdm is not None:
+        total_blocks = ((n_rows + block_size - 1) // block_size) * ((n_cols + block_size - 1) // block_size)
+        iterator = tqdm(iterator, total=total_blocks, desc="Processing blocks", unit="block")
+    elif progress_bar:
+        print("Warning: tqdm not installed, progress bar disabled.")
+
+    # Statistics for performance monitoring
+    blocks_processed = 0
+    blocks_skipped = 0
+    total_corrected_pixels = 0
+
+    # 7. Block-wise processing with early filtering
+    for window in iterator:
+        # Step 1: Read mask FIRST to enable early exit for pure land blocks (mask mode only)
+        if use_mask:
+            mask_block = src_mask.read(1, window=window).astype(np.int8)
+            water_mask = (mask_block == 1)
+        else:
+            mask_block = None
+            water_mask = np.ones((window.height, window.width), dtype=bool)
+        n_water = np.count_nonzero(water_mask)
+
+        if n_water == 0:
+            # Fast path: pure land/no-water block - minimal processing (only when mask is used)
+            blocks_skipped += 1
+            # Still need to read hyperspectral for Rw_brdf output, but skip angles
+            rw_block = src_rw.read(window=window).astype(np.float32)
+            
+            output_blocks = {}
+            for var_name in output_vars:
+                if var_name == 'Rw_brdf':
+                    output_blocks[var_name] = rw_block.copy()
+                else:
+                    output_blocks[var_name] = np.full_like(rw_block, np.nan, dtype=np.float32)
+            
+            for var_name, block_data in output_blocks.items():
+                output_dsts[var_name].write(block_data, window=window)
+            
+            del rw_block, water_mask
+            gc.collect()
+            continue
+
+        # Step 2: Has water pixels - full processing
+        blocks_processed += 1
+        total_corrected_pixels += n_water
+
+        # Read remaining data only for blocks that need processing
+        rw_block = src_rw.read(window=window).astype(np.float32)
+        sza_block = src_ang.read(band_sza, window=window).astype(np.float32)
+        vza_block = src_ang.read(band_vza, window=window).astype(np.float32)
+        saa_block = src_ang.read(band_saa, window=window).astype(np.float32)
+        vaa_block = src_ang.read(band_vaa, window=window).astype(np.float32)
+        raa_block = relative_azimuth_olci(saa_block, vaa_block)
+
+        water_idx = np.where(water_mask)
+
+        # Prepare output blocks for all variables
+        output_blocks = {}
+        for var_name in output_vars:
+            if var_name == 'Rw_brdf':
+                output_blocks[var_name] = rw_block.copy()
+            else:
+                output_blocks[var_name] = np.full_like(rw_block, np.nan, dtype=np.float32)
+
+        # Process water pixels
+        rw_water = rw_block[:, water_idx[0], water_idx[1]].T  # (n_water, bands)
+        sza_water = sza_block[water_idx]
+        vza_water = vza_block[water_idx]
+        raa_water = raa_block[water_idx]
+
+        pixel_ids = np.arange(n_water)
+        ds_block = xr.Dataset(
+            {
+                'Rw': xr.DataArray(rw_water, dims=('n', 'bands'), coords={'n': pixel_ids, 'bands': bands_xr}),
+                'sza': xr.DataArray(sza_water, dims='n', coords={'n': pixel_ids}),
+                'vza': xr.DataArray(vza_water, dims='n', coords={'n': pixel_ids}),
+                'raa': xr.DataArray(raa_water, dims='n', coords={'n': pixel_ids}),
+            }
+        )
+
+        ds_corr = brdf_prototype(
+            ds_block,
+            brdf_model=brdf_model,
+            compute_uncertainty=compute_uncertainty,
+        )
+
+        # Write each requested variable to the water pixels
+        for var_name in output_vars:
+            field_name = var_map.get(var_name, 'rho_ex_w')
+            if var_name == 'Rw_brdf':
+                field_name = 'rho_ex_w'
+            if field_name in ds_corr:
+                values = ds_corr[field_name].values.T
+                output_blocks[var_name][:, water_idx[0], water_idx[1]] = values
+
+        # Write all variables for this block
+        for var_name, block_data in output_blocks.items():
+            output_dsts[var_name].write(block_data, window=window)
+
+        del (rw_block, sza_block, vza_block, saa_block, vaa_block, raa_block,
+             output_blocks, water_mask, water_idx, rw_water, sza_water, vza_water,
+             raa_water, pixel_ids, ds_block, ds_corr)
+        if mask_block is not None:
+            del mask_block
+        gc.collect()
+
+    # Close all files
+    src_rw.close()
+    src_ang.close()
+    if src_mask is not None:
+        src_mask.close()
+    for dst in output_dsts.values():
+        dst.close()
+
+    # Print performance statistics
+    total_blocks = blocks_processed + blocks_skipped
+    if total_blocks > 0:
+        skip_ratio = blocks_skipped / total_blocks * 100
+        print(f"Block processing completed:")
+        print(f"  - Total blocks: {total_blocks}")
+        print(f"  - Processed blocks (corrected): {blocks_processed}")
+        print(f"  - Skipped blocks (no mask pixels): {blocks_skipped} ({skip_ratio:.1f}%)")
+        print(f"  - Total pixels BRDF-corrected: {total_corrected_pixels:,}")
+
+    return list(output_paths.values())
+
+
+# ---------------------------------------------------------------------------
+# Top-level pipeline
+# ---------------------------------------------------------------------------
+def brdf_uncertainty(ds, adf=None):
+    ''' Compute uncertainty of BRDF factor and propagate to nrrs '''
+
+    # Read LUT
+    if adf is None:
+        adf = ADF_OCP
+
+    unc_lut_path = adf % 'UNC'
+    if unc_lut_path not in UNC_LUT_CACHE:
+        UNC_LUT_CACHE[unc_lut_path] = xr.open_dataset(unc_lut_path, engine='netcdf4')
+    LUT = UNC_LUT_CACHE[unc_lut_path]
+
+    # Interpolate relative uncertainty
+    unc = LUT['unc'].interp(lambda_unc=ds.bands, theta_s_unc=ds.sza, theta_v_unc=ds.vza,
+                            delta_phi_unc=ds.raa)
+
+    # Compute absolute uncertainty of factor
+    ds['brdf_unc'] = unc * ds['C_brdf']
+
+    # Flag BRDF_unc where NaN and set to 0
+    ds['brdf_unc_fail'] = np.isnan(ds['brdf_unc'])
+    ds['brdf_unc'] = xr.where(ds['brdf_unc_fail'], 0, ds['brdf_unc'])
+
+    # Propagate to nrrs
+    Rwex_unc2 = ds['brdf_unc'] * ds['brdf_unc'] * ds['Rw'] * ds['Rw']
+    if 'Rw_unc' in ds:
+        Rwex_unc2 += ds['C_brdf'] * ds['C_brdf'] * ds['Rw_unc'] * ds['Rw_unc']
+    ds['nrrs_unc'] = np.sqrt(Rwex_unc2)/np.pi
+
+    return ds
+
+
+def brdf_prototype(ds, adf=None, brdf_model='L11', compute_uncertainty=True):
+    # 测试 BRDF 模型在 GUI 中不支持：强制覆盖
+    # brdf_model = 'M02SeaDAS'
+
+    # 压缩单一维度（例如，投放次数、提取次数等）以避免插值问题
+    ds, squeezedDims = squeeze_trivial_dims(ds)
+
+    # Initialise model
+    if brdf_model == 'M02':
+        BRDF_model = M02(bands=ds.bands, aot=ds.aot, wind=ds.wind, adf=None)  # Don't use brdf_py.ADF context
+    elif brdf_model == 'M02SeaDAS':
+        BRDF_model = M02SeaDAS(bands=ds.bands, adf=None)  # Don't use brdf_py.ADF context
+    elif brdf_model == 'L11':
+        BRDF_model = L11(bands=ds.bands, adf=None)  # Don't use brdf_py.ADF context
+    elif brdf_model == 'O25':
+        BRDF_model = O25(bands=ds.bands, adf=None)  # Don't use brdf_py.ADF context
+    else:
+        print("BRDF model %s not supported" % brdf_model)
+        sys.exit(1)
+
+    # Init pixel
+    BRDF_model.init_pixels(ds['sza'], ds['vza'], ds['raa'])
+
+    # Compute IOP and normalize by iterating
+    ds['nrrs'] = ds['Rw'] / np.pi
+
+    ds['convergeFlag'] = (0 * ds['sza']).astype(bool)
+    ds['C_brdf'] = 0 * ds['nrrs'] + 1
+
+    for iter_brdf in range(int(BRDF_model.niter)):
+
+        # M02: Initialise chl_iter
+        if brdf_model in ['M02', 'M02SeaDAS'] and (iter_brdf == 0):
+            chl_iter = {}
+            ds['log10_chl'] = 0 * ds['sza'] + float(np.log10(BRDF_model.OC4MEchl0))
+            chl_iter[-1]    = 0 * ds['sza'] + float(BRDF_model.OC4MEchl0)
+
+        ds = BRDF_model.backward(ds, iter_brdf)
+
+        # M02: Check convergence (dummy for M02SeaDAS for the moment... epsilon set to 0)
+        if brdf_model in ['M02', 'M02SeaDAS']:
+            chl_iter[iter_brdf] = 10 ** ds['log10_chl']
+            #  Check if convergence is reached |chl_old-chl_new| < epsilon * chl_new
+            ds['convergeFlag'] = (ds['convergeFlag']) | (
+                (np.abs(chl_iter[iter_brdf - 1] - chl_iter[iter_brdf]) < float(BRDF_model.OC4MEepsilon) * chl_iter[
+                    iter_brdf]))
+
+        # Apply forward model in both geometries
+        forward_mod = BRDF_model.forward(ds)
+        forward_mod0 = BRDF_model.forward(ds, normalized=True)
+
+        ratio = forward_mod0 / forward_mod
+
+        # Drop remnant coordinates to avoid ambiguities in the update of the BRDF factor.
+        for coord in ratio.coords:
+            if coord not in ds['C_brdf'].coords:
+                ratio = ratio.drop(coord)
+
+        # Normalize reflectance
+        ds['C_brdf'] = xr.where(ds['convergeFlag'], ds['C_brdf'], ratio)
+        ds['nrrs'] = ds['Rw'] / np.pi * ds['C_brdf']
+
+    # Flag BRDF where NaN and set to 1 (no correction applied).
+    ds['C_brdf_fail'] = np.isnan(ds['C_brdf'])
+    ds['C_brdf'] = xr.where(ds['C_brdf_fail'], 1, ds['C_brdf'])
+    ds['nrrs'] = xr.where(ds['C_brdf_fail'], ds['Rw'] / np.pi, ds['nrrs'])
+
+    # If QAA_fail is raised, raise C_brdf_fail (but still apply C_brdf).
+    if 'QAA_fail' in ds:
+        ds['C_brdf_fail'] = (ds['C_brdf_fail']) | (ds['QAA_fail'])
+
+    # Compute uncertainty only if requested
+    if compute_uncertainty:
+        ds = brdf_uncertainty(ds)
+
+    # Compute flag
+    ds['flags_level2'] = ds['Rw'] * 0  # TODO
+
+    # Convert to reflectance unit
+    ds['rho_ex_w'] = ds['nrrs'] * np.pi
+
+    # Expand squeezed trivial dimensions
+    for dim,d0 in squeezedDims.items():
+        ds = ds.expand_dims(dim,axis=d0)
+
+    return ds
+
+
+def run_brdf_correction(
+    hyperspectral_file,
+    angle_file,
+    mask_file=None,
+    wavelengths=None,
+    brdf_model='L11',
+    chunk_size=4096,
+    output_var=None,          # if None, keep full output (original behaviour)
+    compute_uncertainty=True, # compute uncertainty only if requested
+    progress_bar=True,        # 是否显示进度条
+):
+    """BRDF correction pipeline for water pixels in a hyperspectral scene.
+
+    ``angle_file`` must be a 10-band ENVI geometry cube (see :func:`read_bip_angles`);
+    relative azimuth is computed from to-sun and to-sensor azimuths.
+
+    If ``mask_file`` is None, a mask of all ones is used (full-scene correction).
+    """
+    # ------------------------------------------------------------------
+    # 1. Read all inputs
+    # ------------------------------------------------------------------
+    rw_data, wvl, _img_meta = read_bsq(hyperspectral_file, wavelengths, dtype=np.float32)
+    angles = read_bip_angles(angle_file, dtype=np.float32)
+
+    n_bands, n_rows, n_cols = rw_data.shape
+
+    if mask_file is None:
+        water_mask = np.ones((n_rows, n_cols), dtype=np.int8)
+    else:
+        water_mask = read_water_mask(mask_file)
+
+    # Sanity-check spatial dimensions
+    if angles['sza'].shape != (n_rows, n_cols):
+        raise ValueError(
+            f"Angle file spatial shape {angles['sza'].shape} does not match "
+            f"hyperspectral image shape ({n_rows}, {n_cols})."
+        )
+    if water_mask.shape != (n_rows, n_cols):
+        raise ValueError(
+            f"Mask shape {water_mask.shape} does not match "
+            f"hyperspectral image shape ({n_rows}, {n_cols})."
+        )
+
+    # ------------------------------------------------------------------
+    # 2. Identify water pixels
+    # ------------------------------------------------------------------
+    water_rows, water_cols = np.where(water_mask == 1)
+    n_water = water_rows.size
+    if n_water == 0:
+        raise RuntimeError("No water pixels found in the mask (no pixels == 1).")
+
+    # ------------------------------------------------------------------
+    # 3. Determine which output variable(s) we need to keep
+    # ------------------------------------------------------------------
+    # If output_var is None, we keep everything (original behaviour)
+    keep_all = (output_var is None)
+
+    # Map output_var to the field in ds_corr
+    var_map = {
+        'Rw_brdf': 'rho_ex_w',
+        'rho_ex_w': 'rho_ex_w',
+        'nrrs': 'nrrs',
+        'C_brdf': 'C_brdf'
+    }
+    if not keep_all:
+        if output_var not in var_map:
+            raise ValueError(f"Unknown output_var '{output_var}'. Allowed: {list(var_map.keys())}")
+        field_name = var_map[output_var]
+    else:
+        field_name = None
+
+    # ------------------------------------------------------------------
+    # 4. Allocate output arrays (only those needed)
+    # ------------------------------------------------------------------
+    shape3 = (n_bands, n_rows, n_cols)
+    shape2 = (n_rows, n_cols)
+
+    if keep_all:
+        # Original behaviour: allocate all arrays
+        nrrs_full        = np.full(shape3, np.nan, dtype=np.float32)
+        rho_ex_w_full    = np.full(shape3, np.nan, dtype=np.float32)
+        C_brdf_full      = np.full(shape3, np.nan, dtype=np.float32)
+        brdf_unc_full    = np.full(shape3, np.nan, dtype=np.float32)
+        nrrs_unc_full    = np.full(shape3, np.nan, dtype=np.float32)
+        C_brdf_fail_full = np.zeros(shape2, dtype=bool)
+        rw_brdf_full     = rw_data.copy()
+    else:
+        # Only allocate the requested output array
+        output_full = np.full(shape3, np.nan, dtype=np.float32)
+        # We may also need C_brdf_fail if requested? Not for output_var, but could be used internally? Skip.
+        # For simplicity, we do not allocate extra arrays.
+
+    # ------------------------------------------------------------------
+    # 5–6. Chunked BRDF correction
+    # ------------------------------------------------------------------
+    n_chunks = int(np.ceil(n_water / chunk_size))
+
+    # 初始化进度条
+    iterator = range(n_chunks)
+    if progress_bar and tqdm is not None:
+        iterator = tqdm(iterator, desc="Processing chunks", unit="chunk")
+    elif progress_bar:
+        print("Warning: tqdm not installed, progress bar disabled.")
+
+    for chunk_idx in iterator:
+        lo = chunk_idx * chunk_size
+        hi = min(lo + chunk_size, n_water)
+
+        rows_c = water_rows[lo:hi]
+        cols_c = water_cols[lo:hi]
+
+        # Extract batch vectors
+        rw_c  = rw_data[:, rows_c, cols_c].T          # (chunk, bands)
+        sza_c = angles['sza'][rows_c, cols_c]          # (chunk,)
+        vza_c = angles['vza'][rows_c, cols_c]
+        raa_c = angles['raa'][rows_c, cols_c]
+
+        pixel_idx_c = np.arange(hi - lo)
+        ds_chunk = xr.Dataset(
+            {
+                'Rw':  xr.DataArray(rw_c,  dims=('n', 'bands'),
+                                    coords={'n': pixel_idx_c, 'bands': wvl}),
+                'sza': xr.DataArray(sza_c, dims='n',
+                                    coords={'n': pixel_idx_c}),
+                'vza': xr.DataArray(vza_c, dims='n',
+                                    coords={'n': pixel_idx_c}),
+                'raa': xr.DataArray(raa_c, dims='n',
+                                    coords={'n': pixel_idx_c}),
+            }
+        )
+
+        # Run BRDF model with uncertainty control
+        ds_corr = brdf_prototype(ds_chunk, brdf_model=brdf_model, compute_uncertainty=compute_uncertainty)
+
+        if keep_all:
+            # Write all results back
+            nrrs_full[:, rows_c, cols_c]        = ds_corr['nrrs'].values.T
+            rho_ex_w_full[:, rows_c, cols_c]    = ds_corr['rho_ex_w'].values.T
+            C_brdf_full[:, rows_c, cols_c]      = ds_corr['C_brdf'].values.T
+            if compute_uncertainty:
+                brdf_unc_full[:, rows_c, cols_c]    = ds_corr['brdf_unc'].values.T
+                nrrs_unc_full[:, rows_c, cols_c]    = ds_corr['nrrs_unc'].values.T
+            C_brdf_fail_full[rows_c, cols_c]    = ds_corr['C_brdf_fail'].values
+            rw_brdf_full[:, rows_c, cols_c]     = ds_corr['rho_ex_w'].values.T
+        else:
+            # Only write the requested output field
+            values = ds_corr[field_name].values.T
+            output_full[:, rows_c, cols_c] = values
+
+        # Release batch objects so GC can reclaim memory
+        del rw_c, sza_c, vza_c, raa_c, pixel_idx_c, ds_chunk, ds_corr
+        gc.collect()
+
+    # ------------------------------------------------------------------
+    # 7. Assemble full-scene output Dataset
+    # ------------------------------------------------------------------
+    full_coords = {'bands': wvl, 'y': np.arange(n_rows), 'x': np.arange(n_cols)}
+    dims3 = ('bands', 'y', 'x')
+    dims2 = ('y', 'x')
+
+    if keep_all:
+        # Original behaviour: return all variables
+        ds_out = xr.Dataset(
+            {
+                'Rw':          xr.DataArray(rw_data,          dims=dims3, coords=full_coords),
+                'Rw_brdf':     xr.DataArray(rw_brdf_full,     dims=dims3, coords=full_coords),
+                'nrrs':        xr.DataArray(nrrs_full,        dims=dims3, coords=full_coords),
+                'rho_ex_w':    xr.DataArray(rho_ex_w_full,    dims=dims3, coords=full_coords),
+                'C_brdf':      xr.DataArray(C_brdf_full,      dims=dims3, coords=full_coords),
+                'C_brdf_fail': xr.DataArray(C_brdf_fail_full, dims=dims2),
+                'water_mask':  xr.DataArray(water_mask,       dims=dims2),
+                'sza':         xr.DataArray(angles['sza'],    dims=dims2),
+                'saa':         xr.DataArray(angles['saa'],    dims=dims2),
+                'vza':         xr.DataArray(angles['vza'],    dims=dims2),
+                'vaa':         xr.DataArray(angles['vaa'],    dims=dims2),
+                'raa':         xr.DataArray(angles['raa'],    dims=dims2),
+            }
+        )
+        if compute_uncertainty:
+            ds_out['brdf_unc'] = xr.DataArray(brdf_unc_full, dims=dims3, coords=full_coords)
+            ds_out['nrrs_unc'] = xr.DataArray(nrrs_unc_full, dims=dims3, coords=full_coords)
+    else:
+        # Return only the requested output variable + geometry + mask
+        ds_out = xr.Dataset(
+            {
+                output_var:    xr.DataArray(output_full,      dims=dims3, coords=full_coords),
+                'Rw':          xr.DataArray(rw_data,          dims=dims3, coords=full_coords),
+                'water_mask':  xr.DataArray(water_mask,       dims=dims2),
+                'sza':         xr.DataArray(angles['sza'],    dims=dims2),
+                'saa':         xr.DataArray(angles['saa'],    dims=dims2),
+                'vza':         xr.DataArray(angles['vza'],    dims=dims2),
+                'vaa':         xr.DataArray(angles['vaa'],    dims=dims2),
+                'raa':         xr.DataArray(angles['raa'],    dims=dims2),
+            }
+        )
+        # If uncertainties were computed, we could add them optionally, but user didn't request them.
+        # For memory efficiency, we skip them here.
+
+    return ds_out
+
+
+def process_brdf_files(
+    hyperspectral_file,
+    angle_file,
+    output_file,
+    mask_file=None,
+    wavelengths=None,
+    brdf_model="L11",
+    output_vars=None,
+    output_format="ENVI",
+    chunk_size=4096,
+    block_size=512,
+    compute_uncertainty=False,
+    progress_bar=True,
+):
+    """Run BRDF correction from input files and save multiple variables to separate files.
+
+    The geometry file must contain 10 bands as documented for :func:`read_bip_angles`.
+    If ``mask_file`` is None, the full scene is corrected (same as :func:`run_brdf_correction_block`).
+    """
+    if output_format.upper() != "ENVI":
+        raise NotImplementedError("Only ENVI output format is supported in block-wise version.")
+
+    if output_vars is None:
+        output_vars = ["Rw_brdf"]
+    if isinstance(output_vars, str):
+        output_vars = [output_vars]
+
+    print(f"Computing BRDF correction with variables: {output_vars}")
+    print(f"Output base name: {output_file}")
+
+    saved_paths = run_brdf_correction_block(
+        hyperspectral_file=hyperspectral_file,
+        angle_file=angle_file,
+        mask_file=mask_file,
+        output_file=output_file,
+        wavelengths=wavelengths,
+        brdf_model=brdf_model,
+        output_vars=output_vars,
+        block_size=block_size,
+        compute_uncertainty=compute_uncertainty,
+        progress_bar=progress_bar,
+    )
+    return None, saved_paths
+
+
+def main():
+    import argparse
+
+    parser = argparse.ArgumentParser(description="Water-region BRDF correction for hyperspectral data.")
+    parser.add_argument("hyperspectral_file", help="Input ENVI BSQ hyperspectral file or .hdr path.")
+    parser.add_argument(
+        "angle_file",
+        help="10-band ENVI geometry file or .hdr: path length, VAA, VZA, SAA, SZA, "
+        "phase, slope, aspect, cos i, UTC (RAA computed from VAA and SAA).",
+    )
+    parser.add_argument("output_file", help="Output path prefix. ENVI writes .hdr/.img, NetCDF writes .nc.")
+    parser.add_argument(
+        "--mask",
+        dest="mask_file",
+        default=None,
+        metavar="PATH",
+        help="Optional water mask GeoTIFF; pixels == 1 are corrected. "
+        "If omitted, BRDF correction is applied to the full image.",
+    )
+    parser.add_argument("--brdf-model", default="L11", choices=["L11", "M02", "M02SeaDAS", "O25"])
+    parser.add_argument("--output-var", nargs='+', default=["Rw_brdf"],
+                        help="Output variables to save (one or more). "
+                             "Supported: Rw_brdf, rho_ex_w, nrrs, C_brdf, brdf_unc, nrrs_unc. "
+                             "Example: --output-var Rw_brdf nrrs brdf_unc")
+    parser.add_argument("--output-format", default="ENVI", choices=["ENVI"])
+    parser.add_argument("--chunk-size", type=int, default=4096,
+                        help="Number of water pixels per processing batch (default: 4096).")
+    parser.add_argument("--block-size", type=int, default=512,
+                        help="Spatial block size for block-wise processing (default: 512).")
+    parser.add_argument("--compute-uncertainty", action="store_true",
+                        help="Compute BRDF uncertainty (may increase runtime and output size).")
+    parser.add_argument("--no-progress-bar", action="store_true",
+                        help="Disable progress bar (useful for non-interactive environments).")
+    args = parser.parse_args()
+
+    progress_bar = not args.no_progress_bar
+
+    _ds_out, saved_paths = process_brdf_files(
+        hyperspectral_file=args.hyperspectral_file,
+        angle_file=args.angle_file,
+        mask_file=args.mask_file,
+        output_file=args.output_file,
+        brdf_model=args.brdf_model,
+        output_vars=args.output_var,
+        output_format=args.output_format,
+        chunk_size=args.chunk_size,
+        block_size=args.block_size,
+        compute_uncertainty=args.compute_uncertainty,
+        progress_bar=progress_bar,
+    )
+
+    if isinstance(saved_paths, list):
+        print(f"BRDF correction finished. Saved {len(saved_paths)} files:")
+        for path in saved_paths:
+            print(f"  - {path}")
+    else:
+        print(f"BRDF correction finished. Saved to: {saved_paths}")
+
+if __name__ == "__main__":
+    main()

ocbrdf/brdf_model_L11.py

@@ -0,0 +1,191 @@
+import numpy as np
+import xarray as xr
+
+from brdf_utils import ADF_OCP, solve_2nd_order_poly, drop_unused_coords
+from Raman import Raman
+
+
+''' Lee et al. (2011) BRDF correction
+    With QAA v6 inversion
+
+'''
+
+# Init Raman class
+Raman = Raman()
+
+""" Class for L11 coefficients """
+class Coeffs():
+    def __init__(self,Gw0,Gw1,Gp0,Gp1):
+        self.Gw0 = Gw0
+        self.Gw1 = Gw1
+        self.Gp0 = Gp0
+        self.Gp1 = Gp1
+
+""" Class for L11 BRDF model """
+class L11:
+
+    """ Initialise L11 model: BRDF LUT, coeffs, QAA parameters, water IOPs LUT
+        Note: bands are fixed and defined at class initilization, but could be initialized in init_pixels if needed
+    """
+    def __init__(self, bands, adf=None):
+        if adf is None:
+            adf = ADF_OCP
+
+        # Check required bands are existing, within a 10 nm threshold
+        self.bands = bands
+        threshold = 10.
+        bands_required = [442, 490, 560, 665]
+        bands_ref = bands.sel(bands=bands_required, method='nearest')
+        for band_ref, band_required in zip(bands_ref, bands_required):
+            assert abs(band_ref - band_required) < threshold, 'Band %d nm missing or too far'% band_ref
+        self.b442, self.b490, self.b560, self.b665 = bands_ref
+
+        # Read BRDF LUT and compute default coeffs
+        LUT_OCP = xr.open_dataset(adf % 'L11',engine='netcdf4')
+        self.LUT = xr.Dataset()
+        self.LUT['Gw0'] = LUT_OCP.Gw0
+        self.LUT['Gw1'] = LUT_OCP.Gw1
+        self.LUT['Gp0'] = LUT_OCP.Gp0
+        self.LUT['Gp1'] = LUT_OCP.Gp1
+
+        self.coeffs0 = self.interp(0.,0.,0.)
+        self.coeffs = Coeffs(np.nan,np.nan,np.nan,np.nan)
+
+        # Read IOPs of pure water (store in LUT for further spectral interpolation)
+        self.awLUT = LUT_OCP.aw.rename({'IOP_wl':'bands'})
+        self.bbwLUT = LUT_OCP.bbw.rename({'IOP_wl':'bands'})
+
+        # Read QAA parameters
+        self.a0G = LUT_OCP.a0G.values
+        self.a0R = LUT_OCP.a0R.values
+        self.gamma = LUT_OCP.gamma.values
+        self.niter = LUT_OCP.niter.values
+              
+    """ Initialize pixel: coefficient at current geometry and water IOP at current bands """
+    def init_pixels(self, theta_s, theta_v, delta_phi):
+        self.coeffs = self.interp(theta_s, theta_v, delta_phi)
+
+        # Compute IOPs at current bands
+        self.aw = self.awLUT.interp(bands = self.bands, kwargs={'fill_value':'extrapolate'})
+        self.bbw = self.bbwLUT.interp(bands = self.bands, kwargs={'fill_value':'extrapolate'})
+
+    """ Interpolate coefficients """
+    def interp(self, theta_s, theta_v, delta_phi):
+        Gw0 = self.LUT.Gw0.interp(theta_s=theta_s,theta_v=theta_v,delta_phi=delta_phi)
+        Gw1 = self.LUT.Gw1.interp(theta_s=theta_s,theta_v=theta_v,delta_phi=delta_phi)
+        Gp0 = self.LUT.Gp0.interp(theta_s=theta_s,theta_v=theta_v,delta_phi=delta_phi)
+        Gp1 = self.LUT.Gp1.interp(theta_s=theta_s,theta_v=theta_v,delta_phi=delta_phi)
+        return Coeffs(Gw0,Gw1,Gp0,Gp1)
+
+    """ Compute remote-sensing reflectance, without Raman effect (vanish in the normalization factor) """
+    def forward(self, ds, normalized=False):
+        omega_b = ds['omega_b']
+        eta_b = ds['eta_b']
+
+        if normalized:
+            coeffs = self.coeffs0
+        else:
+            coeffs = self.coeffs
+        mod_Rrs = (coeffs.Gw0+coeffs.Gw1*omega_b*eta_b)*omega_b*eta_b + (coeffs.Gp0+coeffs.Gp1*omega_b*(1-eta_b))*omega_b*(1-eta_b)
+
+        return mod_Rrs
+
+    """ Apply QAA to retrieve IOP (omega_b, eta_b) from Rrs """
+    def backward(self, ds, iter_brdf):
+
+        Rrs = ds['nrrs']
+
+        # Select G coeff according to iteration
+        if iter_brdf == 0:
+            coeffs = self.coeffs
+        else:
+            coeffs = self.coeffs0
+
+        # Apply Raman correction
+        Rrs = Raman.correct(Rrs)
+       
+        # Local renaming of bands 
+        b442, b490, b560, b665 = self.b442, self.b490, self.b560, self.b665
+
+        # Apply upper and lower limits to Rrs(665) #TODO check if not finite or missing?
+        Rrs442 = Rrs.sel(bands=b442)
+        Rrs490 = Rrs.sel(bands=b490)
+        Rrs560 = Rrs.sel(bands=b560)
+        Rrs665 = Rrs.sel(bands=b665)
+        mask= ((Rrs665 > 20*np.power(Rrs560,1.5)) | (Rrs665 < 0.9*np.power(Rrs560, 1.7)))
+        if np.any(mask):
+            Rrs665_ = 1.27*np.power(Rrs560, 1.47) + 0.00018*np.power(Rrs490/Rrs560,-3.19)
+            # Redefine Rrs665 and Rrs[bands=b665] (both important for computations below)
+            Rrs665 = xr.where(mask, Rrs665_, Rrs665)
+            Rrs.loc[dict(bands=b665)] = Rrs665
+
+        # Calculate rrs below water for absorption computation
+        rrs = Rrs / (0.52 + 1.7*Rrs)
+
+        # Define reference band band0 according to Rrs at 665 nm
+        # and compute total absorption
+        mask = Rrs.sel(bands=b665) < 0.0015
+        band0 = xr.where(mask, b560, b665)
+        aw0 = xr.where(mask, self.aw.sel(bands=b560), self.aw.sel(bands=b665))
+        bbw0 = xr.where(mask, self.bbw.sel(bands=b560), self.bbw.sel(bands=b665))
+        Rrs0 = xr.where(mask, Rrs.sel(bands=b560), Rrs.sel(bands=b665))
+        # Compute a0 when band0 = b560
+        rrs442 = rrs.sel(bands=b442)
+        rrs490 = rrs.sel(bands=b490)
+        rrs560 = rrs.sel(bands=b560)
+        rrs665 = rrs.sel(bands=b665)
+        chi = np.log10((rrs442 + rrs490) / (rrs560 + 5.0 * rrs665*rrs665 / rrs490))
+        poly = np.polynomial.polynomial.polyval(chi, self.a0G)
+        a0_560 = aw0 + np.power(10., poly)
+        # Compute a0 when band0 = b665
+        a0_665 = aw0 + self.a0R[0] * np.power(Rrs665 / (Rrs442 + Rrs490), self.a0R[1])
+        # Compute a0 for all pixels
+        a0 = xr.where(mask, a0_560, a0_665)
+
+        # Compute bbp at band0 by 2nd order polynomial inversion
+        k0 = a0 + bbw0
+        cA = coeffs.Gp0 + coeffs.Gp1 - Rrs0
+        cB = coeffs.Gw0 * bbw0 + (coeffs.Gp0 -2*Rrs0) *k0
+        cC = (coeffs.Gw0 * bbw0 - Rrs0 * k0) * k0 + coeffs.Gw1 * bbw0 * bbw0
+        bbp0 = solve_2nd_order_poly(cA, cB, cC)
+
+        # Assume bbp0 = 0 if solve_2nd_order_poly fails to retrieve non-negative numbers
+        # In this case activate bbp0_fail flag (which in turn activates QAA_fail)
+        bbp0_fail = (bbp0 < 0) | (np.isinf(bbp0)) | (np.isnan(bbp0))
+        bbp0 = xr.where(bbp0_fail, 0, bbp0)
+
+        # Compute bbp slope and extrapolate at all bands
+        gamma = self.gamma[0] * (1.0 - self.gamma[1] * np.exp(-self.gamma[2] * (rrs442 / rrs560)))
+        bbp = bbp0 * np.power(band0 / self.bands, gamma)
+
+        # Compute total bb
+        bb = self.bbw + bbp
+
+        # Compute quasi-diffuse attenuation coefficient k at each band
+        # by 2nd order polynomial inversion
+        cA = Rrs
+        cB = - (coeffs.Gw0 * self.bbw + coeffs.Gp0 * bbp)
+        cC = - (coeffs.Gw1 * self.bbw *self.bbw + coeffs.Gp1 * bbp * bbp)
+        k = solve_2nd_order_poly(cA, cB, cC)
+
+        # If k is not a positive value, then
+        #   i) k --> bbw + aw
+        #   ii) k_fail flag is activated
+        k_fail = (k <= 0) | (np.isinf(k)) | (np.isnan(k))
+        k = xr.where(k_fail, self.aw + self.bbw, k)
+
+        # Drop unused coords to avoid issues
+        bb        = drop_unused_coords(bb)
+        k         = drop_unused_coords(k)
+        k_fail    = drop_unused_coords(k_fail)
+        bbp0_fail = drop_unused_coords(bbp0_fail)
+
+        # Set QAA_fail is either bbp0_fail or k_fail are activated
+        ds['QAA_fail'] = (bbp0_fail) | (k_fail)
+
+        # Compute final IOPs
+        ds['omega_b'] = bb / k
+        ds['eta_b'] = self.bbw / bb
+
+        return ds
+

ocbrdf/brdf_model_M02.py

@@ -0,0 +1,182 @@
+import numpy as np
+import xarray as xr
+
+from brdf_utils import ADF_OCP, solve_2nd_order_poly, drop_unused_coords
+
+''' Morel et al. (2002) BRDF 校正
+    包含 R gothic
+    f/Q 查找表，与 OLCI 业务处理器中的一致（参见 OLCI "OCP" ADF）
+    使用 OC4ME 叶绿素反演，进行两次迭代，未应用叶绿素差值收敛标准（1%）
+'''
+
+""" M02 系数的类 """
+
+
+class Coeffs():
+    def __init__(self, Rgoth, foq):
+        self.Rgoth = Rgoth
+        self.foq = foq
+
+
+""" Class for M02 BRDF model """
+class M02:
+    """ Initialise M02 model: BRDF LUT, coeffs, OC4ME parameters, water IOPs LUT
+        Note: bands are fixed and defined at class initilization, but could be initialized in init_pixels if needed
+    """
+
+    def __init__(self, bands, aot, wind, adf=None):
+        if adf is None:
+            adf = ADF_OCP
+
+        # Check required bands are existing, within a 25 nm threshold
+        self.bands = bands
+        threshold = 25.
+        bands_required = [442.5, 490, 510, 560]
+        bands_ref = bands.sel(bands=bands_required, method='nearest')
+        for band_ref, band_required in zip(bands_ref, bands_required):
+            assert abs(band_ref - band_required) < threshold, 'Band %d nm missing or too far' % band_ref
+        self.b442, self.b490, self.b510, self.b560 = bands_ref
+
+        # Read BRDF LUT and compute default coeffs
+        LUT_OCP = xr.open_dataset(adf % 'M02',engine='netcdf4')
+        self.LUT = xr.Dataset()
+
+        # Homogeneise naming convention with other methods... (PZA --> OZA transformation comes below...)
+        self.LUT['Rgoth'] = LUT_OCP.r_goth_LUT.rename({'PZA_r_goth':'theta_v_Rgoth',
+                                                       'wind_speeds_r_goth':'wind_Rgoth'})
+        self.LUT['foq'] = LUT_OCP.f_over_q_LUT.rename({'SZA_FOQ':'theta_s',
+                                                       'PZA_FOQ':'theta_v',
+                                                       'RAA_FOQ':'delta_phi',
+                                                       'wind_speeds_FOQ':'wind_foq',
+                                                       'tau_a_FOQ':'aot_foq',
+                                                       'log_chl_FOQ':'log_chl_foq'})
+
+        # f0 factor to convert Rrs to Irradiance Reflectance (R) --> in order to apply OC4ME
+        self.LUT['f0'] = LUT_OCP.f0_LUT
+
+        # Index of refraction
+        self.n_w = float(LUT_OCP.water_refraction_index.data)
+
+        # Remove trivial aot indexation
+        self.LUT['foq'] = self.LUT['foq'].squeeze()
+
+        # Interpolate for wind speed.
+        self.LUT['Rgoth'] = self.LUT['Rgoth'].interp(wind_Rgoth=np.clip(wind,0,16))
+        self.LUT['foq'  ] = self.LUT['foq'  ].interp(wind_foq  =np.clip(wind,0,10))
+
+        #
+        self.coeffs0 = self.interp_geometries(0., 0., 0.)
+        self.coeffs = Coeffs(np.nan, np.nan)
+
+        # Parameters for the OC4ME chl retrieval
+        self.OC4MEcoeff = LUT_OCP.log10_coeff_LUT.values
+        self.OC4MEepsilon = LUT_OCP.oc4me_epsilon
+        self.OC4MEchl0 = float(LUT_OCP.oc4me_chl0.values)
+        self.niter = LUT_OCP.oc4me_niter
+
+    """ Initialize pixel: coefficient at current geometry and water IOP at current bands """
+
+    def init_pixels(self, theta_s, theta_v, delta_phi):
+        self.coeffs = self.interp_geometries(theta_s, theta_v, delta_phi)
+
+    """ Interpolate coefficients """
+    def interp_geometries(self, theta_s, theta_v, delta_phi):
+        # Transform PZA to VZA (Snell's refraction Law)
+        theta_v = np.rad2deg(np.arcsin(np.sin(np.deg2rad(theta_v)) / self.n_w))
+        theta_v_Rgoth_0 = np.clip(theta_v,
+                                   float(np.min(self.LUT.theta_v_Rgoth)),
+                                   float(np.max(self.LUT.theta_v_Rgoth)))
+
+        theta_v_0 = np.clip(theta_v,
+                                   float(np.min(self.LUT.theta_v)),
+                                   float(np.max(self.LUT.theta_v)))
+
+        Rgoth = self.LUT.Rgoth.interp(theta_v_Rgoth=theta_v_Rgoth_0)
+        foq = self.LUT.foq.interp(theta_s=theta_s, theta_v=theta_v_0, delta_phi=delta_phi)
+
+        return Coeffs(Rgoth, foq)
+
+    """ Compute remote-sensing reflectance"""
+
+    def forward(self, ds, normalized=False):
+
+        if normalized:
+            coeffs = self.coeffs0
+        else:
+            coeffs = self.coeffs
+
+        wave_foq = np.clip(ds['bands'],
+                           float(np.min(coeffs.foq.wavelengths_FOQ)),
+                           float(np.max(coeffs.foq.wavelengths_FOQ)))
+
+        log10_chl_foq = np.clip(ds['log10_chl'],
+                           float(np.min(coeffs.foq.log_chl_foq)/np.log(10)),
+                           float(np.max(coeffs.foq.log_chl_foq)/np.log(10)))
+
+        # f/Q LUT indexed with ln(CHL), i.e. log_e(CHL)
+        log_chl_foq = log10_chl_foq * np.log(10)
+
+        forward_mod = coeffs.Rgoth * coeffs.foq.interp(wavelengths_FOQ=wave_foq).interp(log_chl_foq=log_chl_foq)
+
+        return forward_mod
+
+    """ Apply QAA to retrieve IOP (omega_b, eta_b) from Rrs """
+
+    def backward(self, ds, iter_brdf):
+
+        Rrs = ds['nrrs']
+
+        # Local renaming of bands
+        b442, b490, b510, b560 = self.b442, self.b490, self.b510, self.b560
+
+        # Convert to scalar if np.array of 1 value to avoid issues
+        try:
+            b442=b442.item()
+            b490=b490.item()
+            b510=b510.item()
+            b560=b560.item()
+        except:
+            pass
+
+        # Clip to f0 log10(CHL) values
+        log10_chl_f0 = np.clip(ds['log10_chl'],
+                           float(np.min(self.LUT['f0'].log_chl_f0)/np.log(10)),
+                           float(np.max(self.LUT['f0'].log_chl_f0)/np.log(10)))
+
+        # f/Q LUT indexed with ln(CHL), i.e. log_e(CHL)
+        log_chl_f0 = log10_chl_f0 * np.log(10)
+        f0_chl = self.LUT['f0'].interp(log_chl_f0=log_chl_f0)
+
+        fQ_chl = self.coeffs.foq.interp(log_chl_foq=log10_chl_f0)
+
+        # Drop unused coordinates to avoid ambiguities in indexation...
+        Rrs = drop_unused_coords(Rrs)
+        f0_chl = drop_unused_coords(f0_chl)
+        fQ_chl = drop_unused_coords(fQ_chl)
+
+        # Get Rrs at relevant bands for OC4ME and convert to Irradiance Reflectance (R)
+        R442 = np.pi * Rrs.sel(bands=b442) * f0_chl.sel(bands_f0=2.0) / fQ_chl.interp(wavelengths_FOQ=b442)
+        R490 = np.pi * Rrs.sel(bands=b490) * f0_chl.sel(bands_f0=3.0) / fQ_chl.interp(wavelengths_FOQ=b490)
+        R510 = np.pi * Rrs.sel(bands=b510) * f0_chl.sel(bands_f0=4.0) / fQ_chl.interp(wavelengths_FOQ=b510)
+        R560 = np.pi * Rrs.sel(bands=b560) * f0_chl.sel(bands_f0=5.0) / fQ_chl.interp(wavelengths_FOQ=b560)
+
+        R442 = drop_unused_coords(R442)
+        R490 = drop_unused_coords(R490)
+        R510 = drop_unused_coords(R510)
+        R560 = drop_unused_coords(R560)
+
+        # Compute the OC4ME "R"
+        # NB: Could be the case that R[442-560]<0,
+        #   then ds['log10_chl_OC4ME_Ratio'] will be NaN.
+        #   then ds['log10_chl'] will be NaN.
+        #   then forward_mod0/forward_mod (in brdf_prototype, ocbrdf/main.py) will be NaN.
+        #   then ds['C_BRDF'] = NaN.
+        #   In this case, C_BRDF will be set to 1 and C_BRDF_flag will be raised (see brdf_prototype, ocbrdf/main.py).
+        ds['log10_chl_OC4ME_Ratio'] = np.log10(np.max([R442, R490, R510], axis=0) / R560)
+
+        ds['log10_chl'] = 0 * ds['log10_chl_OC4ME_Ratio']
+        # Apply OC4ME 5-degree polynomial
+        for k, Ak in enumerate(self.OC4MEcoeff):
+            ds['log10_chl'] += Ak * (ds['log10_chl_OC4ME_Ratio'] ** k)
+
+        return ds

ocbrdf/brdf_model_M02SeaDAS.py

@@ -0,0 +1,157 @@
+import numpy as np
+import xarray as xr
+
+from brdf_utils import ADF_OCP, solve_2nd_order_poly, drop_unused_coords
+''' Morel et al. (2002) BRDF 校正
+    未包含 R gothic
+    f/Q 查找表与 NASA SeaDAS 中的一致（详见 BRDF_M02SeaDAS.nc 属性）
+    使用 OC4ME 叶绿素反演，进行一次迭代，未应用收敛标准
+'''
+
+""" M02 系数的类 """
+
+
+class Coeffs():
+    def __init__(self, foq):
+        self.foq = foq
+
+
+""" Class for M02 BRDF model """
+class M02SeaDAS:
+    """ Initialise M02 model: BRDF LUT, coeffs, OC4ME parameters, water IOPs LUT
+        Note: bands are fixed and defined at class initilization, but could be initialized in init_pixels if needed
+    """
+
+    def __init__(self, bands, adf=None):
+        if adf is None:
+            adf = ADF_OCP
+
+        # Check required bands are existing, within a 25 nm threshold
+        self.bands = bands
+        threshold = 25.
+        bands_required = [442.5, 490, 510, 560]
+        bands_ref = bands.sel(bands=bands_required, method='nearest')
+        for band_ref, band_required in zip(bands_ref, bands_required):
+            assert abs(band_ref - band_required) < threshold, 'Band %d nm missing or too far' % band_ref
+        self.b442 = float(bands_ref[0].item())
+        self.b490 = float(bands_ref[1].item())
+        self.b510 = float(bands_ref[2].item())
+        self.b560 = float(bands_ref[3].item())
+
+        # Read BRDF LUT and compute default coeffs
+        LUT_OCP = xr.open_dataset(adf % 'M02SeaDAS',engine='netcdf4')
+        self.LUT = xr.Dataset()
+
+        # Homogeneise naming convention with other methods... (PZA --> OZA transformation comes below...)
+        self.LUT['foq'] = LUT_OCP.f_over_q_LUT.rename({'SZA_FOQ':'theta_s',
+                                                       'PZA_FOQ':'theta_v',
+                                                       'RAA_FOQ':'delta_phi',
+                                                       'log_chl_FOQ':'log_chl_foq'})
+
+        # Index of refraction
+        self.n_w = float(LUT_OCP.water_refraction_index.data)
+
+        # Remove trivial aot indexation
+        self.LUT['foq'] = self.LUT['foq'].squeeze()
+
+        #
+        self.coeffs0 = self.interp_geometries(0., 0., 0.)
+        self.coeffs = Coeffs(np.nan)
+
+        # Parameters for the OC4ME chl retrieval
+        self.OC4MEcoeff = LUT_OCP.log10_coeff_LUT.values
+        self.OC4MEepsilon = LUT_OCP.oc4me_epsilon
+        self.OC4MEchl0 = float(LUT_OCP.oc4me_chl0.values)
+        self.niter = LUT_OCP.oc4me_niter
+
+    """ Initialize pixel: coefficient at current geometry and water IOP at current bands """
+
+    def init_pixels(self, theta_s, theta_v, delta_phi):
+        self.coeffs = self.interp_geometries(theta_s, theta_v, delta_phi)
+        
+        # Cache wavelength and chlorophyll bounds for fast access in forward()
+        self.min_wavelength_FOQ = float(np.min(self.coeffs.foq.wavelengths_FOQ))
+        self.max_wavelength_FOQ = float(np.max(self.coeffs.foq.wavelengths_FOQ))
+        self.min_log_chl_foq = float(np.min(self.coeffs.foq.log_chl_foq) / np.log(10))
+        self.max_log_chl_foq = float(np.max(self.coeffs.foq.log_chl_foq) / np.log(10))
+
+    """ Interpolate coefficients """
+    def interp_geometries(self, theta_s, theta_v, delta_phi):
+        # Transform PZA to VZA (Snell's refraction Law)
+        theta_v = np.rad2deg(np.arcsin(np.sin(np.deg2rad(theta_v)) / self.n_w))
+        theta_v_0 = np.clip(theta_v,
+                                   float(np.min(self.LUT.theta_v)),
+                                   float(np.max(self.LUT.theta_v)))
+
+        foq = self.LUT.foq.interp(theta_s=theta_s, theta_v=theta_v_0, delta_phi=delta_phi)
+
+        return Coeffs(foq)
+
+    """ Compute remote-sensing reflectance"""
+
+    def forward(self, ds, normalized=False):
+
+        if normalized:
+            coeffs = self.coeffs0
+            min_wl = float(np.min(coeffs.foq.wavelengths_FOQ))
+            max_wl = float(np.max(coeffs.foq.wavelengths_FOQ))
+            min_chl = float(np.min(coeffs.foq.log_chl_foq) / np.log(10))
+            max_chl = float(np.max(coeffs.foq.log_chl_foq) / np.log(10))
+        else:
+            coeffs = self.coeffs
+            min_wl = self.min_wavelength_FOQ
+            max_wl = self.max_wavelength_FOQ
+            min_chl = self.min_log_chl_foq
+            max_chl = self.max_log_chl_foq
+
+        wave_foq = np.clip(ds['bands'], min_wl, max_wl)
+        log10_chl_foq = np.clip(ds['log10_chl'], min_chl, max_chl)
+        
+        # f/Q LUT indexed with ln(CHL), i.e. log_e(CHL)
+        log_chl_foq = log10_chl_foq * np.log(10)
+
+        # Merged dual interpolation into single operation for better performance
+        forward_mod = coeffs.foq.interp(
+            wavelengths_FOQ=wave_foq,
+            log_chl_foq=log_chl_foq,
+            method='linear'
+        )
+
+        return forward_mod
+
+    """ Apply QAA to retrieve IOP (omega_b, eta_b) from Rrs """
+
+    def backward(self, ds, iter_brdf):
+
+        Rrs = ds['nrrs']
+
+        # Local renaming of bands (pre-converted to float in __init__)
+        b442 = self.b442
+        b490 = self.b490
+        b510 = self.b510
+        b560 = self.b560
+
+        # Apply upper and lower limits to Rrs(665) #TODO check if not finite or missing?
+        Rrs442 = Rrs.sel(bands=b442)
+        Rrs490 = Rrs.sel(bands=b490)
+        Rrs510 = Rrs.sel(bands=b510)
+        Rrs560 = Rrs.sel(bands=b560)
+        
+        # Drop unused coordinates to reduce overhead
+        Rrs442 = drop_unused_coords(Rrs442)
+        Rrs490 = drop_unused_coords(Rrs490)
+        Rrs510 = drop_unused_coords(Rrs510)
+        Rrs560 = drop_unused_coords(Rrs560)
+
+        # Compute the OC4ME "R"
+        ds['log10_chl_OC4ME_Ratio'] = np.log10(np.max([Rrs442, Rrs490, Rrs510], axis=0) / Rrs560)
+
+        ds['log10_chl'] = 0 * ds['log10_chl_OC4ME_Ratio']
+        # Apply OC4ME 5-degree polynomial using numpy for better performance
+        poly = np.polynomial.polynomial.polyval(
+            ds['log10_chl_OC4ME_Ratio'].values,
+            self.OC4MEcoeff
+        )
+        ds['log10_chl'].values = poly
+
+        return ds

ocbrdf/brdf_model_O25.py

@@ -0,0 +1,170 @@
+import numpy as np
+import xarray as xr
+
+from brdf_utils import ADF_OCP, solve_2nd_order_poly, drop_unused_coords
+from Raman import Raman
+
+
+''' O25 BRDF correction from EUMETSAT BRDF4OLCI study
+    Ref. ATBD_EUM-CO-21-4600002626-JIG, 29/03/2024
+'''
+
+# Init Raman class
+Raman = Raman()
+
+""" Class for O25 coefficients """
+class Coeffs():
+    def __init__(self,Gw0,Gw1,Gp0,Gp1):
+        self.Gw0 = Gw0
+        self.Gw1 = Gw1
+        self.Gp0 = Gp0
+        self.Gp1 = Gp1
+
+""" Class for O25 BRDF model """
+class O25:
+
+    """ Initialise O25 model: BRDF LUT, coeffs, QAA parameters, water IOPs LUT
+        Note: bands are fixed and defined at class initilization, but could be initialized in init_pixels if needed
+    """
+    def __init__(self, bands, adf=None):
+        if adf is None:
+            adf = ADF_OCP
+
+        # Check required bands are existing, within a 10 nm threshold
+        self.bands = bands
+        threshold = 10.
+        bands_required = [442, 490, 560, 665]
+        bands_ref = bands.sel(bands=bands_required, method='nearest')
+        for band_ref, band_required in zip(bands_ref, bands_required):
+            assert abs(band_ref - band_required) < threshold, 'Band %d nm missing or too far'%band_ref
+        self.b442, self.b490, self.b560, self.b665 = bands_ref
+
+        # Read BRDF LUT and compute default coeffs
+        LUT_OCP = xr.open_dataset(adf % 'O25',engine='netcdf4')
+        self.LUT = xr.Dataset()
+        self.LUT['Gw0'] = LUT_OCP.Gw0
+        self.LUT['Gw1'] = LUT_OCP.Gw1 
+        self.LUT['Gp0'] = LUT_OCP.Gp0 
+        self.LUT['Gp1'] = LUT_OCP.Gp1 
+
+        self.coeffs0 = self.interp(0.,0.,0.)
+        self.coeffs = Coeffs(np.nan,np.nan,np.nan,np.nan)
+
+        # Read IOPs of pure water (store in LUT for further spectral interpolation)
+        self.awLUT = LUT_OCP.aw.rename({'IOP_wl':'bands'})
+        self.bbwLUT = LUT_OCP.bbw.rename({'IOP_wl':'bands'})
+
+        # Read QAA parameters
+        self.a0 = LUT_OCP.a0.values
+        self.gamma = LUT_OCP.gamma.values
+        self.niter = LUT_OCP.niter.values
+              
+    """ Initialize pixel: coefficient at current geometry and water IOP at current bands """
+    def init_pixels(self, theta_s, theta_v, delta_phi):
+        self.coeffs = self.interp(theta_s, theta_v, delta_phi)
+
+        # Compute IOPs at current bands
+        self.aw = self.awLUT.interp(bands = self.bands, kwargs={'fill_value':'extrapolate'})
+        self.bbw = self.bbwLUT.interp(bands = self.bands, kwargs={'fill_value':'extrapolate'})
+
+    """ Interpolate coefficients """
+    def interp(self, theta_s, theta_v, delta_phi):
+        Gw0 = self.LUT.Gw0.interp(theta_s=theta_s,theta_v=theta_v,delta_phi=delta_phi)
+        Gw1 = self.LUT.Gw1.interp(theta_s=theta_s,theta_v=theta_v,delta_phi=delta_phi)
+        Gp0 = self.LUT.Gp0.interp(theta_s=theta_s,theta_v=theta_v,delta_phi=delta_phi)
+        Gp1 = self.LUT.Gp1.interp(theta_s=theta_s,theta_v=theta_v,delta_phi=delta_phi)
+        return Coeffs(Gw0,Gw1,Gp0,Gp1)
+
+    """ Compute remote-sensing reflectance, without Raman effect (vanish in the normalization factor) """
+    def forward(self, ds, normalized=False):
+
+        omega_b = ds['omega_b']
+        eta_b = ds['eta_b']
+
+        if normalized:
+            coeffs = self.coeffs0
+        else:
+            coeffs = self.coeffs
+        mod_Rrs = (coeffs.Gw0+coeffs.Gw1*omega_b*eta_b)*omega_b*eta_b + (coeffs.Gp0+coeffs.Gp1*omega_b*(1-eta_b))*omega_b*(1-eta_b)
+
+        return mod_Rrs
+
+    """ Apply QAA to retrieve IOP (omega_b, eta_b) from rrs """
+    def backward(self, ds, iter_brdf):
+
+        Rrs = ds['nrrs']
+
+        # Select G coeff according to iteration
+        if iter_brdf == 0:
+            coeffs = self.coeffs
+        else:
+            coeffs = self.coeffs0
+
+        # Apply Raman correction
+        Rrs = Raman.correct(Rrs)
+       
+        # Local renaming of bands 
+        b442, b490, b560, b665 = self.b442, self.b490, self.b560, self.b665
+
+        # Define reference band band0 at 560 nm
+        # and compute total absorption
+        Rrs0 = Rrs.sel(bands=b560)
+        band0 = xr.zeros_like(Rrs0) + b560
+        aw0 = xr.zeros_like(Rrs0) + self.aw.sel(bands=b560)
+        bbw0 = xr.zeros_like(Rrs0) + self.bbw.sel(bands=b560)
+        # Compute a0 when band0 = b560
+        Rrs442 = Rrs.sel(bands=b442)
+        Rrs490 = Rrs.sel(bands=b490)
+        Rrs560 = Rrs.sel(bands=b560)
+        Rrs665 = Rrs.sel(bands=b665)
+        chi = np.log10((Rrs442 + Rrs490) / (Rrs560 + 5.0 * Rrs665*Rrs665 / Rrs490))
+        poly = np.polynomial.polynomial.polyval(chi, self.a0)
+        a0 = aw0 + np.power(10., poly)
+
+        # Compute bbp at band0 by 2nd order polynomial inversion
+        k0 = a0 + bbw0
+        cA = coeffs.Gp0 + coeffs.Gp1 - Rrs0
+        cB = coeffs.Gw0 * bbw0 + (coeffs.Gp0 -2*Rrs0) *k0
+        cC = (coeffs.Gw0 * bbw0 - Rrs0 * k0) * k0 + coeffs.Gw1 * bbw0 * bbw0
+        bbp0 = solve_2nd_order_poly(cA, cB, cC)
+
+        # Assume bbp0 = 0 if solve_2nd_order_poly fails to retrieve non-negative numbers
+        # In this case activate bbp0_fail flag (which in turn activates QAA_fail)
+        bbp0_fail = (bbp0 < 0) | (np.isinf(bbp0)) | (np.isnan(bbp0))
+        bbp0 = xr.where(bbp0_fail, 0, bbp0)
+
+        # Compute bbp slope and extrapolate at all bands
+        gamma = self.gamma[0] * (1.0 - self.gamma[1] * np.power(Rrs442 / Rrs560, -self.gamma[2]))
+        bbp = bbp0 * np.power(band0 / self.bands, gamma)
+
+        # Compute total bb
+        bb = self.bbw + bbp
+
+        # Compute quasi-diffuse attenuation coefficient k at each band
+        # by 2nd order polynomial inversion
+        cA = Rrs
+        cB = - (coeffs.Gw0 * self.bbw + coeffs.Gp0 * bbp)
+        cC = - (coeffs.Gw1 * self.bbw *self.bbw + coeffs.Gp1 * bbp * bbp)
+        k = solve_2nd_order_poly(cA, cB, cC)
+
+        # If k is not a positive value, then
+        #   i) k --> bbw + aw
+        #   ii) k_fail flag is activated
+        k_fail = (k <= 0) | (np.isinf(k)) | (np.isnan(k))
+        k = xr.where(k_fail, self.aw + self.bbw, k)
+
+        # Drop unused coords to avoid issues
+        bb        = drop_unused_coords(bb)
+        k         = drop_unused_coords(k)
+        k_fail    = drop_unused_coords(k_fail)
+        bbp0_fail = drop_unused_coords(bbp0_fail)
+
+        # Set QAA_fail is either bbp0_fail or k_fail are activated
+        ds['QAA_fail'] = (bbp0_fail) | (k_fail)
+
+        # Compute final IOPs
+        ds['omega_b'] = bb / k
+        ds['eta_b'] = self.bbw / bb
+
+        return ds
+

ocbrdf/brdf_utils.py

@@ -0,0 +1,46 @@
+import numpy as np
+import xarray as xr
+import os
+
+# Define default auxiliary data file (OLCI OCP ADF)
+ref_path = os.path.dirname(os.path.realpath(__file__))
+# ADF_OCP = os.path.join(ref_path, '..', 'AuxiliaryData/OCP/S3A_OL_2_OCP_AX_20160216T000000_20991231T235959_20240327T100000___________________EUM_O_AL_008.SEN3/OL_2_OCP_AX.nc')
+ADF_OCP = os.path.join(ref_path, 'BRDF_LUTs','BRDF_%s.nc')
+
+def solve_2nd_order_poly(A, B, C):
+    """ Solve 2nd order polynomial inversion 
+    where coefficients are xr dataArray
+    Take only positive solution, otherwise provide 0.
+    """
+
+    # Compute solution according to sign of delta
+    delta = B*B - 4*A*C
+    mask = delta > 0
+    
+    # By default (and when delta < 0), take value at extremum
+    x = - B / (2*A)
+    
+    # When delta > 0, take biggest solutions
+    x_1 = (-B.where(mask) - np.sqrt(delta.where(mask))) / (2*A.where(mask))
+    x_2 = (-B.where(mask) + np.sqrt(delta.where(mask))) / (2*A.where(mask))
+    x_sol = xr.where(x_1 > x_2, x_1, x_2)
+    x = xr.where(delta > 0, x_sol, x)
+    
+    # Take only positive value
+    x = xr.where(x > 0, x, 0.)
+
+    return x
+
+def drop_unused_coords(var):
+    for coord in var.coords:
+        if coord not in var.dims:
+            var = var.drop(coord)
+    return var
+
+def squeeze_trivial_dims(ds,dimOrderAccordingTo='Rw'):
+    squeezedDims = {}
+    for d0,dim in enumerate(ds[dimOrderAccordingTo].dims):
+        if ds.sizes[dim] == 1:
+            squeezedDims[dim] = d0
+
+    return ds.squeeze(), squeezedDims

ocbrdf/ocbrdf_main.py

@@ -0,0 +1,577 @@
+import gc
+import os
+import numpy as np
+import sys
+import xarray as xr
+import rasterio
+import spectral
+
+from brdf_model_M02 import M02
+from brdf_model_M02SeaDAS import M02SeaDAS
+from brdf_model_L11 import L11
+from brdf_model_O25 import O25
+from brdf_utils import ADF_OCP, squeeze_trivial_dims
+
+
+UNC_LUT_CACHE = {}
+
+"""
+主 BRDF 校正模块
+输入为 xarray 数据集
+所需的光谱维度为 "bands"，其他维度自由
+输入数据集中必填的字段：
+Rw: 方向性海洋反射率
+sza: 太阳天顶角
+vza: 观测天顶角
+raa: 相对方位角（当太阳和观测在同一侧时，raa=0）
+输入数据集中可选的字段：
+Rw_unc: Rw 的不确定度（若缺失，设为零）
+输出数据集中的字段：
+nrrs: 完全归一化的遥感反射率
+rho_ex_w: nrrs * PI
+omega_b: bb/(a+bb)
+eta_b: bbw/bb
+C_brdf: BRDF 校正因子
+brdf_unc: C_brdf 的不确定度
+nrrs_unc: nrrs 的不确定度
+
+使用和适配 brdf_hypercp 模块需保留的信息：
+Brdf_hypercp 是 EUMETSAT 研究项目的一部分
+"BRDF correction of S3 OLCI water reflectance products"（S3 OLCI 水体反射率产品的 BRDF 校正），
+合同号：RB_EUM-CO-21-4600002626-JIG。
+研究团队成员：Davide D'Alimonte (davide.dalimonte@aequora.org)，
+Tamito Kajiyama (tamito.kajiyama@aequora.org)，
+Jaime Pitarch (jaime.pitarchportero@artov.ismar.cnr.it)，
+Vittorio Brando (vittorio.brando@cnr.it)，
+Marco Talone (talone@icm.csic.es) 以及
+Constant Mazeran (constant.mazeran@solvo.fr)。
+BRDF 查找表中的相对方位角遵循 OLCI 约定。详见 https://www.eumetsat.int/media/50720，图 6。
+"""
+
+
+# ---------------------------------------------------------------------------
+# File I/O helpers
+# ---------------------------------------------------------------------------
+
+def _resolve_envi_header(filepath):
+    """Return the ENVI header path accepted by spectral."""
+    if filepath.lower().endswith(".hdr"):
+        return filepath
+
+    hdr_path = filepath + ".hdr"
+    if os.path.exists(hdr_path):
+        return hdr_path
+
+    base, _ext = os.path.splitext(filepath)
+    hdr_path = base + ".hdr"
+    if os.path.exists(hdr_path):
+        return hdr_path
+
+    return filepath
+
+
+def read_bsq(filepath, wavelengths=None, dtype=np.float32):
+    """Read an ENVI BSQ hyperspectral image.
+
+    Parameters
+    ----------
+    filepath : str
+        Path to the ENVI header (.hdr) or image file.
+    wavelengths : array-like, optional
+        Band centre wavelengths in nm.  When *None*, wavelengths are read
+        from the ENVI header if available.
+
+    Returns
+    -------
+    data : ndarray, shape (bands, rows, cols), float64
+    wavelengths : ndarray, shape (bands,)
+    img : spectral image object (metadata / transform access)
+    """
+    img = spectral.open_image(_resolve_envi_header(filepath))
+    data = np.asarray(img.open_memmap(interleave="source"), dtype=dtype)
+    # data = np.transpose(data, (2, 0, 1)) # -> (bands, rows, cols)
+
+    if wavelengths is None:
+        if hasattr(img, 'bands') and img.bands.centers is not None:
+            wavelengths = np.array(img.bands.centers, dtype=np.float32)
+        else:
+            n_bands = data.shape[0]
+            wavelengths = np.arange(1, n_bands + 1, dtype=np.float32)
+
+    return data, np.asarray(wavelengths, dtype=dtype), img
+
+
+def read_bip_angles(filepath, dtype=np.float32):
+    """Read an ENVI BIP angle file and extract the five geometry bands.
+
+    Band layout (1-indexed):
+        band  7  -> SZA  Solar Zenith Angle
+        band  8  -> SAA  Solar Azimuth Angle
+        band  9  -> VZA  Sensor Zenith Angle
+        band 10  -> VAA  Sensor Azimuth Angle
+        band 11  -> RAA  Relative Azimuth Angle
+
+    Parameters
+    ----------
+    filepath : str
+        Path to the ENVI BIP file (header or image file).
+
+    Returns
+    -------
+    dict with keys 'sza', 'saa', 'vza', 'vaa', 'raa'.
+    Each value is a 2-D ndarray of shape (rows, cols), float64.
+    """
+    img = spectral.open_image(_resolve_envi_header(filepath))
+    data = np.asarray(img.open_memmap(interleave="source"), dtype=dtype)
+
+    n_bands = data.shape[2]
+    if n_bands < 11:
+        raise ValueError(
+            f"Angle file has only {n_bands} bands; at least 11 are required "
+            "(bands 7–11 = SZA, SAA, VZA, VAA, RAA)."
+        )
+
+    return {
+        'sza': data[:, :, 6].copy(),
+        'saa': data[:, :, 7].copy(),
+        'vza': data[:, :, 8].copy(),
+        'vaa': data[:, :, 9].copy(),
+        'raa': data[:, :, 10].copy(),
+    }
+
+
+def read_water_mask(filepath):
+    """Read a single-band GeoTIFF water mask.
+
+    Parameters
+    ----------
+    filepath : str
+        Path to the GeoTIFF file.
+
+    Returns
+    -------
+    mask : ndarray, shape (rows, cols), int8
+        Pixels with value 1 indicate water; all other values are non-water.
+    """
+    with rasterio.open(filepath) as src:
+        mask = src.read(1).astype(np.int8)
+    return mask
+
+
+def save_brdf_result(
+    ds_out,
+    output_file,
+    source_file=None,
+    output_var="Rw_brdf",
+    output_format="ENVI",
+    dtype=np.float32,
+):
+    """Save the BRDF result to ENVI or NetCDF."""
+    output_format = output_format.upper()
+
+    if output_var not in ds_out:
+        raise KeyError(f"Output variable '{output_var}' not found in result dataset.")
+
+    if output_format == "NETCDF":
+        if not output_file.lower().endswith(".nc"):
+            output_file = output_file + ".nc"
+        ds_out.to_netcdf(output_file)
+        return output_file
+
+    if output_format != "ENVI":
+        raise ValueError("output_format must be 'ENVI' or 'NETCDF'.")
+
+    cube = np.asarray(ds_out[output_var].values, dtype=dtype)
+    cube = np.transpose(cube, (1, 2, 0))
+
+    if output_file.lower().endswith(".hdr"):
+        hdr_path = output_file
+    else:
+        hdr_path = output_file + ".hdr"
+
+    metadata = {
+        "interleave": "bsq",
+        "description": f"BRDF corrected result: {output_var}",
+        "bands": cube.shape[2],
+        "lines": cube.shape[0],
+        "samples": cube.shape[1],
+    }
+    if "bands" in ds_out.coords:
+        metadata["wavelength"] = [str(v) for v in ds_out.coords["bands"].values.tolist()]
+
+    if source_file is not None:
+        metadata["source_file"] = source_file
+
+    spectral.envi.save_image(
+        hdr_path,
+        cube,
+        dtype=dtype,
+        interleave="bsq",
+        metadata=metadata,
+        force=True,
+    )
+    return hdr_path
+
+
+# ---------------------------------------------------------------------------
+# Top-level pipeline
+# ---------------------------------------------------------------------------
+
+def run_brdf_correction(
+    hyperspectral_file,
+    angle_file,
+    mask_file,
+    wavelengths=None,
+    brdf_model='L11',
+    chunk_size=4096,
+):
+    """BRDF correction pipeline for water pixels in a hyperspectral scene.
+
+    Reads the three input files, identifies water pixels via the mask,
+    runs ``brdf_prototype`` only on those pixels, and returns a full-scene
+    xarray Dataset where water pixels carry corrected values and non-water
+    pixels retain the original Rw (BRDF-derived fields are NaN outside
+    the water mask).
+
+    Parameters
+    ----------
+    hyperspectral_file : str
+        Path to the ENVI BSQ hyperspectral image (.hdr sidecar required).
+    angle_file : str
+        Path to the ENVI BIP angle file.  Bands 7–11 must be
+        SZA, SAA, VZA, VAA, RAA respectively.
+    mask_file : str
+        Path to the GeoTIFF water mask.  Pixels == 1 are treated as water.
+    wavelengths : array-like, optional
+        Band centre wavelengths in nm.  Overrides the header value when given.
+    brdf_model : str
+        BRDF model identifier: ``'L11'`` (default), ``'M02'``,
+        ``'M02SeaDAS'``, or ``'O25'``.
+    chunk_size : int, optional
+        Number of water pixels processed per batch (default 4096).
+        Smaller values reduce peak memory at the cost of more loop
+        iterations.
+
+    Returns
+    -------
+    xr.Dataset
+        Full-scene dataset with dimensions ``(bands, y, x)``.
+        Variables:
+
+        * ``Rw``          – original directional water reflectance
+        * ``nrrs``        – normalised remote-sensing reflectance (water only)
+        * ``rho_ex_w``    – nrrs × π (water only)
+        * ``C_brdf``      – BRDF correction factor (water only)
+        * ``brdf_unc``    – uncertainty of C_brdf (water only)
+        * ``nrrs_unc``    – uncertainty of nrrs (water only)
+        * ``C_brdf_fail`` – flag: True where correction failed (water only)
+        * ``water_mask``  – the input mask (1 = water)
+        * ``sza``, ``saa``, ``vza``, ``vaa``, ``raa`` – viewing geometry
+    """
+    # ------------------------------------------------------------------
+    # 1. Read all inputs
+    # ------------------------------------------------------------------
+    rw_data, wvl, _img_meta = read_bsq(hyperspectral_file, wavelengths, dtype=np.float32)
+    angles = read_bip_angles(angle_file, dtype=np.float32)
+    water_mask = read_water_mask(mask_file)
+
+    n_bands, n_rows, n_cols = rw_data.shape
+
+    # Sanity-check spatial dimensions
+    if angles['sza'].shape != (n_rows, n_cols):
+        raise ValueError(
+            f"Angle file spatial shape {angles['sza'].shape} does not match "
+            f"hyperspectral image shape ({n_rows}, {n_cols})."
+        )
+    if water_mask.shape != (n_rows, n_cols):
+        raise ValueError(
+            f"Mask shape {water_mask.shape} does not match "
+            f"hyperspectral image shape ({n_rows}, {n_cols})."
+        )
+
+    # ------------------------------------------------------------------
+    # 2. Identify water pixels
+    # ------------------------------------------------------------------
+    water_rows, water_cols = np.where(water_mask == 1)
+    n_water = water_rows.size
+    if n_water == 0:
+        raise RuntimeError("No water pixels found in the mask (no pixels == 1).")
+
+    # ------------------------------------------------------------------
+    # 3. Pre-allocate full-scene output arrays
+    #    Water pixels will be filled in chunks; non-water pixels stay NaN.
+    # ------------------------------------------------------------------
+    shape3 = (n_bands, n_rows, n_cols)
+    shape2 = (n_rows, n_cols)
+
+    nrrs_full        = np.full(shape3, np.nan, dtype=np.float32)
+    rho_ex_w_full    = np.full(shape3, np.nan, dtype=np.float32)
+    C_brdf_full      = np.full(shape3, np.nan, dtype=np.float32)
+    brdf_unc_full    = np.full(shape3, np.nan, dtype=np.float32)
+    nrrs_unc_full    = np.full(shape3, np.nan, dtype=np.float32)
+    C_brdf_fail_full = np.zeros(shape2, dtype=bool)
+    rw_brdf_full     = rw_data.copy()
+
+    # ------------------------------------------------------------------
+    # 4–6. Chunked BRDF correction
+    #      Process water pixels in batches of `chunk_size` to bound
+    #      peak memory.  Each batch builds a minimal xarray Dataset,
+    #      calls brdf_prototype, writes results back to the full-scene
+    #      arrays, then releases the batch objects.
+    # ------------------------------------------------------------------
+    n_chunks = int(np.ceil(n_water / chunk_size))
+    for chunk_idx in range(n_chunks):
+        lo = chunk_idx * chunk_size
+        hi = min(lo + chunk_size, n_water)
+
+        rows_c = water_rows[lo:hi]
+        cols_c = water_cols[lo:hi]
+
+        # Extract batch vectors
+        rw_c  = rw_data[:, rows_c, cols_c].T          # (chunk, bands)
+        sza_c = angles['sza'][rows_c, cols_c]          # (chunk,)
+        vza_c = angles['vza'][rows_c, cols_c]
+        raa_c = angles['raa'][rows_c, cols_c]
+
+        pixel_idx_c = np.arange(hi - lo)
+        ds_chunk = xr.Dataset(
+            {
+                'Rw':  xr.DataArray(rw_c,  dims=('n', 'bands'),
+                                    coords={'n': pixel_idx_c, 'bands': wvl}),
+                'sza': xr.DataArray(sza_c, dims='n',
+                                    coords={'n': pixel_idx_c}),
+                'vza': xr.DataArray(vza_c, dims='n',
+                                    coords={'n': pixel_idx_c}),
+                'raa': xr.DataArray(raa_c, dims='n',
+                                    coords={'n': pixel_idx_c}),
+            }
+        )
+
+        ds_corr = brdf_prototype(ds_chunk, brdf_model=brdf_model)
+
+        # Write results back; ds_corr arrays have dims (n, bands) → transpose
+        nrrs_full[:, rows_c, cols_c]        = ds_corr['nrrs'].values.T
+        rho_ex_w_full[:, rows_c, cols_c]    = ds_corr['rho_ex_w'].values.T
+        C_brdf_full[:, rows_c, cols_c]      = ds_corr['C_brdf'].values.T
+        brdf_unc_full[:, rows_c, cols_c]    = ds_corr['brdf_unc'].values.T
+        nrrs_unc_full[:, rows_c, cols_c]    = ds_corr['nrrs_unc'].values.T
+        C_brdf_fail_full[rows_c, cols_c]    = ds_corr['C_brdf_fail'].values
+        rw_brdf_full[:, rows_c, cols_c]     = ds_corr['rho_ex_w'].values.T
+
+        # Release batch objects so GC can reclaim memory
+        del rw_c, sza_c, vza_c, raa_c, pixel_idx_c, ds_chunk, ds_corr
+        gc.collect()
+
+    # ------------------------------------------------------------------
+    # 7. Assemble full-scene output Dataset
+    # ------------------------------------------------------------------
+    full_coords = {'bands': wvl, 'y': np.arange(n_rows), 'x': np.arange(n_cols)}
+    dims3 = ('bands', 'y', 'x')
+    dims2 = ('y', 'x')
+
+    ds_out = xr.Dataset(
+        {
+            'Rw':          xr.DataArray(rw_data,          dims=dims3, coords=full_coords),
+            'Rw_brdf':     xr.DataArray(rw_brdf_full,     dims=dims3, coords=full_coords),
+            'nrrs':        xr.DataArray(nrrs_full,        dims=dims3, coords=full_coords),
+            'rho_ex_w':    xr.DataArray(rho_ex_w_full,    dims=dims3, coords=full_coords),
+            'C_brdf':      xr.DataArray(C_brdf_full,      dims=dims3, coords=full_coords),
+            'brdf_unc':    xr.DataArray(brdf_unc_full,    dims=dims3, coords=full_coords),
+            'nrrs_unc':    xr.DataArray(nrrs_unc_full,    dims=dims3, coords=full_coords),
+            'C_brdf_fail': xr.DataArray(C_brdf_fail_full, dims=dims2),
+            'water_mask':  xr.DataArray(water_mask,       dims=dims2),
+            'sza':         xr.DataArray(angles['sza'],    dims=dims2),
+            'saa':         xr.DataArray(angles['saa'],    dims=dims2),
+            'vza':         xr.DataArray(angles['vza'],    dims=dims2),
+            'vaa':         xr.DataArray(angles['vaa'],    dims=dims2),
+            'raa':         xr.DataArray(angles['raa'],    dims=dims2),
+        }
+    )
+
+    return ds_out
+
+
+def process_brdf_files(
+    hyperspectral_file,
+    angle_file,
+    mask_file,
+    output_file,
+    wavelengths=None,
+    brdf_model="L11",
+    output_var="Rw_brdf",
+    output_format="ENVI",
+    chunk_size=4096,
+):
+    """Run BRDF correction from three input files and save the result."""
+    ds_out = run_brdf_correction(
+        hyperspectral_file=hyperspectral_file,
+        angle_file=angle_file,
+        mask_file=mask_file,
+        wavelengths=wavelengths,
+        brdf_model=brdf_model,
+        chunk_size=chunk_size,
+    )
+    saved_path = save_brdf_result(
+        ds_out=ds_out,
+        output_file=output_file,
+        source_file=hyperspectral_file,
+        output_var=output_var,
+        output_format=output_format,
+    )
+    return ds_out, saved_path
+
+
+# ---------------------------------------------------------------------------
+# Core BRDF correction (unchanged from original)
+# ---------------------------------------------------------------------------
+
+def brdf_prototype(ds, adf=None, brdf_model='L11'):
+    # 测试 BRDF 模型在 GUI 中不支持：强制覆盖
+    # brdf_model = 'M02SeaDAS'
+
+    # 压缩单一维度（例如，投放次数、提取次数等）以避免插值问题
+    ds, squeezedDims = squeeze_trivial_dims(ds)
+
+    # Initialise model
+    if brdf_model == 'M02':
+        BRDF_model = M02(bands=ds.bands, aot=ds.aot, wind=ds.wind, adf=None)  # Don't use brdf_py.ADF context
+    elif brdf_model == 'M02SeaDAS':
+        BRDF_model = M02SeaDAS(bands=ds.bands, adf=None)  # Don't use brdf_py.ADF context
+    elif brdf_model == 'L11':
+        BRDF_model = L11(bands=ds.bands, adf=None)  # Don't use brdf_py.ADF context
+    elif brdf_model == 'O25':
+        BRDF_model = O25(bands=ds.bands, adf=None)  # Don't use brdf_py.ADF context
+    else:
+        print("BRDF model %s not supported" % brdf_model)
+        sys.exit(1)
+
+    # Init pixel
+    BRDF_model.init_pixels(ds['sza'], ds['vza'], ds['raa'])
+
+    # Compute IOP and normalize by iterating
+    ds['nrrs'] = ds['Rw'] / np.pi
+
+    ds['convergeFlag'] = (0 * ds['sza']).astype(bool)
+    ds['C_brdf'] = 0 * ds['nrrs'] + 1
+
+    for iter_brdf in range(int(BRDF_model.niter)):
+
+        # M02: Initialise chl_iter
+        if brdf_model in ['M02', 'M02SeaDAS'] and (iter_brdf == 0):
+            chl_iter = {}
+            ds['log10_chl'] = 0 * ds['sza'] + float(np.log10(BRDF_model.OC4MEchl0))
+            chl_iter[-1]    = 0 * ds['sza'] + float(BRDF_model.OC4MEchl0)
+
+        ds = BRDF_model.backward(ds, iter_brdf)
+
+        # M02: Check convergence (dummy for M02SeaDAS for the moment... epsilon set to 0)
+        if brdf_model in ['M02', 'M02SeaDAS']:
+            chl_iter[iter_brdf] = 10 ** ds['log10_chl']
+            #  Check if convergence is reached |chl_old-chl_new| < epsilon * chl_new
+            ds['convergeFlag'] = (ds['convergeFlag']) | (
+                (np.abs(chl_iter[iter_brdf - 1] - chl_iter[iter_brdf]) < float(BRDF_model.OC4MEepsilon) * chl_iter[
+                    iter_brdf]))
+
+        # Apply forward model in both geometries
+        # forward_mod = BRDF_model.forward(ds).transpose('n', 'bands')
+        # forward_mod0 = BRDF_model.forward(ds, normalized=True).transpose('n', 'bands')
+        forward_mod = BRDF_model.forward(ds)
+        forward_mod0 = BRDF_model.forward(ds, normalized=True)
+
+        ratio = forward_mod0 / forward_mod
+
+        # Drop remnant coordinates to avoid ambiguities in the update of the BRDF factor.
+        for coord in ratio.coords:
+            if coord not in ds['C_brdf'].coords:
+                ratio = ratio.drop(coord)
+
+        # Normalize reflectance
+        ds['C_brdf'] = xr.where(ds['convergeFlag'], ds['C_brdf'], ratio)
+        ds['nrrs'] = ds['Rw'] / np.pi * ds['C_brdf']
+
+    # Flag BRDF where NaN and set to 1 (no correction applied).
+    ds['C_brdf_fail'] = np.isnan(ds['C_brdf'])
+    ds['C_brdf'] = xr.where(ds['C_brdf_fail'], 1, ds['C_brdf'])
+    ds['nrrs'] = xr.where(ds['C_brdf_fail'], ds['Rw'] / np.pi, ds['nrrs'])
+
+    # If QAA_fail is raised, raise C_brdf_fail (but still apply C_brdf).
+    if 'QAA_fail' in ds:
+        ds['C_brdf_fail'] = (ds['C_brdf_fail']) | (ds['QAA_fail'])
+
+    # Compute uncertainty
+    ds = brdf_uncertainty(ds)
+
+    # Compute flag
+    ds['flags_level2'] = ds['Rw'] * 0  # TODO
+
+    # Convert to reflectance unit
+    ds['rho_ex_w'] = ds['nrrs'] * np.pi
+
+    # Expand squeezed trivial dimensions
+    for dim,d0 in squeezedDims.items():
+        ds = ds.expand_dims(dim,axis=d0)
+
+    return ds
+
+def brdf_uncertainty(ds, adf=None):
+    ''' Compute uncertainty of BRDF factor and propagate to nrrs '''
+
+    # Read LUT
+    if adf is None:
+        adf = ADF_OCP
+
+    unc_lut_path = adf % 'UNC'
+    if unc_lut_path not in UNC_LUT_CACHE:
+        UNC_LUT_CACHE[unc_lut_path] = xr.open_dataset(unc_lut_path, engine='netcdf4')
+    LUT = UNC_LUT_CACHE[unc_lut_path]
+
+    # Interpolate relative uncertainty
+    unc = LUT['unc'].interp(lambda_unc=ds.bands, theta_s_unc=ds.sza, theta_v_unc=ds.vza,
+                            delta_phi_unc=ds.raa)
+
+    # Compute absolute uncertainty of factor
+    ds['brdf_unc'] = unc * ds['C_brdf']
+
+    # Flag BRDF_unc where NaN and set to 0
+    ds['brdf_unc_fail'] = np.isnan(ds['brdf_unc'])
+    ds['brdf_unc'] = xr.where(ds['brdf_unc_fail'], 0, ds['brdf_unc'])
+
+    # Propagate to nrrs
+    Rwex_unc2 = ds['brdf_unc'] * ds['brdf_unc'] * ds['Rw'] * ds['Rw']
+    if 'Rw_unc' in ds:
+        Rwex_unc2 += ds['C_brdf'] * ds['C_brdf'] * ds['Rw_unc'] * ds['Rw_unc']
+    ds['nrrs_unc'] = np.sqrt(Rwex_unc2)/np.pi
+
+    return ds
+
+
+def main():
+    import argparse
+
+    parser = argparse.ArgumentParser(description="Water-region BRDF correction for hyperspectral data.")
+    parser.add_argument("hyperspectral_file", help="Input ENVI BSQ hyperspectral file or .hdr path.")
+    parser.add_argument("angle_file", help="Input ENVI BIP angle file or .hdr path.")
+    parser.add_argument("mask_file", help="Input water mask GeoTIFF; pixels == 1 are corrected.")
+    parser.add_argument("output_file", help="Output path prefix. ENVI writes .hdr/.img, NetCDF writes .nc.")
+    parser.add_argument("--brdf-model", default="L11", choices=["L11", "M02", "M02SeaDAS", "O25"])
+    parser.add_argument("--output-var", default="Rw_brdf", choices=["Rw_brdf", "rho_ex_w", "nrrs", "C_brdf"])
+    parser.add_argument("--output-format", default="ENVI", choices=["ENVI", "NETCDF"])
+    parser.add_argument("--chunk-size", type=int, default=4096,
+                        help="Number of water pixels per processing batch (default: 4096).")
+    args = parser.parse_args()
+
+    _ds_out, saved_path = process_brdf_files(
+        hyperspectral_file=args.hyperspectral_file,
+        angle_file=args.angle_file,
+        mask_file=args.mask_file,
+        output_file=args.output_file,
+        brdf_model=args.brdf_model,
+        output_var=args.output_var,
+        output_format=args.output_format,
+        chunk_size=args.chunk_size,
+    )
+    print(f"BRDF correction finished. Saved to: {saved_path}")
+
+
+if __name__ == "__main__":
+    main()

ocbrdf/ocbrdf_main2.py

@@ -0,0 +1,586 @@
+import gc
+import os
+import numpy as np
+import sys
+import xarray as xr
+import rasterio
+import spectral
+
+from brdf_model_M02 import M02
+from brdf_model_M02SeaDAS import M02SeaDAS
+from brdf_model_L11 import L11
+from brdf_model_O25 import O25
+from brdf_utils import ADF_OCP, squeeze_trivial_dims
+
+
+UNC_LUT_CACHE = {}
+
+"""
+主 BRDF 校正模块
+输入为 xarray 数据集
+所需的光谱维度为 "bands"，其他维度自由
+输入数据集中必填的字段：
+Rw: 方向性海洋反射率
+sza: 太阳天顶角
+vza: 观测天顶角
+raa: 相对方位角（当太阳和观测在同一侧时，raa=0）
+输入数据集中可选的字段：
+Rw_unc: Rw 的不确定度（若缺失，设为零）
+输出数据集中的字段：
+nrrs: 完全归一化的遥感反射率
+rho_ex_w: nrrs * PI
+omega_b: bb/(a+bb)
+eta_b: bbw/bb
+C_brdf: BRDF 校正因子
+brdf_unc: C_brdf 的不确定度
+nrrs_unc: nrrs 的不确定度
+
+使用和适配 brdf_hypercp 模块需保留的信息：
+Brdf_hypercp 是 EUMETSAT 研究项目的一部分
+"BRDF correction of S3 OLCI water reflectance products"（S3 OLCI 水体反射率产品的 BRDF 校正），
+合同号：RB_EUM-CO-21-4600002626-JIG。
+研究团队成员：Davide D'Alimonte (davide.dalimonte@aequora.org)，
+Tamito Kajiyama (tamito.kajiyama@aequora.org)，
+Jaime Pitarch (jaime.pitarchportero@artov.ismar.cnr.it)，
+Vittorio Brando (vittorio.brando@cnr.it)，
+Marco Talone (talone@icm.csic.es) 以及
+Constant Mazeran (constant.mazeran@solvo.fr)。
+BRDF 查找表中的相对方位角遵循 OLCI 约定。详见 https://www.eumetsat.int/media/50720，图 6。
+"""
+
+
+# ---------------------------------------------------------------------------
+# File I/O helpers
+# ---------------------------------------------------------------------------
+
+def _resolve_envi_header(filepath):
+    """Return the ENVI header path accepted by spectral."""
+    if filepath.lower().endswith(".hdr"):
+        return filepath
+
+    hdr_path = filepath + ".hdr"
+    if os.path.exists(hdr_path):
+        return hdr_path
+
+    base, _ext = os.path.splitext(filepath)
+    hdr_path = base + ".hdr"
+    if os.path.exists(hdr_path):
+        return hdr_path
+
+    return filepath
+
+
+def read_bsq(filepath, wavelengths=None, dtype=np.float32):
+    """Read an ENVI BSQ hyperspectral image.
+
+    Parameters
+    ----------
+    filepath : str
+        Path to the ENVI header (.hdr) or image file.
+    wavelengths : array-like, optional
+        Band centre wavelengths in nm.  When *None*, wavelengths are read
+        from the ENVI header if available.
+
+    Returns
+    -------
+    data : ndarray, shape (bands, rows, cols), float64
+    wavelengths : ndarray, shape (bands,)
+    img : spectral image object (metadata / transform access)
+    """
+    img = spectral.open_image(_resolve_envi_header(filepath))
+    data = np.asarray(img.open_memmap(interleave="source"), dtype=dtype)
+    # data = np.transpose(data, (2, 0, 1)) # -> (bands, rows, cols)
+
+    if wavelengths is None:
+        if hasattr(img, 'bands') and img.bands.centers is not None:
+            wavelengths = np.array(img.bands.centers, dtype=np.float32)
+        else:
+            n_bands = data.shape[0]
+            wavelengths = np.arange(1, n_bands + 1, dtype=np.float32)
+
+    return data, np.asarray(wavelengths, dtype=dtype), img
+
+
+def read_bip_angles(filepath, dtype=np.float32):
+    """Read an ENVI BIP angle file and extract the five geometry bands.
+
+    Band layout (1-indexed):
+        band  7  -> SZA  Solar Zenith Angle
+        band  8  -> SAA  Solar Azimuth Angle
+        band  9  -> VZA  Sensor Zenith Angle
+        band 10  -> VAA  Sensor Azimuth Angle
+        band 11  -> RAA  Relative Azimuth Angle
+
+    Parameters
+    ----------
+    filepath : str
+        Path to the ENVI BIP file (header or image file).
+
+    Returns
+    -------
+    dict with keys 'sza', 'saa', 'vza', 'vaa', 'raa'.
+    Each value is a 2-D ndarray of shape (rows, cols), float64.
+    """
+    img = spectral.open_image(_resolve_envi_header(filepath))
+    data = np.asarray(img.open_memmap(interleave="source"), dtype=dtype)
+
+    n_bands = data.shape[2]
+    if n_bands < 11:
+        raise ValueError(
+            f"Angle file has only {n_bands} bands; at least 11 are required "
+            "(bands 7–11 = SZA, SAA, VZA, VAA, RAA)."
+        )
+
+    return {
+        'sza': data[:, :, 6].copy(),
+        'saa': data[:, :, 7].copy(),
+        'vza': data[:, :, 8].copy(),
+        'vaa': data[:, :, 9].copy(),
+        'raa': data[:, :, 10].copy(),
+    }
+
+
+def read_water_mask(filepath):
+    """Read a single-band GeoTIFF water mask.
+
+    Parameters
+    ----------
+    filepath : str
+        Path to the GeoTIFF file.
+
+    Returns
+    -------
+    mask : ndarray, shape (rows, cols), int8
+        Pixels with value 1 indicate water; all other values are non-water.
+    """
+    with rasterio.open(filepath) as src:
+        mask = src.read(1).astype(np.int8)
+    return mask
+
+
+def save_brdf_result(
+    ds_out,
+    output_file,
+    source_file=None,
+    output_var="Rw_brdf",
+    output_format="ENVI",
+    dtype=np.float32,
+):
+    """Save the BRDF result to ENVI or NetCDF."""
+    output_format = output_format.upper()
+
+    if output_var not in ds_out:
+        raise KeyError(f"Output variable '{output_var}' not found in result dataset.")
+
+    if output_format == "NETCDF":
+        if not output_file.lower().endswith(".nc"):
+            output_file = output_file + ".nc"
+        ds_out.to_netcdf(output_file)
+        return output_file
+
+    if output_format != "ENVI":
+        raise ValueError("output_format must be 'ENVI' or 'NETCDF'.")
+
+    cube = np.asarray(ds_out[output_var].values, dtype=dtype)
+    cube = np.transpose(cube, (1, 2, 0))
+
+    if output_file.lower().endswith(".hdr"):
+        hdr_path = output_file
+    else:
+        hdr_path = output_file + ".hdr"
+
+    metadata = {
+        "interleave": "bsq",
+        "description": f"BRDF corrected result: {output_var}",
+        "bands": cube.shape[2],
+        "lines": cube.shape[0],
+        "samples": cube.shape[1],
+    }
+    if "bands" in ds_out.coords:
+        metadata["wavelength"] = [str(v) for v in ds_out.coords["bands"].values.tolist()]
+
+    if source_file is not None:
+        metadata["source_file"] = source_file
+
+    spectral.envi.save_image(
+        hdr_path,
+        cube,
+        dtype=dtype,
+        interleave="bsq",
+        metadata=metadata,
+        force=True,
+    )
+    return hdr_path
+
+
+# ---------------------------------------------------------------------------
+# Top-level pipeline
+# ---------------------------------------------------------------------------
+def brdf_uncertainty(ds, adf=None):
+    ''' Compute uncertainty of BRDF factor and propagate to nrrs '''
+
+    # Read LUT
+    if adf is None:
+        adf = ADF_OCP
+
+    unc_lut_path = adf % 'UNC'
+    if unc_lut_path not in UNC_LUT_CACHE:
+        UNC_LUT_CACHE[unc_lut_path] = xr.open_dataset(unc_lut_path, engine='netcdf4')
+    LUT = UNC_LUT_CACHE[unc_lut_path]
+
+    # Interpolate relative uncertainty
+    unc = LUT['unc'].interp(lambda_unc=ds.bands, theta_s_unc=ds.sza, theta_v_unc=ds.vza,
+                            delta_phi_unc=ds.raa)
+
+    # Compute absolute uncertainty of factor
+    ds['brdf_unc'] = unc * ds['C_brdf']
+
+    # Flag BRDF_unc where NaN and set to 0
+    ds['brdf_unc_fail'] = np.isnan(ds['brdf_unc'])
+    ds['brdf_unc'] = xr.where(ds['brdf_unc_fail'], 0, ds['brdf_unc'])
+
+    # Propagate to nrrs
+    Rwex_unc2 = ds['brdf_unc'] * ds['brdf_unc'] * ds['Rw'] * ds['Rw']
+    if 'Rw_unc' in ds:
+        Rwex_unc2 += ds['C_brdf'] * ds['C_brdf'] * ds['Rw_unc'] * ds['Rw_unc']
+    ds['nrrs_unc'] = np.sqrt(Rwex_unc2)/np.pi
+
+    return ds
+def brdf_prototype(ds, adf=None, brdf_model='L11', compute_uncertainty=True):
+    # 测试 BRDF 模型在 GUI 中不支持：强制覆盖
+    # brdf_model = 'M02SeaDAS'
+
+    # 压缩单一维度（例如，投放次数、提取次数等）以避免插值问题
+    ds, squeezedDims = squeeze_trivial_dims(ds)
+
+    # Initialise model
+    if brdf_model == 'M02':
+        BRDF_model = M02(bands=ds.bands, aot=ds.aot, wind=ds.wind, adf=None)  # Don't use brdf_py.ADF context
+    elif brdf_model == 'M02SeaDAS':
+        BRDF_model = M02SeaDAS(bands=ds.bands, adf=None)  # Don't use brdf_py.ADF context
+    elif brdf_model == 'L11':
+        BRDF_model = L11(bands=ds.bands, adf=None)  # Don't use brdf_py.ADF context
+    elif brdf_model == 'O25':
+        BRDF_model = O25(bands=ds.bands, adf=None)  # Don't use brdf_py.ADF context
+    else:
+        print("BRDF model %s not supported" % brdf_model)
+        sys.exit(1)
+
+    # Init pixel
+    BRDF_model.init_pixels(ds['sza'], ds['vza'], ds['raa'])
+
+    # Compute IOP and normalize by iterating
+    ds['nrrs'] = ds['Rw'] / np.pi
+
+    ds['convergeFlag'] = (0 * ds['sza']).astype(bool)
+    ds['C_brdf'] = 0 * ds['nrrs'] + 1
+
+    for iter_brdf in range(int(BRDF_model.niter)):
+
+        # M02: Initialise chl_iter
+        if brdf_model in ['M02', 'M02SeaDAS'] and (iter_brdf == 0):
+            chl_iter = {}
+            ds['log10_chl'] = 0 * ds['sza'] + float(np.log10(BRDF_model.OC4MEchl0))
+            chl_iter[-1]    = 0 * ds['sza'] + float(BRDF_model.OC4MEchl0)
+
+        ds = BRDF_model.backward(ds, iter_brdf)
+
+        # M02: Check convergence (dummy for M02SeaDAS for the moment... epsilon set to 0)
+        if brdf_model in ['M02', 'M02SeaDAS']:
+            chl_iter[iter_brdf] = 10 ** ds['log10_chl']
+            #  Check if convergence is reached |chl_old-chl_new| < epsilon * chl_new
+            ds['convergeFlag'] = (ds['convergeFlag']) | (
+                (np.abs(chl_iter[iter_brdf - 1] - chl_iter[iter_brdf]) < float(BRDF_model.OC4MEepsilon) * chl_iter[
+                    iter_brdf]))
+
+        # Apply forward model in both geometries
+        forward_mod = BRDF_model.forward(ds)
+        forward_mod0 = BRDF_model.forward(ds, normalized=True)
+
+        ratio = forward_mod0 / forward_mod
+
+        # Drop remnant coordinates to avoid ambiguities in the update of the BRDF factor.
+        for coord in ratio.coords:
+            if coord not in ds['C_brdf'].coords:
+                ratio = ratio.drop(coord)
+
+        # Normalize reflectance
+        ds['C_brdf'] = xr.where(ds['convergeFlag'], ds['C_brdf'], ratio)
+        ds['nrrs'] = ds['Rw'] / np.pi * ds['C_brdf']
+
+    # Flag BRDF where NaN and set to 1 (no correction applied).
+    ds['C_brdf_fail'] = np.isnan(ds['C_brdf'])
+    ds['C_brdf'] = xr.where(ds['C_brdf_fail'], 1, ds['C_brdf'])
+    ds['nrrs'] = xr.where(ds['C_brdf_fail'], ds['Rw'] / np.pi, ds['nrrs'])
+
+    # If QAA_fail is raised, raise C_brdf_fail (but still apply C_brdf).
+    if 'QAA_fail' in ds:
+        ds['C_brdf_fail'] = (ds['C_brdf_fail']) | (ds['QAA_fail'])
+
+    # Compute uncertainty only if requested
+    if compute_uncertainty:
+        ds = brdf_uncertainty(ds)
+
+    # Compute flag
+    ds['flags_level2'] = ds['Rw'] * 0  # TODO
+
+    # Convert to reflectance unit
+    ds['rho_ex_w'] = ds['nrrs'] * np.pi
+
+    # Expand squeezed trivial dimensions
+    for dim,d0 in squeezedDims.items():
+        ds = ds.expand_dims(dim,axis=d0)
+
+    return ds
+
+
+def run_brdf_correction(
+    hyperspectral_file,
+    angle_file,
+    mask_file,
+    wavelengths=None,
+    brdf_model='L11',
+    chunk_size=4096,
+    output_var=None,          # if None, keep full output (original behaviour)
+    compute_uncertainty=True, # compute uncertainty only if requested
+):
+    """BRDF correction pipeline for water pixels in a hyperspectral scene.
+
+    ... (docstring unchanged, update description) ...
+    """
+    # ------------------------------------------------------------------
+    # 1. Read all inputs
+    # ------------------------------------------------------------------
+    rw_data, wvl, _img_meta = read_bsq(hyperspectral_file, wavelengths, dtype=np.float32)
+    angles = read_bip_angles(angle_file, dtype=np.float32)
+    water_mask = read_water_mask(mask_file)
+
+    n_bands, n_rows, n_cols = rw_data.shape
+
+    # Sanity-check spatial dimensions
+    if angles['sza'].shape != (n_rows, n_cols):
+        raise ValueError(
+            f"Angle file spatial shape {angles['sza'].shape} does not match "
+            f"hyperspectral image shape ({n_rows}, {n_cols})."
+        )
+    if water_mask.shape != (n_rows, n_cols):
+        raise ValueError(
+            f"Mask shape {water_mask.shape} does not match "
+            f"hyperspectral image shape ({n_rows}, {n_cols})."
+        )
+
+    # ------------------------------------------------------------------
+    # 2. Identify water pixels
+    # ------------------------------------------------------------------
+    water_rows, water_cols = np.where(water_mask == 1)
+    n_water = water_rows.size
+    if n_water == 0:
+        raise RuntimeError("No water pixels found in the mask (no pixels == 1).")
+
+    # ------------------------------------------------------------------
+    # 3. Determine which output variable(s) we need to keep
+    # ------------------------------------------------------------------
+    # If output_var is None, we keep everything (original behaviour)
+    keep_all = (output_var is None)
+
+    # Map output_var to the field in ds_corr
+    var_map = {
+        'Rw_brdf': 'rho_ex_w',
+        'rho_ex_w': 'rho_ex_w',
+        'nrrs': 'nrrs',
+        'C_brdf': 'C_brdf'
+    }
+    if not keep_all:
+        if output_var not in var_map:
+            raise ValueError(f"Unknown output_var '{output_var}'. Allowed: {list(var_map.keys())}")
+        field_name = var_map[output_var]
+    else:
+        field_name = None
+
+    # ------------------------------------------------------------------
+    # 4. Allocate output arrays (only those needed)
+    # ------------------------------------------------------------------
+    shape3 = (n_bands, n_rows, n_cols)
+    shape2 = (n_rows, n_cols)
+
+    if keep_all:
+        # Original behaviour: allocate all arrays
+        nrrs_full        = np.full(shape3, np.nan, dtype=np.float32)
+        rho_ex_w_full    = np.full(shape3, np.nan, dtype=np.float32)
+        C_brdf_full      = np.full(shape3, np.nan, dtype=np.float32)
+        brdf_unc_full    = np.full(shape3, np.nan, dtype=np.float32)
+        nrrs_unc_full    = np.full(shape3, np.nan, dtype=np.float32)
+        C_brdf_fail_full = np.zeros(shape2, dtype=bool)
+        rw_brdf_full     = rw_data.copy()
+    else:
+        # Only allocate the requested output array
+        output_full = np.full(shape3, np.nan, dtype=np.float32)
+        # We may also need C_brdf_fail if requested? Not for output_var, but could be used internally? Skip.
+        # For simplicity, we do not allocate extra arrays.
+
+    # ------------------------------------------------------------------
+    # 5–6. Chunked BRDF correction
+    # ------------------------------------------------------------------
+    n_chunks = int(np.ceil(n_water / chunk_size))
+    for chunk_idx in range(n_chunks):
+        lo = chunk_idx * chunk_size
+        hi = min(lo + chunk_size, n_water)
+
+        rows_c = water_rows[lo:hi]
+        cols_c = water_cols[lo:hi]
+
+        # Extract batch vectors
+        rw_c  = rw_data[:, rows_c, cols_c].T          # (chunk, bands)
+        sza_c = angles['sza'][rows_c, cols_c]          # (chunk,)
+        vza_c = angles['vza'][rows_c, cols_c]
+        raa_c = angles['raa'][rows_c, cols_c]
+
+        pixel_idx_c = np.arange(hi - lo)
+        ds_chunk = xr.Dataset(
+            {
+                'Rw':  xr.DataArray(rw_c,  dims=('n', 'bands'),
+                                    coords={'n': pixel_idx_c, 'bands': wvl}),
+                'sza': xr.DataArray(sza_c, dims='n',
+                                    coords={'n': pixel_idx_c}),
+                'vza': xr.DataArray(vza_c, dims='n',
+                                    coords={'n': pixel_idx_c}),
+                'raa': xr.DataArray(raa_c, dims='n',
+                                    coords={'n': pixel_idx_c}),
+            }
+        )
+
+        # Run BRDF model with uncertainty control
+        ds_corr = brdf_prototype(ds_chunk, brdf_model=brdf_model, compute_uncertainty=compute_uncertainty)
+
+        if keep_all:
+            # Write all results back
+            nrrs_full[:, rows_c, cols_c]        = ds_corr['nrrs'].values.T
+            rho_ex_w_full[:, rows_c, cols_c]    = ds_corr['rho_ex_w'].values.T
+            C_brdf_full[:, rows_c, cols_c]      = ds_corr['C_brdf'].values.T
+            if compute_uncertainty:
+                brdf_unc_full[:, rows_c, cols_c]    = ds_corr['brdf_unc'].values.T
+                nrrs_unc_full[:, rows_c, cols_c]    = ds_corr['nrrs_unc'].values.T
+            C_brdf_fail_full[rows_c, cols_c]    = ds_corr['C_brdf_fail'].values
+            rw_brdf_full[:, rows_c, cols_c]     = ds_corr['rho_ex_w'].values.T
+        else:
+            # Only write the requested output field
+            values = ds_corr[field_name].values.T
+            output_full[:, rows_c, cols_c] = values
+
+        # Release batch objects so GC can reclaim memory
+        del rw_c, sza_c, vza_c, raa_c, pixel_idx_c, ds_chunk, ds_corr
+        gc.collect()
+
+    # ------------------------------------------------------------------
+    # 7. Assemble full-scene output Dataset
+    # ------------------------------------------------------------------
+    full_coords = {'bands': wvl, 'y': np.arange(n_rows), 'x': np.arange(n_cols)}
+    dims3 = ('bands', 'y', 'x')
+    dims2 = ('y', 'x')
+
+    if keep_all:
+        # Original behaviour: return all variables
+        ds_out = xr.Dataset(
+            {
+                'Rw':          xr.DataArray(rw_data,          dims=dims3, coords=full_coords),
+                'Rw_brdf':     xr.DataArray(rw_brdf_full,     dims=dims3, coords=full_coords),
+                'nrrs':        xr.DataArray(nrrs_full,        dims=dims3, coords=full_coords),
+                'rho_ex_w':    xr.DataArray(rho_ex_w_full,    dims=dims3, coords=full_coords),
+                'C_brdf':      xr.DataArray(C_brdf_full,      dims=dims3, coords=full_coords),
+                'C_brdf_fail': xr.DataArray(C_brdf_fail_full, dims=dims2),
+                'water_mask':  xr.DataArray(water_mask,       dims=dims2),
+                'sza':         xr.DataArray(angles['sza'],    dims=dims2),
+                'saa':         xr.DataArray(angles['saa'],    dims=dims2),
+                'vza':         xr.DataArray(angles['vza'],    dims=dims2),
+                'vaa':         xr.DataArray(angles['vaa'],    dims=dims2),
+                'raa':         xr.DataArray(angles['raa'],    dims=dims2),
+            }
+        )
+        if compute_uncertainty:
+            ds_out['brdf_unc'] = xr.DataArray(brdf_unc_full, dims=dims3, coords=full_coords)
+            ds_out['nrrs_unc'] = xr.DataArray(nrrs_unc_full, dims=dims3, coords=full_coords)
+    else:
+        # Return only the requested output variable + geometry + mask
+        ds_out = xr.Dataset(
+            {
+                output_var:    xr.DataArray(output_full,      dims=dims3, coords=full_coords),
+                'Rw':          xr.DataArray(rw_data,          dims=dims3, coords=full_coords),
+                'water_mask':  xr.DataArray(water_mask,       dims=dims2),
+                'sza':         xr.DataArray(angles['sza'],    dims=dims2),
+                'saa':         xr.DataArray(angles['saa'],    dims=dims2),
+                'vza':         xr.DataArray(angles['vza'],    dims=dims2),
+                'vaa':         xr.DataArray(angles['vaa'],    dims=dims2),
+                'raa':         xr.DataArray(angles['raa'],    dims=dims2),
+            }
+        )
+        # If uncertainties were computed, we could add them optionally, but user didn't request them.
+        # For memory efficiency, we skip them here.
+
+    return ds_out
+
+
+def process_brdf_files(
+    hyperspectral_file,
+    angle_file,
+    mask_file,
+    output_file,
+    wavelengths=None,
+    brdf_model="L11",
+    output_var="Rw_brdf",
+    output_format="ENVI",
+    chunk_size=4096,
+    compute_uncertainty=False,   # new parameter, default False to save time
+):
+    """Run BRDF correction from three input files and save the result."""
+    ds_out = run_brdf_correction(
+        hyperspectral_file=hyperspectral_file,
+        angle_file=angle_file,
+        mask_file=mask_file,
+        wavelengths=wavelengths,
+        brdf_model=brdf_model,
+        chunk_size=chunk_size,
+        output_var=output_var,                # pass output_var to limit output
+        compute_uncertainty=compute_uncertainty,
+    )
+    saved_path = save_brdf_result(
+        ds_out=ds_out,
+        output_file=output_file,
+        source_file=hyperspectral_file,
+        output_var=output_var,
+        output_format=output_format,
+    )
+    return ds_out, saved_path
+
+
+def main():
+    import argparse
+
+    parser = argparse.ArgumentParser(description="Water-region BRDF correction for hyperspectral data.")
+    parser.add_argument("hyperspectral_file", help="Input ENVI BSQ hyperspectral file or .hdr path.")
+    parser.add_argument("angle_file", help="Input ENVI BIP angle file or .hdr path.")
+    parser.add_argument("mask_file", help="Input water mask GeoTIFF; pixels == 1 are corrected.")
+    parser.add_argument("output_file", help="Output path prefix. ENVI writes .hdr/.img, NetCDF writes .nc.")
+    parser.add_argument("--brdf-model", default="L11", choices=["L11", "M02", "M02SeaDAS", "O25"])
+    parser.add_argument("--output-var", default="Rw_brdf", choices=["Rw_brdf", "rho_ex_w", "nrrs", "C_brdf"])
+    parser.add_argument("--output-format", default="ENVI", choices=["ENVI", "NETCDF"])
+    parser.add_argument("--chunk-size", type=int, default=4096,
+                        help="Number of water pixels per processing batch (default: 4096).")
+    parser.add_argument("--compute-uncertainty", action="store_true",
+                        help="Compute BRDF uncertainty (may increase runtime and output size).")
+    args = parser.parse_args()
+
+    _ds_out, saved_path = process_brdf_files(
+        hyperspectral_file=args.hyperspectral_file,
+        angle_file=args.angle_file,
+        mask_file=args.mask_file,
+        output_file=args.output_file,
+        brdf_model=args.brdf_model,
+        output_var=args.output_var,
+        output_format=args.output_format,
+        chunk_size=args.chunk_size,
+        compute_uncertainty=args.compute_uncertainty,
+    )
+    print(f"BRDF correction finished. Saved to: {saved_path}")
+
+if __name__ == "__main__":
+    main()

ocbrdf/ocbrdf_main3.py

@@ -0,0 +1,908 @@
+import gc
+import os
+import numpy as np
+import sys
+import xarray as xr
+import rasterio
+from rasterio.windows import Window
+import spectral
+
+from brdf_model_M02 import M02
+from brdf_model_M02SeaDAS import M02SeaDAS
+from brdf_model_L11 import L11
+from brdf_model_O25 import O25
+from brdf_utils import ADF_OCP, squeeze_trivial_dims
+
+# 尝试导入 tqdm 进度条库，如果不可用则设为 None
+try:
+    from tqdm import tqdm
+except ImportError:
+    tqdm = None
+
+UNC_LUT_CACHE = {}
+
+"""
+主 BRDF 校正模块
+输入为 xarray 数据集
+所需的光谱维度为 "bands"，其他维度自由
+输入数据集中必填的字段：
+Rw: 方向性海洋反射率
+sza: 太阳天顶角
+vza: 观测天顶角
+raa: 相对方位角（当太阳和观测在同一侧时，raa=0）
+输入数据集中可选的字段：
+Rw_unc: Rw 的不确定度（若缺失，设为零）
+输出数据集中的字段：
+nrrs: 完全归一化的遥感反射率
+rho_ex_w: nrrs * PI
+omega_b: bb/(a+bb)
+eta_b: bbw/bb
+C_brdf: BRDF 校正因子
+brdf_unc: C_brdf 的不确定度
+nrrs_unc: nrrs 的不确定度
+
+使用和适配 brdf_hypercp 模块需保留的信息：
+Brdf_hypercp 是 EUMETSAT 研究项目的一部分
+"BRDF correction of S3 OLCI water reflectance products"（S3 OLCI 水体反射率产品的 BRDF 校正），
+合同号：RB_EUM-CO-21-4600002626-JIG。
+研究团队成员：Davide D'Alimonte (davide.dalimonte@aequora.org)，
+Tamito Kajiyama (tamito.kajiyama@aequora.org)，
+Jaime Pitarch (jaime.pitarchportero@artov.ismar.cnr.it)，
+Vittorio Brando (vittorio.brando@cnr.it)，
+Marco Talone (talone@icm.csic.es) 以及
+Constant Mazeran (constant.mazeran@solvo.fr)。
+BRDF 查找表中的相对方位角遵循 OLCI 约定。详见 https://www.eumetsat.int/media/50720，图 6。
+"""
+
+
+# ---------------------------------------------------------------------------
+# File I/O helpers
+# ---------------------------------------------------------------------------
+
+def _resolve_envi_header(filepath):
+    """Return the ENVI header path accepted by spectral."""
+    if filepath.lower().endswith(".hdr"):
+        return filepath
+
+    hdr_path = filepath + ".hdr"
+    if os.path.exists(hdr_path):
+        return hdr_path
+
+    base, _ext = os.path.splitext(filepath)
+    hdr_path = base + ".hdr"
+    if os.path.exists(hdr_path):
+        return hdr_path
+
+    return filepath
+
+
+def read_bsq(filepath, wavelengths=None, dtype=np.float32):
+    """Read an ENVI BSQ hyperspectral image.
+
+    Parameters
+    ----------
+    filepath : str
+        Path to the ENVI header (.hdr) or image file.
+    wavelengths : array-like, optional
+        Band centre wavelengths in nm.  When *None*, wavelengths are read
+        from the ENVI header if available.
+
+    Returns
+    -------
+    data : ndarray, shape (bands, rows, cols), float64
+    wavelengths : ndarray, shape (bands,)
+    img : spectral image object (metadata / transform access)
+    """
+    img = spectral.open_image(_resolve_envi_header(filepath))
+    data = np.asarray(img.open_memmap(interleave="source"), dtype=dtype)
+    # data = np.transpose(data, (2, 0, 1)) # -> (bands, rows, cols)
+
+    if wavelengths is None:
+        if hasattr(img, 'bands') and img.bands.centers is not None:
+            wavelengths = np.array(img.bands.centers, dtype=np.float32)
+        else:
+            n_bands = data.shape[0]
+            wavelengths = np.arange(1, n_bands + 1, dtype=np.float32)
+
+    return data, np.asarray(wavelengths, dtype=dtype), img
+
+
+def read_bip_angles(filepath, dtype=np.float32):
+    """Read an ENVI BIP angle file and extract the five geometry bands.
+
+    Band layout (1-indexed):
+        band  7  -> SZA  Solar Zenith Angle
+        band  8  -> SAA  Solar Azimuth Angle
+        band  9  -> VZA  Sensor Zenith Angle
+        band 10  -> VAA  Sensor Azimuth Angle
+        band 11  -> RAA  Relative Azimuth Angle
+
+    Parameters
+    ----------
+    filepath : str
+        Path to the ENVI BIP file (header or image file).
+
+    Returns
+    -------
+    dict with keys 'sza', 'saa', 'vza', 'vaa', 'raa'.
+    Each value is a 2-D ndarray of shape (rows, cols), float64.
+    """
+    img = spectral.open_image(_resolve_envi_header(filepath))
+    data = np.asarray(img.open_memmap(interleave="source"), dtype=dtype)
+
+    n_bands = data.shape[2]
+    if n_bands < 11:
+        raise ValueError(
+            f"Angle file has only {n_bands} bands; at least 11 are required "
+            "(bands 7–11 = SZA, SAA, VZA, VAA, RAA)."
+        )
+
+    return {
+        'sza': data[:, :, 6].copy(),
+        'saa': data[:, :, 7].copy(),
+        'vza': data[:, :, 8].copy(),
+        'vaa': data[:, :, 9].copy(),
+        'raa': data[:, :, 10].copy(),
+    }
+
+
+def read_water_mask(filepath):
+    """Read a single-band GeoTIFF water mask.
+
+    Parameters
+    ----------
+    filepath : str
+        Path to the GeoTIFF file.
+
+    Returns
+    -------
+    mask : ndarray, shape (rows, cols), int8
+        Pixels with value 1 indicate water; all other values are non-water.
+    """
+    with rasterio.open(filepath) as src:
+        mask = src.read(1).astype(np.int8)
+    return mask
+
+
+def _block_windows(src, block_size=512):
+    """Yield raster windows for block-wise processing."""
+    rows, cols = src.height, src.width
+    for row in range(0, rows, block_size):
+        row_end = min(row + block_size, rows)
+        for col in range(0, cols, block_size):
+            col_end = min(col + block_size, cols)
+            window = Window(col, row, col_end - col, row_end - row)
+            yield window
+
+
+def _resolve_output_envi_path(output_file):
+    """Return the ENVI image path used for rasterio writes."""
+    if output_file.lower().endswith(".img"):
+        return output_file
+    if output_file.lower().endswith(".hdr"):
+        return output_file[:-4] + ".img"
+    return output_file + ".img"
+
+
+def _read_wavelengths_from_header(hyperspectral_file):
+    """Read wavelength metadata from ENVI header if available."""
+    img = spectral.open_image(_resolve_envi_header(hyperspectral_file))
+    if hasattr(img, 'bands') and img.bands.centers is not None:
+        wavelengths = img.bands.centers
+    else:
+        wavelengths = img.metadata.get('wavelength')
+    if wavelengths is None:
+        return None
+    return np.asarray(wavelengths, dtype=np.float32)
+
+
+def save_brdf_result(
+    ds_out,
+    output_file,
+    source_file=None,
+    output_var="Rw_brdf",
+    output_format="ENVI",
+    dtype=np.float32,
+):
+    """Save the BRDF result to ENVI or NetCDF."""
+    output_format = output_format.upper()
+
+    if output_var not in ds_out:
+        raise KeyError(f"Output variable '{output_var}' not found in result dataset.")
+
+    if output_format == "NETCDF":
+        if not output_file.lower().endswith(".nc"):
+            output_file = output_file + ".nc"
+        ds_out.to_netcdf(output_file)
+        return output_file
+
+    if output_format != "ENVI":
+        raise ValueError("output_format must be 'ENVI' or 'NETCDF'.")
+
+    cube = np.asarray(ds_out[output_var].values, dtype=dtype)
+    cube = np.transpose(cube, (1, 2, 0))
+
+    if output_file.lower().endswith(".hdr"):
+        hdr_path = output_file
+    else:
+        hdr_path = output_file + ".hdr"
+
+    metadata = {
+        "interleave": "bsq",
+        "description": f"BRDF corrected result: {output_var}",
+        "bands": cube.shape[2],
+        "lines": cube.shape[0],
+        "samples": cube.shape[1],
+    }
+    if "bands" in ds_out.coords:
+        metadata["wavelength"] = [str(v) for v in ds_out.coords["bands"].values.tolist()]
+
+    if source_file is not None:
+        metadata["source_file"] = source_file
+
+    spectral.envi.save_image(
+        hdr_path,
+        cube,
+        dtype=dtype,
+        interleave="bsq",
+        metadata=metadata,
+        force=True,
+    )
+    return hdr_path
+
+
+def _get_output_filename(base_output_file: str, var_name: str) -> str:
+    """Generate output filename with variable suffix.
+    
+    Examples:
+        input.hdr + Rw_brdf     -> input_Rw_brdf.hdr
+        input.hdr + brdf_unc    -> input_brdf_unc.hdr
+    """
+    base, ext = os.path.splitext(base_output_file)
+    if ext.lower() in ('.hdr', '.img'):
+        base = base  # keep base without extension
+    suffix = f"_{var_name}" if var_name.lower() != "rw_brdf" else ""
+    return f"{base}{suffix}{ext}"
+
+
+def save_multiple_brdf_results(
+    ds_out,
+    output_file,
+    source_file=None,
+    output_vars=None,
+    output_format="ENVI",
+    dtype=np.float32,
+):
+    """Save multiple output variables to separate files."""
+    if output_vars is None or len(output_vars) == 0:
+        output_vars = ["Rw_brdf"]
+    
+    saved_files = []
+    for var in output_vars:
+        out_file = _get_output_filename(output_file, var)
+        try:
+            saved_path = save_brdf_result(
+                ds_out=ds_out,
+                output_file=out_file,
+                source_file=source_file,
+                output_var=var,
+                output_format=output_format,
+                dtype=dtype,
+            )
+            saved_files.append((var, saved_path))
+            print(f"  Saved {var} to: {saved_path}")
+        except KeyError as e:
+            print(f"  Warning: {e}. Skipping {var}.")
+    
+    return saved_files
+
+
+def run_brdf_correction_block(
+    hyperspectral_file,
+    angle_file,
+    mask_file,
+    output_file,
+    wavelengths=None,
+    brdf_model='L11',
+    output_vars=None,           # 支持列表
+    block_size=512,
+    compute_uncertainty=False,
+    progress_bar=True,
+):
+    """
+    Run BRDF correction block-wise and write multiple variables to separate ENVI files.
+    """
+    if output_vars is None:
+        output_vars = ["Rw_brdf"]
+    if isinstance(output_vars, str):
+        output_vars = [output_vars]
+
+    # 1. Open all inputs
+    src_rw = rasterio.open(hyperspectral_file)
+    n_bands = src_rw.count
+    n_rows = src_rw.height
+    n_cols = src_rw.width
+
+    src_ang = rasterio.open(angle_file)
+    if src_ang.height != n_rows or src_ang.width != n_cols:
+        raise ValueError("Angle file dimensions do not match hyperspectral file.")
+
+    src_mask = rasterio.open(mask_file)
+    if src_mask.height != n_rows or src_mask.width != n_cols:
+        raise ValueError("Mask file dimensions do not match hyperspectral file.")
+
+    # 2. Wavelengths
+    if wavelengths is None:
+        wavelengths = _read_wavelengths_from_header(hyperspectral_file)
+    if wavelengths is None:
+        wavelengths = np.arange(1, n_bands + 1, dtype=np.float32)
+    bands_xr = xr.DataArray(wavelengths, dims=('bands',), coords={'bands': wavelengths}, name='bands')
+
+    # 3. Prepare multiple output files (one per variable)
+    out_meta = src_rw.meta.copy()
+    out_meta.update({
+        'driver': 'ENVI',
+        'dtype': 'float32',
+        'count': n_bands,
+        'interleave': 'bsq',
+    })
+
+    # Open one output file for each variable
+    output_dsts = {}
+    output_paths = {}
+    for var_name in output_vars:
+        var_file = _get_output_filename(output_file, var_name)
+        out_path = _resolve_output_envi_path(var_file)
+        output_dsts[var_name] = rasterio.open(out_path, 'w', **out_meta)
+        output_paths[var_name] = out_path
+        print(f"  Writing {var_name} to {out_path}")
+
+    # 4. Band indices in angle file (1-based)
+    band_sza = 7
+    band_vza = 9
+    band_raa = 11
+
+    # 5. Map output_var to ds_corr field (support uncertainty variables)
+    var_map = {
+        'Rw_brdf': 'rho_ex_w',
+        'rho_ex_w': 'rho_ex_w',
+        'nrrs': 'nrrs',
+        'C_brdf': 'C_brdf',
+        'brdf_unc': 'brdf_unc',
+        'nrrs_unc': 'nrrs_unc',
+    }
+
+    # Validate all requested output variables
+    for v in output_vars:
+        if v not in var_map and v not in ['Rw_brdf']:
+            raise ValueError(
+                f"Unknown output_var '{v}'. Allowed: {list(var_map.keys())}"
+            )
+
+    # 6. Progress bar setup
+    iterator = _block_windows(src_rw, block_size=block_size)
+    if progress_bar and tqdm is not None:
+        total_blocks = ((n_rows + block_size - 1) // block_size) * ((n_cols + block_size - 1) // block_size)
+        iterator = tqdm(iterator, total=total_blocks, desc="Processing blocks", unit="block")
+    elif progress_bar:
+        print("Warning: tqdm not installed, progress bar disabled.")
+
+    # Statistics for performance monitoring
+    blocks_processed = 0
+    blocks_skipped = 0
+    total_water_pixels = 0
+
+    # 7. Block-wise processing with early filtering
+    for window in iterator:
+        # Step 1: Read mask FIRST to enable early exit for pure land blocks
+        mask_block = src_mask.read(1, window=window).astype(np.int8)
+        water_mask = (mask_block == 1)
+        n_water = np.count_nonzero(water_mask)
+
+        if n_water == 0:
+            # Fast path: pure land/no-water block - minimal processing
+            blocks_skipped += 1
+            # Still need to read hyperspectral for Rw_brdf output, but skip angles
+            rw_block = src_rw.read(window=window).astype(np.float32)
+            
+            output_blocks = {}
+            for var_name in output_vars:
+                if var_name == 'Rw_brdf':
+                    output_blocks[var_name] = rw_block.copy()
+                else:
+                    output_blocks[var_name] = np.full_like(rw_block, np.nan, dtype=np.float32)
+            
+            for var_name, block_data in output_blocks.items():
+                output_dsts[var_name].write(block_data, window=window)
+            
+            del rw_block, mask_block, water_mask
+            gc.collect()
+            continue
+
+        # Step 2: Has water pixels - full processing
+        blocks_processed += 1
+        total_water_pixels += n_water
+
+        # Read remaining data only for blocks that need processing
+        rw_block = src_rw.read(window=window).astype(np.float32)
+        sza_block = src_ang.read(band_sza, window=window).astype(np.float32)
+        vza_block = src_ang.read(band_vza, window=window).astype(np.float32)
+        raa_block = src_ang.read(band_raa, window=window).astype(np.float32)
+
+        water_idx = np.where(water_mask)
+
+        # Prepare output blocks for all variables
+        output_blocks = {}
+        for var_name in output_vars:
+            if var_name == 'Rw_brdf':
+                output_blocks[var_name] = rw_block.copy()
+            else:
+                output_blocks[var_name] = np.full_like(rw_block, np.nan, dtype=np.float32)
+
+        # Process water pixels
+        rw_water = rw_block[:, water_idx[0], water_idx[1]].T  # (n_water, bands)
+        sza_water = sza_block[water_idx]
+        vza_water = vza_block[water_idx]
+        raa_water = raa_block[water_idx]
+
+        pixel_ids = np.arange(n_water)
+        ds_block = xr.Dataset(
+            {
+                'Rw': xr.DataArray(rw_water, dims=('n', 'bands'), coords={'n': pixel_ids, 'bands': bands_xr}),
+                'sza': xr.DataArray(sza_water, dims='n', coords={'n': pixel_ids}),
+                'vza': xr.DataArray(vza_water, dims='n', coords={'n': pixel_ids}),
+                'raa': xr.DataArray(raa_water, dims='n', coords={'n': pixel_ids}),
+            }
+        )
+
+        ds_corr = brdf_prototype(
+            ds_block,
+            brdf_model=brdf_model,
+            compute_uncertainty=compute_uncertainty,
+        )
+
+        # Write each requested variable to the water pixels
+        for var_name in output_vars:
+            field_name = var_map.get(var_name, 'rho_ex_w')
+            if var_name == 'Rw_brdf':
+                field_name = 'rho_ex_w'
+            if field_name in ds_corr:
+                values = ds_corr[field_name].values.T
+                output_blocks[var_name][:, water_idx[0], water_idx[1]] = values
+
+        # Write all variables for this block
+        for var_name, block_data in output_blocks.items():
+            output_dsts[var_name].write(block_data, window=window)
+
+        del (rw_block, sza_block, vza_block, raa_block, mask_block, output_blocks, 
+             water_mask, water_idx, rw_water, sza_water, vza_water, raa_water, 
+             pixel_ids, ds_block, ds_corr)
+        gc.collect()
+
+    # Close all files
+    src_rw.close()
+    src_ang.close()
+    src_mask.close()
+    for dst in output_dsts.values():
+        dst.close()
+
+    # Print performance statistics
+    total_blocks = blocks_processed + blocks_skipped
+    if total_blocks > 0:
+        skip_ratio = blocks_skipped / total_blocks * 100
+        print(f"Block processing completed:")
+        print(f"  - Total blocks: {total_blocks}")
+        print(f"  - Processed blocks (with water): {blocks_processed}")
+        print(f"  - Skipped blocks (no water): {blocks_skipped} ({skip_ratio:.1f}%)")
+        print(f"  - Total water pixels processed: {total_water_pixels:,}")
+    
+    return list(output_paths.values())
+
+
+# ---------------------------------------------------------------------------
+# Top-level pipeline
+# ---------------------------------------------------------------------------
+def brdf_uncertainty(ds, adf=None):
+    ''' Compute uncertainty of BRDF factor and propagate to nrrs '''
+
+    # Read LUT
+    if adf is None:
+        adf = ADF_OCP
+
+    unc_lut_path = adf % 'UNC'
+    if unc_lut_path not in UNC_LUT_CACHE:
+        UNC_LUT_CACHE[unc_lut_path] = xr.open_dataset(unc_lut_path, engine='netcdf4')
+    LUT = UNC_LUT_CACHE[unc_lut_path]
+
+    # Interpolate relative uncertainty
+    unc = LUT['unc'].interp(lambda_unc=ds.bands, theta_s_unc=ds.sza, theta_v_unc=ds.vza,
+                            delta_phi_unc=ds.raa)
+
+    # Compute absolute uncertainty of factor
+    ds['brdf_unc'] = unc * ds['C_brdf']
+
+    # Flag BRDF_unc where NaN and set to 0
+    ds['brdf_unc_fail'] = np.isnan(ds['brdf_unc'])
+    ds['brdf_unc'] = xr.where(ds['brdf_unc_fail'], 0, ds['brdf_unc'])
+
+    # Propagate to nrrs
+    Rwex_unc2 = ds['brdf_unc'] * ds['brdf_unc'] * ds['Rw'] * ds['Rw']
+    if 'Rw_unc' in ds:
+        Rwex_unc2 += ds['C_brdf'] * ds['C_brdf'] * ds['Rw_unc'] * ds['Rw_unc']
+    ds['nrrs_unc'] = np.sqrt(Rwex_unc2)/np.pi
+
+    return ds
+
+
+def brdf_prototype(ds, adf=None, brdf_model='L11', compute_uncertainty=True):
+    # 测试 BRDF 模型在 GUI 中不支持：强制覆盖
+    # brdf_model = 'M02SeaDAS'
+
+    # 压缩单一维度（例如，投放次数、提取次数等）以避免插值问题
+    ds, squeezedDims = squeeze_trivial_dims(ds)
+
+    # Initialise model
+    if brdf_model == 'M02':
+        BRDF_model = M02(bands=ds.bands, aot=ds.aot, wind=ds.wind, adf=None)  # Don't use brdf_py.ADF context
+    elif brdf_model == 'M02SeaDAS':
+        BRDF_model = M02SeaDAS(bands=ds.bands, adf=None)  # Don't use brdf_py.ADF context
+    elif brdf_model == 'L11':
+        BRDF_model = L11(bands=ds.bands, adf=None)  # Don't use brdf_py.ADF context
+    elif brdf_model == 'O25':
+        BRDF_model = O25(bands=ds.bands, adf=None)  # Don't use brdf_py.ADF context
+    else:
+        print("BRDF model %s not supported" % brdf_model)
+        sys.exit(1)
+
+    # Init pixel
+    BRDF_model.init_pixels(ds['sza'], ds['vza'], ds['raa'])
+
+    # Compute IOP and normalize by iterating
+    ds['nrrs'] = ds['Rw'] / np.pi
+
+    ds['convergeFlag'] = (0 * ds['sza']).astype(bool)
+    ds['C_brdf'] = 0 * ds['nrrs'] + 1
+
+    for iter_brdf in range(int(BRDF_model.niter)):
+
+        # M02: Initialise chl_iter
+        if brdf_model in ['M02', 'M02SeaDAS'] and (iter_brdf == 0):
+            chl_iter = {}
+            ds['log10_chl'] = 0 * ds['sza'] + float(np.log10(BRDF_model.OC4MEchl0))
+            chl_iter[-1]    = 0 * ds['sza'] + float(BRDF_model.OC4MEchl0)
+
+        ds = BRDF_model.backward(ds, iter_brdf)
+
+        # M02: Check convergence (dummy for M02SeaDAS for the moment... epsilon set to 0)
+        if brdf_model in ['M02', 'M02SeaDAS']:
+            chl_iter[iter_brdf] = 10 ** ds['log10_chl']
+            #  Check if convergence is reached |chl_old-chl_new| < epsilon * chl_new
+            ds['convergeFlag'] = (ds['convergeFlag']) | (
+                (np.abs(chl_iter[iter_brdf - 1] - chl_iter[iter_brdf]) < float(BRDF_model.OC4MEepsilon) * chl_iter[
+                    iter_brdf]))
+
+        # Apply forward model in both geometries
+        forward_mod = BRDF_model.forward(ds)
+        forward_mod0 = BRDF_model.forward(ds, normalized=True)
+
+        ratio = forward_mod0 / forward_mod
+
+        # Drop remnant coordinates to avoid ambiguities in the update of the BRDF factor.
+        for coord in ratio.coords:
+            if coord not in ds['C_brdf'].coords:
+                ratio = ratio.drop(coord)
+
+        # Normalize reflectance
+        ds['C_brdf'] = xr.where(ds['convergeFlag'], ds['C_brdf'], ratio)
+        ds['nrrs'] = ds['Rw'] / np.pi * ds['C_brdf']
+
+    # Flag BRDF where NaN and set to 1 (no correction applied).
+    ds['C_brdf_fail'] = np.isnan(ds['C_brdf'])
+    ds['C_brdf'] = xr.where(ds['C_brdf_fail'], 1, ds['C_brdf'])
+    ds['nrrs'] = xr.where(ds['C_brdf_fail'], ds['Rw'] / np.pi, ds['nrrs'])
+
+    # If QAA_fail is raised, raise C_brdf_fail (but still apply C_brdf).
+    if 'QAA_fail' in ds:
+        ds['C_brdf_fail'] = (ds['C_brdf_fail']) | (ds['QAA_fail'])
+
+    # Compute uncertainty only if requested
+    if compute_uncertainty:
+        ds = brdf_uncertainty(ds)
+
+    # Compute flag
+    ds['flags_level2'] = ds['Rw'] * 0  # TODO
+
+    # Convert to reflectance unit
+    ds['rho_ex_w'] = ds['nrrs'] * np.pi
+
+    # Expand squeezed trivial dimensions
+    for dim,d0 in squeezedDims.items():
+        ds = ds.expand_dims(dim,axis=d0)
+
+    return ds
+
+
+def run_brdf_correction(
+    hyperspectral_file,
+    angle_file,
+    mask_file,
+    wavelengths=None,
+    brdf_model='L11',
+    chunk_size=4096,
+    output_var=None,          # if None, keep full output (original behaviour)
+    compute_uncertainty=True, # compute uncertainty only if requested
+    progress_bar=True,        # 是否显示进度条
+):
+    """BRDF correction pipeline for water pixels in a hyperspectral scene.
+
+    ... (docstring unchanged, update description) ...
+    """
+    # ------------------------------------------------------------------
+    # 1. Read all inputs
+    # ------------------------------------------------------------------
+    rw_data, wvl, _img_meta = read_bsq(hyperspectral_file, wavelengths, dtype=np.float32)
+    angles = read_bip_angles(angle_file, dtype=np.float32)
+    water_mask = read_water_mask(mask_file)
+
+    n_bands, n_rows, n_cols = rw_data.shape
+
+    # Sanity-check spatial dimensions
+    if angles['sza'].shape != (n_rows, n_cols):
+        raise ValueError(
+            f"Angle file spatial shape {angles['sza'].shape} does not match "
+            f"hyperspectral image shape ({n_rows}, {n_cols})."
+        )
+    if water_mask.shape != (n_rows, n_cols):
+        raise ValueError(
+            f"Mask shape {water_mask.shape} does not match "
+            f"hyperspectral image shape ({n_rows}, {n_cols})."
+        )
+
+    # ------------------------------------------------------------------
+    # 2. Identify water pixels
+    # ------------------------------------------------------------------
+    water_rows, water_cols = np.where(water_mask == 1)
+    n_water = water_rows.size
+    if n_water == 0:
+        raise RuntimeError("No water pixels found in the mask (no pixels == 1).")
+
+    # ------------------------------------------------------------------
+    # 3. Determine which output variable(s) we need to keep
+    # ------------------------------------------------------------------
+    # If output_var is None, we keep everything (original behaviour)
+    keep_all = (output_var is None)
+
+    # Map output_var to the field in ds_corr
+    var_map = {
+        'Rw_brdf': 'rho_ex_w',
+        'rho_ex_w': 'rho_ex_w',
+        'nrrs': 'nrrs',
+        'C_brdf': 'C_brdf'
+    }
+    if not keep_all:
+        if output_var not in var_map:
+            raise ValueError(f"Unknown output_var '{output_var}'. Allowed: {list(var_map.keys())}")
+        field_name = var_map[output_var]
+    else:
+        field_name = None
+
+    # ------------------------------------------------------------------
+    # 4. Allocate output arrays (only those needed)
+    # ------------------------------------------------------------------
+    shape3 = (n_bands, n_rows, n_cols)
+    shape2 = (n_rows, n_cols)
+
+    if keep_all:
+        # Original behaviour: allocate all arrays
+        nrrs_full        = np.full(shape3, np.nan, dtype=np.float32)
+        rho_ex_w_full    = np.full(shape3, np.nan, dtype=np.float32)
+        C_brdf_full      = np.full(shape3, np.nan, dtype=np.float32)
+        brdf_unc_full    = np.full(shape3, np.nan, dtype=np.float32)
+        nrrs_unc_full    = np.full(shape3, np.nan, dtype=np.float32)
+        C_brdf_fail_full = np.zeros(shape2, dtype=bool)
+        rw_brdf_full     = rw_data.copy()
+    else:
+        # Only allocate the requested output array
+        output_full = np.full(shape3, np.nan, dtype=np.float32)
+        # We may also need C_brdf_fail if requested? Not for output_var, but could be used internally? Skip.
+        # For simplicity, we do not allocate extra arrays.
+
+    # ------------------------------------------------------------------
+    # 5–6. Chunked BRDF correction
+    # ------------------------------------------------------------------
+    n_chunks = int(np.ceil(n_water / chunk_size))
+
+    # 初始化进度条
+    iterator = range(n_chunks)
+    if progress_bar and tqdm is not None:
+        iterator = tqdm(iterator, desc="Processing chunks", unit="chunk")
+    elif progress_bar:
+        print("Warning: tqdm not installed, progress bar disabled.")
+
+    for chunk_idx in iterator:
+        lo = chunk_idx * chunk_size
+        hi = min(lo + chunk_size, n_water)
+
+        rows_c = water_rows[lo:hi]
+        cols_c = water_cols[lo:hi]
+
+        # Extract batch vectors
+        rw_c  = rw_data[:, rows_c, cols_c].T          # (chunk, bands)
+        sza_c = angles['sza'][rows_c, cols_c]          # (chunk,)
+        vza_c = angles['vza'][rows_c, cols_c]
+        raa_c = angles['raa'][rows_c, cols_c]
+
+        pixel_idx_c = np.arange(hi - lo)
+        ds_chunk = xr.Dataset(
+            {
+                'Rw':  xr.DataArray(rw_c,  dims=('n', 'bands'),
+                                    coords={'n': pixel_idx_c, 'bands': wvl}),
+                'sza': xr.DataArray(sza_c, dims='n',
+                                    coords={'n': pixel_idx_c}),
+                'vza': xr.DataArray(vza_c, dims='n',
+                                    coords={'n': pixel_idx_c}),
+                'raa': xr.DataArray(raa_c, dims='n',
+                                    coords={'n': pixel_idx_c}),
+            }
+        )
+
+        # Run BRDF model with uncertainty control
+        ds_corr = brdf_prototype(ds_chunk, brdf_model=brdf_model, compute_uncertainty=compute_uncertainty)
+
+        if keep_all:
+            # Write all results back
+            nrrs_full[:, rows_c, cols_c]        = ds_corr['nrrs'].values.T
+            rho_ex_w_full[:, rows_c, cols_c]    = ds_corr['rho_ex_w'].values.T
+            C_brdf_full[:, rows_c, cols_c]      = ds_corr['C_brdf'].values.T
+            if compute_uncertainty:
+                brdf_unc_full[:, rows_c, cols_c]    = ds_corr['brdf_unc'].values.T
+                nrrs_unc_full[:, rows_c, cols_c]    = ds_corr['nrrs_unc'].values.T
+            C_brdf_fail_full[rows_c, cols_c]    = ds_corr['C_brdf_fail'].values
+            rw_brdf_full[:, rows_c, cols_c]     = ds_corr['rho_ex_w'].values.T
+        else:
+            # Only write the requested output field
+            values = ds_corr[field_name].values.T
+            output_full[:, rows_c, cols_c] = values
+
+        # Release batch objects so GC can reclaim memory
+        del rw_c, sza_c, vza_c, raa_c, pixel_idx_c, ds_chunk, ds_corr
+        gc.collect()
+
+    # ------------------------------------------------------------------
+    # 7. Assemble full-scene output Dataset
+    # ------------------------------------------------------------------
+    full_coords = {'bands': wvl, 'y': np.arange(n_rows), 'x': np.arange(n_cols)}
+    dims3 = ('bands', 'y', 'x')
+    dims2 = ('y', 'x')
+
+    if keep_all:
+        # Original behaviour: return all variables
+        ds_out = xr.Dataset(
+            {
+                'Rw':          xr.DataArray(rw_data,          dims=dims3, coords=full_coords),
+                'Rw_brdf':     xr.DataArray(rw_brdf_full,     dims=dims3, coords=full_coords),
+                'nrrs':        xr.DataArray(nrrs_full,        dims=dims3, coords=full_coords),
+                'rho_ex_w':    xr.DataArray(rho_ex_w_full,    dims=dims3, coords=full_coords),
+                'C_brdf':      xr.DataArray(C_brdf_full,      dims=dims3, coords=full_coords),
+                'C_brdf_fail': xr.DataArray(C_brdf_fail_full, dims=dims2),
+                'water_mask':  xr.DataArray(water_mask,       dims=dims2),
+                'sza':         xr.DataArray(angles['sza'],    dims=dims2),
+                'saa':         xr.DataArray(angles['saa'],    dims=dims2),
+                'vza':         xr.DataArray(angles['vza'],    dims=dims2),
+                'vaa':         xr.DataArray(angles['vaa'],    dims=dims2),
+                'raa':         xr.DataArray(angles['raa'],    dims=dims2),
+            }
+        )
+        if compute_uncertainty:
+            ds_out['brdf_unc'] = xr.DataArray(brdf_unc_full, dims=dims3, coords=full_coords)
+            ds_out['nrrs_unc'] = xr.DataArray(nrrs_unc_full, dims=dims3, coords=full_coords)
+    else:
+        # Return only the requested output variable + geometry + mask
+        ds_out = xr.Dataset(
+            {
+                output_var:    xr.DataArray(output_full,      dims=dims3, coords=full_coords),
+                'Rw':          xr.DataArray(rw_data,          dims=dims3, coords=full_coords),
+                'water_mask':  xr.DataArray(water_mask,       dims=dims2),
+                'sza':         xr.DataArray(angles['sza'],    dims=dims2),
+                'saa':         xr.DataArray(angles['saa'],    dims=dims2),
+                'vza':         xr.DataArray(angles['vza'],    dims=dims2),
+                'vaa':         xr.DataArray(angles['vaa'],    dims=dims2),
+                'raa':         xr.DataArray(angles['raa'],    dims=dims2),
+            }
+        )
+        # If uncertainties were computed, we could add them optionally, but user didn't request them.
+        # For memory efficiency, we skip them here.
+
+    return ds_out
+
+
+def process_brdf_files(
+    hyperspectral_file,
+    angle_file,
+    mask_file,
+    output_file,
+    wavelengths=None,
+    brdf_model="L11",
+    output_vars=None,
+    output_format="ENVI",
+    chunk_size=4096,
+    block_size=512,
+    compute_uncertainty=False,
+    progress_bar=True,
+):
+    """Run BRDF correction from three input files and save multiple variables to separate files."""
+    if output_format.upper() != "ENVI":
+        raise NotImplementedError("Only ENVI output format is supported in block-wise version.")
+
+    if output_vars is None:
+        output_vars = ["Rw_brdf"]
+    if isinstance(output_vars, str):
+        output_vars = [output_vars]
+
+    print(f"Computing BRDF correction with variables: {output_vars}")
+    print(f"Output base name: {output_file}")
+
+    saved_paths = run_brdf_correction_block(
+        hyperspectral_file=hyperspectral_file,
+        angle_file=angle_file,
+        mask_file=mask_file,
+        output_file=output_file,
+        wavelengths=wavelengths,
+        brdf_model=brdf_model,
+        output_vars=output_vars,
+        block_size=block_size,
+        compute_uncertainty=compute_uncertainty,
+        progress_bar=progress_bar,
+    )
+    return None, saved_paths
+
+
+def main():
+    import argparse
+
+    parser = argparse.ArgumentParser(description="Water-region BRDF correction for hyperspectral data.")
+    parser.add_argument("hyperspectral_file", help="Input ENVI BSQ hyperspectral file or .hdr path.")
+    parser.add_argument("angle_file", help="Input ENVI BIP angle file or .hdr path.")
+    parser.add_argument("mask_file", help="Input water mask GeoTIFF; pixels == 1 are corrected.")
+    parser.add_argument("output_file", help="Output path prefix. ENVI writes .hdr/.img, NetCDF writes .nc.")
+    parser.add_argument("--brdf-model", default="L11", choices=["L11", "M02", "M02SeaDAS", "O25"])
+    parser.add_argument("--output-var", nargs='+', default=["Rw_brdf"],
+                        help="Output variables to save (one or more). "
+                             "Supported: Rw_brdf, rho_ex_w, nrrs, C_brdf, brdf_unc, nrrs_unc. "
+                             "Example: --output-var Rw_brdf nrrs brdf_unc")
+    parser.add_argument("--output-format", default="ENVI", choices=["ENVI"])
+    parser.add_argument("--chunk-size", type=int, default=4096,
+                        help="Number of water pixels per processing batch (default: 4096).")
+    parser.add_argument("--block-size", type=int, default=512,
+                        help="Spatial block size for block-wise processing (default: 512).")
+    parser.add_argument("--compute-uncertainty", action="store_true",
+                        help="Compute BRDF uncertainty (may increase runtime and output size).")
+    parser.add_argument("--no-progress-bar", action="store_true",
+                        help="Disable progress bar (useful for non-interactive environments).")
+    args = parser.parse_args()
+
+    progress_bar = not args.no_progress_bar
+
+    _ds_out, saved_paths = process_brdf_files(
+        hyperspectral_file=args.hyperspectral_file,
+        angle_file=args.angle_file,
+        mask_file=args.mask_file,
+        output_file=args.output_file,
+        brdf_model=args.brdf_model,
+        output_vars=args.output_var,
+        output_format=args.output_format,
+        chunk_size=args.chunk_size,
+        block_size=args.block_size,
+        compute_uncertainty=args.compute_uncertainty,
+        progress_bar=progress_bar,
+    )
+
+    if isinstance(saved_paths, list):
+        print(f"BRDF correction finished. Saved {len(saved_paths)} files:")
+        for path in saved_paths:
+            print(f"  - {path}")
+    else:
+        print(f"BRDF correction finished. Saved to: {saved_paths}")
+
+if __name__ == "__main__":
+    main()

ocbrdf/src/change_bsq.py

@@ -0,0 +1,115 @@
+import numpy as np
+import os
+import struct
+
+def read_envi_header(header_path):
+    """
+    读取 ENVI 头文件，返回元数据字典。
+    """
+    metadata = {}
+    with open(header_path, 'r') as f:
+        lines = f.readlines()
+    for line in lines:
+        line = line.strip()
+        if not line or line.startswith(';'):
+            continue
+        if '=' in line:
+            key, value = line.split('=', 1)
+            key = key.strip().lower()
+            value = value.strip().strip('{}')
+            # 尝试转换数值
+            if value.isdigit():
+                metadata[key] = int(value)
+            elif value.replace('.', '', 1).isdigit():
+                metadata[key] = float(value)
+            else:
+                metadata[key] = value
+    return metadata
+
+def write_envi_header(output_header, metadata):
+    """
+    将元数据写入 ENVI 头文件。
+    """
+    with open(output_header, 'w') as f:
+        f.write('ENVI\n')
+        for key, val in metadata.items():
+            if isinstance(val, (int, float)):
+                f.write(f'{key} = {val}\n')
+            else:
+                f.write(f'{key} = {{{val}}}\n')
+
+def convert_bsq(input_bsq, output_bsq):
+    """
+    主转换函数：逐波段读取，缩放并转换数据类型。
+    """
+    # 获取头文件路径
+    input_hdr = input_bsq + '.hdr' if not input_bsq.endswith('.hdr') else input_bsq.replace('.hdr', '')
+    if not os.path.exists(input_hdr):
+        base, ext = os.path.splitext(input_bsq)
+        input_hdr = base + '.hdr'
+        if not os.path.exists(input_hdr):
+            raise FileNotFoundError(f"头文件 {input_hdr} 不存在")
+
+    # 读取头文件信息
+    meta = read_envi_header(input_hdr)
+    samples = meta.get('samples')
+    lines = meta.get('lines')
+    bands = meta.get('bands')
+    data_type = meta.get('data type')
+    interleave = meta.get('interleave', 'bsq').lower()
+
+    if interleave != 'bsq':
+        raise ValueError("本脚本仅支持 BSQ 格式，当前 interleave = {}".format(interleave))
+
+    # 根据 data_type 确定输入数据类型
+    dtype_map = {1: np.uint8, 2: np.int16, 3: np.int32, 4: np.float32, 5: np.float64, 12: np.uint16}
+    if data_type not in dtype_map:
+        raise ValueError(f"不支持的 data type: {data_type}")
+    in_dtype = dtype_map[data_type]
+
+    print(f"输入文件尺寸: 波段数={bands}, 行数={lines}, 列数={samples}")
+    print(f"原始 data type = {data_type} ({in_dtype}), 将缩放为 uint16 (type 12)")
+
+    # 计算每个波段的字节数
+    band_size = lines * samples * np.dtype(in_dtype).itemsize
+
+    # 打开输入文件，输出文件
+    with open(input_bsq, 'rb') as fin, open(output_bsq, 'wb') as fout:
+        for band_idx in range(bands):
+            # 定位到当前波段的起始位置
+            fin.seek(band_idx * band_size)
+            # 读取整个波段
+            band_data = np.fromfile(fin, dtype=in_dtype, count=lines * samples)
+            # 重塑为 (lines, samples)
+            band_data = band_data.reshape((lines, samples))
+
+            # 缩放并转换为 uint16
+            scaled = np.round(band_data * 10000).astype(np.uint16)
+
+            # 写入输出文件（BSQ 顺序）
+            fout.write(scaled.tobytes())
+
+            # 进度提示
+            if (band_idx + 1) % 10 == 0 or band_idx == bands - 1:
+                print(f"已处理波段 {band_idx+1}/{bands}")
+
+    # 生成新的头文件
+    output_hdr = output_bsq + '.hdr'
+    new_meta = meta.copy()
+    new_meta['data type'] = 12          # uint16
+    new_meta['interleave'] = 'bsq'
+    # 可选：修改描述
+    if 'description' in new_meta:
+        new_meta['description'] = new_meta['description'] + ' (scaled by 10000)'
+    else:
+        new_meta['description'] = 'Scaled to uint16 (factor 10000)'
+    write_envi_header(output_hdr, new_meta)
+
+    print(f"转换完成：{output_bsq} 已生成，data type = 12 (uint16)")
+
+if __name__ == '__main__':
+    import sys
+    if len(sys.argv) != 3:
+        print("用法: python script.py input.bsq output.bsq")
+        sys.exit(1)
+    convert_bsq(sys.argv[1], sys.argv[2])

ocbrdf/src/cli_mask.py

@@ -0,0 +1,95 @@
+#!/usr/bin/env python3
+"""
+Mask alignment tool: crop and align a TIFF mask to a hyperspectral BSQ file,
+then set areas outside the hyperspectral valid data mask to 0.
+"""
+
+import argparse
+import numpy as np
+import rasterio
+from rasterio.warp import reproject, Resampling
+from rasterio.enums import Resampling as ResamplingEnum
+import os
+
+def get_valid_mask_all_bands(hyperspectral_file):
+    with rasterio.open(hyperspectral_file) as src:
+        # Read all bands
+        data = src.read()  # shape (bands, height, width)
+        # Mask where all bands are zero
+        all_zero = np.all(data == 0, axis=0)
+        valid = ~all_zero
+        return valid, src.profile
+
+def align_mask_to_hyperspectral(mask_file, hyperspectral_profile, resampling=Resampling.bilinear):
+    """
+    Reproject/resample the mask to match the hyperspectral spatial grid.
+    Returns the aligned mask as a 2D array.
+    """
+    with rasterio.open(mask_file) as mask_src:
+        # Create an empty array for the output
+        out_shape = (hyperspectral_profile['height'], hyperspectral_profile['width'])
+        out_array = np.empty(out_shape, dtype=mask_src.dtypes[0])
+        # Reproject the first band
+        reproject(
+            source=rasterio.band(mask_src, 1),
+            destination=out_array,
+            src_transform=mask_src.transform,
+            src_crs=mask_src.crs,
+            dst_transform=hyperspectral_profile['transform'],
+            dst_crs=hyperspectral_profile['crs'],
+            resampling=resampling
+        )
+        return out_array
+
+def main():
+    parser = argparse.ArgumentParser(
+        description="Crop TIFF mask to hyperspectral extent and zero-out outside valid data mask."
+    )
+    parser.add_argument('hyperspectral', help='Input hyperspectral BSQ file (with .hdr)')
+    parser.add_argument('mask', help='Input TIFF mask file')
+    parser.add_argument('output', help='Output TIFF file')
+    parser.add_argument('--band', type=int, default=0,
+                        help='Band index (0-based) to use for valid data mask (default: 0)')
+    parser.add_argument('--nodata', type=float, default=None,
+                        help='NoData value in hyperspectral data; if not given, uses file metadata')
+    parser.add_argument('--resample', choices=['nearest', 'bilinear', 'cubic'], default='bilinear',
+                        help='Resampling method for mask alignment (default: bilinear)')
+    args = parser.parse_args()
+
+    # Mapping of resampling names to rasterio enums
+    resample_map = {
+        'nearest': Resampling.nearest,
+        'bilinear': Resampling.bilinear,
+        'cubic': Resampling.cubic
+    }
+    resampling = resample_map[args.resample]
+
+    # Get valid data mask from hyperspectral file
+    print(f"Reading valid mask from {args.hyperspectral} (band {args.band})...")
+    valid_mask, hs_profile = get_valid_mask_all_bands(args.hyperspectral)
+    print(f"Hyperspectral dimensions: {hs_profile['height']} x {hs_profile['width']}")
+
+    # Reproject mask to hyperspectral grid
+    print(f"Aligning mask {args.mask} to hyperspectral grid...")
+    aligned_mask = align_mask_to_hyperspectral(args.mask, hs_profile, resampling)
+
+    # Zero-out areas outside valid_mask
+    print("Applying valid mask...")
+    result = np.where(valid_mask, aligned_mask, 0).astype(aligned_mask.dtype)
+
+    # Save result as TIFF
+    print(f"Saving to {args.output}...")
+    out_profile = hs_profile.copy()
+    out_profile.update({
+        'driver': 'GTiff',
+        'count': 1,
+        'dtype': result.dtype,
+        'compress': 'lzw'          # optional compression
+    })
+    with rasterio.open(args.output, 'w', **out_profile) as dst:
+        dst.write(result, 1)
+
+    print("Done.")
+
+if __name__ == '__main__':
+    main()

ocbrdf/src/sha2dat.py

@@ -0,0 +1,72 @@
+import numpy as np
+import rasterio
+from rasterio import features
+import geopandas as gpd
+from shapely.geometry import mapping
+
+def shp_to_mask(shp_path, reference_image_path, output_mask_path, burn_value=1):
+    """
+    将 SHP 多边形文件转换为与参考图像对齐的掩膜 ENVI 文件（.dat + .hdr）。
+
+    参数
+    ----------
+    shp_path : str
+        输入 SHP 文件路径（多边形要素）
+    reference_image_path : str
+        参考高光谱图像路径（用于获取空间范围、分辨率、投影）
+    output_mask_path : str
+        输出掩膜文件路径（应包含 .dat 扩展名）
+    burn_value : int, optional
+        掩膜内部像素值（默认为1，外部为0）
+    """
+    # 1. 读取 SHP 文件
+    gdf = gpd.read_file(shp_path)
+    if not gdf.geometry.geom_type.isin(['Polygon', 'MultiPolygon']).all():
+        raise ValueError("SHP 必须包含 Polygon 或 MultiPolygon 几何类型。")
+
+    # 2. 打开参考图像获取元数据
+    with rasterio.open(reference_image_path) as src:
+        ref_crs = src.crs
+        ref_transform = src.transform
+        ref_width = src.width
+        ref_height = src.height
+
+        # 3. 坐标转换（如果需要）
+        if ref_crs is not None and gdf.crs != ref_crs:
+            gdf = gdf.to_crs(ref_crs)
+
+        # 4. 生成掩膜数组
+        geometries = [mapping(geom) for geom in gdf.geometry]
+        mask = features.geometry_mask(
+            geometries,
+            out_shape=(ref_height, ref_width),
+            transform=ref_transform,
+            all_touched=False,      # 仅中心点在多边形内
+            invert=True            # True=内部, False=外部（这里取反）
+        )
+        # 转换布尔为整数掩膜（True -> burn_value, False -> 0）
+        mask = mask.astype(np.uint8) * burn_value
+
+        # 5. 写入 ENVI 格式（.dat + .hdr）
+        # 继承参考图像的元数据，修改为单波段 ENVI 格式
+        profile = src.profile
+        profile.update(
+            driver='ENVI',          # 输出格式
+            count=1,
+            dtype='uint8',
+            compress=None          # ENVI 不支持 LZW 压缩
+        )
+        # 确保输出文件扩展名为 .dat
+        if not output_mask_path.lower().endswith('.dat'):
+            output_mask_path += '.dat'
+        with rasterio.open(output_mask_path, 'w', **profile) as dst:
+            dst.write(mask, 1)
+
+    print(f"掩膜已保存至: {output_mask_path} (及同名的 .hdr 头文件)")
+
+# ==================== 使用示例 ====================
+if __name__ == "__main__":
+    shp_file = r"E:\is2\yaopu\925\roi.shp"
+    reference_img = r"E:\is2\yaopu\output\before.dat"   # 参考图像（可以是任意一幅高光谱图像）
+    output_mask = r"E:\is2\yaopu\925\roi_mask.dat"       # 输出的掩膜文件
+    shp_to_mask(shp_file, reference_img, output_mask, burn_value=1)

ocbrdf/src/跨轨变化规律.py

@@ -0,0 +1,217 @@
+import os
+import numpy as np
+import matplotlib.pyplot as plt
+import rasterio
+from rasterio.windows import from_bounds
+import warnings
+warnings.filterwarnings("ignore", category=RuntimeWarning)
+
+
+def extract_mask_values(image_path, mask_src, band_index=0):
+    """
+    从图像中提取与掩膜重叠区域内掩膜值为1的像素值，按行主序返回。
+
+    参数
+    ----------
+    image_path : str
+        图像文件路径（需包含 .hdr 头文件）
+    mask_src : rasterio.DatasetReader
+        已打开的掩膜文件（用于读取掩膜数组和变换）
+    band_index : int
+        波段索引（0‑based）
+
+    返回
+    -------
+    list of float
+        按行主序排列的像素值
+    """
+    with rasterio.open(image_path) as img_src:
+        # 获取图像与掩膜的交集窗口（以图像坐标系为准）
+        # 计算重叠区域的边界（在掩膜投影坐标系下）
+        img_bounds = img_src.bounds
+        mask_bounds = mask_src.bounds
+        overlap_left = max(img_bounds.left, mask_bounds.left)
+        overlap_bottom = max(img_bounds.bottom, mask_bounds.bottom)
+        overlap_right = min(img_bounds.right, mask_bounds.right)
+        overlap_top = min(img_bounds.top, mask_bounds.top)
+        if overlap_left >= overlap_right or overlap_bottom >= overlap_top:
+            return []   # 无重叠
+
+        # 将重叠边界转换到图像像素窗口
+        img_window = from_bounds(overlap_left, overlap_bottom, overlap_right, overlap_top, img_src.transform)
+        # 将重叠边界转换到掩膜像素窗口
+        mask_window = from_bounds(overlap_left, overlap_bottom, overlap_right, overlap_top, mask_src.transform)
+
+        # 读取图像的重叠区域波段数据
+        img_data = img_src.read(band_index + 1, window=img_window)
+        # 读取掩膜的重叠区域（单波段）
+        mask_data = mask_src.read(1, window=mask_window)
+
+        # 提取掩膜值为1的像素位置（在重叠窗口内）
+        rows, cols = np.where(mask_data == 1)
+        if len(rows) == 0:
+            return []
+
+        # 按行主序排序（先按行，再按列）
+        idx = np.lexsort((cols, rows))
+        sorted_rows = rows[idx]
+        sorted_cols = cols[idx]
+
+        # 从图像数据中获取对应像素值
+        values = img_data[sorted_rows, sorted_cols].tolist()
+        return values
+
+
+def fit_segment(x, y, degree):
+    """对一段数据执行多项式拟合，返回多项式函数和拟合点（用于绘图）"""
+    if len(x) < degree + 1:
+        # 有效点数不足，降低阶数
+        degree = max(1, len(x) - 1)
+        if degree < 1:
+            return None, None, None
+    coeffs = np.polyfit(x, y, degree)
+    poly = np.poly1d(coeffs)
+    x_fine = np.linspace(x[0], x[-1], 200)
+    y_fine = poly(x_fine)
+    return poly, x_fine, y_fine
+
+
+def main():
+    # ==================== 用户可修改的参数 ====================
+    folder = r"D:\BaiduNetdiskDownload\20250902\_3_52_52\316\R\test"   # 存放BSQ图像的文件夹
+    mask_path = r"E:\is2\yaopu\925\roi_mask.dat"                 # 单波段掩膜 DAT 文件（值为1的区域为ROI）
+    band_idx = 10                                               # 要提取的波段索引（0‑based）
+    poly_degree = 2                                               # 多项式拟合阶数（全局，各航带独立拟合）
+    # ========================================================
+
+    # 打开掩膜文件（保持打开状态以便重复读取）
+    if not os.path.exists(mask_path):
+        print(f"掩膜文件不存在: {mask_path}")
+        return
+    mask_src = rasterio.open(mask_path)
+
+    # 获取文件夹内所有高光谱文件（支持 .dat, .img, .bsq）
+    files = [os.path.join(folder, f) for f in os.listdir(folder)
+             if f.lower().endswith(('.dat', '.img', '.bsq'))]
+    files.sort()
+
+    # 存储全局信息
+    all_filtered_values = []      # 所有有效像素值（按顺序）
+    segment_info = []             # 每个航带的 (全局起始索引, 全局结束索引, 有效像素值数组, 拟合结果)
+    cum_len = 0                   # 累积的有效像素数（全局x坐标）
+
+    for i, f in enumerate(files):
+        print(f"Processing {f}...")
+        values = extract_mask_values(f, mask_src, band_idx)
+        if values:
+            # 去除零值，保留顺序
+            seg_vals = [v for v in values if v != 0]
+            n_seg = len(seg_vals)
+            if n_seg > 0:
+                start_idx = cum_len + 1
+                end_idx = cum_len + n_seg
+                cum_len += n_seg
+                all_filtered_values.extend(seg_vals)
+
+                # 对该段进行拟合（局部坐标1..n_seg）
+                x_local = np.arange(1, n_seg + 1)
+                y_local = np.array(seg_vals)
+                # 调用拟合函数
+                poly, x_fine_local, y_fine_local = fit_segment(x_local, y_local, poly_degree)
+                # 将局部拟合坐标转换为全局坐标
+                if x_fine_local is not None:
+                    x_fine_global = x_fine_local + start_idx - 1
+                    segment_info.append({
+                        'start_idx': start_idx,
+                        'end_idx': end_idx,
+                        'seg_vals': seg_vals,
+                        'poly': poly,
+                        'x_fine_global': x_fine_global,
+                        'y_fine': y_fine_local,
+                        'n_vals': n_seg
+                    })
+                else:
+                    # 拟合失败（点数太少），仍记录段信息但无拟合曲线
+                    segment_info.append({
+                        'start_idx': start_idx,
+                        'end_idx': end_idx,
+                        'seg_vals': seg_vals,
+                        'poly': None,
+                        'x_fine_global': None,
+                        'y_fine': None,
+                        'n_vals': n_seg
+                    })
+            else:
+                # 该航带无有效像素
+                print(f"  -> 航带 {i+1} 无有效像素（所有值为0）")
+                segment_info.append(None)   # 占位，便于后续处理（可选）
+        else:
+            print(f"  -> 无重叠或掩膜内无像素")
+            segment_info.append(None)
+
+    mask_src.close()
+
+    if not all_filtered_values:
+        print("没有有效的采样值。")
+        return
+
+    n_total = len(all_filtered_values)
+    x_global = np.arange(1, n_total + 1)
+
+    # 绘制图形
+    plt.figure(figsize=(14, 6))
+
+    # 绘制原始数据折线（所有有效像素）
+    plt.plot(x_global, all_filtered_values, 'b-', linewidth=1, label='Measured values (zeros removed)')
+
+    # 绘制各航带拟合曲线（不同颜色）
+    colors = plt.cm.tab10(np.linspace(0, 1, len([seg for seg in segment_info if seg is not None])))
+    color_idx = 0
+    for seg in segment_info:
+        if seg is None or seg['poly'] is None:
+            continue
+        # 使用不同颜色绘制拟合曲线
+        plt.plot(seg['x_fine_global'], seg['y_fine'], '--', linewidth=2, color=colors[color_idx],
+                 label=f'Segment {color_idx+1} fit (deg={poly_degree})')
+        color_idx += 1
+
+    # 添加图例（自动合并重复项，这里不合并）
+    plt.legend()
+
+    plt.xlabel("Pixel index after removing zeros (row-major order, cumulative across scenes)")
+    plt.ylabel(f"Band {band_idx + 1} value")
+    plt.title("Spectral variation across water pixels in ROI (per‑segment polynomial fit)")
+    plt.grid(True, linestyle='--', alpha=0.5)
+
+    # 标注每个航带的分界线（基于有效像素累积）
+    ylim = plt.ylim()
+    for seg in segment_info:
+        if seg is None:
+            continue
+        end_idx = seg['end_idx']
+        n_vals = seg['n_vals']
+        if end_idx < n_total:   # 非最后一个段画竖线
+            plt.axvline(x=end_idx, color='gray', linestyle='--', alpha=0.7)
+            plt.text(end_idx, ylim[1] * 0.95, f"{n_vals}", ha='center', fontsize=8)
+        else:
+            # 最后一个段：在右边界标注数量
+            plt.text(end_idx - 0.5, ylim[1] * 0.95, f"{n_vals}", ha='right', fontsize=8)
+
+    plt.tight_layout()
+    out_path = r"E:\is2\yaopu\ocbrdf\test\roi_variation_original_10band.png"
+    plt.savefig(out_path)
+    print(f"图形已保存至: {out_path}")
+
+    # 可选：打印各段拟合系数
+    print("\n各航带拟合多项式系数（从最高次到常数项）：")
+    for i, seg in enumerate(segment_info):
+        if seg is not None and seg['poly'] is not None:
+            coeffs = seg['poly'].coeffs
+            print(f"航带 {i+1}（有效像素数 {seg['n_vals']}）：")
+            for j, c in enumerate(coeffs):
+                print(f"  x^{len(coeffs)-1-j}: {c:.6f}")
+        elif seg is not None:
+            print(f"航带 {i+1}：无有效像素或拟合失败")
+
+if __name__ == "__main__":
+    main()

ocbrdf/src/跨轨变化规律多波段.py

@@ -0,0 +1,250 @@
+import os
+import numpy as np
+import matplotlib.pyplot as plt
+import rasterio
+from rasterio.windows import from_bounds
+import warnings
+warnings.filterwarnings("ignore", category=RuntimeWarning)
+
+
+def extract_mask_values_multi_band(image_path, mask_src, band_indices):
+    """
+    从图像中提取与掩膜重叠区域内掩膜值为1的像素值（多个波段），按行主序返回。
+    返回值包含零值（保留原始顺序）。
+    """
+    with rasterio.open(image_path) as img_src:
+        img_bounds = img_src.bounds
+        mask_bounds = mask_src.bounds
+        overlap_left = max(img_bounds.left, mask_bounds.left)
+        overlap_bottom = max(img_bounds.bottom, mask_bounds.bottom)
+        overlap_right = min(img_bounds.right, mask_bounds.right)
+        overlap_top = min(img_bounds.top, mask_bounds.top)
+        if overlap_left >= overlap_right or overlap_bottom >= overlap_top:
+            return {}
+
+        img_window = from_bounds(overlap_left, overlap_bottom, overlap_right, overlap_top, img_src.transform)
+        mask_window = from_bounds(overlap_left, overlap_bottom, overlap_right, overlap_top, mask_src.transform)
+
+        img_data = img_src.read(band_indices, window=img_window)
+        mask_data = mask_src.read(1, window=mask_window)
+
+        rows, cols = np.where(mask_data == 1)
+        if len(rows) == 0:
+            return {}
+
+        idx = np.lexsort((cols, rows))
+        sorted_rows = rows[idx]
+        sorted_cols = cols[idx]
+
+        result = {}
+        for i, band in enumerate(band_indices):
+            band_data = img_data[i]
+            values = band_data[sorted_rows, sorted_cols].tolist()
+            result[band] = values
+        return result
+
+
+def fit_segment(x, y, degree):
+    """对一段数据执行多项式拟合，返回多项式函数和拟合点（用于绘图）"""
+    if len(x) < degree + 1:
+        degree = max(1, len(x) - 1)
+        if degree < 1:
+            return None, None, None
+    coeffs = np.polyfit(x, y, degree)
+    poly = np.poly1d(coeffs)
+    x_fine = np.linspace(x[0], x[-1], 200)
+    y_fine = poly(x_fine)
+    return poly, x_fine, y_fine
+
+
+def plot_single_band(band, x_global, y_original, seg_info_list, segment_sizes_nonzero,
+                     out_dir, poly_degree):
+    """绘制单个波段的原始数据和分段拟合曲线，并保存图片。"""
+    plt.figure(figsize=(14, 6))
+
+    # 原始数据
+    plt.plot(x_global, y_original, 'b-', linewidth=1.5, alpha=0.7, label='Original data')
+
+    # 分段拟合曲线，只给第一个有效段添加图例
+    labeled = False
+    for seg in seg_info_list:
+        if seg is not None:
+            label = 'Fitted curve' if not labeled else ""
+            plt.plot(seg['x_fine_global'], seg['y_fine'], 'r--', linewidth=2, label=label)
+            labeled = True
+
+    plt.xlabel("Pixel index after removing zeros (row-major order, cumulative across scenes)")
+    plt.ylabel(f"Band {band+1} value")
+    plt.title(f"Band {band+1}: Original vs. Fitted (per‑segment polynomial, deg={poly_degree})")
+    plt.grid(True, linestyle='--', alpha=0.5)
+    plt.legend()
+
+    # 标注航带分界线
+    ylim = plt.ylim()
+    cum_nz = 0
+    total_nz = len(y_original)
+    for i, nz_cnt in enumerate(segment_sizes_nonzero):
+        if nz_cnt == 0:
+            continue
+        cum_nz += nz_cnt
+        if cum_nz < total_nz:
+            plt.axvline(x=cum_nz, color='gray', linestyle='--', alpha=0.7)
+            plt.text(cum_nz, ylim[1] * 0.95, f"{nz_cnt}", ha='center', fontsize=8)
+        else:
+            plt.text(cum_nz - 0.5, ylim[1] * 0.95, f"{nz_cnt}", ha='right', fontsize=8)
+
+    plt.tight_layout()
+    out_path = os.path.join(out_dir, f"band_{band+1}_original_and_fitted.png")
+    plt.savefig(out_path)
+    plt.close()
+    print(f"  -> 波段 {band+1} 图已保存至: {out_path}")
+
+
+def main():
+    # ==================== 用户可修改的参数 ====================
+    folder = r"D:\BaiduNetdiskDownload\20250902\_3_52_52\316\R\test"
+    mask_path = r"E:\is2\yaopu\925\roi_mask.dat"
+    band_indices = [35,36,37]          # 波段 9,10,11 (0‑based)
+    poly_degree = 2
+    out_dir = r"E:\is2\yaopu\ocbrdf\test"
+    # ========================================================
+
+    os.makedirs(out_dir, exist_ok=True)
+
+    if not os.path.exists(mask_path):
+        print(f"掩膜文件不存在: {mask_path}")
+        return
+    mask_src = rasterio.open(mask_path)
+
+    files = [os.path.join(folder, f) for f in os.listdir(folder)
+             if f.lower().endswith(('.dat', '.img', '.bsq'))]
+    files.sort()
+
+    band_values_with_zeros = {band: [] for band in band_indices}
+    segment_sizes_total = []
+
+    for f in files:
+        print(f"Processing {f}...")
+        extracted = extract_mask_values_multi_band(f, mask_src, band_indices)
+        if extracted:
+            seg_len = len(extracted[band_indices[0]])
+            segment_sizes_total.append(seg_len)
+            for band in band_indices:
+                band_values_with_zeros[band].extend(extracted[band])
+        else:
+            print(f"  -> 无重叠或掩膜内无像素")
+            segment_sizes_total.append(0)
+
+    mask_src.close()
+
+    if not band_values_with_zeros[band_indices[0]]:
+        print("没有有效的采样值。")
+        return
+
+    # 去除零值
+    band_values_filtered = {}
+    nonzero_masks = {}
+    for band in band_indices:
+        vals = np.array(band_values_with_zeros[band])
+        mask_nonzero = (vals != 0)
+        band_values_filtered[band] = vals[mask_nonzero].tolist()
+        nonzero_masks[band] = mask_nonzero
+
+    # 计算每个航带有效像素数
+    segment_sizes_nonzero = []
+    cum = 0
+    for sz in segment_sizes_total:
+        seg_mask = nonzero_masks[band_indices[0]][cum:cum+sz]
+        nz_count = np.sum(seg_mask)
+        segment_sizes_nonzero.append(nz_count)
+        cum += sz
+
+    # 构建每个波段的分段拟合信息
+    segment_info_per_band = {}
+    for band in band_indices:
+        seg_info = []
+        cum_nz = 0
+        for nz_cnt in segment_sizes_nonzero:
+            if nz_cnt == 0:
+                seg_info.append(None)
+                continue
+            start_nz = cum_nz
+            end_nz = cum_nz + nz_cnt
+            seg_vals = band_values_filtered[band][start_nz:end_nz]
+            x_local = np.arange(1, nz_cnt + 1)
+            y_local = np.array(seg_vals)
+            poly, x_fine_local, y_fine_local = fit_segment(x_local, y_local, poly_degree)
+            if x_fine_local is not None:
+                x_fine_global = x_fine_local + cum_nz
+                seg_info.append({
+                    'x_fine_global': x_fine_global,
+                    'y_fine': y_fine_local,
+                    'nz_count': nz_cnt,
+                    'start_idx': cum_nz + 1,
+                    'end_idx': cum_nz + nz_cnt
+                })
+            else:
+                seg_info.append(None)
+            cum_nz += nz_cnt
+        segment_info_per_band[band] = seg_info
+
+    # 1. 绘制三个波段拟合曲线对比图
+    plt.figure(figsize=(14, 6))
+    colors = ['blue', 'green', 'red']
+    labels = ['Band 9', 'Band 10', 'Band 11']
+    for idx, band in enumerate(band_indices):
+        seg_info = segment_info_per_band[band]
+        labeled = False
+        for seg in seg_info:
+            if seg is not None:
+                label = labels[idx] if not labeled else ""
+                plt.plot(seg['x_fine_global'], seg['y_fine'],
+                         color=colors[idx], linestyle='--', linewidth=2, label=label)
+                labeled = True
+    plt.xlabel("Pixel index after removing zeros (row-major order, cumulative across scenes)")
+    plt.ylabel("Pixel value")
+    plt.title(f"Fitted variation for bands 9, 10, 11 (per-segment polynomial fit, deg={poly_degree})")
+    plt.grid(True, linestyle='--', alpha=0.5)
+    plt.legend()
+    ylim = plt.ylim()
+    cum_nz = 0
+    total_nz = len(band_values_filtered[band_indices[0]])
+    for i, nz_cnt in enumerate(segment_sizes_nonzero):
+        if nz_cnt == 0:
+            continue
+        cum_nz += nz_cnt
+        if cum_nz < total_nz:
+            plt.axvline(x=cum_nz, color='gray', linestyle='--', alpha=0.7)
+            plt.text(cum_nz, ylim[1] * 0.95, f"{nz_cnt}", ha='center', fontsize=8)
+        else:
+            plt.text(cum_nz - 0.5, ylim[1] * 0.95, f"{nz_cnt}", ha='right', fontsize=8)
+    plt.tight_layout()
+    combined_path = os.path.join(out_dir, "roi_fitted_bands_9_10_11_zeros_removed.png")
+    plt.savefig(combined_path)
+    plt.close()
+    print(f"组合对比图已保存至: {combined_path}")
+
+    # 2. 为每个波段单独绘制原始数据+拟合曲线图
+    for band in band_indices:
+        x_global = np.arange(1, len(band_values_filtered[band]) + 1)
+        y_original = np.array(band_values_filtered[band])
+        seg_info = segment_info_per_band[band]
+        plot_single_band(band, x_global, y_original, seg_info,
+                         segment_sizes_nonzero, out_dir, poly_degree)
+
+    # 打印拟合系数摘要
+    print("\n各航带拟合多项式系数（仅显示每个波段第一个有拟合的航带示例）")
+    for band in band_indices:
+        print(f"\n波段 {band+1}:")
+        for i, seg in enumerate(segment_info_per_band[band]):
+            if seg is not None:
+                nz_cnt = seg['nz_count']
+                seg_vals = band_values_filtered[band][seg['start_idx']-1:seg['end_idx']]
+                x_local = np.arange(1, nz_cnt+1)
+                coeffs = np.polyfit(x_local, seg_vals, poly_degree)
+                print(f"  航带 {i+1}: {coeffs}")
+                break
+
+
+if __name__ == "__main__":
+    main()

ocbrdf/test3.py

@@ -0,0 +1,609 @@
+import gc
+import os
+import numpy as np
+import sys
+import xarray as xr
+import rasterio
+import spectral
+
+from brdf_model_M02 import M02
+from brdf_model_M02SeaDAS import M02SeaDAS
+from brdf_model_L11 import L11
+from brdf_model_O25 import O25
+from brdf_utils import ADF_OCP, squeeze_trivial_dims
+
+# 尝试导入 tqdm 进度条库，如果不可用则设为 None
+try:
+    from tqdm import tqdm
+except ImportError:
+    tqdm = None
+
+UNC_LUT_CACHE = {}
+
+"""
+主 BRDF 校正模块
+输入为 xarray 数据集
+所需的光谱维度为 "bands"，其他维度自由
+输入数据集中必填的字段：
+Rw: 方向性海洋反射率
+sza: 太阳天顶角
+vza: 观测天顶角
+raa: 相对方位角（当太阳和观测在同一侧时，raa=0）
+输入数据集中可选的字段：
+Rw_unc: Rw 的不确定度（若缺失，设为零）
+输出数据集中的字段：
+nrrs: 完全归一化的遥感反射率
+rho_ex_w: nrrs * PI
+omega_b: bb/(a+bb)
+eta_b: bbw/bb
+C_brdf: BRDF 校正因子
+brdf_unc: C_brdf 的不确定度
+nrrs_unc: nrrs 的不确定度
+
+使用和适配 brdf_hypercp 模块需保留的信息：
+Brdf_hypercp 是 EUMETSAT 研究项目的一部分
+"BRDF correction of S3 OLCI water reflectance products"（S3 OLCI 水体反射率产品的 BRDF 校正），
+合同号：RB_EUM-CO-21-4600002626-JIG。
+研究团队成员：Davide D'Alimonte (davide.dalimonte@aequora.org)，
+Tamito Kajiyama (tamito.kajiyama@aequora.org)，
+Jaime Pitarch (jaime.pitarchportero@artov.ismar.cnr.it)，
+Vittorio Brando (vittorio.brando@cnr.it)，
+Marco Talone (talone@icm.csic.es) 以及
+Constant Mazeran (constant.mazeran@solvo.fr)。
+BRDF 查找表中的相对方位角遵循 OLCI 约定。详见 https://www.eumetsat.int/media/50720，图 6。
+"""
+
+
+# ---------------------------------------------------------------------------
+# File I/O helpers
+# ---------------------------------------------------------------------------
+
+def _resolve_envi_header(filepath):
+    """Return the ENVI header path accepted by spectral."""
+    if filepath.lower().endswith(".hdr"):
+        return filepath
+
+    hdr_path = filepath + ".hdr"
+    if os.path.exists(hdr_path):
+        return hdr_path
+
+    base, _ext = os.path.splitext(filepath)
+    hdr_path = base + ".hdr"
+    if os.path.exists(hdr_path):
+        return hdr_path
+
+    return filepath
+
+
+def read_bsq(filepath, wavelengths=None, dtype=np.float32):
+    """Read an ENVI BSQ hyperspectral image.
+
+    Parameters
+    ----------
+    filepath : str
+        Path to the ENVI header (.hdr) or image file.
+    wavelengths : array-like, optional
+        Band centre wavelengths in nm.  When *None*, wavelengths are read
+        from the ENVI header if available.
+
+    Returns
+    -------
+    data : ndarray, shape (bands, rows, cols), float64
+    wavelengths : ndarray, shape (bands,)
+    img : spectral image object (metadata / transform access)
+    """
+    img = spectral.open_image(_resolve_envi_header(filepath))
+    data = np.asarray(img.open_memmap(interleave="source"), dtype=dtype)
+    # data = np.transpose(data, (2, 0, 1)) # -> (bands, rows, cols)
+
+    if wavelengths is None:
+        if hasattr(img, 'bands') and img.bands.centers is not None:
+            wavelengths = np.array(img.bands.centers, dtype=np.float32)
+        else:
+            n_bands = data.shape[0]
+            wavelengths = np.arange(1, n_bands + 1, dtype=np.float32)
+
+    return data, np.asarray(wavelengths, dtype=dtype), img
+
+
+def read_bip_angles(filepath, dtype=np.float32):
+    """Read an ENVI BIP angle file and extract the five geometry bands.
+
+    Band layout (1-indexed):
+        band  7  -> SZA  Solar Zenith Angle
+        band  8  -> SAA  Solar Azimuth Angle
+        band  9  -> VZA  Sensor Zenith Angle
+        band 10  -> VAA  Sensor Azimuth Angle
+        band 11  -> RAA  Relative Azimuth Angle
+
+    Parameters
+    ----------
+    filepath : str
+        Path to the ENVI BIP file (header or image file).
+
+    Returns
+    -------
+    dict with keys 'sza', 'saa', 'vza', 'vaa', 'raa'.
+    Each value is a 2-D ndarray of shape (rows, cols), float64.
+    """
+    img = spectral.open_image(_resolve_envi_header(filepath))
+    data = np.asarray(img.open_memmap(interleave="source"), dtype=dtype)
+
+    n_bands = data.shape[2]
+    if n_bands < 11:
+        raise ValueError(
+            f"Angle file has only {n_bands} bands; at least 11 are required "
+            "(bands 7–11 = SZA, SAA, VZA, VAA, RAA)."
+        )
+
+    return {
+        'sza': data[:, :, 6].copy(),
+        'saa': data[:, :, 7].copy(),
+        'vza': data[:, :, 8].copy(),
+        'vaa': data[:, :, 9].copy(),
+        'raa': data[:, :, 10].copy(),
+    }
+
+
+def read_water_mask(filepath):
+    """Read a single-band GeoTIFF water mask.
+
+    Parameters
+    ----------
+    filepath : str
+        Path to the GeoTIFF file.
+
+    Returns
+    -------
+    mask : ndarray, shape (rows, cols), int8
+        Pixels with value 1 indicate water; all other values are non-water.
+    """
+    with rasterio.open(filepath) as src:
+        mask = src.read(1).astype(np.int8)
+    return mask
+
+
+def save_brdf_result(
+    ds_out,
+    output_file,
+    source_file=None,
+    output_var="Rw_brdf",
+    output_format="ENVI",
+    dtype=np.float32,
+):
+    """Save the BRDF result to ENVI or NetCDF."""
+    output_format = output_format.upper()
+
+    if output_var not in ds_out:
+        raise KeyError(f"Output variable '{output_var}' not found in result dataset.")
+
+    if output_format == "NETCDF":
+        if not output_file.lower().endswith(".nc"):
+            output_file = output_file + ".nc"
+        ds_out.to_netcdf(output_file)
+        return output_file
+
+    if output_format != "ENVI":
+        raise ValueError("output_format must be 'ENVI' or 'NETCDF'.")
+
+    cube = np.asarray(ds_out[output_var].values, dtype=dtype)
+    cube = np.transpose(cube, (1, 2, 0))
+
+    if output_file.lower().endswith(".hdr"):
+        hdr_path = output_file
+    else:
+        hdr_path = output_file + ".hdr"
+
+    metadata = {
+        "interleave": "bsq",
+        "description": f"BRDF corrected result: {output_var}",
+        "bands": cube.shape[2],
+        "lines": cube.shape[0],
+        "samples": cube.shape[1],
+    }
+    if "bands" in ds_out.coords:
+        metadata["wavelength"] = [str(v) for v in ds_out.coords["bands"].values.tolist()]
+
+    if source_file is not None:
+        metadata["source_file"] = source_file
+
+    spectral.envi.save_image(
+        hdr_path,
+        cube,
+        dtype=dtype,
+        interleave="bsq",
+        metadata=metadata,
+        force=True,
+    )
+    return hdr_path
+
+
+# ---------------------------------------------------------------------------
+# Top-level pipeline
+# ---------------------------------------------------------------------------
+def brdf_uncertainty(ds, adf=None):
+    ''' Compute uncertainty of BRDF factor and propagate to nrrs '''
+
+    # Read LUT
+    if adf is None:
+        adf = ADF_OCP
+
+    unc_lut_path = adf % 'UNC'
+    if unc_lut_path not in UNC_LUT_CACHE:
+        UNC_LUT_CACHE[unc_lut_path] = xr.open_dataset(unc_lut_path, engine='netcdf4')
+    LUT = UNC_LUT_CACHE[unc_lut_path]
+
+    # Interpolate relative uncertainty
+    unc = LUT['unc'].interp(lambda_unc=ds.bands, theta_s_unc=ds.sza, theta_v_unc=ds.vza,
+                            delta_phi_unc=ds.raa)
+
+    # Compute absolute uncertainty of factor
+    ds['brdf_unc'] = unc * ds['C_brdf']
+
+    # Flag BRDF_unc where NaN and set to 0
+    ds['brdf_unc_fail'] = np.isnan(ds['brdf_unc'])
+    ds['brdf_unc'] = xr.where(ds['brdf_unc_fail'], 0, ds['brdf_unc'])
+
+    # Propagate to nrrs
+    Rwex_unc2 = ds['brdf_unc'] * ds['brdf_unc'] * ds['Rw'] * ds['Rw']
+    if 'Rw_unc' in ds:
+        Rwex_unc2 += ds['C_brdf'] * ds['C_brdf'] * ds['Rw_unc'] * ds['Rw_unc']
+    ds['nrrs_unc'] = np.sqrt(Rwex_unc2)/np.pi
+
+    return ds
+
+
+def brdf_prototype(ds, adf=None, brdf_model='L11', compute_uncertainty=True):
+    # 测试 BRDF 模型在 GUI 中不支持：强制覆盖
+    # brdf_model = 'M02SeaDAS'
+
+    # 压缩单一维度（例如，投放次数、提取次数等）以避免插值问题
+    ds, squeezedDims = squeeze_trivial_dims(ds)
+
+    # Initialise model
+    if brdf_model == 'M02':
+        BRDF_model = M02(bands=ds.bands, aot=ds.aot, wind=ds.wind, adf=None)  # Don't use brdf_py.ADF context
+    elif brdf_model == 'M02SeaDAS':
+        BRDF_model = M02SeaDAS(bands=ds.bands, adf=None)  # Don't use brdf_py.ADF context
+    elif brdf_model == 'L11':
+        BRDF_model = L11(bands=ds.bands, adf=None)  # Don't use brdf_py.ADF context
+    elif brdf_model == 'O25':
+        BRDF_model = O25(bands=ds.bands, adf=None)  # Don't use brdf_py.ADF context
+    else:
+        print("BRDF model %s not supported" % brdf_model)
+        sys.exit(1)
+
+    # Init pixel
+    BRDF_model.init_pixels(ds['sza'], ds['vza'], ds['raa'])
+
+    # Compute IOP and normalize by iterating
+    ds['nrrs'] = ds['Rw'] / np.pi
+
+    ds['convergeFlag'] = (0 * ds['sza']).astype(bool)
+    ds['C_brdf'] = 0 * ds['nrrs'] + 1
+
+    for iter_brdf in range(int(BRDF_model.niter)):
+
+        # M02: Initialise chl_iter
+        if brdf_model in ['M02', 'M02SeaDAS'] and (iter_brdf == 0):
+            chl_iter = {}
+            ds['log10_chl'] = 0 * ds['sza'] + float(np.log10(BRDF_model.OC4MEchl0))
+            chl_iter[-1]    = 0 * ds['sza'] + float(BRDF_model.OC4MEchl0)
+
+        ds = BRDF_model.backward(ds, iter_brdf)
+
+        # M02: Check convergence (dummy for M02SeaDAS for the moment... epsilon set to 0)
+        if brdf_model in ['M02', 'M02SeaDAS']:
+            chl_iter[iter_brdf] = 10 ** ds['log10_chl']
+            #  Check if convergence is reached |chl_old-chl_new| < epsilon * chl_new
+            ds['convergeFlag'] = (ds['convergeFlag']) | (
+                (np.abs(chl_iter[iter_brdf - 1] - chl_iter[iter_brdf]) < float(BRDF_model.OC4MEepsilon) * chl_iter[
+                    iter_brdf]))
+
+        # Apply forward model in both geometries
+        forward_mod = BRDF_model.forward(ds)
+        forward_mod0 = BRDF_model.forward(ds, normalized=True)
+
+        ratio = forward_mod0 / forward_mod
+
+        # Drop remnant coordinates to avoid ambiguities in the update of the BRDF factor.
+        for coord in ratio.coords:
+            if coord not in ds['C_brdf'].coords:
+                ratio = ratio.drop(coord)
+
+        # Normalize reflectance
+        ds['C_brdf'] = xr.where(ds['convergeFlag'], ds['C_brdf'], ratio)
+        ds['nrrs'] = ds['Rw'] / np.pi * ds['C_brdf']
+
+    # Flag BRDF where NaN and set to 1 (no correction applied).
+    ds['C_brdf_fail'] = np.isnan(ds['C_brdf'])
+    ds['C_brdf'] = xr.where(ds['C_brdf_fail'], 1, ds['C_brdf'])
+    ds['nrrs'] = xr.where(ds['C_brdf_fail'], ds['Rw'] / np.pi, ds['nrrs'])
+
+    # If QAA_fail is raised, raise C_brdf_fail (but still apply C_brdf).
+    if 'QAA_fail' in ds:
+        ds['C_brdf_fail'] = (ds['C_brdf_fail']) | (ds['QAA_fail'])
+
+    # Compute uncertainty only if requested
+    if compute_uncertainty:
+        ds = brdf_uncertainty(ds)
+
+    # Compute flag
+    ds['flags_level2'] = ds['Rw'] * 0  # TODO
+
+    # Convert to reflectance unit
+    ds['rho_ex_w'] = ds['nrrs'] * np.pi
+
+    # Expand squeezed trivial dimensions
+    for dim,d0 in squeezedDims.items():
+        ds = ds.expand_dims(dim,axis=d0)
+
+    return ds
+
+
+def run_brdf_correction(
+    hyperspectral_file,
+    angle_file,
+    mask_file,
+    wavelengths=None,
+    brdf_model='L11',
+    chunk_size=4096,
+    output_var=None,          # if None, keep full output (original behaviour)
+    compute_uncertainty=True, # compute uncertainty only if requested
+    progress_bar=True,        # 是否显示进度条
+):
+    """BRDF correction pipeline for water pixels in a hyperspectral scene.
+
+    ... (docstring unchanged, update description) ...
+    """
+    # ------------------------------------------------------------------
+    # 1. Read all inputs
+    # ------------------------------------------------------------------
+    rw_data, wvl, _img_meta = read_bsq(hyperspectral_file, wavelengths, dtype=np.float32)
+    angles = read_bip_angles(angle_file, dtype=np.float32)
+    water_mask = read_water_mask(mask_file)
+
+    n_bands, n_rows, n_cols = rw_data.shape
+
+    # Sanity-check spatial dimensions
+    if angles['sza'].shape != (n_rows, n_cols):
+        raise ValueError(
+            f"Angle file spatial shape {angles['sza'].shape} does not match "
+            f"hyperspectral image shape ({n_rows}, {n_cols})."
+        )
+    if water_mask.shape != (n_rows, n_cols):
+        raise ValueError(
+            f"Mask shape {water_mask.shape} does not match "
+            f"hyperspectral image shape ({n_rows}, {n_cols})."
+        )
+
+    # ------------------------------------------------------------------
+    # 2. Identify water pixels
+    # ------------------------------------------------------------------
+    water_rows, water_cols = np.where(water_mask == 1)
+    n_water = water_rows.size
+    if n_water == 0:
+        raise RuntimeError("No water pixels found in the mask (no pixels == 1).")
+
+    # ------------------------------------------------------------------
+    # 3. Determine which output variable(s) we need to keep
+    # ------------------------------------------------------------------
+    # If output_var is None, we keep everything (original behaviour)
+    keep_all = (output_var is None)
+
+    # Map output_var to the field in ds_corr
+    var_map = {
+        'Rw_brdf': 'rho_ex_w',
+        'rho_ex_w': 'rho_ex_w',
+        'nrrs': 'nrrs',
+        'C_brdf': 'C_brdf'
+    }
+    if not keep_all:
+        if output_var not in var_map:
+            raise ValueError(f"Unknown output_var '{output_var}'. Allowed: {list(var_map.keys())}")
+        field_name = var_map[output_var]
+    else:
+        field_name = None
+
+    # ------------------------------------------------------------------
+    # 4. Allocate output arrays (only those needed)
+    # ------------------------------------------------------------------
+    shape3 = (n_bands, n_rows, n_cols)
+    shape2 = (n_rows, n_cols)
+
+    if keep_all:
+        # Original behaviour: allocate all arrays
+        nrrs_full        = np.full(shape3, np.nan, dtype=np.float32)
+        rho_ex_w_full    = np.full(shape3, np.nan, dtype=np.float32)
+        C_brdf_full      = np.full(shape3, np.nan, dtype=np.float32)
+        brdf_unc_full    = np.full(shape3, np.nan, dtype=np.float32)
+        nrrs_unc_full    = np.full(shape3, np.nan, dtype=np.float32)
+        C_brdf_fail_full = np.zeros(shape2, dtype=bool)
+        rw_brdf_full     = rw_data.copy()
+    else:
+        # Only allocate the requested output array
+        output_full = np.full(shape3, np.nan, dtype=np.float32)
+        # We may also need C_brdf_fail if requested? Not for output_var, but could be used internally? Skip.
+        # For simplicity, we do not allocate extra arrays.
+
+    # ------------------------------------------------------------------
+    # 5–6. Chunked BRDF correction
+    # ------------------------------------------------------------------
+    n_chunks = int(np.ceil(n_water / chunk_size))
+
+    # 初始化进度条
+    iterator = range(n_chunks)
+    if progress_bar and tqdm is not None:
+        iterator = tqdm(iterator, desc="Processing chunks", unit="chunk")
+    elif progress_bar:
+        print("Warning: tqdm not installed, progress bar disabled.")
+
+    for chunk_idx in iterator:
+        lo = chunk_idx * chunk_size
+        hi = min(lo + chunk_size, n_water)
+
+        rows_c = water_rows[lo:hi]
+        cols_c = water_cols[lo:hi]
+
+        # Extract batch vectors
+        rw_c  = rw_data[:, rows_c, cols_c].T          # (chunk, bands)
+        sza_c = angles['sza'][rows_c, cols_c]          # (chunk,)
+        vza_c = angles['vza'][rows_c, cols_c]
+        raa_c = angles['raa'][rows_c, cols_c]
+
+        pixel_idx_c = np.arange(hi - lo)
+        ds_chunk = xr.Dataset(
+            {
+                'Rw':  xr.DataArray(rw_c,  dims=('n', 'bands'),
+                                    coords={'n': pixel_idx_c, 'bands': wvl}),
+                'sza': xr.DataArray(sza_c, dims='n',
+                                    coords={'n': pixel_idx_c}),
+                'vza': xr.DataArray(vza_c, dims='n',
+                                    coords={'n': pixel_idx_c}),
+                'raa': xr.DataArray(raa_c, dims='n',
+                                    coords={'n': pixel_idx_c}),
+            }
+        )
+
+        # Run BRDF model with uncertainty control
+        ds_corr = brdf_prototype(ds_chunk, brdf_model=brdf_model, compute_uncertainty=compute_uncertainty)
+
+        if keep_all:
+            # Write all results back
+            nrrs_full[:, rows_c, cols_c]        = ds_corr['nrrs'].values.T
+            rho_ex_w_full[:, rows_c, cols_c]    = ds_corr['rho_ex_w'].values.T
+            C_brdf_full[:, rows_c, cols_c]      = ds_corr['C_brdf'].values.T
+            if compute_uncertainty:
+                brdf_unc_full[:, rows_c, cols_c]    = ds_corr['brdf_unc'].values.T
+                nrrs_unc_full[:, rows_c, cols_c]    = ds_corr['nrrs_unc'].values.T
+            C_brdf_fail_full[rows_c, cols_c]    = ds_corr['C_brdf_fail'].values
+            rw_brdf_full[:, rows_c, cols_c]     = ds_corr['rho_ex_w'].values.T
+        else:
+            # Only write the requested output field
+            values = ds_corr[field_name].values.T
+            output_full[:, rows_c, cols_c] = values
+
+        # Release batch objects so GC can reclaim memory
+        del rw_c, sza_c, vza_c, raa_c, pixel_idx_c, ds_chunk, ds_corr
+        gc.collect()
+
+    # ------------------------------------------------------------------
+    # 7. Assemble full-scene output Dataset
+    # ------------------------------------------------------------------
+    full_coords = {'bands': wvl, 'y': np.arange(n_rows), 'x': np.arange(n_cols)}
+    dims3 = ('bands', 'y', 'x')
+    dims2 = ('y', 'x')
+
+    if keep_all:
+        # Original behaviour: return all variables
+        ds_out = xr.Dataset(
+            {
+                'Rw':          xr.DataArray(rw_data,          dims=dims3, coords=full_coords),
+                'Rw_brdf':     xr.DataArray(rw_brdf_full,     dims=dims3, coords=full_coords),
+                'nrrs':        xr.DataArray(nrrs_full,        dims=dims3, coords=full_coords),
+                'rho_ex_w':    xr.DataArray(rho_ex_w_full,    dims=dims3, coords=full_coords),
+                'C_brdf':      xr.DataArray(C_brdf_full,      dims=dims3, coords=full_coords),
+                'C_brdf_fail': xr.DataArray(C_brdf_fail_full, dims=dims2),
+                'water_mask':  xr.DataArray(water_mask,       dims=dims2),
+                'sza':         xr.DataArray(angles['sza'],    dims=dims2),
+                'saa':         xr.DataArray(angles['saa'],    dims=dims2),
+                'vza':         xr.DataArray(angles['vza'],    dims=dims2),
+                'vaa':         xr.DataArray(angles['vaa'],    dims=dims2),
+                'raa':         xr.DataArray(angles['raa'],    dims=dims2),
+            }
+        )
+        if compute_uncertainty:
+            ds_out['brdf_unc'] = xr.DataArray(brdf_unc_full, dims=dims3, coords=full_coords)
+            ds_out['nrrs_unc'] = xr.DataArray(nrrs_unc_full, dims=dims3, coords=full_coords)
+    else:
+        # Return only the requested output variable + geometry + mask
+        ds_out = xr.Dataset(
+            {
+                output_var:    xr.DataArray(output_full,      dims=dims3, coords=full_coords),
+                'Rw':          xr.DataArray(rw_data,          dims=dims3, coords=full_coords),
+                'water_mask':  xr.DataArray(water_mask,       dims=dims2),
+                'sza':         xr.DataArray(angles['sza'],    dims=dims2),
+                'saa':         xr.DataArray(angles['saa'],    dims=dims2),
+                'vza':         xr.DataArray(angles['vza'],    dims=dims2),
+                'vaa':         xr.DataArray(angles['vaa'],    dims=dims2),
+                'raa':         xr.DataArray(angles['raa'],    dims=dims2),
+            }
+        )
+        # If uncertainties were computed, we could add them optionally, but user didn't request them.
+        # For memory efficiency, we skip them here.
+
+    return ds_out
+
+
+def process_brdf_files(
+    hyperspectral_file,
+    angle_file,
+    mask_file,
+    output_file,
+    wavelengths=None,
+    brdf_model="L11",
+    output_var="Rw_brdf",
+    output_format="ENVI",
+    chunk_size=4096,
+    compute_uncertainty=False,   # new parameter, default False to save time
+    progress_bar=True,            # 是否显示进度条
+):
+    """Run BRDF correction from three input files and save the result."""
+    ds_out = run_brdf_correction(
+        hyperspectral_file=hyperspectral_file,
+        angle_file=angle_file,
+        mask_file=mask_file,
+        wavelengths=wavelengths,
+        brdf_model=brdf_model,
+        chunk_size=chunk_size,
+        output_var=output_var,                # pass output_var to limit output
+        compute_uncertainty=compute_uncertainty,
+        progress_bar=progress_bar,            # 传递进度条标志
+    )
+    saved_path = save_brdf_result(
+        ds_out=ds_out,
+        output_file=output_file,
+        source_file=hyperspectral_file,
+        output_var=output_var,
+        output_format=output_format,
+    )
+    return ds_out, saved_path
+
+
+def main():
+    import argparse
+
+    parser = argparse.ArgumentParser(description="Water-region BRDF correction for hyperspectral data.")
+    parser.add_argument("hyperspectral_file", help="Input ENVI BSQ hyperspectral file or .hdr path.")
+    parser.add_argument("angle_file", help="Input ENVI BIP angle file or .hdr path.")
+    parser.add_argument("mask_file", help="Input water mask GeoTIFF; pixels == 1 are corrected.")
+    parser.add_argument("output_file", help="Output path prefix. ENVI writes .hdr/.img, NetCDF writes .nc.")
+    parser.add_argument("--brdf-model", default="L11", choices=["L11", "M02", "M02SeaDAS", "O25"])
+    parser.add_argument("--output-var", default="Rw_brdf", choices=["Rw_brdf", "rho_ex_w", "nrrs", "C_brdf"])
+    parser.add_argument("--output-format", default="ENVI", choices=["ENVI", "NETCDF"])
+    parser.add_argument("--chunk-size", type=int, default=4096,
+                        help="Number of water pixels per processing batch (default: 4096).")
+    parser.add_argument("--compute-uncertainty", action="store_true",
+                        help="Compute BRDF uncertainty (may increase runtime and output size).")
+    parser.add_argument("--no-progress-bar", action="store_true",
+                        help="Disable progress bar (useful for non-interactive environments).")
+    args = parser.parse_args()
+
+    progress_bar = not args.no_progress_bar
+
+    _ds_out, saved_path = process_brdf_files(
+        hyperspectral_file=args.hyperspectral_file,
+        angle_file=args.angle_file,
+        mask_file=args.mask_file,
+        output_file=args.output_file,
+        brdf_model=args.brdf_model,
+        output_var=args.output_var,
+        output_format=args.output_format,
+        chunk_size=args.chunk_size,
+        compute_uncertainty=args.compute_uncertainty,
+        progress_bar=progress_bar,
+    )
+    print(f"BRDF correction finished. Saved to: {saved_path}")
+
+if __name__ == "__main__":
+    main()