diff --git a/src/joseki/accessor.py b/src/joseki/accessor.py index 6ffac0ed..922723d1 100644 --- a/src/joseki/accessor.py +++ b/src/joseki/accessor.py @@ -27,6 +27,46 @@ def molecular_mass(m: str) -> pint.Quantity: # type: ignore[type-arg] return ureg.Quantity(float(molmass.Formula(m).mass), "dalton") +def _scaling_factor( + initial_amount: pint.Quantity, # type: ignore[type-arg] + target_amount: pint.Quantity, # type: ignore[type-arg] +) -> float: + """ + Compute scaling factor given initial and target amounts. + + Parameters + ---------- + initial_amount: :class:`~pint.Quantity` + Initial amount. + + target_amount: :class:`~pint.Quantity` + Target amount. + + Raises + ------ + ValueError + when the initial amount has a zero magnitude and the target amount + has a non-zero magnitude. + + Returns + ------- + float + Scaling factor + """ + if initial_amount.m == 0.0: + if target_amount.m == 0.0: + return 0.0 + else: + raise ValueError( + f"Cannot compute scaling factor when initial amount has " + f"magnitude of zero and target amount has a non-zero magnitude " + f"(got {target_amount})." + ) + else: + factor = (target_amount / initial_amount).m_as(ureg.dimensionless) + return float(factor) + + @xr.register_dataset_accessor("joseki") # type: ignore[no-untyped-call] class JosekiAccessor: # pragma: no cover """Joseki accessor.""" @@ -149,7 +189,7 @@ def number_density_at_sea_level( def mass_density_at_sea_level( self, ) -> t.Dict[str, pint.Quantity]: # type: ignore[type-arg] - """Compute mass mass density at sea level. + """Compute mass density at sea level. Returns ------- @@ -162,6 +202,80 @@ def mass_density_at_sea_level( for m in self._obj.m.values } + @property + def volume_mixing_fraction_at_sea_level( + self, + ) -> t.Dict[str, pint.Quantity]: # type: ignore[type-arg] + """Compute volume mixing fraction at sea level. + + Returns + ------- + dict: + A mapping of molecule and volume mixing fraction at sea level. + """ + ds = self._obj + return {m: to_quantity(ds.x.sel(m=m).isel(z=0)) for m in ds.m.values} + + def scaling_factors( + self, target: t.MutableMapping[str, pint.Quantity] # type: ignore[type-arg] + ) -> t.MutableMapping[str, float]: + """ + Compute scaling factor(s) to reach specific target amount(s). + + Seealso + ------- + :meth:`rescale` + + Parameters + ---------- + target + Mapping of molecule and target amount. + + Raises + ------ + ValueError + If a target amount has dimensions that are not supported. + + Returns + ------- + dict + Mapping of molecule and scaling factors. + + Notes + ----- + For each molecule in the ``target`` mapping, the target amount is + interpreted, depending on its dimensions (indicated in square brackets), + as: + + * a column number density [``length^-2``], + * a column mass density [``mass * length^-2``], + * a number densitx at sea level [``length^-3``], + * a mass density at sea level [``mass * length^-3``], + * a volume mixing fraction at sea level [``dimensionless``] + + The scaling factor is then evaluated as the ratio of the target amount + with the original amount, for each molecule. + """ + compute_initial_amount = { + "[length]^-2": self.column_number_density, + "[mass] * [length]^-2": self.column_mass_density, + "[length]^-3": self.number_density_at_sea_level, + "[mass] * [length]^-3": self.mass_density_at_sea_level, + "": self.volume_mixing_fraction_at_sea_level, + } + factors = {} + for m, target_amount in target.items(): + initial_amount = None + for dim in compute_initial_amount.keys(): + if target_amount.check(dim): + initial_amount = compute_initial_amount[dim][m] + if initial_amount is None: + raise ValueError + factors[m] = _scaling_factor( + initial_amount=initial_amount, target_amount=target_amount + ) + return factors + def rescale(self, factors: t.MutableMapping[str, float]) -> None: """Rescale molecules concentration in atmospheric profile. diff --git a/tests/test_accessor.py b/tests/test_accessor.py index 452ee5ff..2e495d4c 100644 --- a/tests/test_accessor.py +++ b/tests/test_accessor.py @@ -4,6 +4,7 @@ import xarray as xr import joseki +from joseki.accessor import _scaling_factor from joseki.units import to_quantity from joseki.units import ureg @@ -82,6 +83,46 @@ def test_mass_density_at_sea_level(test_dataset: xr.Dataset) -> None: ) +def test_volume_mixing_fraction_at_sea_level(test_dataset: xr.Dataset) -> None: + """CO2 volume mixing fraction at sea level in us_standard is 0.000333.""" + assert ( + test_dataset.joseki.volume_mixing_fraction_at_sea_level["CO2"] + == 0.000330 * ureg.dimensionless + ) + + +def test_scaling_factor() -> None: + """Returns 2.0 when target amount is twice larger than initial amount.""" + initial = 5 * ureg.m + target = 10 * ureg.m + assert _scaling_factor(initial_amount=initial, target_amount=target) == 2.0 + + +def test_scaling_factor_zero() -> None: + """Returns 0.0 when target amount and initial amount are both zero.""" + initial = 0.0 * ureg.m + target = 0.0 * ureg.m + assert _scaling_factor(initial_amount=initial, target_amount=target) == 0.0 + + +def test_scaling_factor_raises() -> None: + """Raises when the initial amount is zero but not the target amount.""" + initial = 0 * ureg.m + target = 10 * ureg.m + with pytest.raises(ValueError): + _scaling_factor(initial_amount=initial, target_amount=target) + + +def test_scaling_factors(test_dataset: xr.Dataset) -> None: + """Scaling factors keys match target amounts keys.""" + target = { + "H2O": 20.0 * ureg.kg * ureg.m**-2, + "O3": 350.0 * ureg.dobson_unit, + } + factors = test_dataset.joseki.scaling_factors(target=target) + assert all([k1 == k2 for k1, k2 in zip(target.keys(), factors.keys())]) + + def test_rescale(test_dataset: xr.Dataset) -> None: """Scaling factors are applied.""" factors = dict(