MultiScene blending with weights

86 views
Skip to first unread message

lobsiger...@gmail.com

unread,
Jan 31, 2023, 4:25:44 PM1/31/23
to pytroll
Dear developers,


now that PR 2275 has been merged to main I wonder whether the
new functionality can be used for my LEO multi pass scripts.
TBH I have no clue how I get the weights list. My pseudo code:


Find the LEO files that possibly contribute to 'area' ascending
or descending at day X. Sort them to N passes from East to West.

from satpy import Scene, MultiScene
...
scenes = []
for N passes:
    scenes.append(Scene(filenames=passfiles, reader=reader)
mscn = MultiScene(scenes)
mscn.load(composites, generate=generate)
new_mscn = mscn.resample(area, **resampler_kwargs)
blended_scn = new_mscn.blend()
blended_scn.save_dataset(...)


I understand I need a weights list with N (xarrays?) entries:
weights = [ 1./(C + satellite_zenith_angle), ...] with C >= 1.
This should give the pixels at SSP (nadir) the biggest weight.


There are different stages of difficulty:

1) Multipass of one satellite, all channels have the same resolution (avhrr).
2) Multipass of one satellite, different sorts of channels (MXX, DNB, viirs)
3) Multipass of two identical satellites (MetopB, MetopC, MetopB, ..., avhrr)
4) Multipass of two identical satellites (NOAA20, SNPP, NOAA20, ..., viirs)


Now it seems to me that for 1) weighted I need something like:


from satpy import Scene, MultiScene
from functools import partial
...
scenes = []
weights = []
for N passes:
    scenes.append(Scene(filenames=passfiles, reader=reader)
    weights.append(1./(C + satellite_zenith_angle)
mscn = MultiScene(scenes)
mscn.load(composites, generate=generate)
new_mscn = mscn.resample(area, **resampler_kwargs)
stack_with_weights = partial(stack, weights=weights)
blended_scn = new_mscn.blend(blend_function=stack_with_weights)
blended_scene.save_dataset(...)

For viirs I might have to group channels? But let's
take the most simple (??) case 1) first. How do I
get the weights list? What's wrong with the rest?


Best regards,
Ernst

Aqua-20230130-DAY-1247-overview-eurol-multi.jpg
NOAAX-20230130-NIG-0229-night_overview-westminster-multi.jpg

Adam Dybbroe

unread,
Feb 1, 2023, 3:17:51 AM2/1/23
to pyt...@googlegroups.com, adam.d...@smhi.se

Ernst,

Some quick comments inline below.

You did of course consult the documentation, right?

https://satpy.readthedocs.io/en/latest/multiscene.html#stacking-scenes-using-weights

:-)

On 31/01/2023 22:25, lobsiger...@gmail.com wrote:
Dear developers,


now that PR 2275 has been merged to main I wonder whether the
new functionality can be used for my LEO multi pass scripts.
TBH I have no clue how I get the weights list. My pseudo code:


Find the LEO files that possibly contribute to 'area' ascending
or descending at day X. Sort them to N passes from East to West.

from satpy import Scene, MultiScene
...
scenes = []
for N passes:
    scenes.append(Scene(filenames=passfiles, reader=reader)
mscn = MultiScene(scenes)
mscn.load(composites, generate=generate)
new_mscn = mscn.resample(area, **resampler_kwargs)
blended_scn = new_mscn.blend()
blended_scn.save_dataset(...)

You are missing the grouping part, see documentation.

I understand I need a weights list with N (xarrays?) entries:
weights = [ 1./(C + satellite_zenith_angle), ...] with C >= 1.
This should give the pixels at SSP (nadir) the biggest weight.

Yes, it will be a list of XArrays with weights, each with the same dimension NxM as the resampled scenes. So the weights should already be remapped to the area.

Why do you need the "C"?


There are different stages of difficulty:

1) Multipass of one satellite, all channels have the same resolution (avhrr).
2) Multipass of one satellite, different sorts of channels (MXX, DNB, viirs)
3) Multipass of two identical satellites (MetopB, MetopC, MetopB, ..., avhrr)
4) Multipass of two identical satellites (NOAA20, SNPP, NOAA20, ..., viirs)

These are interesting and useful use cases, so I am interested in seeing your progress!


Now it seems to me that for 1) weighted I need something like:


from satpy import Scene, MultiScene
from functools import partial
...
scenes = []
weights = []
for N passes:
    scenes.append(Scene(filenames=passfiles, reader=reader)
    weights.append(1./(C + satellite_zenith_angle)
mscn = MultiScene(scenes)
mscn.load(composites, generate=generate)
new_mscn = mscn.resample(area, **resampler_kwargs)
stack_with_weights = partial(stack, weights=weights)
blended_scn = new_mscn.blend(blend_function=stack_with_weights)
blended_scene.save_dataset(...)

Again look at the example at RTD, you need to tell Satpy which datasets to group together.


For viirs I might have to group channels? But let's
take the most simple (??) case 1) first. How do I
get the weights list? What's wrong with the rest?


I haven't tried other cases than blending level 2 cloud parameters (integers) so have no idea how it works for radiances/Tbs, but I suppose and hope it does.

-Adam


Best regards,
Ernst

--
You received this message because you are subscribed to the Google Groups "pytroll" group.
To unsubscribe from this group and stop receiving emails from it, send an email to pytroll+u...@googlegroups.com.
To view this discussion on the web, visit https://groups.google.com/d/msgid/pytroll/cf566d4c-2306-4d22-8968-928d6ab37e6an%40googlegroups.com.
-- 
Adam Dybbroe,
Satellite Remote Sensing Scientist
Research lead - Nowcasting and Remote Sensing
Meteorological Research unit, Research and Development
Swedish Meteorological and Hydrological Institute (SMHI)
www.pytroll.org
nwcsaf.smhi.se
www.smhi.se

lobsiger...@gmail.com

unread,
Feb 1, 2023, 4:19:45 AM2/1/23
to pytroll
Adam,

Thanks for your quick reply. I'll wait doing real
coding until 2275 is available in a conda package.

Yes I have RTFM but this currently goes over my head.

The upper part of my post shows the pseudo code of what
I use now. This code is working (see examples attached).
The lower part is first ideas for the weighting pseudo code.

Yes I intend to blend radiances. If I'd start with an
instrument that has the same resolution in all channels
simply I assumed I could skip the grouping part ;-).
I will always blend data from the the same (or kind of)
satellite only. Of course there is a time delay of 100
minutes between passes or 50 minutes for combined sats.

