Add low level moisture composite

[gerritholl]

Add a composite to visualise low level moisture using the ratio between 0.91 µm and 0.86 µm. Experimental.

Greyscale version:

Test image produced from FCI test data for 2017-09-20 11:20, stretched between 0.35 and 0.86.

In the greyscale version, clouds are white. Areas with high low level moisture are dark.

The fine structure will be more visible in a colorised version.

Based on work by Hans-Peter Rösli and Pieter Groenemeijer.

• [ ] Closes #xxxx
• [ ] Tests added
• [ ] Fully documented
• [ ] Add your name to AUTHORS.md if not there already

[codecov[bot]]

Codecov Report

Merging #2164 (f341742) into main (5ff1611) will increase coverage by 0.00%. The diff coverage is 100.00%.

@@           Coverage Diff           @@
##             main    #2164   +/-   ##
=======================================
  Coverage   94.19%   94.20%           
=======================================
  Files         295      297    +2     
  Lines       45376    45438   +62     
=======================================
+ Hits        42743    42805   +62     
  Misses       2633     2633           

Flag	Coverage Δ
behaviourtests	4.67% <15.87%> (+0.01%)	:arrow_up:
unittests	94.85% <100.00%> (+<0.01%)	:arrow_up:

Flags with carried forward coverage won't be shown. Click here to find out more.

Impacted Files	Coverage Δ
satpy/enhancements/atmosphere.py	100.00% <100.00%> (ø)
satpy/tests/enhancement_tests/test_atmosphere.py	100.00% <100.00%> (ø)
satpy/tests/test_scene.py	99.49% <100.00%> (ø)

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

[coveralls]

Coverage remained the same at 94.801% when pulling f34174276040a5f08478ab2b6346fa317e6e2799 on gerritholl:low-moisture into 5ff16116d824ef38b809f448643b007a90b1f887 on pytroll:main.

[simonrp84]

Is there any documentation available for this? A link to a paper or even to a powerpoint deck or something would be nice to explain what this composite is showing, I suspect ESSL has some training docs, for example?

[gerritholl]

Is there any documentation available for this? A link to a paper or even to a powerpoint deck or something would be nice to explain what this composite is showing, I suspect ESSL has some training docs, for example?

I don't know if anything is public. I'm in contact with ESSL for the colour version and will ask about public documentation.

[gerritholl]

The colorized version is still wrong:

Investigating what could be the problem.

[gerritholl]

Should the colorised version be done in the compositor or in an enhancement? At first I tried it in the compositor:

    def __call__(self, projectables, others=None, **info):
        """Generate the ESSL low level moisture composite."""
        (nir_086, nir_091) = projectables
        with xr.set_options(keep_attrs=True):
            ratio = nir_091 / nir_086
            ratio = self._scale_and_clip(ratio)
            red = self._calc_red(ratio)
            green = self._calc_green(ratio)
            blue = self._calc_blue(ratio)
        return super().__call__([red, green, blue], **info)

where red, green, and blue are clipped to [0, 1]. This worked only if I passed enhance=False, dtype=uint8 to save_datasets, which is not ideal. Then I thought the colorisation should be in a enhancement instead:

    ratio = img.data

    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")
    return data

but I'm not sure if that's the correct way to create an enhancement (it results in a RGB image with bands 'L', 'L', 'L').

Looking elsewhere in Satpy, it would seem an enhancement is the better place to define and apply a colormap, but I'm not sure.

From a conversation on slack, I learned that it's possible to define an enhancement that does nothing, and that this will suppress the application of any other enhancements. I will try that next.

[djhoese]

It depends what you consider this operation. If you consider it similar to "applying a colormap" then it makes sense in the enhancements. If you consider it "generating an RGB from a single band" then I'd say it is a composite. It is a fine line.

In the enhancement case, the end result has to be within a normalized 0 to 1 range. The data is then scaled during saving to the output format (ex. geotiff) to the range of the data type (ex. 0 to 255 for uint8) and then clipped. Your bands dimension coordinates being ['L', 'L', 'L'] may be confusing the XRImage class, but I'm not sure. Even if it isn't, just for correctness you should overwrite the bands coordinates to be ['R', 'G', 'B']. This is just a consequence of combining DataArrays and not knowing the underlying "meaning" of the bands dimension.

[gerritholl]

I finally managed to implement it as an enhancement. The trick was to modify img.data in-place; the return value from an enhancement function is unused.

What remains is probably not a Satpy issue but rather related to tuning the parameters of the enhancement. For MODIS I'm getting close, but when I try the same for FCI, the results look totally off. It may be that the tuning parameters are quite sensitive to the exact boundaries, for example, for the scaling. Those may need to be different for FCI than for MODIS.

For this test I downloaded Terra MODIS data for the same date for which the FCI test data were simulated, and then generated the low level moisture composites from this PR as well as the input channels using either the MODIS data or the FCI test data. The simulated FCI 0.865 µm channel appears much brighter than the MODIS 0.8585 µm channel. The simulated FCI 0.914 µm channel is also brighter than the MODIS 0.905 µm channel, but the difference is smaller here. The greyscale ratio image is much brighter for MODIS than for simulated FCI. The colorised MODIS image looks OK, but the colorised FCI image looks completely off.

