`AreaDefinition.get_area_slices` doesn't work for cropped geostationary areas like SEVIRI RSS

[ameraner]

There seems to be a bug in the computation of the boundaries and intersections between a cropped geostationary projection, like SEVIRI RSS, and a lon-lat area.

To reproduce:

from pyresample.geometry import AreaDefinition
from satpy.area import get_area_def

rss_area = get_area_def("msg_seviri_rss_3km")
# the saved areadef in satpy is for full-disc, so lets use the actual area_extent from real data
rss_area_loc = rss_area.copy(area_extent=(-5570248.477339745, 1393687.2705221176, 5567248.074173927, 5570248.477339745))

ll_bbox = [0, 40, 1, 41]  # fails (South Spain)
# ll_bbox = [0, 50, 1, 51]  # works (South UK)
crop_area = AreaDefinition( 
    "crop_area", "crop_area", "crop_latlong",
    {"proj": "latlong"}, 100, 100, ll_bbox)


s = rss_area_loc.get_area_slices(crop_area)

which results in

  File "/tcenas/home/andream/.config/JetBrains/PyCharm2025.3/scratches/scratch_18.py", line 13, in <module>
    s = rss_area_loc.get_area_slices(crop_area)
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/tcenas/home/andream/code/satpy_latest/pyresample/pyresample/geometry.py", line 2648, in get_area_slices
    return get_area_slices(self, area_to_cover, shape_divisible_by)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/tcenas/home/andream/code/satpy_latest/pyresample/pyresample/_caching.py", line 49, in __call__
    return self._callable(*args)
           ^^^^^^^^^^^^^^^^^^^^^
  File "/tcenas/home/andream/code/satpy_latest/pyresample/pyresample/future/geometry/_subset.py", line 55, in get_area_slices
    raise NotImplementedError
NotImplementedError

as it fails to find an intersection here

pyresample/pyresample/future/geometry/_subset.py

Lines 48 to 55 in a7e1f25

	data_boundary = _get_area_boundary(src_area)
	area_boundary = _get_area_boundary(area_to_cover)
	intersection = data_boundary.contour_poly.intersection(
	area_boundary.contour_poly)
	if intersection is None:
	logger.debug('Cannot determine appropriate slicing. '
	"Data and projection area do not overlap.")
	raise NotImplementedError

The lon-lat bbox used is definitely in the area (it's somewhere in Spain). Using a bbox south of UK works again (see code example).

Additional context

This of course makes also the satpy Scene.crop fail accordingly, to reproduce:

import satpy

path = '/tcenas/fbf/MSG/in/MSG4/OPE4/SEVI-MSG15/MSG4-SEVI-MSG15-0100-NA-20251230083420.370000000Z-NA.nat'
scene = satpy.Scene(reader='seviri_l1b_native', filenames=[path], reader_kwargs={'fill_disk': False})

# load channel
channel = 'VIS008'
scene.load([channel], upper_right_corner='NE')

# crop on ll_box: (lon_min, lat_min, lon_max, lat_max)
ll_box = (0, 40, 1, 41)
scene_sub = scene.crop(ll_bbox=ll_box)

here, the workaround is to use 'fill_disk': True

Additional context II

@simonrp84 reported this already on Slack in 2023 https://pytroll.slack.com/archives/C0LNH7LMB/p1686061927840339 and then also for Himawary when passing a single segment, see https://pytroll.slack.com/archives/C0LNH7LMB/p1716292859837229 .

Additional context III

According to my colleague, it doesn't happen with pyresample v1.25.1

[djhoese]

Unfortunately various bugs being fixed in this portion of pyresample's code caused some previously working cases that worked "by chance" to stop working. At least that was my experience and understanding of things as they changed.

For the record, I'm guessing you're using the current newest pyresample version 1.35.0?

If you convert the lon/lat bbox to an x/y bbox by using pyproj to transform the coordinates manually, do you get a valid solution?

[ameraner]

Yes, the current pyresample main.

I guess you mean something like this, replicating what Scene.crop does for xy_bbox: this works!

from pyresample.geometry import AreaDefinition
from satpy.area import get_area_def

rss_area = get_area_def("msg_seviri_rss_3km")
rss_area_loc = rss_area.copy(area_extent=(-5570248.477339745, 1393687.2705221176, 5567248.074173927, 5570248.477339745))

ll_bbox = (0, 40, 1, 41)
xy_bbox = (
    rss_area_loc.get_projection_coordinates_from_lonlat(ll_bbox[0], ll_bbox[1])[0],
    rss_area_loc.get_projection_coordinates_from_lonlat(ll_bbox[0], ll_bbox[1])[1],
    rss_area_loc.get_projection_coordinates_from_lonlat(ll_bbox[2], ll_bbox[3])[0],
    rss_area_loc.get_projection_coordinates_from_lonlat(ll_bbox[2], ll_bbox[3])[1],
)

crop_area = AreaDefinition(
    "crop_area", "crop_area", "crop_xy",
    rss_area_loc.crs, rss_area_loc.width, rss_area_loc.height,
    xy_bbox)

s = rss_area_loc.get_area_slices(crop_area)

[djhoese]

I'm realizing now that this isn't technically the same operation unless the ll bbox is a rectangle in the xy space. Otherwise the fancy polygon math is needed. Plus I think there's a short cut in the get_area_slices logic if the 2 areas being the same CRS.

Anyway, I don't have time to track down all the issues and PRs and github discussions but there are existing things talking about:

1. Fixing some of the issues related to these problems as best we can. The problem is some of the PRs doing this are in a complex state where too much happened in one PR and was never fully cleaned up. Some of the fixes also seemed to be fixed most easily by causing backwards incompatibility or reinventing the interfaces. I can't pretend to understand all of this.
2. Switch to an entirely different third-party library to handle this polygon math. I found the discussion on their repository: Replacement for pyresample's spherical polygons benbovy/spherely#102

Why this particular case of yours is such an easy problem in human-space but so hard for the current algorithm I don't know.

[mraspaud]

Side question: @djhoese were you happy with the answers you got from the spherely maintainers, and thus should we attempt to plat a migration?

[djhoese]

I think the project is still experimental but if anyone has time then it is definitely worth trying a proof of concept for some of pyresample's use cases.