py2sess

py2sess is a Python implementation of the optimized 2S-ESS radiative-transfer model. It supports solar and thermal forward calculations with NumPy and optional torch backends. It does not call the original Fortran code.

Install

python3 -m pip install py2sess

For local development:

python3 -m pip install -e ".[torch,dev]"

The optional native backend uses the installed PyTorch shared libraries. Install PyTorch first when building native wheels from source or when using a backend that needs torch tensors.

Build

py2sess uses CMake through scikit-build-core. By default, source builds prepare the Python package without compiling the optional native backend.

cmake -S . -B build
cmake --build build
python3 -m build

To build a local wheel with the native backend, install PyTorch and build without PEP 517 isolation so CMake can find Torch:

python3 -m pip install build scikit-build-core setuptools-scm torch
python3 -m build --wheel --no-isolation -Ccmake.define.PY2SESS_BUILD_NATIVE=ON

Native wheels link against the PyTorch shared libraries supplied by the installed torch package.

Quick Start

Solar:

import numpy as np
from py2sess import TwoStreamEss, TwoStreamEssOptions

solver = TwoStreamEss(TwoStreamEssOptions(nlyr=3, mode="solar"))
result = solver.forward(
    tau=np.full(3, 0.02),
    ssa=np.full(3, 0.2),
    g=np.full(3, 0.1),
    z=np.array([3.0, 2.0, 1.0, 0.0]),
    angles=[30.0, 20.0, 0.0],  # sza, vza, relative azimuth
    albedo=0.3,
)
print(result.radiance)

Thermal:

solver = TwoStreamEss(TwoStreamEssOptions(nlyr=3, mode="thermal"))
result = solver.forward(
    tau=np.full(3, 0.1),
    ssa=np.zeros(3),
    g=np.zeros(3),
    z=np.array([3.0, 2.0, 1.0, 0.0]),
    angles=20.0,
    planck=np.array([1.0, 1.1, 1.2, 1.3]),
    surface_planck=1.4,
    emissivity=1.0,
)

Batched wavelengths use leading dimensions:

solver = TwoStreamEss(TwoStreamEssOptions(nlyr=3, mode="thermal"))
tau = np.full((100, 3), 0.02)
result = solver.forward(
    tau=tau,
    ssa=np.zeros_like(tau),
    g=np.zeros_like(tau),
    z=np.array([3.0, 2.0, 1.0, 0.0]),
    angles=20.0,
    planck=np.ones((100, 4)),
    surface_planck=np.ones(100),
    emissivity=np.ones(100),
)
print(result.radiance.shape)  # (100,)

Level fluxes use the final axis for TOA-to-BOA levels. This clear absorbing solar case has an analytic Beer-Lambert flux solution:

import numpy as np
from py2sess import TwoStreamEss, TwoStreamEssOptions

sza = 30.0
mu0 = np.cos(np.deg2rad(sza))
fbeam = 1.0
tau = np.array([0.1, 0.2])
z = np.array([2.0, 1.0, 0.0])

solver = TwoStreamEss(
    TwoStreamEssOptions(
        nlyr=tau.size,
        mode="solar",
        plane_parallel=True,
        delta_scaling=False,
        downwelling=True,
        output_levels=True,
        output_fluxes=True,
        fo_flux_n_mu=8,
    )
)
result = solver.forward(
    tau=tau,
    ssa=np.zeros_like(tau),  # pure absorption
    g=np.zeros_like(tau),
    z=z,
    angles=[sza, 0.0, 0.0],
    fbeam=fbeam,
    albedo=0.0,  # black surface: no upward reflected flux
    delta_m_truncation_factor=np.zeros_like(tau),
    include_fo=True,
)

level_tau = np.concatenate(([0.0], np.cumsum(tau)))
analytic_down = fbeam * mu0 * np.exp(-level_tau / mu0)

np.testing.assert_allclose(result.flux_down[0], analytic_down, atol=1.0e-9)
np.testing.assert_allclose(result.flux_up[0], 0.0, atol=1.0e-8)
np.testing.assert_allclose(result.flux_net, result.flux_up - result.flux_down)

print(result.flux_down[0])

Torch CPU float64:

solver = TwoStreamEss(
    TwoStreamEssOptions(nlyr=3, mode="solar", backend="torch", torch_dtype="float64")
)

API Notes

Core inputs are tau, ssa, g, z, angles, and the surface/source terms needed by the selected mode. Solar angles are [sza, vza, raz] in degrees; thermal angles are viewing zenith angles. Heights are in km, ordered top to bottom.

See docs/api_arguments.md for the full argument table and conventions. Level-flux conventions are summarized in docs/level_fluxes.md.

Examples

python3 examples/level_flux_beer_lambert.py
python3 examples/build_thermal_source_from_temperature.py
python3 examples/retrieve_synthetic_spectra.py --case uv --noise-level 0

Scene/profile runs:

from py2sess.scene import load_scene

scene = load_scene(profile="profile.txt", config="scene.yaml")
result = scene.forward(backend="numpy", include_fo=True)

Full-spectrum benchmark details are in docs/full_spectrum_benchmarks.md.

Tests

python3 -m unittest discover -s tests -v
python3 -m ruff check .
python3 -m ruff format --check .

Full-spectrum benchmarks use profile text plus scene YAML inputs and Python optical preprocessing. Keep large local cross-section tables, benchmark bundles, and generated outputs out of git.

References

• Fortran repository: vnatraj1/2S-ESS
• Paper: Natraj et al., 2022, JQSRT