In night-time geo_color, mid-level clouds give false impression of clear-sky conditions

[gerritholl]

Describe the bug

For the night part, the geo_color composite uses a blending of geo_color_high_clouds and geo_color_low_clouds. For the high clouds, brightness temperature is shown as transparency using HighCloudCompositor / CloudCompositor, where 300 K are fully transparent (not a high cloud) and a latitude-dependent threshold (either 210 K or 230 K) are fully opaque / white.

When mid-level clouds cover low-level clouds, those mid-level clouds are shown mostly transparent. However, they are actually opaque, and it is unknowable whether low-level clouds exist underneath the mid-level clouds. By displaying opaque clouds as semi- or even nearly fully transparent, the user may be misled to believe there are clear-sky conditions.

An example can be seen on 2024-10-15 over Central Europe. At 04:20, low clouds cover valleys in German mountain valleys:

By 06:00, advection of mid-level clouds looks just like the fog is dissipating. Note that the mid-level clouds themselves are nearly invisible in this visualisation, but the city lights underneath are clearly shown.

After sunrise, at 07:20, we can see that actually, the low-level clouds didn't go anywhere, and are partially covered by higher clouds, now clearly visible in the true colour:

For additional views or videos, see eumetview.

To Reproduce

import hdf5plugin
import os
from satpy import Scene
from glob import glob
from sattools.io import plotdir
fci_files = {
    "good": glob("/media/nas/x23352/MTG/FCI/L1c-cases/20241015-hidden-fog/04/W_XX*BODY*C_0027*nc"),
    "bad": glob("/media/nas/x23352/MTG/FCI/L1c-cases/20241015-hidden-fog/06/W_XX*BODY*C_0037*nc"),
    "day": glob("/media/nas/x23352/MTG/FCI/L1c-cases/20241015-hidden-fog/08/W_XX*BODY*C_0052*nc")}
for fns in fci_files.values():
    sc = Scene(filenames={"fci_l1c_nc": fns})
    sc.load(["geo_color", "geo_color_high_clouds", "geo_color_low_clouds", "dwd_ir105"])
    ls = sc.resample("nqgerm1000m")
    ls.save_datasets(
        filename=os.fspath(plotdir() / "{start_time:%Y%m%d%H%M}-{area.area_id}-{name}.tif"))

Expected behaviour

I expect that any (image) product clearly shows the difference between clear sky conditions and uncertain conditions. The mid-level clouds are not undetected; in fact, looking at infrared imagery at 06:00, they are clearly detectable:

So this is a visualisation problem, raising the question if it really is a good idea to scale high clouds by transparency in this case.

Actual results

See above.

Environment Info:

• Satpy Version: 0.52.0

Additional context

I realise this problem is not easy to solve, but I think it is dangerous to present the product to forecasters in this way. One can of course warn about this problem in training, but that will get lost over time (forecasters have a lot of products to keep in mind), and in particular the city lights give the impression of clear-sky conditions where there are not (users need to remember that the city lights are a static background image). I expect users might falsely conclude that fog conditions have lifted.

[strandgren]

I definitely agree that we can improve the transparency of the high cloud layer to avoid having fully or almost fully transparent clouds. Something similar to the IR image you show would maybe be better. The same comment has been raised by Ivan who is using EUMETView, where the same recipe is used. so we could try to converge on an optimal transparency together with him, since he has a lot of experience with RGBs and also quite some interaction with users.

However, I'm actually not convinced that the masked low-level clouds is (only) because of the high-level clouds. I rather think that the IR3.8-IR10.5 test used to detect low-level clouds at night struggles at twilight, when we start to get solar contribution.. Perhaps we should try using the IR3.8 BT data corrected for the solar contribution, because we have something like that in Satpy, right?

I'm not sure if I agree with this being a bug, but yeah certainly room for improvement.

[gerritholl]

This is geo_color_night with various settings for transition_min_limits, while keeping transition_max at 300 and ɣ=1. The image shows Germany and surroundings.

210/230 (this is the default):

230/250:

260/260:

270/270:

