Reduce the number of lines in flux calculations

atmtools · Nov 6, 2024 · 5dd0bf2 · 5dd0bf2
1 parent 4a21d4d
commit 5dd0bf2
Show file tree

Hide file tree

Showing 7 changed files with 200 additions and 8 deletions.
diff --git a/examples/recipes/AtmosphericFlux/atmospheric_flux.ipynb b/examples/recipes/AtmosphericFlux/atmospheric_flux.ipynb
@@ -37,7 +37,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "fop = pyarts.recipe.AtmosphericFlux()"
+    "fop = pyarts.recipe.AtmosphericFlux(species=[\"H2O-161\", \"O2-66\", \"N2-44\", \"CO2-626\", \"O3-XFIT\"],\n",
+    "                                    remove_lines_percentile={\"H2O\": 70})"
    ]
   },
   {
@@ -120,7 +121,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3 (ipykernel)",
+   "display_name": "Python 3",
    "language": "python",
    "name": "python3"
   },
@@ -134,7 +135,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.12.6"
+   "version": "3.12.7"
   }
  },
  "nbformat": 4,

diff --git a/examples/recipes/SpectralAtmosphericFlux/spectral_atmospheric_flux.ipynb b/examples/recipes/SpectralAtmosphericFlux/spectral_atmospheric_flux.ipynb
@@ -38,7 +38,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "fop = pyarts.recipe.SpectralAtmosphericFlux()"
+    "fop = pyarts.recipe.SpectralAtmosphericFlux(species=[\"H2O-161\", \"O2-66\", \"N2-44\", \"CO2-626\", \"O3-XFIT\"],\n",
+    "                                            remove_lines_percentile={\"H2O\": 70})"
    ]
   },
   {
@@ -130,7 +131,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3 (ipykernel)",
+   "display_name": "Python 3",
    "language": "python",
    "name": "python3"
   },
@@ -144,7 +145,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.12.6"
+   "version": "3.12.7"
   }
  },
  "nbformat": 4,

diff --git a/python/src/pyarts/recipe/AtmosphericFlux.py b/python/src/pyarts/recipe/AtmosphericFlux.py
@@ -33,6 +33,7 @@ def __init__(
         solar_latitude: float = 0.0,
         solar_longitude: float = 0.0,
         species: list = ["H2O-161", "O2-66", "N2-44", "CO2-626", "O3-XFIT"],
+        remove_lines_percentile: dict[pyarts.arts.SpeciesEnum, float] | float | None = None,
     ):
         """Compute the total flux for a given atmospheric profile and surface temperature
 
@@ -50,6 +51,7 @@ def __init__(
             solar_latitude (float, optional): Latitude of sun [degrees]. Defaults to 0.0.
             solar_longitude (float, optional): Longitude of sun [degrees]. Defaults to 0.0.
             species (list, optional): The list of absorption species. Defaults to [ "H2O-161", "O2-66", "N2-44", "CO2-626", "O3-XFIT", ].
+            remove_lines_percentile (dict | float | None, optional): The percentile of lines to remove [0, 100].  Per species if dict. Defaults to None.
         """
 
         self.visible_surface_reflectivity = visible_surface_reflectivity
