EUMETSAT SEVIRI Level 1.5 data - Cannot generate composites with the HRV channel

[cbincus]

Describe the bug Satpy fails to generate SEVIRI composites with the HRV channel (e.g., HRV Clouds). Composites that do not require this channel are generated correctly.

To Reproduce

from satpy import Scene

filenames = ['MSG2-SEVI-MSG15-0100-NA-20081005135741.556000000Z-NA.nat']
global_scene = Scene(reader='seviri_l1b_native', filenames=filenames)

composite = 'hrv_clouds'
global_scene.load([composite], upper_right_corner='NE')

euro_scn = global_scene.resample('euro4', resampler='nearest')
euro_scn.show(composite)

Link to the data file: http://usb.iphoto.md/zpJuc

Expected behavior The composite image is generated correctly.

Actual results

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
File ~/miniconda3/envs/sat/lib/python3.10/site-packages/dask/highlevelgraph.py:550, in HighLevelGraph.get_all_external_keys(self)
    549 try:
--> 550     return self._all_external_keys
    551 except AttributeError:

AttributeError: 'HighLevelGraph' object has no attribute '_all_external_keys'

During handling of the above exception, another exception occurred:

ValueError                                Traceback (most recent call last)
File ~/miniconda3/envs/sat/lib/python3.10/site-packages/IPython/core/formatters.py:344, in BaseFormatter.__call__(self, obj)
    342     method = get_real_method(obj, self.print_method)
    343     if method is not None:
--> 344         return method()
    345     return None
    346 else:

File ~/miniconda3/envs/sat/lib/python3.10/site-packages/trollimage/xrimage.py:1469, in XRImage._repr_png_(self)
   1467 import io
   1468 b = io.BytesIO()
-> 1469 self.pil_image().save(b, format='png')
   1470 return b.getvalue()

File ~/miniconda3/envs/sat/lib/python3.10/site-packages/trollimage/xrimage.py:923, in XRImage.pil_image(self, fill_value, compute)
    921 img = dask.delayed(PILImage.fromarray)(np.squeeze(res.data), mode)
    922 if compute:
--> 923     img = img.compute()
    924 return img

File ~/miniconda3/envs/sat/lib/python3.10/site-packages/dask/base.py:314, in DaskMethodsMixin.compute(self, **kwargs)
    290 def compute(self, **kwargs):
    291     """Compute this dask collection
    292 
    293     This turns a lazy Dask collection into its in-memory equivalent.
   (...)
    312     dask.compute
    313     """
--> 314     (result,) = compute(self, traverse=False, **kwargs)
    315     return result

File ~/miniconda3/envs/sat/lib/python3.10/site-packages/dask/base.py:593, in compute(traverse, optimize_graph, scheduler, get, *args, **kwargs)
    585     return args
    587 schedule = get_scheduler(
    588     scheduler=scheduler,
    589     collections=collections,
    590     get=get,
    591 )
--> 593 dsk = collections_to_dsk(collections, optimize_graph, **kwargs)
    594 keys, postcomputes = [], []
    595 for x in collections:

File ~/miniconda3/envs/sat/lib/python3.10/site-packages/dask/base.py:366, in collections_to_dsk(collections, optimize_graph, optimizations, **kwargs)
    364 for opt, val in groups.items():
    365     dsk, keys = _extract_graph_and_keys(val)
--> 366     dsk = opt(dsk, keys, **kwargs)
    368     for opt_inner in optimizations:
    369         dsk = opt_inner(dsk, keys, **kwargs)

File ~/miniconda3/envs/sat/lib/python3.10/site-packages/dask/delayed.py:514, in optimize(dsk, keys, **kwargs)
    512 if not isinstance(dsk, HighLevelGraph):
    513     dsk = HighLevelGraph.from_collections(id(dsk), dsk, dependencies=())
--> 514 dsk = dsk.cull(set(flatten(keys)))
    515 return dsk

File ~/miniconda3/envs/sat/lib/python3.10/site-packages/dask/highlevelgraph.py:707, in HighLevelGraph.cull(self, keys)
    703 from dask.layers import Blockwise
    705 keys_set = set(flatten(keys))
--> 707 all_ext_keys = self.get_all_external_keys()
    708 ret_layers: dict = {}
    709 ret_key_deps: dict = {}

File ~/miniconda3/envs/sat/lib/python3.10/site-packages/dask/highlevelgraph.py:557, in HighLevelGraph.get_all_external_keys(self)
    552 keys: set = set()
    553 for layer in self.layers.values():
    554     # Note: don't use `keys |= ...`, because the RHS is a
    555     # collections.abc.Set rather than a real set, and this will
    556     # cause a whole new set to be constructed.
--> 557     keys.update(layer.get_output_keys())
    558 self._all_external_keys = keys
    559 return keys

File ~/miniconda3/envs/sat/lib/python3.10/site-packages/dask/blockwise.py:486, in Blockwise.get_output_keys(self)
    480     return {(self.output, *p) for p in self.output_blocks}
    482 # Return all possible output keys (no culling)
    483 return {
    484     (self.output, *p)
    485     for p in itertools.product(
--> 486         *[range(self.dims[i]) for i in self.output_indices]
    487     )
    488 }

File ~/miniconda3/envs/sat/lib/python3.10/site-packages/dask/blockwise.py:486, in (.0)
    480     return {(self.output, *p) for p in self.output_blocks}
    482 # Return all possible output keys (no culling)
    483 return {
    484     (self.output, *p)
    485     for p in itertools.product(
--> 486         *[range(self.dims[i]) for i in self.output_indices]
    487     )
    488 }

File ~/miniconda3/envs/sat/lib/python3.10/site-packages/dask/blockwise.py:446, in Blockwise.dims(self)
    442 """Returns a dictionary mapping between each index specified in
    443 `self.indices` and the number of output blocks for that indice.
    444 """
    445 if not hasattr(self, "_dims"):
--> 446     self._dims = _make_dims(self.indices, self.numblocks, self.new_axes)
    447 return self._dims

File ~/miniconda3/envs/sat/lib/python3.10/site-packages/dask/blockwise.py:1479, in _make_dims(indices, numblocks, new_axes)
   1475 def _make_dims(indices, numblocks, new_axes):
   1476     """Returns a dictionary mapping between each index specified in
   1477     `indices` and the number of output blocks for that indice.
   1478     """
-> 1479     dims = broadcast_dimensions(indices, numblocks)
   1480     for k, v in new_axes.items():
   1481         dims[k] = len(v) if isinstance(v, tuple) else 1

File ~/miniconda3/envs/sat/lib/python3.10/site-packages/dask/blockwise.py:1470, in broadcast_dimensions(argpairs, numblocks, sentinels, consolidate)
   1467     return toolz.valmap(consolidate, g2)
   1469 if g2 and not set(map(len, g2.values())) == {1}:
-> 1470     raise ValueError("Shapes do not align %s" % g)
   1472 return toolz.valmap(toolz.first, g2)

ValueError: Shapes do not align {'.1': {2, 12}, '.0': {1}}

Screenshots Not applicable.

Environment Info:

• OS: Ubuntu 20.04.6 LTS
• Python version: 3.10.11
• Satpy Version: 0.42.2
• PyResample Version: 1.26.1

Additional context None.

[mraspaud]

@cbincus thanks for reporting this! I can reproduce the error with the file you provide. I also tried with another data format (HRIT) and that seems to work fine. So the error could be in the file or in the nat format reader in satpy. @sjoro @ameraner do you have any idea of what could be happening?

[ameraner]

Hi @cbincus , the issue is solved by using reader_kwargs={'fill_disk': True} in your Scene initialisation, like this

# ...
global_scene = Scene(reader='seviri_l1b_native', filenames=filenames, reader_kwargs={'fill_disk': True})
#...

background: the HRV comes in two separate windows, and this sometimes causes problems with its usage. fill_disk makes sure that the HRV is padded to full-disc during the reading, making it equivalent to the other channels and avoiding this issue.

Output image looks like this in my own test granule - you may want to add radius_of_influence=30e3 to your resample call to fill the gaps.

[gerritholl]

MCVE to reproduce:

import os
from satpy import Scene
from glob import glob
from pathlib import Path
from satpy.utils import debug_on; debug_on()
fn = "/media/nas/x21308/scratch/data/MSG4-SEVI-MSG15-0100-NA-20220216002743.011000000Z-NA.nat"
sc = Scene(filenames=[fn], reader=["seviri_l1b_native"])
sc.load(["ir_sandwich"])
ls = sc.resample("eurol")
ls["ir_sandwich"].compute()