Multiscene blend with weights

[adybbroe]

This PR aims adding a way to blend several scenes (polar and geo cloud parameters is what we have in mind) using pixel weights. Weights could for instance simply be the satellite zenith angle (or rather 1/satz) so that the pixel with the smallest satellite zenith angle is chosen from all scenes being blended. Weights could also be a comlex mixture of several parameters, like time, satellite zenith, retrieval uncertainty (if available) and so on.

• Closes #xxxx
• Tests added
• Fully documented

Example/pseudo code creating a multi scene of one geo (SEVIRI) and two polar (here AVHRR/3 from N18 and N19) remapping to a pre-defined area, and blending the NWCSAF cloudtype datasets by stacking with weighs. Here the weights are the inverse of the satellite-zenith angles:

    mscn = MultiScene([geo, polar1, polar2])
    groups = {DataQuery(name='CTY_group'): ['ct']}
    mscn.group(groups)

    resampled = mscn.resample(areaid, reduce_data=False)

    weights = [1./geo_satz_angles, 1./polar1_satz_angles, 1./polar2_satz_angles]
    blended = resampled.blend(blend_function=stack_weighted, weights=weights)

    # blended.show('CTY_group')
    blended.save_dataset('CTY_group', filename='./blended_stack_weighted_geo_polar_2polars.nc')

[adybbroe]

@djhoese Just a question. The stack function retains the attributes from the first dataset, but wouldn't it more intuitive to store the attributes from the last one in the stack, the on on top? I realize that can be debated, so if you had a particular reason for keeping the first I am absolutely fine. Just need to be clear what we do when creating/adapting the unittests. Don't think this aspect was tested before...

[codecov]

Codecov Report

Merging #2275 (8baca79) into main (9a496eb) will increase coverage by 0.03%.
The diff coverage is 100.00%.

@@            Coverage Diff             @@
##             main    #2275      +/-   ##
==========================================
+ Coverage   94.64%   94.67%   +0.03%     
==========================================
  Files         324      328       +4     
  Lines       48339    48544     +205     
==========================================
+ Hits        45749    45960     +211     
+ Misses       2590     2584       -6     

Flag	Coverage Δ
behaviourtests	4.43% <1.42%> (-0.02%)	⬇️
unittests	95.30% <100.00%> (+0.03%)	⬆️

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

Impacted Files	Coverage Δ
satpy/modifiers/angles.py	99.16% <ø> (ø)
satpy/multiscene.py	94.17% <100.00%> (+2.27%)	⬆️
satpy/tests/multiscene_tests/__init__.py	100.00% <100.00%> (ø)
satpy/tests/multiscene_tests/test_blend.py	100.00% <100.00%> (ø)
satpy/tests/multiscene_tests/test_misc.py	100.00% <100.00%> (ø)
...atpy/tests/multiscene_tests/test_save_animation.py	99.09% <100.00%> (ø)
satpy/tests/multiscene_tests/test_utils.py	100.00% <100.00%> (ø)

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

[simonrp84]

@adybbroe This old PR of mine might be relevant here: #1157

The purpose of the PR was to compute the pixel size to enable blending of various datasets based on whichever had the smallest pixel size.

I chose that method as, for example, in some cases GOES/ABI has a larger zenith angle but smaller pixel size than SEVIRI as ABI has a higher native pixel resolution.

[coveralls]

Coverage: 95.263% (+0.03%) from 95.231% when pulling 8baca79 on adybbroe:multiscene_blend_with_weights into 9a496eb on pytroll:main.

[adybbroe]

@adybbroe This old PR of mine might be relevant here: #1157

The purpose of the PR was to compute the pixel size to enable blending of various datasets based on whichever had the smallest pixel size.

I chose that method as, for example, in some cases GOES/ABI has a larger zenith angle but smaller pixel size than SEVIRI as ABI has a higher native pixel resolution.

Oh, that's interesting! Didn't quite thought of that, might be useful also when blending AVHRR, MODIS and VIIRS (and eventually MetImage) for instance. I will have a look, but that PR was never merged, so work is "half" finished right?

[djhoese]

A couple questions inline. Otherwise, regarding the attrs from the first DataArray being used...no probably no real reason that I can think of. I'm not sure I'm even the one who added that function originally. I could see there being a chance of performance differences in that I think the first Scene is always generated first (if Scenes were created via a generator), but I think once this blend function is called it is already a list (maybe?) so it might not matter. I'm open to changing it.