@@ -69,6 +71,9 @@ def __init__(
             self.ws.absorption_bands[band].cutoff = "ByLine"
             self.ws.absorption_bands[band].cutoff_value = 750e9
 
+        if remove_lines_percentile is not None:
+            self.ws.absorption_bands.keep_hitran_s(remove_lines_percentile)
+
         self.ws.propagation_matrix_agendaAuto()
 
         self.ws.surface_fieldSetPlanetEllipsoid(option="Earth")

diff --git a/python/src/pyarts/recipe/SpectralAtmosphericFlux.py b/python/src/pyarts/recipe/SpectralAtmosphericFlux.py
@@ -32,6 +32,7 @@ def __init__(
         solar_latitude: float = 0.0,
         solar_longitude: float = 0.0,
         species=["H2O-161", "O2-66", "N2-44", "CO2-626", "O3-XFIT"],
+        remove_lines_percentile: dict[pyarts.arts.SpeciesEnum, float] | float | None = None,
     ):
         """Compute the total flux for a given atmospheric profile and surface temperature
 
@@ -49,6 +50,7 @@ def __init__(
             solar_latitude (float, optional): Latitude of sun [degrees]. Defaults to 0.0.
             solar_longitude (float, optional): Longitude of sun [degrees]. Defaults to 0.0.
             species (list, optional): The list of absorption species. Defaults to [ "H2O-161", "O2-66", "N2-44", "CO2-626", "O3-XFIT", ].
+            remove_lines_percentile (dict | float | None, optional): The percentile of lines to remove [0, 100].  Per species if dict. Defaults to None.
         """
 
         self.visible_surface_reflectivity = visible_surface_reflectivity
@@ -68,6 +70,9 @@ def __init__(
             self.ws.absorption_bands[band].cutoff = "ByLine"
             self.ws.absorption_bands[band].cutoff_value = 750e9
 
+        if remove_lines_percentile is not None:
+            self.ws.absorption_bands.keep_hitran_s(remove_lines_percentile)
+
         self.ws.propagation_matrix_agendaAuto()
 
         self.ws.surface_fieldSetPlanetEllipsoid(option="Earth")

diff --git a/src/core/lbl/lbl_data.cpp b/src/core/lbl/lbl_data.cpp
@@ -237,4 +237,75 @@ std::unordered_set<SpeciesEnum> species_in_bands(
   }
   return out;
 }
+
+void keep_hitran_s(std::unordered_map<QuantumIdentifier, band_data>& bands,
+                   const std::unordered_map<SpeciesEnum, Numeric>& keep,
+                   const Numeric T0) {
+  for (auto& [key, data] : bands) {
+    const auto ptr = keep.find(key.Species());
+    if (ptr != keep.end()) {
+      std::erase_if(data.lines, [&key, &T0, &min_s = ptr->second](line& line) {
+        return line.hitran_s(key.Isotopologue(), T0) < min_s;
+      });
+    }
+  }
+}
+
+std::unordered_map<SpeciesEnum, Numeric> percentile_hitran_s(
+    const std::unordered_map<QuantumIdentifier, band_data>& bands,
+    const Numeric approx_percentile,
+    const Numeric T0) {
+  ARTS_USER_ERROR_IF(approx_percentile < 0 or approx_percentile > 100,
+                     "Approximate percentile must be between 0 and 100");
+
+  std::unordered_map<SpeciesEnum, std::vector<Numeric>> compute;
+  for (auto& [key, data] : bands) {
+    for (auto& line : data.lines) {
+      compute[key.Species()].push_back(line.hitran_s(key.Isotopologue(), T0));
+    }
+  }
+
+  std::unordered_map<SpeciesEnum, Numeric> out;
+  for (auto& [spec, values] : compute) {
+    if (const Size N = values.size(); N != 0) {
+      std::ranges::sort(values);
+      const Size i =
+          static_cast<Size>(static_cast<Numeric>(N) * approx_percentile * 0.01);
+      out[spec] = values[std::clamp<Size>(i, 0, N - 1)];
+    }
+  }
+
+  return out;
+}
+
+std::unordered_map<SpeciesEnum, Numeric> percentile_hitran_s(
+    const std::unordered_map<QuantumIdentifier, band_data>& bands,
+    const std::unordered_map<SpeciesEnum, Numeric>& approx_percentile,
+    const Numeric T0) {
+  ARTS_USER_ERROR_IF(
+      std::ranges::any_of(approx_percentile | std::views::values,
+                          [](auto i) { return i < 0 or i > 100; }),
+      "Approximate percentile must be between 0 and 100");
+
+  std::unordered_map<SpeciesEnum, std::vector<Numeric>> compute;
+  for (auto& [key, data] : bands) {
+    if (approx_percentile.contains(key.Species())) {
+      for (auto& line : data.lines) {
+        compute[key.Species()].push_back(line.hitran_s(key.Isotopologue(), T0));
+      }
+    }
+  }
+
+  std::unordered_map<SpeciesEnum, Numeric> out;
+  for (auto& [spec, values] : compute) {
+    if (const Size N = values.size(); N != 0) {
+      std::ranges::sort(values);
+      const Size i = static_cast<Size>(static_cast<Numeric>(N) *
+                                       approx_percentile.at(spec) * 0.01);
+      out[spec]    = values[std::clamp<Size>(i, 0, N - 1)];
+    }
+  }
+
+  return out;
+}
 }  // namespace lbl
diff --git a/src/core/lbl/lbl_data.h b/src/core/lbl/lbl_data.h
@@ -15,12 +15,10 @@
 #include <limits>
 #include <unordered_set>
 #include <vector>
-#include <unordered_set>
 
 #include "lbl_lineshape_model.h"
 #include "lbl_zeeman.h"
 
