Merge pull request #2164 from gerritholl/low-moisture

doc/source/conf.py

@@ -86,7 +86,8 @@ def __getattr__(cls, name):
 # coming with Sphinx (named 'sphinx.ext.*') or your custom ones.
 extensions = ['sphinx.ext.autodoc', 'sphinx.ext.intersphinx', 'sphinx.ext.todo', 'sphinx.ext.coverage',
               'sphinx.ext.doctest', 'sphinx.ext.napoleon', 'sphinx.ext.autosummary', 'doi_role',
-              'sphinx.ext.viewcode', 'sphinxcontrib.apidoc']
+              'sphinx.ext.viewcode', 'sphinxcontrib.apidoc',
+              'sphinx.ext.mathjax']
 
 # API docs
 apidoc_module_dir = "../../satpy"

satpy/enhancements/atmosphere.py

@@ -0,0 +1,110 @@
+# Copyright (c) 2022- Satpy developers
+#
+# This file is part of satpy.
+#
+# satpy is free software: you can redistribute it and/or modify it under the
+# terms of the GNU General Public License as published by the Free Software
+# Foundation, either version 3 of the License, or (at your option) any later
+# version.
+#
+# satpy is distributed in the hope that it will be useful, but WITHOUT ANY
+# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+# A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License along with
+# satpy.  If not, see <http://www.gnu.org/licenses/>.
+"""Enhancements related to visualising atmospheric phenomena."""
+
+import datetime
+
+import dask.array as da
+import xarray as xr
+
+
+def essl_moisture(img, low=1.1, high=1.6) -> None:
+    r"""Low level moisture by European Severe Storms Laboratory (ESSL).
+
+    Expects a mode L image with data corresponding to the ratio of the
+    calibrated reflectances for the 0.86 µm and 0.906 µm channel.
+
+    This composite and its colorisation were developed by ESSL.
+
+    Ratio values are scaled from the range ``[low, high]``, which is by default
+    between 1.1 and 1.6, but might be tuned based on region or sensor,
+    to ``[0, 1]``.  Values outside this range are clipped.  Color values
+    for red, green, and blue are calculated as follows, where ``x`` is the
+    ratio between the 0.86 µm and 0.905 µm channels:
+
+    .. math::
+
+        R = \max(1.375 - 2.67 x, -0.75 + x) \\
+        G = 1 - \frac{8x}{7} \\
+        B = \max(0.75 - 1.5 x, 0.25 - (x - 0.75)^2) \\
+
+    The value of ``img.data`` is modified in-place.
+
+    A color interpretation guide is pending further adjustments to the
+    parameters for current and future sensors.
+
+    Args:
+        img: XRImage containing the relevant composite
+        low: optional, low end for scaling, defaults to 1.1
+        high: optional, high end for scaling, defaults to 1.6
+    """
+    ratio = img.data
+    if _is_fci_test_data(img.data):
+        # Due to a bug in the FCI pre-launch simulated test data,
+        # the 0.86 µm channel is too bright.  To correct for this, its
+        # reflectances should be multiplied by 0.8.
+        ratio *= 0.8
+
+    with xr.set_options(keep_attrs=True):
+        ratio = _scale_and_clip(ratio, low, high)
+        red = _calc_essl_red(ratio)
+        green = _calc_essl_green(ratio)
+        blue = _calc_essl_blue(ratio)
+        data = xr.concat([red, green, blue], dim="bands")
+        data.attrs["mode"] = "RGB"
+        data["bands"] = ["R", "G", "B"]
+    img.data = data
+
+
+def _scale_and_clip(ratio, low, high):
+    """Scale ratio values to [0, 1] and clip values outside this range."""
+    scaled = (ratio - low) / (high - low)
+    scaled.data = da.clip(scaled.data, 0, 1)
+    return scaled
+
+
+def _calc_essl_red(ratio):
+    """Calculate values for red based on scaled and clipped ratio."""
+    red_a = 1.375 - 2.67 * ratio
+    red_b = -0.75 + ratio
+    red = xr.where(red_a > red_b, red_a, red_b)
+    red.data = da.clip(red.data, 0, 1)
+    return red
+
+
+def _calc_essl_green(ratio):
+    """Calculate values for green based on scaled and clipped ratio."""
+    green = 1 - (8/7) * ratio
+    green.data = da.clip(green.data, 0, 1)
+    return green
+
+
+def _calc_essl_blue(ratio):
+    """Calculate values for blue based on scaled and clipped ratio."""
+    blue_a = 0.75 - 1.5 * ratio
+    blue_b = 0.25 - (ratio - 0.75)**2
+    blue = xr.where(blue_a > blue_b, blue_a, blue_b)
+    blue.data = da.clip(blue.data, 0, 1)
+    return blue
+
+
+def _is_fci_test_data(data):
+    """Check if we are working with FCI test data."""
+    return ("sensor" in data.attrs and
+            "start_time" in data.attrs and
+            data.attrs["sensor"] == "fci" and
+            isinstance(data.attrs["start_time"], datetime.datetime) and
+            data.attrs["start_time"] < datetime.datetime(2022, 11, 30))

satpy/etc/composites/visir.yaml

@@ -505,3 +505,47 @@ composites:
       - wavelength: 0.64
       - wavelength: 1.61
     standard_name: cimss_cloud_type
