Error loading composites with seviri_l1b_native

[Marty Sullivan]

Hello,

I'm downloading some SEVIRI l1b data in the .nat format. I can load some composites like "colorized_ir_clouds" but I get a ValueError when I try to load hi res composites like "realistic_colors" or "overview". I'd appreciate any help debugging this problem so that I can generate these composites, thanks!

Here is the code I'm running:

from satpy import Scene

scene = Scene(filenames=['MSG4-SEVI-MSG15-0100-NA-20210311041243.393000000Z-NA.nat'], reader='seviri_l1b_native')

scene.load(['realistic_colors'])

And here is the traceback I get: 

Traceback (most recent call last):

  File "<stdin>", line 1, in <module>

  File "/usr/local/lib/python3.8/site-packages/satpy/scene.py", line 1182, in load

    self.generate_possible_composites(generate, unload)

  File "/usr/local/lib/python3.8/site-packages/satpy/scene.py", line 1239, in generate_possible_composites

    keepables = self._generate_composites_from_loaded_datasets()

  File "/usr/local/lib/python3.8/site-packages/satpy/scene.py", line 1251, in _generate_composites_from_loaded_datasets

    return self._generate_composites_nodes_from_loaded_datasets(nodes)

  File "/usr/local/lib/python3.8/site-packages/satpy/scene.py", line 1257, in _generate_composites_nodes_from_loaded_datasets

    self._generate_composite(node, keepables)

  File "/usr/local/lib/python3.8/site-packages/satpy/scene.py", line 1281, in _generate_composite

    prereq_datasets = self._get_prereq_datasets(

  File "/usr/local/lib/python3.8/site-packages/satpy/scene.py", line 1363, in _get_prereq_datasets

    self._generate_composite(prereq_node, keepables)

  File "/usr/local/lib/python3.8/site-packages/satpy/scene.py", line 1315, in _generate_composite

    composite = compositor(prereq_datasets,

  File "/usr/local/lib/python3.8/site-packages/satpy/modifiers/geometry.py", line 67, in __call__

    lons, lats = vis.attrs["area"].get_lonlats(chunks=vis.data.chunks)

  File "/usr/local/lib/python3.8/site-packages/pyresample/geometry.py", line 2349, in get_lonlats

    lons, lats = definition.get_lonlats(nprocs=nprocs, data_slice=(local_row_slice, col_slice),

  File "/usr/local/lib/python3.8/site-packages/pyresample/geometry.py", line 1949, in get_lonlats

    target_x, target_y = self.get_proj_coords(data_slice=data_slice, chunks=chunks, dtype=dtype)

  File "/usr/local/lib/python3.8/site-packages/pyresample/geometry.py", line 1852, in get_proj_coords

    target_x, target_y = self._get_proj_vectors(dtype=dtype, check_rotation=False, chunks=chunks)

  File "/usr/local/lib/python3.8/site-packages/pyresample/geometry.py", line 1795, in _get_proj_vectors

    target_y = arange(self.height, **y_kwargs) * -self.pixel_size_y + self.pixel_upper_left[1]

  File "/usr/local/lib/python3.8/site-packages/dask/array/creation.py", line 385, in arange

    chunks = normalize_chunks(chunks, (num,), dtype=dtype)

  File "/usr/local/lib/python3.8/site-packages/dask/array/core.py", line 2726, in normalize_chunks

    raise ValueError(

ValueError: Chunks do not add up to shape. Got chunks=((3711, 3711, 3714),), shape=(8320,)

[Marty Sullivan]

It seems like this issue might be related to satpy #1261, from when the HRV read was still not dask compatible. It seems like progress was made here but I'm unable to follow the discussion in Slack mentioned in the linked issue on GitHub since I can't join the workspace.

I'd be happy to work on this a bit if there's some code I can contribute to get this reader working, if it's still not quite right.

Best,

Marty

[David Hoese]

Hi Marty,

I believe that issue is fixed, but requires a new version of pyresample.
The related fix *should* have been released already. Do you know what
version of pyresample you are using?

Dave

On 3/11/21 9:51 AM, Marty Sullivan wrote:
> It seems like this issue might be related to satpy #1261

> <https://github.com/pytroll/satpy/issues/1261>, from when the HRV read

> --
> 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/f118fe9d-1c4b-4cef-8fd4-c8c4b6521888n%40googlegroups.com
> <https://groups.google.com/d/msgid/pytroll/f118fe9d-1c4b-4cef-8fd4-c8c4b6521888n%40googlegroups.com?utm_medium=email&utm_source=footer>.

[Marty Sullivan]

I'm using pyresample 1.17.0

[David Hoese]

Ah I think the issue was fixed in this pull request by Andrea:

https://github.com/pytroll/pyresample/pull/313

Looks like, based on the dates on PyPI and the PR, that this was merged
5 days after 1.17.0 was released. Any chance you could try installing
pyresample from the github master branch?

Dave

> <https://groups.google.com/d/msgid/pytroll/f118fe9d-1c4b-4cef-8fd4-c8c4b6521888n%40googlegroups.com?utm_medium=email&utm_source=footer
> <https://groups.google.com/d/msgid/pytroll/f118fe9d-1c4b-4cef-8fd4-c8c4b6521888n%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/44e9e35e-7527-4a75-bf10-4c3b2747ebc9n%40googlegroups.com
> <https://groups.google.com/d/msgid/pytroll/44e9e35e-7527-4a75-bf10-4c3b2747ebc9n%40googlegroups.com?utm_medium=email&utm_source=footer>.

[Marty Sullivan]

I'll test it out but I am using this in production so would rather avoid installing from master. I suppose if this fix works, I can just patch the two updated files in that PR in my install though and that would suffice.

Thanks for the info though, I'll update once I have a chance to try this out.

Marty

[David Hoese]

For what its worth, we hope to have a pyresample release in the next
couple business days (followed by a Satpy release).

Dave

> <https://groups.google.com/d/msgid/pytroll/44e9e35e-7527-4a75-bf10-4c3b2747ebc9n%40googlegroups.com?utm_medium=email&utm_source=footer
> <https://groups.google.com/d/msgid/pytroll/44e9e35e-7527-4a75-bf10-4c3b2747ebc9n%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/7e1de8d7-60a3-45d6-bd5c-0431f3b67cb7n%40googlegroups.com
> <https://groups.google.com/d/msgid/pytroll/7e1de8d7-60a3-45d6-bd5c-0431f3b67cb7n%40googlegroups.com?utm_medium=email&utm_source=footer>.

[Marty Sullivan]

Patching the files does seem to work to get the composite to load. And great to hear a new release is coming.

However, once I Ioad the composite, when I use scene.max_area() it tells me ValueError: Can't compare areas of different types. Am I correct in assuming I'll need to use a resampler other than native for these composites?

Thanks,
Marty

[David Hoese]

I believe so. I'm not a SEVIRI expert but I think because the way HRV is
structured, you can't use `max_area` because it assumes they are all the
same extents (HRV is different) and projection. There are some SEVIRI
related areas that are being released in the next Satpy that you could
resample to:

https://github.com/pytroll/satpy/blob/master/satpy/etc/areas.yaml#L15

They have different names in the current release, but I don't remember
what they are. Sorry.

Dave

> <https://groups.google.com/d/msgid/pytroll/7e1de8d7-60a3-45d6-bd5c-0431f3b67cb7n%40googlegroups.com?utm_medium=email&utm_source=footer
> <https://groups.google.com/d/msgid/pytroll/7e1de8d7-60a3-45d6-bd5c-0431f3b67cb7n%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/0c86477d-7aa6-48b4-87dc-0e5ad6fa427dn%40googlegroups.com
> <https://groups.google.com/d/msgid/pytroll/0c86477d-7aa6-48b4-87dc-0e5ad6fa427dn%40googlegroups.com?utm_medium=email&utm_source=footer>.

[Marty Sullivan]

Ok, that makes sense. 

I'm having a bit of trouble finding the docs on using those predefined areas in satpy. I'm used to using create_area_def but it seems like these are set up to avoid having to create the AreaDefinition from scratch. Can you point me to how I can load these predefined areas?

Thanks,

Marty

[Andrea Meraner]

Hi Marty, 

glad to see that my fix worked and was helpful for you. 

You can still use max_area and the native resampler with HRV data by padding the HRV dataset to the full disk. This is done by passing the fill_disk reader kwarg:
scene = Scene(filenames=['MSG4-SEVI-MSG15-0100-NA-20210311041243.393000000Z-NA.nat'], reader='seviri_l1b_native' , reader_kwargs={'fill_disk': True})

This lets you past that Can't compare areas of different types error.

You can use the predefined areas directly in the resampling call by their name: 

scn_r= scn.resample('seviri_0deg', resampler='nearest')

seviri_0deg is the old name for the current msg_seviri_fes_3km area. In the old releases there was no 1km SEVIRI area (current msg_seviri_fes_1km) .

Hope this is helpful,

Cheers, 
Andrea

[Andrea Meraner]

Oh and if you're interested in loading the AreaDefinition instances for those areas, there's also

from satpy.resample import get_area_def
get_area_def('seviri_0deg')

Andrea

[Marty Sullivan]

Thanks for the guidance, Andrea.

I tried adding in the fill_disk kwarg and got further, but now I get a new ValueError; any insight on the following? Thanks in advance.

Code:

from satpy import Scene

scene = Scene(filenames=['MSG4-SEVI-MSG15-0100-NA-20210311041243.393000000Z-NA.nat'], reader='seviri_l1b_native', reader_kwargs={'fill_disk': True})

resampled = scene.resample(scene.max_area(), resampler='native')

Traceback:

Traceback (most recent call last):

  File "<stdin>", line 1, in <module>

  File "/usr/local/lib/python3.8/site-packages/satpy/scene.py", line 837, in resample

    self._resampled_scene(new_scn, destination, resampler=resampler,

  File "/usr/local/lib/python3.8/site-packages/satpy/scene.py", line 793, in _resampled_scene

    res = resample_dataset(dataset, destination_area, **kwargs)

  File "/usr/local/lib/python3.8/site-packages/satpy/resample.py", line 1321, in resample_dataset

    new_data = resample(source_area, dataset, destination_area, fill_value=fill_value, **kwargs)

  File "/usr/local/lib/python3.8/site-packages/satpy/resample.py", line 1284, in resample

    res = resampler_instance.resample(data, **kwargs)

  File "/usr/local/lib/python3.8/site-packages/satpy/resample.py", line 946, in resample

    return super(NativeResampler, self).resample(data,

  File "/usr/local/lib/python3.8/site-packages/satpy/resample.py", line 424, in resample

    return self.compute(data, cache_id=cache_id, **kwargs)

  File "/usr/local/lib/python3.8/site-packages/satpy/resample.py", line 1036, in compute

    d_arr = self.expand_reduce(data.data, repeats)

  File "/usr/local/lib/python3.8/site-packages/satpy/resample.py", line 1002, in expand_reduce

    return cls.aggregate(d_arr, y_size, x_size)

  File "/usr/local/lib/python3.8/site-packages/satpy/resample.py", line 964, in aggregate

    raise ValueError("Aggregation requires arrays with "

ValueError: Aggregation requires arrays with shapes and chunks divisible by the factor

[David Hoese]

I think we ran into this just the other day. There is something up with
the chunk size (low-ish level dask thing) for the HRV channel where it
has different sized chunks compared to the other channels and that
causes the native resampling to fail.

Although not great for performance you could do `resampler='nearest'`
which is also the default. You could improve performance for future runs
by adding `cache_dir='/path/to/some/directory'`.

Andrea, do you remember who ran into this issue on slack the other day?
Any idea if we have an issue for it?

Dave

> <https://groups.google.com/d/msgid/pytroll/0c86477d-7aa6-48b4-87dc-0e5ad6fa427dn%40googlegroups.com?utm_medium=email&utm_source=footer
> <https://groups.google.com/d/msgid/pytroll/0c86477d-7aa6-48b4-87dc-0e5ad6fa427dn%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/acacecf7-2ec4-4eda-9fff-cab4b6cad1c3n%40googlegroups.com
> <https://groups.google.com/d/msgid/pytroll/acacecf7-2ec4-4eda-9fff-cab4b6cad1c3n%40googlegroups.com?utm_medium=email&utm_source=footer>.

[Andrea Meraner]

Hi both, 
@Dave, see Slack for the discussion, I also created an issue for it: https://github.com/pytroll/satpy/issues/1595

@Marty: just to double check, as this min_area/max_area has caused considerable confusion, plus there's currently a bug that reverts the logic for flipped SEVIRI data:

the error you're getting shows up if you try to downsample the data (to the 3km grid). If you're trying to get the composite in high-res (1km), then min_area is what will do it and is working already.

FYI, min_area and max_area have been recently renamed to coarsest_area and finest_area, and the bug affecting the SEVIRI data will be solved soon (preparing the fix).

Cheers, 
Andrea

[Marty Sullivan]

Thanks for this info! Yes I definitely had not noticed that min_area/max_area was flipped for this and it works great by switching to min_area for seviri. 

With all that, I do have this working now, so thanks for all the guidance!

I did see the documentation update for finest/coarsest area so it sounds like that will make this even easier in the next satpy release. 

Best,
Marty