MolecularAtmosphere.ussa_1976() fix

CHANGELOG.md

@@ -74,6 +74,9 @@ Entries marked with a ︎⚠️ symbol require particular attention during upgra
   {ghpr}`346`).
 * Fixed a bug where {class}`.Layout` constructors would not raise if passed
   invalid azimuth values (*i.e.* outside the [0, 180]° range) ({ghpr}`345`).
+* Added `wavelength_range` optional parameter to `MolecularAtmosphere.ussa_1976()` 
+  constructor to automatically open absorption datasets and restore working default 
+  constructor ({ghpr}`334`).
 
 ### Documentation
 

src/eradiate/radprops/_afgl1986.py

@@ -50,6 +50,18 @@ class AFGL1986RadProfile(RadProfile):
     This class does not support ``mono`` modes.
     """
 
+    absorption_dataset: xr.Dataset = documented(
+        attrs.field(
+            converter=converters.to_dataset(
+                load_from_id=None  # TODO: open dataset (not load)
+            ),
+            validator=attrs.validators.instance_of(xr.Dataset),
+        ),
+        doc="Absorption coefficient dataset.",
+        type="Dataset",
+        init_type="Dataset or :class:`.PathLike`",
+    )
+
     _thermoprops: xr.Dataset = documented(
         attrs.field(
             factory=lambda: afgl_1986.make_profile(),
@@ -88,38 +100,6 @@ class AFGL1986RadProfile(RadProfile):
         default="True",
     )
 
-    absorption_dataset: xr.Dataset | None = documented(
-        attrs.field(
-            default=None,
-            converter=attrs.converters.optional(
-                converters.to_dataset(load_from_id=None)
-            ),  # TODO: open dataset (not load)
-            validator=attrs.validators.optional(
-                attrs.validators.instance_of(xr.Dataset)
-            ),
-        ),
-        doc="Absorption coefficient dataset. If ``None``, the absorption "
-        "coefficient is set to zero.",
-        type="Dataset or None",
-        init_type="Dataset or :class:`.PathLike`",
-        default="None",
-    )
-
-    @absorption_dataset.validator
-    @has_absorption.validator
-    def _check_absorption_dataset(self, attribute, value):
-        if not self.has_absorption and self.absorption_dataset is not None:
-            warnings.warn(
-                "When validating attribute 'absorption_dataset': specified "
-                "absorption dataset will be ignored because absorption is "
-                "disabled."
-            )
-        if self.has_absorption and self.absorption_dataset is None:
-            raise ValueError(
-                "When validating attribute 'absorption_dataset': no absorption "
-                "dataset was specified, absorption will be disabled."
-            )
-
     _zgrid: ZGrid | None = attrs.field(default=None, init=False)
 
     def __attrs_post_init__(self):

src/eradiate/radprops/_us76_approx.py

@@ -1,17 +1,17 @@
 from __future__ import annotations
 
 import typing as t
+import warnings
 
 import attrs
 import numpy as np
 import pint
 import xarray as xr
 
 from ._core import RadProfile, ZGrid, make_dataset
-from ._util_mono import get_us76_u86_4_spectrum_filename
 from .absorption import compute_sigma_a
 from .rayleigh import compute_sigma_s_air
-from .. import converters, data
+from .. import converters
 from ..attrs import documented, parse_docs
 from ..thermoprops import us76
 from ..units import to_quantity
@@ -106,6 +106,15 @@ class US76ApproxRadProfile(RadProfile):
       not include.
     """
 
