Rayleigh correction with elevation?

[yukaribbba]

CREFL uses taur = self._tau * np.exp(-self._height / SCALEHEIGHT) to fix rayleigh optical depth on given altitude. What about reflectance here in pyspectral? Any way to approach it?

[simonrp84]

Hi, I have thought about this in the past and it'd be a really nice addition. But, due to the way we do the Rayleigh correction it'd require regenerating the lookup tables with a new column for surface altitude. Not impossible but quite time consuming.

If I have the time I'd love to implement it though, so let's see :)

[yukaribbba]

require regenerating the lookup tables with a new column for surface altitude

You mean those on Zenodo? How to do that? I know very little about theories but maybe I could help with some repetitive tasks.

[yukaribbba]

I build this script to try to approach the generation of the original LUTs. Like @simonrp84 said it's very time cosuming if you want to complete the whole thing. Only some of the value groups are tested and the differences are around 3%. The result may be acceptable since the original ones were generated 8 years ago. It also includes elevation but for now only 0 altitude available.
If this script got verified we could see the impact of elevation. Then:

1. How do we set aerosol_set_tau_at_wvl according to the altitude?
2. Altitude steps and interpolation?

@simonrp84 @adybbroe @mraspaud When you got some time... Thanks!

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

import os
import math
from concurrent.futures import ThreadPoolExecutor, as_completed
from typing import Optional, Tuple, Any

import h5py
import numpy as np
import subprocess
from pathlib import Path
import time as time_calc

LIBRADTRAN_ROOT = Path("/home/ll/Downloads/libRadtran-2.0.6")
UVSPEC = LIBRADTRAN_ROOT / "bin" / "uvspec"
DATA_FILES_PATH = LIBRADTRAN_ROOT / "data"
SOLAR_FILE = DATA_FILES_PATH / "solar_flux" / "kurudz_full.dat"
OUTPUT_ROOT = Path("/home/ll/Downloads/LUTs")

AEROSOL_TYPES_TAU = {
    "continental_clean": 0.052,
    "continental_average": 0.12,
    "continental_polluted": 0.263,
    "urban": 0.536,
    "maritime_clean": 0.163,
    "maritime_polluted": 0.171,
    "maritime_tropical": 0.075,
    "desert": 0.276,
    "antarctic": 0.032,
    "rural": 0.162,
    "rayleigh_only": 0
}
ATM_FILES = {
    "midlatitude_summer": DATA_FILES_PATH / "atmmod" / "afglms.dat",
    "midlatitude_winter": DATA_FILES_PATH / "atmmod" / "afglmw.dat",
    "subarctic_summer": DATA_FILES_PATH / "atmmod" / "afglss.dat",
    "subarctic_winter": DATA_FILES_PATH / "atmmod" / "afglsw.dat",
    "tropical": DATA_FILES_PATH / "atmmod" / "afglt.dat",
    "us-standard": DATA_FILES_PATH / "atmmod" / "afglus.dat",
}


WAVELENGTH_TARGETS = np.arange(400, 801, 5, dtype=np.float64)

SUN_SEC = np.round(np.arange(1.0, 24.75 + 1e-9, 0.25, dtype=np.float64), 2)
SAT_SEC = np.round(np.arange(1.0, 3.0 + 1e-9, 0.25, dtype=np.float64), 2)
AZI_DIFF = np.arange(0, 181, 10, dtype=np.float64)
ALTITUDE_KM = np.round(np.arange(0, 0.05, 0.05, dtype=np.float64), 2)
# ALTITUDE_KM = np.round(np.arange(0, 8.85, 0.05, dtype=np.float64), 2)

AEROSOL_SPECIES_LIBRARY = "OPAC"
NUMBER_OF_STREAMS = 16
ALBEDO = 0.0

MAX_WORKERS = 8
REPORT_EVERY = 500

def load_solar_wavelengths(solar_file: Path) -> np.ndarray:
    waves = []
    with solar_file.open("r") as f:
        for line in f:
            s = line.strip()
            if not s or s.startswith("#"):
                continue
            parts = s.split()
            try:
                wl = float(parts[0])
                waves.append(wl)
            except Exception:
                pass
    if not waves:
        raise RuntimeError(f"Failed to load wavelengths from: {solar_file}")
    return np.array(waves, dtype=np.float64)

def map_wavelengths_to_nearest(targets: np.ndarray, available: np.ndarray) -> np.ndarray:
    idx = np.abs(available[:, None] - targets[None, :]).argmin(axis=0)
    return available[idx].astype(np.float64)

def sec_to_deg(sec_val: float) -> float:
    if sec_val < 1.0:
        sec_val = 1.0
    cosv = 1.0 / sec_val
    cosv = min(1.0, max(-1.0, cosv))
    return math.degrees(math.acos(cosv))