+
+  essl_low_level_moisture:
+    description: >
+      Greyscale low level moisture using the ratio between the
+      0.91 µm and the 0.86 µm channels.  Developed by the
+      European Severe Storms Laboratory (ESSL).  For a color version,
+      see essl_colorized_low_level_moisture.
+    compositor: !!python/name:satpy.composites.RatioCompositor
+    prerequisites:
+    - wavelength: 0.905
+    - wavelength: 0.86
+    standard_name: essl_low_level_moisture
+
+  day_essl_low_level_moisture:
+    description: >
+      Daytime only version of essl_low_level_moisture.
+      Nighttime part of the scene will be masked out.
+    compositor: !!python/name:satpy.composites.DayNightCompositor
+    day_night: day_only
+    prerequisites:
+    - name: essl_low_level_moisture
+    standard_name: day_essl_low_level_moisture
+
+  essl_colorized_low_level_moisture:
+    description: >
+      Colorized low level moisture using the ratio between the
+      0.91 µm and the 0.86 µm channels.  Developed by the
+      European Severe Storms Laboratory (ESSL).  The colorization
+      is still under development and may be subject to change.
+    compositor: !!python/name:satpy.composites.RatioCompositor
+    prerequisites:
+    - wavelength: 0.86
+    - wavelength: 0.905
+    standard_name: essl_colorized_low_level_moisture
+
+  day_essl_colorized_low_level_moisture:
+    description: >
+      Daytime only version of essl_colorized_low_level_moisture.
+      Nighttime part of the scene will be masked out.
+    compositor: !!python/name:satpy.composites.DayNightCompositor
+    day_night: day_only
+    prerequisites:
+    - name: essl_colorized_low_level_moisture
+    standard_name: day_essl_colorized_low_level_moisture

satpy/etc/enhancements/generic.yaml

@@ -954,3 +954,24 @@ enhancements:
           stretch: crude
           min_stretch: [ 0,  0,   0]
           max_stretch: [50, 50, 100]
+
+  essl_low_level_moisture:
+    name: essl_low_level_moisture
+    operations:
+      - name: linear_stretch
+        method: !!python/name:satpy.enhancements.stretch
+        kwargs: {stretch: 'crude', min_stretch: 0.35, max_stretch: 0.85}
+
+  day_essl_low_level_moisture:
+    standard_name: day_essl_low_level_moisture
+    operations: []
+
+  essl_colorized_low_level_moisture:
+    name: essl_colorized_low_level_moisture
+    operations:
+      - name: essl_moisture
+        method: !!python/name:satpy.enhancements.atmosphere.essl_moisture
+
+  day_essl_colorized_low_level_moisture:
+    standard_name: day_essl_colorized_low_level_moisture
+    operations: []

satpy/tests/enhancement_tests/test_atmosphere.py

@@ -0,0 +1,61 @@
+# Copyright (c) 2022- Satpy developers
+#
+# This file is part of satpy.
+#
+# satpy is free software: you can redistribute it and/or modify it under the
+# terms of the GNU General Public License as published by the Free Software
+# Foundation, either version 3 of the License, or (at your option) any later
+# version.
+#
+# satpy is distributed in the hope that it will be useful, but WITHOUT ANY
+# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+# A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License along with
+# satpy.  If not, see <http://www.gnu.org/licenses/>.
+"""Tests for enhancements in enhancements/atmosphere.py."""
+
+import datetime
+
+import dask.array as da
+import numpy as np
+import xarray as xr
+from trollimage.xrimage import XRImage
+
+
+def test_essl_moisture():
+    """Test ESSL moisture compositor."""
+    from satpy.enhancements.atmosphere import essl_moisture
+
+    ratio = xr.DataArray(
+        da.linspace(1.0, 1.7, 25, chunks=5).reshape((5, 5)),
+        dims=("y", "x"),
+        attrs={"name": "ratio",
+               "calibration": "reflectance",
+               "units": "%",
+               "mode": "L"})
+    im = XRImage(ratio)
+
+    essl_moisture(im)
+    assert im.data.attrs["mode"] == "RGB"
+    np.testing.assert_array_equal(im.data["bands"], ["R", "G", "B"])
+    assert im.data.sel(bands="R")[0, 0] == 1
+    np.testing.assert_allclose(im.data.sel(bands="R")[2, 2], 0.04, rtol=1e-4)
+    np.testing.assert_allclose(im.data.sel(bands="G")[2, 2], 0.42857, rtol=1e-4)
+    np.testing.assert_allclose(im.data.sel(bands="B")[2, 2], 0.1875, rtol=1e-4)
+
+    # test FCI test data correction
+    ratio = xr.DataArray(
+        da.linspace(1.0, 1.7, 25, chunks=5).reshape((5, 5)),
+        dims=("y", "x"),
+        attrs={"name": "ratio",
+               "calibration": "reflectance",
+               "units": "%",
+               "mode": "L",
+               "sensor": "fci",
+               "start_time": datetime.datetime(1999, 1, 1)})
+    im = XRImage(ratio)
+    essl_moisture(im)
+    np.testing.assert_allclose(im.data.sel(bands="R")[3, 3], 0.7342, rtol=1e-4)
+    np.testing.assert_allclose(im.data.sel(bands="G")[3, 3], 0.7257, rtol=1e-4)
+    np.testing.assert_allclose(im.data.sel(bands="B")[3, 3], 0.39, rtol=1e-4)

satpy/tests/test_scene.py

@@ -833,7 +833,7 @@ def test_all_datasets_one_reader(self):
         num_reader_ds = 21 + 6
         assert len(id_list) == num_reader_ds
         id_list = scene.all_dataset_ids(composites=True)
-        assert len(id_list) == num_reader_ds + 29
+        assert len(id_list) == num_reader_ds + 33
 
     def test_all_datasets_multiple_reader(self):
         """Test all datasets for multiple readers."""