I see from .../satpy/etc/readers that a datasets for
satellite_zenith_angle generally are available. In the
case of the viirs instrument with MXX and DNB channels
and different resolutions there are even two datasets
satellite_zenith_angle and satellite_zenith_angle_dnb.
Also scn.available_dataset_names() (or IDs) show
 that 'satellite_zenith_angle' must be "available".

AGAIN: Above names are just IDs (I guess) and I have
no clue how to get these datasets into list weights=[].

I have used the C>=1 in 1./(C + satellite_zenith_angle)
as a parameter and to avoid a possible division by zero.


Regards,
Ernst

lobsiger...@gmail.com

unread,
Feb 2, 2023, 4:26:24 AM2/2/23
to pytroll

Dear developers,

I can easily add a use case 0) with (hopefully) lowest difficulty:

0) Multipass of polar Metop-B, 'CT', 'CTTH' or 'CMA' only (avhrr)

Below is what I find in the documentation. I still not fully understand
why I should use the grouping code as my satpy script only loads the
'CT', 'CTTH' and 'CMA' files for Metop-B anyway (see attached output).

But that's not the main difficulty I see. My question is again the same:

The weights list as documented uses 'n18_satz' which in my understanding
is an xarry containing the 'satellite_zenith_angle' dataset of NOAA-18.
Where does this 'n18_satz' xarray come from? How is this generated ?

I understand that I must generate such xarrays for all three passes
of the Metop-B multi pass (but only stacked) example I attach below.
Blending could improve the attached image between Iceland and Norway.


Cheers,
Ernst


<cite>
from satpy import Scene, MultiScene,  DataQuery
from functools import partial
from satpy.resample import get_area_def
areaid = get_area_def("myarea")
geo_scene = Scene(filenames=glob('/data/to/nwcsaf/geo/files/*nc'), reader='nwcsaf-geo')
geo_scene.load(['ct'])
polar_scene = Scene(filenames=glob('/data/to/nwcsaf/pps/noaa18/files/*nc'), reader='nwcsaf-pps_nc')
polar_scene.load(['cma', 'ct'])
mscn = MultiScene([geo_scene, polar_scene])
groups = {DataQuery(name='CTY_group'): ['ct']}
mscn.group(groups)
resampled = mscn.resample(areaid, reduce_data=False)
weights = [1./geo_satz, 1./n18_satz]
stack_with_weights = partial(stack, weights=weights)
blended = resampled.blend(blend_function=stack_with_weights)
blended_scene.save_dataset('CTY_group', filename='./blended_stack_weighted_geo_polar.nc')
</cite>
CT_Metop-B_output.txt
cloudtype.png
MetopB-20230201-NIG-2105-cloudtype-eurol-multi.jpg

Adam Dybbroe

unread,
Feb 2, 2023, 4:58:55 AM2/2/23
to pyt...@googlegroups.com, adam.d...@smhi.se

Ernst,

To get the angles I am using this one for SEVIRI:

from satpy.modifiers.angles import get_satellite_zenith_angle

And here is the methods I have in a small PPS Scene loader class:

    def prepare_satz_angles_on(self):
        """Derive the satellite zenith angles and attach data to Satpy scene object."""
        self.scene['satz'] = self._get_satz_angles()
        self.scene['satz'].attrs['area'] = self.scene['cma'].attrs['area']

    def _get_satz_angles(self):
        """Calculate the satellite zenith angles using Pyorbital."""
        self._shape = self.scene['cma'].shape

        satname = self.scene['cma'].attrs['platform_name']
        orb = Orbital(satname)

        time_step = (self.scene['cma'].attrs['end_time'] -
                     self.scene['cma'].attrs['start_time'])/self._shape[0]
        starttime = self.scene['cma'].attrs['start_time']
        obs_time = starttime + time_step * np.arange(self._shape[0])
        obs_time = np.dstack([obs_time]*self._shape[1])[0]

        lon, lat = self.scene['cma'].attrs['area'].get_lonlats()
        LOG.debug('Get satellite elevation via Pyorbital')
        _, elevation = orb.get_observer_look(obs_time, lon, lat, 0.)
        del _
        sat_zen = 90. - elevation
        LOG.debug('Satellite zenith angles derived on data')

        satz = xr.DataArray(data=sat_zen, dims=["y", "x"],
                            coords=dict(
                                lon=(["y", "x"], lon),
                                lat=(["y", "x"], lat)),
                            attrs=dict(
                                description="Satellite zenith angle.",
                                units="deg",
        ),)

        return satz

So, when using the Multiscene I use the methods above to derive the angles and attach them to the scene object like this:

    n18 = PPSCloudProductsLoader(N18_FILES)
    n18.load()
    n18.prepare_satz_angles_on()

And when I do the resampling on the multi-scene I then also have the reampled sat-zenith angles for the polar data:

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

    local_scn = resampled.scenes[0]['ct']
    geo_satz = get_satellite_zenith_angle(local_scn)
    n18_satz = resampled.scenes[1]['satz']

For the Geo part I use the Satpy function above to directly get the angles on the area after the (cloud) scene as been resampled to the area.


Does this make sense?


-Adam

lobsiger...@gmail.com

unread,
Feb 2, 2023, 9:53:04 AM2/2/23
to pytroll
Adam,

thanks a lot for this deep insight. It should be easier for me as 'satellite_zenith_angle' is part of the LEO files I use while for N18 you only have 'longitude' and 'latitude'.
My use of pyorbital is only (rather approximate) for selecting the right segment files for passes over specific areas. I'll come back and report if I make some progress.

Best regards,
Ernst

lobsiger...@gmail.com

unread,
Feb 4, 2023, 3:52:21 AM2/4/23
to pytroll
Hi All,

I have updated my satpy 0.39 manually by just overwriting multiscene.py with the latest version of Adam (PR 2275).
It's my understanding that all other file changes are not relevant for  what I intend to do blending AVHRR radiances.


Everything with my multi pass images seems to work as before. When I want to apply weights I get an error with partial:

  File "/home/eumetcast/SPStools/DEVscripts/LEOstuff.py", line 1048, in multi_single_pass
    stack_with_weights = partial(stack, weights=weights)
NameError: name 'stack' is not defined

Function stack() is in multiscene.py as before.  No more info with debugging on. Any help appreciated.

Best regards,
Ernst

David Hoese

unread,
Feb 4, 2023, 7:08:07 AM2/4/23
to pyt...@googlegroups.com
Ernst,

That error seems to be coming from your own script. You would need to import "stack" from satpy/multiscene.py to use it there. Put another way, that "stack" is a function you need to import, not a string.

Dave
>> blended_scene.save_dataset('CTY_group', filename='./blended_stack_weighted_geo_polar.nc <http://blended_stack_weighted_geo_polar.nc>')
>> https://satpy.readthedocs.io/en/latest/multiscene.html#stacking-scenes-using-weights <https://satpy.readthedocs.io/en/latest/multiscene.html#stacking-scenes-using-weights>
>>
>> :-)
>>
>> On 31/01/2023 22:25, lobsiger...@gmail.com wrote:
>>> Dear developers,
>>>
>>> https://github.com/pytroll/satpy/pull/2275 <https://github.com/pytroll/satpy/pull/2275>
>>> To view this discussion on the web, visit https://groups.google.com/d/msgid/pytroll/cf566d4c-2306-4d22-8968-928d6ab37e6an%40googlegroups.com <https://groups.google.com/d/msgid/pytroll/cf566d4c-2306-4d22-8968-928d6ab37e6an%40googlegroups.com?utm_medium=email&utm_source=footer>.
>>
>> --
>> Adam Dybbroe,
>> Satellite Remote Sensing Scientist
>> Research lead - Nowcasting and Remote Sensing
>> Meteorological Research unit, Research and Development
>> Swedish Meteorological and Hydrological Institute (SMHI)
>> www.pytroll.org <http://www.pytroll.org>
>> nwcsaf.smhi.se <http://nwcsaf.smhi.se>
>> www.smhi.se <http://www.smhi.se>
>>
>> --
>> You received this message because you are subscribed to the Google Groups "pytroll" group.
>> To unsubscribe from this group and stop receiving emails from it, send an email to pytroll+u...@googlegroups.com.
>> To view this discussion on the web, visit https://groups.google.com/d/msgid/pytroll/01ec262d-ae84-46ee-8416-edb8d2d9f669n%40googlegroups.com <https://groups.google.com/d/msgid/pytroll/01ec262d-ae84-46ee-8416-edb8d2d9f669n%40googlegroups.com?utm_medium=email&utm_source=footer>.
>
> --
> Adam Dybbroe,
> Satellite Remote Sensing Scientist
> Research lead - Nowcasting and Remote Sensing
> Meteorological Research unit, Research and Development
> Swedish Meteorological and Hydrological Institute (SMHI)
> www.pytroll.org <http://www.pytroll.org>
> nwcsaf.smhi.se <http://nwcsaf.smhi.se>
> www.smhi.se <http://www.smhi.se>
>
> --
> You received this message because you are subscribed to the Google Groups "pytroll" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to pytroll+u...@googlegroups.com <mailto:pytroll+u...@googlegroups.com>.
> To view this discussion on the web, visit https://groups.google.com/d/msgid/pytroll/e09abc61-2072-44ee-821d-2f85e4c44132n%40googlegroups.com <https://groups.google.com/d/msgid/pytroll/e09abc61-2072-44ee-821d-2f85e4c44132n%40googlegroups.com?utm_medium=email&utm_source=footer>.

lobsiger...@gmail.com

unread,
Feb 4, 2023, 7:25:09 AM2/4/23
to pytroll
Hello Dave,

thanks for that. I figured it out 5 minutes after my last post.
It's not in the documentation. So far so good. I added on top

from satpy.multiscene import stack

but now it seems that something is wrong with my weights list.
I have little knowledge of numpy and no clue regarding xarr ;-).

I have figured out two methods to generate weights lists.
Both generate different errors when applied in partial():

Method 0) stacking works as before, new_mscn is a multi scene
of three Metop-B AVHRR passes, resampled to area 'eurol' with
composites = ['overview', 'satellite_zenith_angle']. It can
be seen in attached images that 's_z_a' is resampled properly.
This even works with both generate=True and generate=False:

        new_scn = new_mscn.blend()

Method 1) see output1.txt attached:

        weights=[]
        for s in new_mscn.scenes:
            weights.append(1. / (1. + s.to_xarray_dataset(['satellite_zenith_angle'])))
        print(weights)
        stack_with_weights = partial(stack, weights=weights)
        new_scn = new_mscn.blend(blend_function=stack_with_weights)

Method 2) see output2.txt attached:

        weights=[]
        for s in new_mscn.scenes:
            weights.append(1. / (1. + s['satellite_zenith_angle']))
        print(weights)
        stack_with_weights = partial(stack, weights=weights)
        new_scn = new_mscn.blend(blend_function=stack_with_weights)


What do I have to change to get an accepted weights list?
I would also prefer 32bit floating point for these weights.


Best regards,
Ernst
output2.txt
output0.txt
MetopB-20230203-DAY-1038-overview-eurol-multi.jpg
MetopB-20230203-DAY-1038-satellite_zenith_angle-eurol-multi.jpg
output1.txt

lobsiger...@gmail.com

unread,
Feb 4, 2023, 11:38:38 AM2/4/23
to pytroll
Adam, Dave and Martin,
Dear core developers,

Method 1) with added ".to_array" (as suggested in output1.txt) gives exactly the same error like Method 2). So I went on with weights from Method 2) only.

I found out that the new weighting works with single channels (single channel "composites" like 'ir108_3d') only including 'satellite_zenith_angle' :-).

Improvement by weights looks not yet really breath taking. Maybe the weighting can still be improved. Script execution times are rather high as expected.

Multi channel composites end with errors (output2.txt). These are now definitely outside my possibilities.

Cheers,
Ernst


    new_scn = new_mscn.blend(blend_function=stack_with_weights)
  File "/home/eumetcast/miniconda3/envs/pytroll/lib/python3.10/site-packages/satpy/multiscene.py", line 423, in blend   new_scn[ds_id] = blend_function(datasets)
  File "/home/eumetcast/miniconda3/envs/pytroll/lib/python3.10/site-packages/satpy/multiscene.py", line 60, in stack   return _stack_weighted(datasets, weights, combine_times)
  File "/home/eumetcast/miniconda3/envs/pytroll/lib/python3.10/site-packages/satpy/multiscene.py", line 84, in _stack_weighted   weighted_array = xr.DataArray(da.choose(indices, datasets), dims=dims, attrs=attrs)
  File "/home/eumetcast/miniconda3/envs/pytroll/lib/python3.10/site-packages/dask/array/routines.py", line 2048, in choose    return elemwise(variadic_choose, a, *choices)
  File "/home/eumetcast/miniconda3/envs/pytroll/lib/python3.10/site-packages/dask/array/core.py", line 4762, in elemwise  broadcast_shapes(*shapes)
  File "/home/eumetcast/miniconda3/envs/pytroll/lib/python3.10/site-packages/dask/array/core.py", line 4690, in broadcast_shapes   raise ValueError(
ValueError: operands could not be broadcast together with shapes (3, 2048) (3, 2048, 2560) (3, 2048, 2560) (3, 2048, 2560)

MetopB-20230204-DAY-1017-1-eurol-multi.jpgMetopB-20230204-DAY-1017-1-eurol-multi-weighted.jpg

MetopB-20230204-DAY-1017-2-eurol-multi.jpgMetopB-20230204-DAY-1017-2-eurol-multi-weighted.jpg

MetopB-20230204-DAY-1017-ir108_3d-eurol-multi.jpgMetopB-20230204-DAY-1017-ir108_3d-eurol-multi-weighted.jpg

MetopB-20230204-DAY-1017-satellite_zenith_angle-eurol-multi.jpgMetopB-20230204-DAY-1017-satellite_zenith_angle-eurol-multi-weighted.jpg

Adam Dybbroe

unread,
Feb 6, 2023, 2:49:07 AM2/6/23
to pyt...@googlegroups.com, adam.d...@smhi.se

Ernst,

Thanks for the report! I have so far only bothered about and tested single array data, so can imagine that multi-channel composites has issues. I am sorry but have very little time this week diving into it. Perhaps I can try have a look on Friday, but can't promise.

-Adam

-- 
Adam Dybbroe,
Satellite Remote Sensing Scientist
Research lead - Nowcasting and Remote Sensing
Meteorological Research unit, Research and Development
Swedish Meteorological and Hydrological Institute (SMHI)
www.pytroll.org

lobsiger...@gmail.com

unread,
Feb 6, 2023, 12:47:55 PM2/6/23
to pytroll
Adam and other core developers,

without really knowing what I do I figured maybe my problem is now a shape issue of the arrays in the weights list:

A)
Assuming I have a Metop-B scene resampled to my area 'westminster' containing "composites " 'satellite_zenith_angle', '1',
then I will have two arrays in this scene of shape (2048, 2048). Weights has also two entries with shape (2048, 2048).

B)
But if I have in the Metop-B scene resampled to 'westminster' the composites 'satellite_zenith_angle', 'natural_color',
then I will get two arrays in this scene, one of shape (2048, 2048) and one of shape (3, 2048, 2048). But my weights
still has two entries with identical shape (2048, 2048). This ends in error giving some hints to shapes ...

QUESTIONS)
The weights list for B) should now contain two arrays, one with shape (2048, 2048) and one with shape (3, 2048, 2048)?
As these arrays are calculated from the 'satellite_zenith_angle' I would have to expand_dims in the case of RGB composites?



Any hints welcome,
Ernst

lobsiger...@gmail.com

unread,
Feb 7, 2023, 3:25:41 PM2/7/23
to pytroll
Hi All,

my last post was BS. Every scene in resampled multi scene just
takes on 2-D array. The weights list has one entry per sat pass.
I can image a multi scene with six Sen3A/Sen3B passes over area
'westminster' with 21 OLCI channels each with six weight arrays.
The problem apparently is in  multiscene.py  where below code
(used so far with 'ct', 'cma',  'ctth') cannot treat RGB composites:


def _stack_weighted(datasets, weights, combine_times):
    """Stack datasets using weights."""
    weights = set_weights_to_zero_where_invalid(datasets, weights)

    indices = da.argmax(da.dstack(weights), axis=-1)
    attrs = combine_metadata(*[x.attrs for x in datasets])

    if combine_times:
        if 'start_time' in attrs and 'end_time' in attrs:
            attrs['start_time'], attrs['end_time'] = _get_combined_start_end_times(*[x.attrs for x in datasets])

    dims = datasets[0].dims

    weighted_array = xr.DataArray(da.choose(indices, datasets), dims=dims, attrs=attrs)
    return weighted_array


that's where the shape error comes from. Well over my head :-(.

Regards,
Ernst

lobsiger...@gmail.com

unread,
Feb 9, 2023, 4:42:12 PM2/9/23
to pytroll
Adam,
numpy, xarr, dask specialists,

that's what I finally understood (but wasn't able to fix):

Assuming I have 5 satellite passes over area 'eurol' and have
a list of composites (including 'satellite_zenith_angle'), then
I must prepare a list weights with 5 xarray datasets in it. As
'satellite_zenith_angle' is resampled to 'eurol' like any other
composite it's easy to setup weights for passes over 'eurol'.


def _stack_weighted(datasets, weights, combine_times):

is then called for every of my composites with a list of
5 datasets corresponding to the 5 passes. The problem is
that xarrays in weights have always 'eurol' shape (2048, 2560)
 like 'satellite_zenith_angle' itself or any single channel
"composite". But for true composites we get 5  RGB datasets
that  instead of ('y', 'x') come with ('bands', 'y', 'x') and shape
 (3, 2048, 2560) for the 3 bands entries 'R', 'G', 'B'. That said the
function above must check on entry whether it has to treat an
RGB composite and has to redim and reshape (?) all weights
xarrays from (2048, 2560) to (3, 2048, 2560) for this to work.
For single channel "composites" it works out of the box.

My two cents,
Ernst

P.S.:
Once the code works for (identical) LEOs as expected it will also work for
combining MSG4 (0° FES) with MSG2 (IODC) and GOES16 with GOES18

lobsiger...@gmail.com

unread,
Feb 12, 2023, 1:57:33 PM2/12/23
to pytroll
Dear all

After days of try and error I finally understood that the
subject of this monologue, as taken from the PR 2275 title,
is not exactly what I'm after. I should have known from the
beginning as the respective function name and doc string

def _stack_weighted(datasets, weights, combine_times):
    """Stack datasets using weights."""

does not sound like what I think my blending should be.

This function was designed for GEO and LEO data up north.
It selects (chooses) either GEO or LEO pixels with smallest
'satellite_zenith_angle'. In my application this choice is
between two or more overlapping passes of a single type of LEO
satellite. This works without grouping as all the channels are
the same. The weighting datasets list can be generated with any
function that decreases  when 'satellite_zenith_angle' increases.

weight = 1./(1. + satellite_zenith_angle)
weight = 1./(100. + satellite_zenith_angle)
weight = cos(Pi/180.*satellite_zenith_angle)
weight = 180. - satellite_zenith_angle

The final 'weighted_array' selected will always be the same.
This is certainly interesting as well as we always get the
pixels from different LEO passes where the instrument has the
best resolution. But my idea of blending LEO passes is this:

A)
I have a weighting function with values across track from
1.0 where 'satellite_zenith_angle' is zero (SSP) to 0.0 where
'satellite_zenith_angle' is maximum (left and right of scan).
My favorite (as adapted from the German "AzB") is below:


            weights=[]
            for s in new_mscn.scenes:
                w = s['satellite_zenith_angle']
                w = w / w.max(dim = ['x', 'y'])
                weights.append(1. - w*w)

B)
Every satellite pass has an overlapping zone with the next
pass. In this zone I'd like to truly blend so that any edge
of scan should disappear in the resulting image. I'am aware
that this will introduce some "blurr" because of the time
difference of the blended passes (101 or just 50.5 minutes).
It will not work well for OLCI that has an asymmetric look.
   
C)
To avoid some intensity variation all the different weights
from different passes applied to pixels in the ovelapping
or non overlapping zones must be normalized to sum up to 1.


These are my two questions:
*****************************

1)
From my list of "AzB" weighted passes I can display every
single pass: How do I combine them into one final image?
I tried to add them up (even without step C)) to no avail.

2)
I have almost the same problem with RGB composites as before:
How do I weight (multiply) a list of RGB datasets each of shape
(3, 2048, 2560) with a list of weights each of shape (2048, 2560)?


Best regards,
Ernst


Attached below:
Three "AzB" weighted passes of Metop-B over area 'eurol'.
AVHRR single channel '2'. Shape of raw image (2048, 2560).
MetopB-20230210-DAY-0953-2-eurol-multi-weighted-pass_1.jpg
MetopB-20230210-DAY-0953-2-eurol-multi-weighted-pass_0.jpg
MetopB-20230210-DAY-0953-2-eurol-multi-weighted-pass_2.jpg

lobsiger...@gmail.com

unread,
Feb 13, 2023, 9:03:43 AM2/13/23
to pytroll
tldw.

MetopB-20230210-DAY-0953-natural_color-eurol-multi-weighted.jpg

Ernst

Raspaud Martin

unread,
Feb 13, 2023, 11:06:41 AM2/13/23
to pytroll
It look great Ernst, well done!

Best regards,
Martin
________________________________________
From: pyt...@googlegroups.com <pyt...@googlegroups.com> on behalf of lobsiger...@gmail.com <lobsiger...@gmail.com>
Sent: 13 February 2023 15:03:43
To: pytroll
Subject: Re: [pytroll] MultiScene blending with weights

tldw.

[MetopB-20230210-DAY-0953-natural_color-eurol-multi-weighted.jpg]
[MetopB-20230204-DAY-1017-1-eurol-multi.jpg][MetopB-20230204-DAY-1017-1-eurol-multi-weighted.jpg]

[MetopB-20230204-DAY-1017-2-eurol-multi.jpg][MetopB-20230204-DAY-1017-2-eurol-multi-weighted.jpg]

[MetopB-20230204-DAY-1017-ir108_3d-eurol-multi.jpg][MetopB-20230204-DAY-1017-ir108_3d-eurol-multi-weighted.jpg]

[MetopB-20230204-DAY-1017-satellite_zenith_angle-eurol-multi.jpg][MetopB-20230204-DAY-1017-satellite_zenith_angle-eurol-multi-weighted.jpg]
>> blended_scene.save_dataset('CTY_group', filename='./blended_stack_weighted_geo_polar.nc<http://blended_stack_weighted_geo_polar.nc> <http://blended_stack_weighted_geo_polar.nc>')
>> www.pytroll.org<http://www.pytroll.org> <http://www.pytroll.org>
>> nwcsaf.smhi.se<http://nwcsaf.smhi.se> <http://nwcsaf.smhi.se>
>> www.smhi.se<http://www.smhi.se> <http://www.smhi.se>
>>
>> --
>> You received this message because you are subscribed to the Google Groups "pytroll" group.
>> To unsubscribe from this group and stop receiving emails from it, send an email to pytroll+u...@googlegroups.com.
>> To view this discussion on the web, visit https://groups.google.com/d/msgid/pytroll/01ec262d-ae84-46ee-8416-edb8d2d9f669n%40googlegroups.com <https://groups.google.com/d/msgid/pytroll/01ec262d-ae84-46ee-8416-edb8d2d9f669n%40googlegroups.com?utm_medium=email&utm_source=footer>.
>
> --
> Adam Dybbroe,
> Satellite Remote Sensing Scientist
> Research lead - Nowcasting and Remote Sensing
> Meteorological Research unit, Research and Development
> Swedish Meteorological and Hydrological Institute (SMHI)
> www.pytroll.org<http://www.pytroll.org> <http://www.pytroll.org>
> nwcsaf.smhi.se<http://nwcsaf.smhi.se> <http://nwcsaf.smhi.se>
> www.smhi.se<http://www.smhi.se> <http://www.smhi.se>
>
> --
> You received this message because you are subscribed to the Google Groups "pytroll" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to pytroll+u...@googlegroups.com <mailto:pytroll+u...@googlegroups.com>.
> To view this discussion on the web, visit https://groups.google.com/d/msgid/pytroll/e09abc61-2072-44ee-821d-2f85e4c44132n%40googlegroups.com <https://groups.google.com/d/msgid/pytroll/e09abc61-2072-44ee-821d-2f85e4c44132n%40googlegroups.com?utm_medium=email&utm_source=footer>.
--
You received this message because you are subscribed to the Google Groups "pytroll" group.
To unsubscribe from this group and stop receiving emails from it, send an email to pytroll+u...@googlegroups.com.
To view this discussion on the web, visit https://groups.google.com/d/msgid/pytroll/9d613f38-2038-4592-a8a1-d7d5adaa8a88n%40googlegroups.com<https://groups.google.com/d/msgid/pytroll/9d613f38-2038-4592-a8a1-d7d5adaa8a88n%40googlegroups.com?utm_medium=email&utm_source=footer>.

--
Adam Dybbroe,
Satellite Remote Sensing Scientist
Research lead - Nowcasting and Remote Sensing
Meteorological Research unit, Research and Development
Swedish Meteorological and Hydrological Institute (SMHI)
www.pytroll.org<http://www.pytroll.org>


nwcsaf.smhi.se<http://nwcsaf.smhi.se>
www.smhi.se<http://www.smhi.se>

--
You received this message because you are subscribed to the Google Groups "pytroll" group.
To unsubscribe from this group and stop receiving emails from it, send an email to pytroll+u...@googlegroups.com<mailto:pytroll+u...@googlegroups.com>.
To view this discussion on the web, visit https://groups.google.com/d/msgid/pytroll/c3c005ce-9f90-4bce-8308-917c720ee676n%40googlegroups.com<https://groups.google.com/d/msgid/pytroll/c3c005ce-9f90-4bce-8308-917c720ee676n%40googlegroups.com?utm_medium=email&utm_source=footer>.

lobsiger...@gmail.com

unread,
Feb 13, 2023, 12:25:31 PM2/13/23
to pytroll
Hi Martin,

the last image was yet without normalization. Africa looks somewhat dim
and we have dim edges left and right. Attached two images of my final code
version as *.webp. Below is my code. Sorry, I'm completely new to this array
stuff. There is probably a lot that can be improved. I have no clue how
how to to adapt Adam's satpy code for composites though. It's maybe
the "indices" array that must get 'R', 'G', 'B' on an additional axis?

Best regards,
Ernst



def _stack_weighted(datasets, weights, combine_times):
    """Truly blend datasets using "AzB" type weights."""


    attrs = combine_metadata(*[x.attrs for x in datasets])

    if combine_times:
        if 'start_time' in attrs and 'end_time' in attrs:
            attrs['start_time'], attrs['end_time'] = _get_combined_start_end_times(*[x.attrs for x in datasets])

    dims = datasets[0].dims

    # Normalization where total = 0 ??
    total = weights[0].fillna(0).copy()
    for n in range(1, len(weights)):
        total += weights[n].fillna(0)

    datasets0=[]
    for n in range(0, len(datasets)):
        weights[n] /= total
        datasets0.append(datasets[n].fillna(0))
        # RGB composite
        if dims[0] == 'bands':
            for b in [0, 1, 2]:
               datasets0[n][b] *= weights[n].fillna(0)
        # Single channel
        else:
            datasets0[n] *= weights[n].fillna(0)

    foobar = datasets0[0].copy()
    for n in range(1, len(datasets)):
        foobar += datasets0[n]

    weighted_array = xr.DataArray(foobar, dims=dims, attrs=attrs)

    return weighted_array
MetopB-20230210-DAY-0953-natural_color-eurol10-multi-weighted.webp
MetopB-20230210-DAY-0953-2-eurol10-multi-weighted.webp

lobsiger...@gmail.com

unread,
Feb 14, 2023, 8:07:07 AM2/14/23
to pytroll

Hi,

O.K. I composite enhanced Adam's maximum elevation selecting satpy code too.
You can clearly see a cut in the Sahara between pass_0 and pass_1 and above
Maspalomas between pass_1 and pass_2. This cut is also seen between Iceland
and Norway. Processing times are extremely high ("to say the least") though.

Cheers,
Ernst
MetopB-20230210-DAY-0953-natural_color-eurol10-multi-selected.webp

lobsiger...@gmail.com

unread,
Feb 17, 2023, 11:55:16 AM2/17/23
to pytroll
To finally end this monologue I made this:


Best regards,
Ernst

lobsiger...@gmail.com

unread,
Mar 5, 2023, 12:59:10 PM3/5/23
to pytroll
Dear All,

although I kind of managed to enhance (real) multiscene blending for RGB-composites with PR2394 I still cannot reproduce what Adam
did for NWCSAF 'cloudtype'. This is due to my missing understanding of OOP (classes). The satpy multiscene documentation does not
explain how to get the satellite_zenith_angle ('n18_satz') that is, as far as I can see, not part of the polar NWCSAF cloud data.

Adam was so kind to explain how he finds 'n18_satz' here:


While I basically understand what the code does I still do not get how I can apply it to my problem. As I do only use
data of the same kind of LEO sats (Metop-B/C) I think I can simplify and speed up the methods used but still don't see
the whole loader class and how and where it comes into play. Here are the methods I think will do for my Metop-B/C
multi passes (maybe I could even combine that to one method only):

def prepare_weights_on(self):
    """Derive very simple weights and attach them to Satpy scene object."""
    self.scene['pixw'] = self._get_weights()
    self.scene['pixw'].attrs['area'] = self.scene['cma'].attrs['area']

def _get_weights(self):
    """Calculate weights from index only."""
    i0, j0 = self.scene['cma'].shape

    lon, lat = self.scene['cma'].attrs['area'].get_lonlats()
    # Treat one scan line
    warr = np.zeros((i0))
    for i in range(i0):
        # Special German 'AzB' weights
        w = (2 * i + 1 - i0) / (i0 - 1)
        warr[i] = 1 - w * w
    # Broadcast along first axis
    warr = np.ones((j0, i0)) * warr

    pixw = xr.DataArray(data=warr, dims=["y", "x"],

                        coords=dict(
                            lon=(["y", "x"], lon),
                            lat=(["y", "x"], lat)),
                        attrs=dict(
                            description="Simple weights from index.",
                            units="[1]",
    ),)

    return pixw


My scripts pseudo code for LEO NWCSAF data works as follows (this is supposed to run already in current satpy v0.40):

# Import what is needed, maybe numpy as np and xarray as xr


from satpy import Scene, MultiScene
from functools import partial
from satpy.multiscene import stack

# Make sure 'cloud_mask' (data set ID='cma') is in composites

composites = ['cloudtype', 'cloud_mask', ...]

# Find the files for consecutive LEO passes 0, 1, 2,..
# over some given 'area' and produce a list of scenes:

scenes=[]
for pass in passes:
    scenes.append(Scene(filenames=passfiles, reader=reader))

# I usually get 2 to 4 passes depending on the 'area' at hand


mscn = MultiScene(scenes)
mscn.load(composites, generate=generate)

# Now here is probably my last chance to add the weights
# to the different scenes in mscn. BUT HOW DO I DO THAT?

for s in mscn.scenes:
    ????????????????????????????????????

# Resample mscn with the different weights attached

new_mscn = mscn.resample(area, **resampler_kwargs)

# Get a list of my weights now resampled to 'area'


weights=[]
for s in new_mscn.scenes:
    weights = s['pixw']


stack_with_weights = partial(stack, weights=weights)
new_scn = new_mscn.blend(blend_function=stack_with_weights)


In a nutshell:
How can I complete and apply what Adam calls 'a small Scene loader class' with my simplified methods above to attach the weights as datasets with ID 'pixw'?


Best regards,
Ernst

David Hoese

unread,
Mar 6, 2023, 10:22:44 AM3/6/23
to pyt...@googlegroups.com
Ernst,

I'm sorry we're having such trouble with this. It is a complex problem and your conversation with Adam moved so fast that the rest of us had trouble keeping up.

Anyway, after reading Adam's answer that you linked to, I think you should ignore the idea of this being a class. I think that was easiest for what Adam had in mind and it made sense to him. If putting this functionality in a class doesn't make sense to you, then don't do it. I think I'd advise adding/using the functionality you've coded up into the body of the "for pass in passes:" loop. There you could create your Scene, load the datasets/composites you need (Scene.load), then create and add your weights "pixw" to each Scene. Once you have your list of Scenes with all of the weights (`scenes` in your current code) you can go back to creating the MultiScene object and straight into resampling the MultiScene (you will have already loaded the DataArrays) and so on with the rest of your script.

At least I think this should work. Let me know how it goes.

Dave
> https://github.com/pytroll/satpy/pull/2394 <https://github.com/pytroll/satpy/pull/2394>
> https://docs.xarray.dev/en/stable/generated/xarray.DataArray.expand_dims.html <https://docs.xarray.dev/en/stable/generated/xarray.DataArray.expand_dims.html>
> > https://raw.githubusercontent.com/pytroll/satpy/ee5a1abc659242b06f01c112b8ec7996f2be287a/satpy/multiscene.py <https://raw.githubusercontent.com/pytroll/satpy/ee5a1abc659242b06f01c112b8ec7996f2be287a/satpy/multiscene.py>
> >> blended_scene.save_dataset('CTY_group', filename='./blended_stack_weighted_geo_polar.nc <http://blended_stack_weighted_geo_polar.nc><http://blended_stack_weighted_geo_polar.nc <http://blended_stack_weighted_geo_polar.nc>> <http://blended_stack_weighted_geo_polar.nc <http://blended_stack_weighted_geo_polar.nc>>')
> >> https://satpy.readthedocs.io/en/latest/multiscene.html#stacking-scenes-using-weights <https://satpy.readthedocs.io/en/latest/multiscene.html#stacking-scenes-using-weights> <https://satpy.readthedocs.io/en/latest/multiscene.html#stacking-scenes-using-weights <https://satpy.readthedocs.io/en/latest/multiscene.html#stacking-scenes-using-weights>>
> >>
> >> :-)
> >>
> >> On 31/01/2023 22:25, lobsiger...@gmail.com wrote:
> >>> Dear developers,
> >>>
> >>> https://github.com/pytroll/satpy/pull/2275 <https://github.com/pytroll/satpy/pull/2275> <https://github.com/pytroll/satpy/pull/2275 <https://github.com/pytroll/satpy/pull/2275>>
> >>> To view this discussion on the web, visit https://groups.google.com/d/msgid/pytroll/cf566d4c-2306-4d22-8968-928d6ab37e6an%40googlegroups.com <https://groups.google.com/d/msgid/pytroll/cf566d4c-2306-4d22-8968-928d6ab37e6an%40googlegroups.com> <https://groups.google.com/d/msgid/pytroll/cf566d4c-2306-4d22-8968-928d6ab37e6an%40googlegroups.com?utm_medium=email&utm_source=footer <https://groups.google.com/d/msgid/pytroll/cf566d4c-2306-4d22-8968-928d6ab37e6an%40googlegroups.com?utm_medium=email&utm_source=footer>>.
> >>
> >> --
> >> Adam Dybbroe,
> >> Satellite Remote Sensing Scientist
> >> Research lead - Nowcasting and Remote Sensing
> >> Meteorological Research unit, Research and Development
> >> Swedish Meteorological and Hydrological Institute (SMHI)
> >> www.pytroll.org <http://www.pytroll.org><http://www.pytroll.org <http://www.pytroll.org>> <http://www.pytroll.org <http://www.pytroll.org>>
> >> nwcsaf.smhi.se <http://nwcsaf.smhi.se><http://nwcsaf.smhi.se <http://nwcsaf.smhi.se>> <http://nwcsaf.smhi.se <http://nwcsaf.smhi.se>>
> >> www.smhi.se <http://www.smhi.se><http://www.smhi.se <http://www.smhi.se>> <http://www.smhi.se <http://www.smhi.se>>
> >>
> >> --
> >> You received this message because you are subscribed to the Google Groups "pytroll" group.
> >> To unsubscribe from this group and stop receiving emails from it, send an email to pytroll+u...@googlegroups.com.
> >> To view this discussion on the web, visit https://groups.google.com/d/msgid/pytroll/01ec262d-ae84-46ee-8416-edb8d2d9f669n%40googlegroups.com <https://groups.google.com/d/msgid/pytroll/01ec262d-ae84-46ee-8416-edb8d2d9f669n%40googlegroups.com> <https://groups.google.com/d/msgid/pytroll/01ec262d-ae84-46ee-8416-edb8d2d9f669n%40googlegroups.com?utm_medium=email&utm_source=footer <https://groups.google.com/d/msgid/pytroll/01ec262d-ae84-46ee-8416-edb8d2d9f669n%40googlegroups.com?utm_medium=email&utm_source=footer>>.
> >
> > --
> > Adam Dybbroe,
> > Satellite Remote Sensing Scientist
> > Research lead - Nowcasting and Remote Sensing
> > Meteorological Research unit, Research and Development
> > Swedish Meteorological and Hydrological Institute (SMHI)
> > www.pytroll.org <http://www.pytroll.org><http://www.pytroll.org <http://www.pytroll.org>> <http://www.pytroll.org <http://www.pytroll.org>>
> > nwcsaf.smhi.se <http://nwcsaf.smhi.se><http://nwcsaf.smhi.se <http://nwcsaf.smhi.se>> <http://nwcsaf.smhi.se <http://nwcsaf.smhi.se>>
> > www.smhi.se <http://www.smhi.se><http://www.smhi.se <http://www.smhi.se>> <http://www.smhi.se <http://www.smhi.se>>
> >
> > --
> > You received this message because you are subscribed to the Google Groups "pytroll" group.
> > To unsubscribe from this group and stop receiving emails from it, send an email to pytroll+u...@googlegroups.com <mailto:pytroll+u...@googlegroups.com>.
> > To view this discussion on the web, visit https://groups.google.com/d/msgid/pytroll/e09abc61-2072-44ee-821d-2f85e4c44132n%40googlegroups.com <https://groups.google.com/d/msgid/pytroll/e09abc61-2072-44ee-821d-2f85e4c44132n%40googlegroups.com> <https://groups.google.com/d/msgid/pytroll/e09abc61-2072-44ee-821d-2f85e4c44132n%40googlegroups.com?utm_medium=email&utm_source=footer <https://groups.google.com/d/msgid/pytroll/e09abc61-2072-44ee-821d-2f85e4c44132n%40googlegroups.com?utm_medium=email&utm_source=footer>>.
> --
> You received this message because you are subscribed to the Google Groups "pytroll" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to pytroll+u...@googlegroups.com.
> To view this discussion on the web, visit https://groups.google.com/d/msgid/pytroll/9d613f38-2038-4592-a8a1-d7d5adaa8a88n%40googlegroups.com <https://groups.google.com/d/msgid/pytroll/9d613f38-2038-4592-a8a1-d7d5adaa8a88n%40googlegroups.com><https://groups.google.com/d/msgid/pytroll/9d613f38-2038-4592-a8a1-d7d5adaa8a88n%40googlegroups.com?utm_medium=email&utm_source=footer <https://groups.google.com/d/msgid/pytroll/9d613f38-2038-4592-a8a1-d7d5adaa8a88n%40googlegroups.com?utm_medium=email&utm_source=footer>>.
>
> --
> Adam Dybbroe,
> Satellite Remote Sensing Scientist
> Research lead - Nowcasting and Remote Sensing
> Meteorological Research unit, Research and Development
> Swedish Meteorological and Hydrological Institute (SMHI)
> www.pytroll.org <http://www.pytroll.org><http://www.pytroll.org <http://www.pytroll.org>>
>
>
> nwcsaf.smhi.se <http://nwcsaf.smhi.se><http://nwcsaf.smhi.se <http://nwcsaf.smhi.se>>
> www.smhi.se <http://www.smhi.se><http://www.smhi.se <http://www.smhi.se>>
>
> --
> You received this message because you are subscribed to the Google Groups "pytroll" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to pytroll+u...@googlegroups.com<mailto:pytroll+u...@googlegroups.com>.
> To view this discussion on the web, visit https://groups.google.com/d/msgid/pytroll/c3c005ce-9f90-4bce-8308-917c720ee676n%40googlegroups.com <https://groups.google.com/d/msgid/pytroll/c3c005ce-9f90-4bce-8308-917c720ee676n%40googlegroups.com><https://groups.google.com/d/msgid/pytroll/c3c005ce-9f90-4bce-8308-917c720ee676n%40googlegroups.com?utm_medium=email&utm_source=footer <https://groups.google.com/d/msgid/pytroll/c3c005ce-9f90-4bce-8308-917c720ee676n%40googlegroups.com?utm_medium=email&utm_source=footer>>.
>
> --
> You received this message because you are subscribed to the Google Groups "pytroll" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to pytroll+u...@googlegroups.com <mailto:pytroll+u...@googlegroups.com>.
> To view this discussion on the web, visit https://groups.google.com/d/msgid/pytroll/c24e755f-e8ef-4c1c-bc31-9d44ef59884cn%40googlegroups.com <https://groups.google.com/d/msgid/pytroll/c24e755f-e8ef-4c1c-bc31-9d44ef59884cn%40googlegroups.com?utm_medium=email&utm_source=footer>.

lobsiger...@gmail.com

unread,
Mar 6, 2023, 11:17:24 AM3/6/23
to pytroll
Dave,

thanks for that. As a matter of fact this moved too fast for me too and the problem Adam wanted to solve initially wasn't even a
priority to me. My LEO scripts use a module LEOstuff.py that makes about everything. It's not straight forward to add all this new stuff.
I also thought of adding the weights for NWCSAF cloud data passes when I make the list of scenes. I have to code a minimal version
of this stuff to experiment with. For all the AVHRR/VIIRS/MODIS/MERSI2 RGB-composites it's very easy as 'satellite_zenith_angle' is
available in the data. I just have to put it into the composites list before I load and resample the multi scene. Maybe I even forget to try
this for cloud data (it will only work for 'select_with_weights' anyway). I have already learnt a lot (but with 70+ forget a lot fast too ;-\).

Thanks again,
Ernst

lobsiger...@gmail.com

unread,
Mar 9, 2023, 10:48:18 AM3/9/23
to pytroll
Dear All,

I made three minimal scripts that use the new MultiScene blend_types for

GOES16+GOES18 (same time slot, ABI)
MetopB (3 consecutive passes, AVHRR)
MetopB (3 consecutive passes, NWCSAF)

The NWCSAF 'select_with_weights' should in theory work in plain satpy 0.40.
Unfortunately it does not. I assume this is either a bug or a "feature" of
the reader for Metop NWCSAF data (or I do not get all data on EUMETCast?).

You can download these scripts including all necessary data on wetransfer.
I hope they help to bring forward my PR#2394, which is quite strong stuff.


App. 1.3GB zip file containing 3 *.tgz.
This link expires in one week from now.


Best regards,
Ernst

I attach the scripts as *.txt and the errors from NWCSAF clouds.txt (.py)
debug.txt
MetopB.txt
GOESXX.txt
clouds.txt

lobsiger...@gmail.com

unread,
Mar 10, 2023, 3:26:29 PM3/10/23
to pytroll
Dear All,

I made some more tests with my MetopB LEO NWCSAF cloud data.

Here are the composite names that should work with palettes:
['cloudmask', 'cloudmask_extended','cloud_top_height',
 'cloud_top_pressure', 'cloud_top_temperature', 'cloudtype']

Here are the corresponding (limited) integer dataset names:
['cma', 'cma_extended', 'ctth_alti',
 'ctth_pres', 'ctth_tempe', 'ct']

I have tested with my code of PR#2394 using 3 blend_types:

From the composites only 'cloudmask' and 'cloudtype' worked
and these *only* worked with blend_type='stack_no_weights'.
Other blend_types threw errors that I assume are compositor
problems related to ISSUE#2392 and maybe fixed with PR#2403.
IIRC all above composites did work one year back from now.

From dataset names all worked (somehow) with all blend_types
*but* 'ct' that only worked exactly like 'cloudtype' above.

My results show that the code to attach simple weights works.
Interestingly the results also depend on generate=True/False.

To summarize my results I made two images (composed with IM).
First image is for generate=True, second for generate=False.

The columns left to right are the 5 datasets named below:
|'cma'|'cma_extended'|'ctth_alti'|'ctth_pres'|'ctth_tempe'|

The rows from top to bottom are the 3 blend_types below:
blend_type = 'stack_no_weights'
blend_type = 'select_with_weights'
blend_type = 'blend_with_weights'


Cheers,
Ernst
MetopB-NWCSAF-generate-True.jpg
MetopB-NWCSAF-generate-False.jpg

lobsiger...@gmail.com

unread,
Mar 11, 2023, 11:35:49 AM3/11/23
to pytroll
Hi,

taking good old satpy version 0.36 the composites basically worked.
The triangle (where all 3 passes overlap) has yet to be understood.

The columns left to right are the 6 composites named below:
|'cloudmask'|'cloudmask_extended'|'cloud_top_height'|'cloud_top_pressure'| 'cloud_top_temperature'|'cloudtype'|
There are color palettes for each composite in the *.nc data.


The rows from top to bottom are the 3 blend_types below:
blend_type = 'stack_no_weights'
blend_type = 'select_with_weights'
blend_type = 'blend_with_weights'

It seems I'll probably have to wait and see what PR#2403 brings.

Regards,
Ernst
MetopB-NWCSAF-COLOR-generate-True.jpg
MetopB-NWCSAF-COLOR-generate-False.jpg
Reply all
Reply to author
Forward
0 new messages