Fix NaN handling in NDVIHybridGreen calculation

[djhoese]

While working on generating FCI true colors @kathys noticed that the night portion of the image was always transparent instead of black like other instruments' true colors. I tracked it down to the NDVI Hybrid Green correction. The "problem" comes down to NDVI having / 0 in some cases which in numpy results in NaNs. This results in the correction/contribution fraction being NaN and the resulting green pixel value being NaN. In this PR I do an additional "clip" on the NDVI to convert NaNs to the minimum allowed NDVI as configured by the user. Given that the current result is different than what the regular hybrid compositor produces (ex. AHI's current default, black night pixels) and how simple the fix is I think this is the right way to go.

I still need to add a test, but wanted to get people's opinion.

Old

New

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

[pnuu]

Apart from the excessive Dask computation, I'm fine with the change. As you said in Slack, there's the DayNightCompositor to make the night-side transparent.

[gerritholl]

The same happens in day_microphysics.

[djhoese]

The same happens in day_microphysics.

I assume that's due to the nir_reflectance on the green channel?

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 96.34%. Comparing base (6c571f2) to head (9928fec).
⚠️ Report is 100 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #3368      +/-   ##
==========================================
- Coverage   96.35%   96.34%   -0.02%     
==========================================
  Files         466      466              
  Lines       59083    59086       +3     
==========================================
- Hits        56931    56926       -5     
- Misses       2152     2160       +8     

Flag	Coverage Δ
behaviourtests	3.58% <0.00%> (-0.01%)	⬇️
unittests	96.43% <100.00%> (-0.02%)	⬇️

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

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[pnuu]

Yep, the calculation of reflected portion generally breaks down before the SZA reaches 90 degrees, and those areas are masked out. The cut-off isn't actually even strictly by any SZA limit, also the signal level matters and the edge is often very ragged.

[strandgren]

Thanks for fixing, I think this change in behavior makes sense! See two comments in-line

[strandgren]

You write in the description that you "do an additional 'clip' on the NDVI to convert NaNs to the minimum allowed NDVI as configured by the user", but for that you'd need to do ndvi.data = np.nan_to_num(ndvi.data, nan=self.ndvi_min), right?

That being said I don't think that's a wanted behavior since we'd then get something neither black nor transparent I think

[djhoese]

Ah yes, I should re-add that. I got some error from dask or numpy about passing nan=self.ndvi_min as a keyword argument so I removed it to get things working then forgot to put it back. For what I was going for I think it makes sense to have it be ndvi_min as we're saying "we want the smallest valid/legal value of NDVI to move on in this processing". If we think about it purely as correcting for the denominator of the NDVI being 0 then maybe we do want it to be ndim_max and not ndvi_min as X / almost-zero would be a large number.

[djhoese]

Thinking about this more, another way to put this is that we're controlling the fraction by controlling the NDVI. It isn't that we'd get something that isn't black and isn't transparent, we'd get something that is a blend of the two bands still and the common cases of NaNs in NDVI are NaNs in the data or a divide by 0 where both NDVI inputs are 0 at a pixel. So the end result will typically be 0 (black).

[strandgren]

yes, that's true so we're probably talking about something that doesn't affect the final output. That being said, setting nan=self.ndvi_min would correspond to theoretical logic of the current behaviour I think 😄

[djhoese]

I can do that. Want thought experiment though, let's say the two NDVI are near zero but one negative and one positive. So p2 is 0.01 and p1 is -0.009. I know reflectances should be negative but they happen in the data. That would make ndvi = (0.019 / 0.001) which gives us 19.0 on a product that should be 0 to 1. If we flip these values then we get ndvi = (-0.019 / 0.001) I think so -19.0. So either ndvi min or max would make sense with this made up example.

I see your documentation says:

    In this example, pixels with NDVI=0.0 will be a weighted average with 15% contribution from the
    near-infrared vis_08 channel and the remaining 85% from the native green vis_05 channel, whereas
    pixels with NDVI=1.0 will be a weighted average with 5% contribution from the near-infrared
    vis_08 channel and the remaining 95% from the native green vis_05 channel. For other values of
    NDVI a linear interpolation between these values will be performed.

Is the 85%/15% case (NDVI min) the desired outcome or is the 95%/5% (NDVI max) what we want?

[strandgren]

I don't get what you mean. Are p1 and p2 the reflectance data used to compute ndvi, or something else? I don't see the connection between So p2 is 0.01 and p1 is -0.009 and ndvi = (0.019 / 0.001)

Is the 85%/15% case (NDVI min) the desired outcome or is the 95%/5% (NDVI max) what we want?

We want to map the lowest valid NDVI value to 85%/15% (i.e. limits[0]) and the highest valid NDVI to 95%/5% (limits[1]). So if the user says that the smallest NDVI should be 0.1 (ndvi_min=0.1), an NDVI of 0.1 should be mapped to 85%/15% (limits[0]) and an NDVI of 1.0 to 95%/5% (limits[1]) as before. Hence, I think it makes sense to set NaN values to ndvi_min. Whatever the final implementation is, I'd still be happy if you can test and verify that the night-side is still black🙂

[djhoese]