def _aerosol_block(aerosol_type: str) -> str:
    if aerosol_type == "rayleigh_only":
        return "\n"
    if aerosol_type == "rural":
        return (
            "aerosol_default\n"
            f"aerosol_haze              1\n"
            f"aerosol_vulcan            0\n"
            f"aerosol_set_tau_at_wvl    550 {AEROSOL_TYPES_TAU[aerosol_type]}\n"
        )
    else:
        return (
            "aerosol_default\n"
            f"aerosol_species_library  {AEROSOL_SPECIES_LIBRARY}\n"
            f"aerosol_species_file     {aerosol_type}\n"
            f"aerosol_set_tau_at_wvl   550 {AEROSOL_TYPES_TAU[aerosol_type]}\n"
        )

def make_uvspec_input(data_files_path: Path,
                      atm_file: Path,
                      solar_file: Path,
                      aerosol_type: str,
                      wavelength_nm: float,
                      sza_deg: float,
                      sat_sec: float,
                      azi_diff_deg: float,
                      altitude: float) -> str:
    umu = 1.0 / sat_sec  # satellite zenith cosine
    aero_txt = _aerosol_block(aerosol_type)
    txt = f"""
data_files_path          {data_files_path}
atmosphere_file          {atm_file}
source                   solar {solar_file}
rte_solver               disort
number_of_streams        {NUMBER_OF_STREAMS}
pseudospherical
wavelength               {wavelength_nm:.10f} {wavelength_nm:.10f}
albedo                   {ALBEDO}
altitude                 {altitude}
{aero_txt}
mol_abs_param            reptran fine
sza                      {sza_deg:.10f}
phi0                     0
umu                      {umu:.10f}
phi                      {azi_diff_deg:.10f}
zout                     TOA
output_quantity          reflectivity
output_user              lambda uu
quiet
"""

    return txt

def run_uvspec(uvspec_path: Path, input_text: str) -> float:
    """return wavelength and reflectance"""
    proc = subprocess.run(
        [str(uvspec_path)],
        input=input_text.encode("utf-8"),
        stdout=subprocess.PIPE,
        stderr=subprocess.PIPE,
        check=False,
    )
    if proc.returncode != 0:
        raise RuntimeError(f"uvspec Failed!（code={proc.returncode}）\n{proc.stderr.decode(errors='ignore')}")
    for line in proc.stdout.decode("utf-8", errors="ignore").splitlines():
        s = line.strip()
        if not s or s.startswith("#"):
            continue
        parts = s.split()
        if len(parts) >= 2:
            try:
                # wavelength = float(parts[0])
                reflectance = float(parts[1])
                return reflectance
            except Exception:
                continue
    raise RuntimeError("Failed to run LibRadTran to get reflectance。")

def write_h5(out_path: Path,
             wavelengths: np.ndarray,
             sun_sec: np.ndarray,
             azi_diff: np.ndarray,
             sat_sec: np.ndarray,
             altitude: np.ndarray,
             reflectance: np.ndarray,
             atmosphere_name: str):
    info_text = (
        "Some results of LibRadTran simulations for the correction of rayleigh scattering\n"
        "and clean marine aerosol composition (Hess et al., 1998) within satellite images \n"
        "of channels in the visible spectral range.\n"
    )
    with h5py.File(out_path, "w") as h5:
        h5.attrs["info"] = np.bytes_(info_text.encode("utf-8"))

        h5.create_dataset("wavelengths", data=wavelengths.astype(np.float64))
        h5.create_dataset("sun_zenith_secant", data=sun_sec.astype(np.float64))
        h5.create_dataset("azimuth_difference", data=azi_diff.astype(np.float64))
        h5.create_dataset("satellite_zenith_secant", data=sat_sec.astype(np.float64))
        if len(ALTITUDE_KM) > 1:
            h5.create_dataset("altitude", data=altitude.astype(np.float64))

        dset = h5.create_dataset(
            "reflectance",
            data=reflectance.astype(np.float64),
            dtype=np.float64,
            compression="gzip",
            compression_opts=4,
            shuffle=True,
            fletcher32=True,
        )
        dset.attrs["1-axis"] = np.bytes_(b"wavelength")
        dset.attrs["2-axis"] = np.bytes_(b"sun zenith secant")
        dset.attrs["3-axis"] = np.bytes_(b"azimuth difference angle")
        dset.attrs["4-axis"] = np.bytes_(b"satellite zenith secant")
        if len(ALTITUDE_KM) > 1:
            dset.attrs["5-axis"] = np.bytes_(b"altitude in km")
        dset.attrs["atmosphere"] = np.bytes_(atmosphere_name.encode("utf-8"))
        dset.attrs["info"] = np.bytes_(b"Rayleigh reflectance as a function of wavelength, "
                                       b"sun-zenith secant, azimuth difference angle and satellite zenith secant")

def task_compute_reflectance(params: Tuple[Any, ...]) -> Tuple[Tuple[int, int, int, int, Optional[int]], float]:
    """Single task for ThreadPoolExecutor."""
    (iw, isun, iazi, isat, idem,
     wl, sza_deg, vsec, dphi, alt_km,
     data_files_path, atm_file, solar_file, aerosol_type) = params

    uvin = make_uvspec_input(
        data_files_path=data_files_path,
        atm_file=atm_file,
        solar_file=solar_file,
        aerosol_type=aerosol_type,
        wavelength_nm=wl,
        sza_deg=sza_deg,
        sat_sec=vsec,
        azi_diff_deg=dphi,
        altitude=alt_km
    )
    ref = run_uvspec(UVSPEC, uvin)
    idx = (iw, isun, iazi, isat, idem) if idem is not None else (iw, isun, iazi, isat, None)
    return idx, ref

def main():
    # File check
    if not UVSPEC.exists():
        raise FileNotFoundError(f"Can't find uvspec: {UVSPEC}")
    if not SOLAR_FILE.exists():
        raise FileNotFoundError(f"Can't find solar file: {SOLAR_FILE}")

    # Map integer wavelengths to the nearest float one in solar file
    solar_waves = load_solar_wavelengths(SOLAR_FILE)
    mapped_waves = map_wavelengths_to_nearest(WAVELENGTH_TARGETS, solar_waves)

    # Dims
    n_w, n_sun, n_azi, n_sat, n_alt = len(mapped_waves), len(SUN_SEC), len(AZI_DIFF), len(SAT_SEC), len(ALTITUDE_KM)
    print(n_w, n_sun, n_azi, n_sat, n_alt)

    for aerosol_type in AEROSOL_TYPES_TAU:
        out_folder = aerosol_type + "_aerosol" if aerosol_type != "rayleigh_only" else aerosol_type
        (OUTPUT_ROOT / out_folder).mkdir(parents=True, exist_ok=True)

        for atm_name, atm_path in ATM_FILES.items():
            if n_alt > 1:
                reflectance = np.full((n_w, n_sun, n_azi, n_sat, n_alt), np.nan, dtype=np.float64)
            else:
                reflectance = np.full((n_w, n_sun, n_azi, n_sat), np.nan, dtype=np.float64)

            total = n_w * n_sun * n_azi * n_sat * n_alt
            completed = 0
            block_completed = 0
            block_t0 = time_calc.time()

            print(f"Start [{atm_name}/{aerosol_type}] total nodes = {total}, threads = {MAX_WORKERS}")

            with ThreadPoolExecutor(max_workers=MAX_WORKERS) as ex:
                futures = []

                for iw, wl in enumerate(mapped_waves):
                    for isun, ssec in enumerate(SUN_SEC):
                        sza_deg = sec_to_deg(ssec)
                        for iazi, dphi in enumerate(AZI_DIFF):
                            for isat, vsec in enumerate(SAT_SEC):

                                if n_alt > 1:
                                    for idem, alt_km in enumerate(ALTITUDE_KM):
                                        params = (iw, isun, iazi, isat, idem,
                                                  wl, sza_deg, vsec, dphi, alt_km,
                                                  DATA_FILES_PATH, atm_path, SOLAR_FILE, aerosol_type)
                                        futures.append(ex.submit(task_compute_reflectance, params))
                                else:
                                    alt_km = ALTITUDE_KM[0] if len(ALTITUDE_KM) > 0 else 0.0
                                    params = (iw, isun, iazi, isat, None,
                                              wl, sza_deg, vsec, dphi, alt_km,
                                              DATA_FILES_PATH, atm_path, SOLAR_FILE, aerosol_type)
                                    futures.append(ex.submit(task_compute_reflectance, params))

                    for fut in as_completed(futures):
                        try:
                            idx, ref = fut.result()
                            if n_alt > 1:
                                iw_i, isun_i, iazi_i, isat_i, idem_i = idx
                                reflectance[iw_i, isun_i, iazi_i, isat_i, idem_i] = ref
                            else:
                                iw_i, isun_i, iazi_i, isat_i, _ = idx
                                reflectance[iw_i, isun_i, iazi_i, isat_i] = ref
                        except Exception:
                            pass
                        finally:
                            completed += 1
                            block_completed += 1
                            if block_completed == REPORT_EVERY:
                                t1 = time_calc.time()
                                print(f"[{atm_name}/{aerosol_type}] "
                                      f"{completed}/{total} ({completed / total:.1%}), "
                                      f"last {REPORT_EVERY} nodes took {t1 - block_t0:.2f}s")
                                block_completed = 0
                                block_t0 = t1
                    futures.clear()

            out_name = f"rayleigh_lut_{atm_name}.h5"
            out_path = OUTPUT_ROOT / out_folder / out_name
            write_h5(out_path, mapped_waves, SUN_SEC, AZI_DIFF, SAT_SEC, ALTITUDE_KM, reflectance, atm_name)
            print(f"Finish write：{out_path}")

if __name__ == "__main__":
    main()