Issues regarding blending datasets in Satpy

[isaac wilmur]

Hello everyone,

I'm currently attempting to combine data from multiple geostationary satellites to create a global composite.

Satpy is already able to work with combinations below:

(1) Himawari 8 + GOES 16 + 17 : works

(2) Meteosat 8 + 11 : works

However Satpy is unable to save generated datasets when combining data between satellites from two sets above.

It is currently possible for me to generate those separately and then combine them with external tools but I'd like to have a more in-line solution if possible.

Below are my current scripts, debug logs, related datasets and results. Most of the codes are inspired by the MultiScene section of the Documentation.

https://drive.google.com/drive/folders/1TX6AeOH4XHISvpLbsvF-NxBnJVLVzh_w?usp=sharing

Additional information:

OS: Gentoo Linux 2.8

Satpy version: 0.35.0

Related packages are installed through pip

Thanks,

Isaac

[David Hoese]

Hi Isaac,

That's a lot of info to sift through and I'm a little short on time so I wasn't able to get through all of it or try it myself, but I looked at your "debug_all.txt" and I don't see anything wrong. The code looks like I would expect it to and things seem to be behaving as expected until that very last call to save_dataset. Could you try printing "print(list(blended.keys()))" and letting me know the output? It seems like something is silently failing in the blending phase.

Also, I notice the wavelength in your group in your separate script has 4 elements (9, 10, 11, 12). This should be invalid and I'm very surprised it doesn't cause a failure. I suppose it might be considering "12" as the units for the wavelength. We might need to add a check on that parameter to make sure it is a string unit. This probably doesn't effect your script (as it works), but it is something I noticed.

Dave

[Isaac Wilmur]

Hi Dave,

Here's the printing output, appended right after the blend() function call in my gen_all script, it seems like that it returns an empty list.

>>> resampled = mscn.resample('worldeqc3km70', reduce_data=False, radius_of_influence=10000)

>>> blended = resampled.blend()

>>> print(list(blended.keys()))
[]
>>>

It behaves differently in my separate generation script, interestingly enough in the Meteosat blend it also lists my inversion composite.

>>> resampled = mscn.resample('worldeqc3km70', reduce_data=False, radius_of_influence=10000)

>>> blended = resampled.blend()

>>> print(list(blended.keys()))
[DataID(name='ir_group0', resolution=2000)]
>>>

>>> resampled = mscn.resample('worldeqc3km70', reduce_data=False, radius_of_influence=10000)

>>> blended = resampled.blend()

>>> print(list(blended.keys()))
[DataID(name='i18_inv', resolution=3000.403165817), DataID(name='ir_group1', resolution=3000.403165817)]

>>>

As for the second issue, I'm informed through the documentation that it should only takes 3 elements for min, central, max wavelength and an extra string for an optional unit, it must have gotten lost in there during my troubleshooting process. Things should still behave as expected when that element is removed.

Hope that was helpful,

Isaac

[David Hoese]

Very interesting. So as far as I know the "group" logic is only supposed
to provide an "alias" (new name) to a copy of the original DataArray. I
suppose in the Meteosat case the DataIDs are the same so it is able to
produce both the group product and the original name product.

Could you may try, in the gen_all.py script, doing `mscn.save_datasets`
to get some geotiffs (or PNGs if you want) for the individual resampled
DataArrays? You might want to pass `datasets=["ir_group"]`. If you get a
geotiff for each input instrument and the images look as expected then
we'll need to dive further into the metadata and the `blend` method to
track down what's going on.

Dave

> --
> You received this message because you are subscribed to a topic in the
> Google Groups "pytroll" group.
> To unsubscribe from this topic, visit
> https://groups.google.com/d/topic/pytroll/XVTTdXFDewk/unsubscribe
> <https://groups.google.com/d/topic/pytroll/XVTTdXFDewk/unsubscribe>.
> To unsubscribe from this group and all its topics, 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/9820d54f-8d8c-46f5-bb69-620639a4b30en%40googlegroups.com
> <https://groups.google.com/d/msgid/pytroll/9820d54f-8d8c-46f5-bb69-620639a4b30en%40googlegroups.com?utm_medium=email&utm_source=footer>.