Yes, sorry. p2 is projectables[2] and p1 is projectables[1] from the compositor code. So plugging in the numbers I gave into the NDVI formula of (p2 - p1) / (p2 + p1) would give us the values I showed if I did my math correctly. Regardless...I'll set it to ndvi_min and make an image.

[strandgren]

Thanks for the clarification. The theoretical range of NDVI is [-1, 1], but indeed, if the reflectance values (p1 and p2) are very small we can get significantly higher/lower values. Below is an example including the distribution of the NDVI values (log-scale):

For the sake of the green band correction I decided to limit the NDVI to [0,1] here and it's mainly clear-sky water which has NDVI < 0.0. Furthermore, we now clip the data to the NDVI limits, so with the default configuration there will never be an NDVI outside the range [0,1]. If any user decides to modify this to allow NDVI values up to e.g. +-19 they will have to deal with the consequences😄

[gerritholl]

Our operational true color imagery is black at night.

MTG_EUROPA_ZENTRAL_true_color_nqcxeur750m_202604060440-geotiff.tif

import hdf5plugin
from satpy import Scene
from glob import glob
fci_files = glob("/media/nas/x23352/MTG/FCI/L1c-cases/202309_10-cwg/10/01/08/W_XX-EUMETSAT-Darmstadt,IMG+SAT,MTI1+FCI-1C-RRAD-*FDHSI*-FD--CHK-BODY---NC4E_C_EUMT_20231001*_IDPFI_VAL_20231001*_20231001*_N__C_0054_*.nc")
sc = Scene(filenames={"fci_l1c_nc": fci_files})
sc.load(["true_color"])
sc.save_datasets(fill_value=0)

In fact, my attempts to set the night part transparent using the DayNightCompositor have failed, see #3003

[gerritholl]

The same happens in day_microphysics.

I assume that's due to the nir_reflectance on the green channel?

Yes. In fact, as a workaround for #3003, I've actually experimented with a DayNightCompositor setting the night part to day_microphysics, which works, but it's ugly and relying on that is no better than relying on spacebar heating… I'm planning to work on a solution for #3003.

[djhoese]

@gerritholl what is your composite and enhancement configuration for your operational true color?

[pnuu]

Our operational true color imagery is black at night.

That's what I'd expect with

sc.save_datasets(fill_value=0)

Right?

[djhoese]

Oh yes, I missed that in @gerritholl's code. That is expected as all the NaNs are converted to 0s (the fill value you specified). This PR and this problem is specifically about wanting the space pixels to be NaN (invalid) and the night pixels to be 0 or some other minimal value or more accurately I want them to be whatever the "real" true color would see (a real red channel, a real green channel, a real blue channel)...so black.

[gerritholl]

Our operational true color imagery is black at night.

That's what I'd expect with

sc.save_datasets(fill_value=0)

Right?

No, that's not what I would expect. Our intention is for the night pixels to become transparent, but satpy "helpfully" sets the night pixels to 1 so that they turn out black, defeating the purpose of them being 0 before (either due to the lack of light, or due to the DayNightCompositor setting them to 0). This is a long-standing bug in our operational production that is on my list of things to fix in 2026.

1. User sets pixels to 0.
2. User passes fill_value=0 to save_datasets, so that those 0 pixels become transparent.

Step 2 fails. But that problem is covered by #3003.

[gerritholl]

@gerritholl what is your composite and enhancement configuration for your operational true color?

We use the built-in satpy composite and enhancement true_color.

[pnuu]

I've always considered fill_value=X to mean "set all missing values to X". Missing data are (in most cases at least) set to NaN by the file handler / reader. So setting them to zero would give black. There might be differences how the writers handle the kwarg and whether they add the metadata to the final file or not, and in the end it's up to the image display software the show the fill values as transparency, like they would be with the alpha band.

[djhoese]

@gerritholl These are two or more separate issues and you have #3003 for the specific DayNightCompositor issue of trying to get the night pixels to be a certain value. For the save_datasets with fill_value=0, at least for the geotiff and simple_image writer, the NaNs will always be made 0 and all of the valid data (including 0) will always be 1-255 for a uint8 geotiff. The entire valid data space is scaled to 1-255. There is no way to specify fill_value=0 to the geotiff writer and have NaNs appear in the final out. That is/was the goal of the fill_value kwarg.

I mentioned this in the meeting yesterday, but I think this PR is correct and we should generally have psuedo-bands try to use a reasonable value where pixels aren't strictly invalid (0-ish in this case). Theoretically the DayNightCompositor can be used to convert night pixels to NaN but #3003 needs to be figured out. Alternatively there is no way for me to get black night pixels and NaN space pixels with the current way this NDVIHybridGreen compositor works or at least I think it is harder.

The question for this PR still remains from the previous comment thread of what should the "fraction" be based on for night pixels where NDVI is invalid? NDVI min or NDVI max or handle the NaNs in a different portion of the compositor?

[djhoese]

With the latest commit using ndvi_min as the replacement for NaNs the image looks like this with default Satpy enhancements:

[djhoese]

And for @kathys here is with the P2G enhancement:

which is just a slightly brighter version and shows more of the cloud tops I think.

[kathys]

@djhoese Looks good to me.

[strandgren]

Thanks for the images, looks good to me!

[djhoese]

Thanks. Any other comments before merging?