+    absorption_dataset: xr.Dataset = documented(
+        attrs.field(
+            converter=converters.to_dataset(load_from_id=None),
+            validator=attrs.validators.instance_of(xr.Dataset),
+        ),
+        doc="Absorption coefficient data set.",
+        type=":class:`xarray.Dataset`",
+    )
+
     _thermoprops: xr.Dataset = documented(
         attrs.field(
             factory=lambda: us76.make_profile(),
@@ -144,22 +153,6 @@ class US76ApproxRadProfile(RadProfile):
         default="True",
     )
 
-    absorption_dataset: xr.Dataset | None = documented(
-        attrs.field(
-            default=None,
-            converter=attrs.converters.optional(
-                converters.to_dataset(load_from_id=None)
-            ),
-            validator=attrs.validators.optional(
-                attrs.validators.instance_of(xr.Dataset)
-            ),
-        ),
-        doc="Absorption coefficient data set. If ``None``, the default "
-        "absorption coefficient data set is opened.",
-        type=":class:`xarray.Dataset` or None",
-        default="None",
-    )
-
     _zgrid: ZGrid | None = attrs.field(default=None, init=False)
 
     def __attrs_post_init__(self):
@@ -222,12 +215,8 @@ def eval_sigma_t_mono(self, w: pint.Quantity, zgrid: ZGrid) -> pint.Quantity:
     def eval_sigma_a_mono(self, w: pint.Quantity, zgrid: ZGrid) -> pint.Quantity:
         profile = self._thermoprops_interp(zgrid)
 
-        # TODO: refactor so that absorption dataset must be provided at init
         if self.has_absorption:
-            if self.absorption_dataset is None:  # ! this is never tested
-                ds = data.open_dataset(get_us76_u86_4_spectrum_filename(w))
-            else:
-                ds = self.absorption_dataset
+            ds = self.absorption_dataset
 
             # Compute scattering coefficient
             result = compute_sigma_a(

src/eradiate/scenes/atmosphere/_molecular_atmosphere.py

@@ -17,7 +17,7 @@
 from ._core import AbstractHeterogeneousAtmosphere
 from ..core import traverse
 from ..phase import PhaseFunction, RayleighPhaseFunction, phase_function_factory
-from ... import converters
+from ... import converters, data
 from ...attrs import documented, parse_docs
 from ...contexts import KernelContext
 from ...quad import Quad
@@ -34,6 +34,45 @@
 from ...units import unit_registry as ureg
 from ...util.misc import summary_repr
 
+DEFAULT_WAVELENGTH_RANGE = [545.0, 555.0] * ureg.nm
+
+
+def open_us76_dataset(wavelength_range: pint.Quantity) -> xr.Dataset:
+    """
+    Open the monochromatic absorption dataset corresponding to the given
+    wavelength range, for the U.S. 1976 Standard Atmosphere.
+    """
+    path = "spectra/absorption/us76_u86_4"
+    wlmin, wlmax = wavelength_range
+    wnmin = (1 / wlmax).m_as("1/cm")
+    wnmax = (1 / wlmin).m_as("1/cm")
+    w = np.concatenate([np.arange(4000, 25711, 1000), [25711]])
+    if wnmin < w[0] or wnmax > w[-1]:
+        raise ValueError(
+            f"Invalid wavelength range: {wavelength_range}. "
+            f"Supported wavelength range: {1e7 / w[-1]:.1f}-"
+            f"{1e7 / w[0]:.1f} nm."
+        )
+
+    wmin = w[wnmin > w][-1]
+    wmax = w[wnmax < w][0]
+    iwmin = w.tolist().index(wmin)
+    iwmax = w.tolist().index(wmax)
+    paths = [
+        f"{path}/us76_u86_4-spectra-{w[i]}_{w[i+1]}.nc" for i in range(iwmin, iwmax)
+    ]
+
+    if len(paths) == 1:
+        absorption_dataset = data.open_dataset(paths[0])
+    else:
+        datasets = [
+            data.open_dataset(path).isel(w=slice(0, -1))  # strip endpoints
+            for path in paths
+        ]
+        absorption_dataset = xr.concat(datasets, dim="w")
+
+    return absorption_dataset
+
 
 @parse_docs
 @attrs.define(eq=False, slots=False)
@@ -153,6 +192,12 @@ def update(self) -> None:
         self.phase.id = self.phase_id
 
         if self.thermoprops.title == "U.S. Standard Atmosphere 1976":
+            # temporary fix
+            if self.absorption_dataset is None and self.has_absorption:
+                self.absorption_dataset = (
+                    "spectra/absorption/us76_u86_4/us76_u86_4-spectra-18000_19000.nc"
+                )
+
             self._radprops_profile = US76ApproxRadProfile(
                 thermoprops=self.thermoprops,
                 has_scattering=self.has_scattering,
@@ -363,16 +408,15 @@ def afgl_1986(
 
         thermoprops = afgl_1986.make_profile(model_id=model)
 
+        # open the absorption dataset corresponding to 'binset' and 'model'
         path = "ckd/absorption"
-        if kwargs.get("has_absorption", True):
-            if binset == "10nm":
-                absorption_dataset = f"{path}/10nm/afgl_1986-{model}-10nm-v3.nc"
-            elif binset == "1nm":
-                absorption_dataset = f"{path}/1nm/afgl_1986-{model}-1nm-v3.nc"
-            else:
-                raise ValueError(f"Invalid binset: {binset}")
+
+        if binset == "10nm":
+            absorption_dataset = f"{path}/10nm/afgl_1986-{model}-10nm-v3.nc"
+        elif binset == "1nm":
+            absorption_dataset = f"{path}/1nm/afgl_1986-{model}-1nm-v3.nc"
         else:
-            absorption_dataset = None
+            raise ValueError(f"Invalid binset: {binset}")
 
         if concentrations is not None:
             factors = compute_scaling_factors(
@@ -399,6 +443,7 @@ def ussa_1976(
         cls,
         levels: pint.Quantity | None = None,
         concentrations: dict[str, pint.Quantity] | None = None,
+        wavelength_range: pint.Quantity | None = None,
         **kwargs: dict[str, t.Any],
     ) -> MolecularAtmosphere:
         """
@@ -413,6 +458,13 @@ def ussa_1976(
         concentrations : dict
             Molecules concentrations as a ``{str: quantity}`` mapping.
 
+        wavelength_range: quantity
+            Wavelength range wherein the atmosphere radiative properties are
+            expected to be computed.
+            This information is used to select the monochromatic absorption
+            dataset(s) to open.
+            If None, defaults to [545.0, 555.0] nm.
+
         **kwargs
             Keyword arguments passed to the :class:`.MolecularAtmosphere`
             constructor.
@@ -422,6 +474,19 @@ def ussa_1976(
         :class:`.MolecularAtmosphere`
             U.S. Standard Atmosphere (1976) molecular atmosphere object.
         """
+        if "absorption_dataset" in kwargs:
+            raise TypeError(
+                "Cannot pass 'absorption_dataset' keyword argument. The "
+                "'ussa_1976' constructor sets the absorption dataset "
+                "automatically."
+            )
+
+        # open the absorption dataset corresponding to 'wavelength_range
+        if wavelength_range is None:
+            wavelength_range = DEFAULT_WAVELENGTH_RANGE
+
+        absorption_dataset = open_us76_dataset(wavelength_range)
+
         if levels is None:
             levels = np.linspace(0, 100, 101) * ureg.km
 
@@ -431,4 +496,8 @@ def ussa_1976(
             factors = compute_scaling_factors(ds=ds, concentration=concentrations)
             ds = rescale_concentration(ds=ds, factors=factors)
 
-        return cls(thermoprops=ds, **kwargs)
+        return cls(
+            thermoprops=ds,
+            absorption_dataset=absorption_dataset,
+            **kwargs,
+        )

tests/01_eradiate/01_unit/experiments/test_atmosphere.py

@@ -12,13 +12,17 @@
 from eradiate.scenes.shapes import RectangleShape
 from eradiate.scenes.spectra import MultiDeltaSpectrum
 from eradiate.test_tools.types import check_scene_element
+from eradiate.test_tools.util import skipif_data_not_found
 
 
 def test_atmosphere_experiment_construct_default(modes_all_double):
     """
     AtmosphereExperiment initializes successfully with default parameters in
     all modes.
     """
+    skipif_data_not_found(
+        "spectra/absorption/us76_u86_4/us76_u86_4-spectra-18000_19000.nc"
+    )
     exp = AtmosphereExperiment()
 
     # Check that the atmosphere's geometry is overridden by the experiment's
@@ -29,6 +33,9 @@ def test_atmosphere_experiment_extra_objects(mode_mono):
     """
     Extra objects can be added to the scene.
     """
+    skipif_data_not_found(
+        "spectra/absorption/us76_u86_4/us76_u86_4-spectra-18000_19000.nc"
+    )
     exp = AtmosphereExperiment(
         extra_objects={
             "reference_surface": {
@@ -49,6 +56,10 @@ def test_atmosphere_experiment_construct_measures(modes_all_double):
     """
     A variety of measure specifications are acceptable.
     """
+    if eradiate.mode().is_mono:
+        skipif_data_not_found(
+            "spectra/absorption/us76_u86_4/us76_u86_4-spectra-18000_19000.nc"
+        )
     # Init with a single measure (not wrapped in a sequence)
     assert AtmosphereExperiment(measures=MultiDistantMeasure())
 
@@ -73,6 +84,9 @@ def test_atmosphere_experiment_construct_normalize_measures(
     """
     Default measure target is set to ground level.
     """
+    skipif_data_not_found(
+        "spectra/absorption/us76_u86_4/us76_u86_4-spectra-18000_19000.nc"
+    )
     exp = AtmosphereExperiment(geometry=geometry)
     assert np.allclose(exp.measures[0].target.xyz, expected)
 
@@ -81,7 +95,9 @@ def test_atmosphere_experiment_kernel_dict(mode_mono):
     """
     Test non-trivial kernel dict generation behaviour.
     """
-
+    skipif_data_not_found(
+        "spectra/absorption/us76_u86_4/us76_u86_4-spectra-18000_19000.nc"
+    )
     # With an atmosphere
     exp = AtmosphereExperiment(
         geometry={"type": "plane_parallel", "width": 42.0, "width_units": "km"},
@@ -128,16 +144,15 @@ def test_atmosphere_experiment_kernel_dict(mode_mono):
 @pytest.mark.slow
 def test_atmosphere_experiment_real_life(mode_mono):
     # Construct with typical parameters
-    test_absorption_data_set = eradiate.data.data_store.fetch(
-        "tests/spectra/absorption/us76_u86_4-spectra-4000_25711.nc"
+    skipif_data_not_found(
+        "spectra/absorption/us76_u86_4/us76_u86_4-spectra-18000_19000.nc"
     )
     exp = AtmosphereExperiment(
         surface={"type": "rpv"},
         atmosphere={
             "type": "heterogeneous",
             "molecular_atmosphere": {
                 "construct": "ussa_1976",
-                "absorption_dataset": test_absorption_data_set,
             },
         },
         illumination={"type": "directional", "zenith": 45.0},
@@ -159,15 +174,9 @@ def test_atmosphere_experiment_real_life(mode_mono):
 @pytest.mark.slow
 def test_atmosphere_experiment_run_basic(modes_all_double):
     if eradiate.mode().is_mono:
-        molecular_atmosphere_params = {
-            "construct": "ussa_1976",
-            "absorption_dataset": "tests/spectra/absorption/us76_u86_4-spectra-4000_25711.nc",
-        }
+        molecular_atmosphere_params = {"construct": "ussa_1976"}
     else:
-        molecular_atmosphere_params = {
-            "construct": "afgl_1986",
-            # "absorption_dataset": "ckd/absorption/10nm/afgl_1986-us_standard-10nm.nc"
-        }
+        molecular_atmosphere_params = {"construct": "afgl_1986"}
 
     exp = AtmosphereExperiment(
         atmosphere={
@@ -193,6 +202,10 @@ def test_atmosphere_experiment_run_detailed(modes_all):
     """
     Test for correctness of the result dataset generated by AtmosphereExperiment.
     """
+    if eradiate.mode().is_mono:
+        skipif_data_not_found(
+            "spectra/absorption/us76_u86_4/us76_u86_4-spectra-18000_19000.nc"
+        )
     # Create simple scene
     exp = AtmosphereExperiment(
         measures=[