[simonrp84]

I will have a look, but that PR was never merged, so work is "half" finished right?

It was pretty much finished in terms of functionality, but I didn't persue it as Dave said it might not be suitable for satpy and I no longer needed the PR itself for my project. It's something I'd be interested in picking up though!

[adybbroe]

Here a use case, with one polar and one geo scene. Looking at the NWCSAF Cloud Type dataset.
First the sat-zenith angles:

And then the individual datasets before blending:

And then the stacked and the stacked using weights. The weights are the inverse of the satellite zenith angle, so 1/satz.
The stacked:

And the weighted:

Maybe a bit hard to see in these images but the difference is in the middle lower part of the image where the satellite zenith angles of the polar scene are higher than the corresponding angles of the Geo scene. The Geo Cloud Type gets weighted in more there. Here the difference between the two blended composites show this more clearly:

[mraspaud]

Looks good!

[adybbroe]

And here two more images, now adding another polar scene, so two polar and one geo. Here the satz angles for the 2nd polar scene (NOAA-19) and the difference between the weighted composite of three scenes and the previous (with only one geo and one polar):

One can see that the new polar scene in the north complements the previous polar adding where it has better viewing compared to previous.

[djhoese]

I'm not sure I like the final stacked result of your first example better, but I see the usefulness. It is kind of amazing how much can be done just by throwing on some weights and/or masks.

[adybbroe]

@djhoese test_multiscene is now a bit long, CodeScene health check complains it is more than 600 lines long. I am sure my additions could be improved, but think we might have to think of splitting it in two. When testing off line all the mp4 tests complaint though I didn't touch them, guess it is some 3rd library I am missing? I ended up running: pytest -s satpy/tests/test_multiscene.py -k 'test_ and not mp4'. Yes, I am yet/still not using Pycharm... ;-)
But could we move the mp4 and timeseries stuff to a separate test file perhaps?

[djhoese]

When testing off line all the mp4 tests complaint though I didn't touch them, guess it is some 3rd library I am missing?

If it is the ones talking about "distributed" then it is probably the dask distributed library. Otherwise it is imageio and imageio-ffmpeg that you are missing.

I am OK splitting the tests if you make it a test_multiscene/ directory and put tests in test_save_animation.py and test_blend.py and test_misc.py or something like that.

[lobsiger]

@adybbroe @djhoese just for the record: I'm very interested in a kind of blending feature as developed and discussed here. In my version 4.0 distribution of PyTROLL/Satpy scripts for EUMETCast I have multi pass scripts for all LEO satellites available and even scripts for combined NOAA20/SNPP, MetopB/MetopC and Sen3A/Sen3B (ERR). That said being able to really blend and not just stack the data would be most welcome. See my latest benchmark here: https://groups.io/g/MSG-1/message/33860

[adybbroe]

When testing off line all the mp4 tests complaint though I didn't touch them, guess it is some 3rd library I am missing?

If it is the ones talking about "distributed" then it is probably the dask distributed library. Otherwise it is imageio and imageio-ffmpeg that you are missing.

I am OK splitting the tests if you make it a test_multiscene/ directory and put tests in test_save_animation.py and test_blend.py and test_misc.py or something like that.

I forgot this, should I try split the test functions as suggested. I can have a go, if that's all which you think is left?

[adybbroe]

@djhoese I have done a bit around the attributes and some minor refactoring and documentation. I have a feeling you will have an opinion. Maybe my new function should be get_combined_start_end_times should be in the metadata.py module?
And actually rather be built into the combine_metadata function?
So, not sure how to finish this nicely? Should I drop it for now perhaps? But it does make sense to make the combination of start and end times as I do I think. I thought it was more intuitive than an average in this case?

Also, I think my section in the doc pages should be improved, but have to run for now.

Otherwise I perhaps only have left to settle for the standard_name as you mentioned above. Just need to make some system testing first.

[djhoese]

Otherwise I perhaps only have left to settle for the standard_name as you mentioned above. Just need to make some system testing first.

Are you saying you're not sure if the standard_name changes broke something?

[djhoese]

Just some small changes. Otherwise this looks good assuming the standard_name thing isn't an issue.

[adybbroe]

@djhoese Perhaps you saw I reverted back on the standard_name? So, unless you have suggestions on the documentation, I think we could be ready to merge this?
Don't know if someone else wants to have a look?

[djhoese]

Merge it and let's see what breaks. I'll let you click the button.