We do have the [nir_emissive] modifier, but for geo_color_night / geo_color_low_clouds it fails with:

  File "/data/gholl/mambaforge/envs/py312/lib/python3.12/site-packages/pyspectral/near_infrared_reflectance.py", line 267, in reflectance_from_tbs
    self.derive_rad39_corr(tb_therm, tbco2)
  File "/data/gholl/mambaforge/envs/py312/lib/python3.12/site-packages/pyspectral/near_infrared_reflectance.py", line 159, in derive_rad39_corr
    self._rad3x_correction = (bt11 - 0.25 * (bt11 - bt13)) ** 4 / bt11 ** 4
                                             ~~~~~^~~~~~
  File "/data/gholl/mambaforge/envs/py312/lib/python3.12/site-packages/dask/array/core.py", line 215, in wrapper
    return f(self, other)
           ^^^^^^^^^^^^^^
  File "/data/gholl/mambaforge/envs/py312/lib/python3.12/site-packages/dask/array/core.py", line 2400, in __sub__
    return elemwise(operator.sub, self, other)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/data/gholl/mambaforge/envs/py312/lib/python3.12/site-packages/dask/array/core.py", line 4823, in elemwise
    broadcast_shapes(*shapes)
  File "/data/gholl/mambaforge/envs/py312/lib/python3.12/site-packages/dask/array/core.py", line 4751, in broadcast_shapes
    raise ValueError(
ValueError: operands could not be broadcast together with shapes (11136, 11136) (5568, 5568)

[gerritholl]

nir_emissive cannot be used when combining FDHSI and HRFI because of #2460. A workaround for #2460 is to use generate=False, but then geo_color runs into #2733.

[gerritholl]

Using only FDHSI to avoid #2460, applying nir_emissive when loading geo_color_low_clouds does not look great.

geo_color:

[strandgren]

Using only FDHSI to avoid #2460, applying nir_emissive when loading geo_color_low_clouds does not look great.

geo_color:

Yeah, not great.. What repeat cycle is this? Looks like we have sun-rise to the east? Because then I don't understand why we don't have any data at night? Seems like it's being masked somehow. Perhaps we need to change the sunz_threshold, but I'm also not very familiar with the modifier.

[strandgren]

This is geo_color_night with various settings for transition_min_limits, while keeping transition_max at 300 and ɣ=1. The image shows Germany and surroundings.

210/230 (this is the default):

230/250:

260/260:

270/270:

We do have the [nir_emissive] modifier, but for geo_color_night / geo_color_low_clouds it fails with:

  File "/data/gholl/mambaforge/envs/py312/lib/python3.12/site-packages/pyspectral/near_infrared_reflectance.py", line 267, in reflectance_from_tbs
    self.derive_rad39_corr(tb_therm, tbco2)
  File "/data/gholl/mambaforge/envs/py312/lib/python3.12/site-packages/pyspectral/near_infrared_reflectance.py", line 159, in derive_rad39_corr
    self._rad3x_correction = (bt11 - 0.25 * (bt11 - bt13)) ** 4 / bt11 ** 4
                                             ~~~~~^~~~~~
  File "/data/gholl/mambaforge/envs/py312/lib/python3.12/site-packages/dask/array/core.py", line 215, in wrapper
    return f(self, other)
           ^^^^^^^^^^^^^^
  File "/data/gholl/mambaforge/envs/py312/lib/python3.12/site-packages/dask/array/core.py", line 2400, in __sub__
    return elemwise(operator.sub, self, other)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/data/gholl/mambaforge/envs/py312/lib/python3.12/site-packages/dask/array/core.py", line 4823, in elemwise
    broadcast_shapes(*shapes)
  File "/data/gholl/mambaforge/envs/py312/lib/python3.12/site-packages/dask/array/core.py", line 4751, in broadcast_shapes
    raise ValueError(
ValueError: operands could not be broadcast together with shapes (11136, 11136) (5568, 5568)

Interesting. Definitely get's better at some point before the low-level clouds start to become grey as well. But we probably need to keep the the use of different values (e.g. 240/260), in order to have good visualization in the tropics as well.

[pnuu]

Using only FDHSI to avoid #2460, applying nir_emissive when loading geo_color_low_clouds does not look great.

I think the strict limit comes from the way the emissive part is computed via first getting the reflected portion and deriving the emissive part from the reflectances. The derivation of the reflectances requires a lot of masking in Pyspectral, one of which is with the SZA.

[gerritholl]

What repeat cycle is this?

2024-10-15 06:00Z (repeat cycle 37). The earlier one 04:20 (RC 27).

Definitely get's better at some point before the low-level clouds start to become grey as well.

Yes, this is the tradeoff.

But we probably need to keep the the use of different values (e.g. 240/260), in order to have good visualization in the tropics as well.

Yes. My experiments were limited to Central Europe.

[gerritholl]

That the nir_emissive modifier messes things up may not be the fault of geo_color_night. This is what the IR 3.8 µm channel looks like with the modifier applied:

[gerritholl]

If the problem is (in part) due to the emissive part of the 3,8 µm band, it might also help to change the lim_low and lim_high parameters passed to DayNightCompositor. Those are normally set to lim_low=78 and lim_high=88.

This is what we get with transition_min_limits = [240, 260], lim_low=86, lim_high=89 (still 2024-10-15 06:00):

Considering how much better the signal is by day as by night, it would seem we could try to use the day part closer to the terminal than SZA=78°.

[gerritholl]

For information: in our internal production, we have now changed the parameters:

Parameter	Old	New
lim_low	78	87
lim_high	88	89
transition_min_limits[0]	210	250
transition_min_limits[1]	230	260

Before:

After:

Note that this was optimised optically for a single Central European case. We have not tested it yet on other cases or regions.

[gerritholl]

In this geo colour example from 2024-10-28, 05:00-07:00, with my adapted settings, low clouds around and between Main and Danube are clearly shown. However, the low clouds in Styria are not only invisible due to higher clouds, but also (still) shown as clear-sky.

For comparison, the version on EUMETVIEW again gives the impression the low clouds are disappearing, before reappearing:

[gerritholl]

This case from 2024-11-11 15:30 shows that with adapted parameters (the same ones as commented three weeks ago), there is a rather dark band near the terminator. Over the North Sea, we can see high clouds that are still sunlit, but the surface is already dark. This puts us in a bit of a double bind. We want to use the solar information for the high clouds, because this has better spatial and spectral information than the thermal channels and the 3.8 µm channel is hard to use due to a mixed signal. But at the surface we would probably like to show city lights, as it's already dark there.

In this case, eumetview might look better: