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
:-)
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.
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)
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
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
To view this discussion on the web, visit https://groups.google.com/d/msgid/pytroll/01ec262d-ae84-46ee-8416-edb8d2d9f669n%40googlegroups.com.
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
To view this discussion on the web, visit https://groups.google.com/d/msgid/pytroll/9d613f38-2038-4592-a8a1-d7d5adaa8a88n%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