This could either be due to inaccuracies in the FCI simulated radiances, or due to actual differences in the spectral response function. If this is due to real differences in the spectral response function, we will need different recipes for this composite between MODIS and FCI.

FCI simulated, 0.865 µm reflectivity

FCI simulated, 0.914 µm reflectivity

MODIS 0.8585 µm reflectivity

MODIS, 0.905 µm reflectivity

MODIS greyscale low level moisture, ESSL scaling

MODIS greyscale low level moisture, automatic scaling (default stretch)

MODIS colorized low level moisture

Simulated FCI greyscale low level moisture, ESSL scaling

Simulated FCI low level moisture, automatic scaling (default stretch)

Simulated FCI colorized low level moisture

[gerritholl]

Are the test failures in test_all_datasets_one_reader and test_get_satpos_from_satname possibly related to this PR or are those failures also occurring elsewhere?

[djhoese]

The satpos test was mentioned by someone else in another PR and had a atol added to it. I don't remember which PR at the moment. The other one is...new. It shouldn't be effect by your changes or any other tests unless you are writing to something inplace from the fake readers...but even then that's odd.

[gerritholl]

In a comment to https://github.com/pytroll/satpy/pull/2181 @ameraner has pointed out that there is a bug in the FCI simulated 0.865 µm reflectivity and that the values should be multiplied by 0.8 prior to applying the composite.

When I do so, the FCI-based colorized composite looks better, but could still do with a cloud mask and tuning:

[gerritholl]

Failing in test_all_datasets_one_reader now...

[gerritholl]

Failing in test_all_datasets_one_reader now...

I have no idea why, but I can actually reproduce this failure locally. I will check.

[gerritholl]

Like for the true color RGB, here too now all tests pass except for the experimental run, which is allowed to fail. The colours will still need tweaking, but this should be done only after we have real measurements.

[gerritholl]

I tried to add a day_ version, but for some reason it looks different even in the day part, even though the enhancements are identical:

  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

and

  essl_low_level_moisture:
    name: essl_low_level_moisture
    operations: &moisture
      - 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:
    name: day_essl_low_level_moisture
    operations: *moisture

but the results look different:

essl_low_level_moisture:

day_essl_low_level_moisture:

I don't understand why they don't look the same. From the logs, both appear to get the configured enhancement:

[DEBUG: 2022-10-27 11:37:02 : satpy.writers] Data for DataID(name='day_essl_low_level_moisture', resolution=1000) will be enhanced with options:
        [{'name': 'linear_stretch', 'method': <function stretch at 0x7f1b63635d80>, 'kwargs': {'stretch': 'crude', 'min_stretch': 0.35, 'max_stretch': 0.85}}]

and

[DEBUG: 2022-10-27 11:37:02 : satpy.writers] Data for DataID(name='essl_low_level_moisture', resolution=1000) will be enhanced with options:
        [{'name': 'linear_stretch', 'method': <function stretch at 0x7f1b63635d80>, 'kwargs': {'stretch': 'crude', 'min_stretch': 0.35, 'max_stretch': 0.85}}]

[gerritholl]

I found out the problem. Satpy was applying the default enhancement. I had to add operations: [] and standard_name: ..., even though for other enhancements I'm just adding name: .... Not sure why, but it works now.

[gerritholl]

ESSL low level moisture:

Day only version:

ESSL colorized low level moisture:

Day only version:

NB: the background colour for the day-only versions is transparent, not black. Conversion to (classic) JPEG and visualisation in GitHub makes it black.

[gerritholl]

The failing tests

FAILED satpy/tests/reader_tests/test_modis_l1b.py::TestModisL1b::test_load_longitude_latitude[modis_l1b_nasa_mod021km_file-True-False-False-1000] - ValueError: Buffer dtype mismatch, expected 'float32_t' but got 'double'
FAILED satpy/tests/reader_tests/test_modis_l1b.py::TestModisL1b::test_load_longitude_latitude[modis_l1b_imapp_1000m_file-True-False-False-1000] - ValueError: Buffer dtype mismatch, expected 'float32_t' but got 'double'
FAILED satpy/tests/reader_tests/test_modis_l1b.py::TestModisL1b::test_load_longitude_latitude[modis_l1b_nasa_1km_mod03_files-True-True-True-250] - ValueError: Buffer dtype mismatch, expected 'float32_t' but got 'double'
FAILED satpy/tests/reader_tests/test_modis_l2.py::TestModisL2::test_load_longitude_latitude[modis_l2_nasa_mod35_file-True-False-False-1000] - ValueError: Buffer dtype mismatch, expected 'float32_t' but got 'double'

pass locally, and don't seem likely to be related to my changes.

[djhoese]

Ugh I must have broken something in the new geotiepoints. I'll take a look.

[gerritholl]

codebeat is confused (again), it says I have introduced problems in modules I haven't touched.