[Isaac Wilmur]

Hi Dave,

I've run through your suggestions and updated my drive repository above with those new articles. In it you'll find:

- a modified gen_all.py script (gen_all_noblend_withgroup.py)
- associated debug log

- generated gtiffs (./result_noblend_withgroup/)

The results are pretty much what I'd expect: Separate images for each instrument/satellite under the "ir_group" alias with roughly the same timestamps.

Here's the link again just in case:

https://drive.google.com/drive/folders/1TX6AeOH4XHISvpLbsvF-NxBnJVLVzh_w?usp=sharing

Hope that was helpful,

Isaac

[David Hoese]

Isaac,

One last thing to confirm something that I think is going on. Could you
print out `print(r_mscn.shared_dataset_ids)` (or "resampled" instead of
"r_mscn" depending on which version of the script you're using).

I think what is happening is that the grouping is naively renaming the
datasets but not removing the other identifying information like
resolution. You can see that in this function:

https://github.com/pytroll/satpy/blob/c20c3bb7042c819baffbb187f1ccaa9efbd74347/satpy/multiscene.py#L72-L94

You can also see how in the blend method of the MultiScene if there are
no datasets that are considered "common" between all of the sub-Scenes
that the for loop won't do anything:

https://github.com/pytroll/satpy/blob/c20c3bb7042c819baffbb187f1ccaa9efbd74347/satpy/multiscene.py#L319-L325

So this is a bug. As a quick workaround you could re-implement the blend
method as your own function:


from satpy.multiscene import stack

def my_blend(mscn, common_datasets, blend_function=stack):
new_scn = Scene()
for ds_id in common_datasets:
datasets = [scn[ds_id] for scn in mscn.scenes if ds_id in scn]
new_scn[ds_id] = blend_function(datasets)

return new_scn

blended = my_blend(resampled, ["ir_group"])


Good luck.

Dave

> <https://groups.google.com/d/msgid/pytroll/9820d54f-8d8c-46f5-bb69-620639a4b30en%40googlegroups.com?utm_medium=email&utm_source=footer
> <https://groups.google.com/d/msgid/pytroll/9820d54f-8d8c-46f5-bb69-620639a4b30en%40googlegroups.com?utm_medium=email&utm_source=footer>>.

>
>
> --
> You received this message because you are subscribed to a topic in the
> Google Groups "pytroll" group.
> To unsubscribe from this topic, visit
> https://groups.google.com/d/topic/pytroll/XVTTdXFDewk/unsubscribe
> <https://groups.google.com/d/topic/pytroll/XVTTdXFDewk/unsubscribe>.
> To unsubscribe from this group and all its topics, 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/627aeb9b-4dfa-449b-bec9-edf7e295a4a9n%40googlegroups.com
> <https://groups.google.com/d/msgid/pytroll/627aeb9b-4dfa-449b-bec9-edf7e295a4a9n%40googlegroups.com?utm_medium=email&utm_source=footer>.

[Isaac Wilmur]

Hello,

The output of 'print(r_mscn.shared_dataset_ids)'  is included below

As for the re-implemented function, it worked and the result is wonderful. I've also shuffled the satellites order around a little, added borders and coastlines for georeferencing. The script used is also included below.

I consider this to be solved on my end for now, but please reach out if you need any extra input.

Thanks,

Isaac

[David Hoese]

FYI to anyone reading this in the future, the bug described in this thread should now be fixed in the `main` branch of the satpy repository. It was fixed in this PR: https://github.com/pytroll/satpy/issues/2089

It should be released in the next release of Satpy (no date yet, but I would expect in the next week or two).

Dave

[Gis Geo]

Thank you David

Best regards

--
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/82d63067-e33e-4841-9810-ac14ee05af37n%40googlegroups.com.