-
 namespace lbl {
 struct line {
   //! Einstein A coefficient
@@ -246,8 +244,51 @@ std::ostream& operator<<(
     std::ostream& os,
     const std::unordered_map<QuantumIdentifier, band_data>& x);
 
+/** Returns all species in the band, including those that are broadening species
+ * 
+ * @param bands The bands to search
+ * @return std::unordered_set<SpeciesEnum> 
+ */
 std::unordered_set<SpeciesEnum> species_in_bands(
     const std::unordered_map<QuantumIdentifier, band_data>& bands);
+
+/** Wraps keep_hitran_s for band_data per species, to remove all lines that are not in the keep map
+ * 
+ * @param bands The bands to use
+ * @param keep A map of species to minimum hitran_s values to keep.  Missing species keep all their lines.
+ * @param T0 The reference temperature.  Defaults to 296.0.
+ */
+void keep_hitran_s(std::unordered_map<QuantumIdentifier, band_data>& bands,
+                   const std::unordered_map<SpeciesEnum, Numeric>& keep,
+                   const Numeric T0 = 296);
+
+/** Compute what lines should be kept.  Meant to be used in conjunction with keep_hitran_s.
+ * 
+ * The same percentile of lines are kept for all species
+ * 
+ * @param bands The bands to use
+ * @param approx_percentile The percentile to keep [0, 100]
+ * @param T0 The reference temperature.  Defaults to 296.0.
+ * @return A map of species to minimum hitran_s values to keep
+ */
+std::unordered_map<SpeciesEnum, Numeric> percentile_hitran_s(
+    const std::unordered_map<QuantumIdentifier, band_data>& bands,
+    const Numeric approx_percentile,
+    const Numeric T0 = 296);
+
+/** Compute what lines should be kept.  Meant to be used in conjunction with keep_hitran_s.
+ * 
+ * Only species in the approx_percentile map are affected.  Otherwise acts like the pure index version.
+ * 
+ * @param bands The bands to use
+ * @param approx_percentile The percentile to keep species: [0, 100]
+ * @param T0 The reference temperature.  Defaults to 296.0.
+ * @return A map of species to minimum hitran_s values to keep
+ */
+std::unordered_map<SpeciesEnum, Numeric> percentile_hitran_s(
+    const std::unordered_map<QuantumIdentifier, band_data>& bands,
+    const std::unordered_map<SpeciesEnum, Numeric>& approx_percentile,
+    const Numeric T0 = 296);
 }  // namespace lbl
 
 //! Support hashing of line keys

diff --git a/src/python_interface/py_lbl.cpp b/src/python_interface/py_lbl.cpp
@@ -6,12 +6,14 @@
 #include <nanobind/stl/shared_ptr.h>
 #include <nanobind/stl/string.h>
 #include <nanobind/stl/unordered_map.h>
+#include <nanobind/stl/variant.h>
 #include <nanobind/stl/vector.h>
 #include <python_interface.h>
 
 #include <iomanip>
 #include <memory>
 #include <stdexcept>
+#include <unordered_map>
 
 #include "hpy_arts.h"
 #include "hpy_numpy.h"
@@ -415,6 +417,72 @@ int : The number of removed variables
       "spec"_a = SpeciesEnum::Bath,
       "Return the total number of lines");
 
+  aoab.def(
+      "remove_hitran_s",
+      &lbl::keep_hitran_s,
+      R"(Removes all lines with a weaker HITRAN-like line strength than those provided by the remove map.
+
+Parameters
+----------
+remove : dict
+    The species to keep with their respective minimum HITRAN-like line strengths
+T0 : float
+    The reference temperature. Defaults to 296.0.
+)",
+      "remove"_a,
+      "T0"_a = 296.0);
+
+  aoab.def(
+      "percentile_hitran_s",
+      [](const AbsorptionBands& self,
+         const std::variant<Numeric, std::unordered_map<SpeciesEnum, Numeric>>&
+             percentile,
+         const Numeric T0) {
+        return std::visit(
+            [&](auto& i) { return lbl::percentile_hitran_s(self, i, T0); },
+            percentile);
+      },
+      R"(Map of HITRAN linestrengths at a given percentile
+
+.. note::
+
+  The percentile is approximated by floating point arithmetic on the sorted HITRAN line strenght values.
+
+Parameters
+----------
+percentile : float or dict
+    The percentile to keep. If a float, the same percentile is used for all species. If a dict, the species are mapped to their respective percentiles.  Values must be [0, 100].
+T0 : float
+    The reference temperature. Defaults to 296.0.
+)",
+      "approximate_percentile"_a,
+      "T0"_a = 296.0);
+
+  aoab.def(
+      "keep_hitran_s",
+      [](AbsorptionBands& self,
+         const std::variant<Numeric, std::unordered_map<SpeciesEnum, Numeric>>&
+             percentile,
+         const Numeric T0) {
+        lbl::keep_hitran_s(
+            self,
+            std::visit(
+                [&](auto& i) { return lbl::percentile_hitran_s(self, i, T0); },
+                percentile),
+            T0);
+      },
+      R"(Wraps calling percentile_hitran_s followed by remove_hitran_s.
+
+Parameters
+----------
+percentile : float or dict
+    See percentile_hitran_s.
+T0 : float
+    The reference temperature. Defaults to 296.0.
+)",
+      "approximate_percentile"_a,
+      "T0"_a = 296.0);
+
   py::class_<LinemixingEcsData> led(m, "LinemixingEcsData");
   workspace_group_interface(led);