Fix get_area_slices under-crop for cropped geostationary cross-CRS areas

[Zaczero]

AreaDefinition.get_area_slices can under-crop cropped geostationary cross-CRS areas by merging the existing boundary-based slice with chunked destination-coverage bounds.

Before (notice cropping at the bottom)

After

---

• Closes Data reduction in resampling yields incorrect result satpy#986
• Tests added
• Tests passed

[djhoese]

Wow! I hope to take a look at all your PRs today or at least start. Since this is the first PR I have to ask: did you have all of these fixes/changes sitting around and waited to submit them or did you use AI to find the problems or...how did you have so many problems and have the fixes for them to make so many pull requests in one weekend?

[Zaczero]

@djhoese I used AI to investigate the problems in my local pipeline, prioritizing correctness bugs and performance bottlenecks. I then reviewed and iterated on fixes (human-in-the-loop). I am not a vibe coder if that's what you're worrying about; otherwise the code would be slop! Plus I am a workaholic.

I am new to this project, so obviously I may have missed some things.

[djhoese]

Thanks for putting this together. This is a very interesting idea and I've made some comments inline. I think my main concerns are:

1. Performance. You have a for loop in Python that is making 2D chunk arrays the size of the destination area.
2. The right solution: We (pyresample maintainers) know that the current bounding logic is not good and has a lot of flaws and we have some alternative implementations (PRs and third-party packages) that we want to investigate. I'm wondering if this workaround in this PR should be included because it gets us closer to the answer we want or if this is just another sign that the current implementation needs to be improved.

@pnuu @mraspaud I'd love to hear what you guys think.

[djhoese]

How was this value decided on?

[Zaczero]

768*768 is just under 600k. 1024x1024 is 1M. I tested several values and picked the one that in my opinion best matches CPU/memory tradeoffs. 1M was minimally faster for bigger images but used much more memory. Anything lower was not saving much memory and performance got lower.

[djhoese]

This concerns me efficiency wise. If I'm not mistaken this is making a full 2D array, right? For real-world areas this could 1000s of rows and 1000s of columns per chunk.

This also assumes that the area_to_cover is best represented by longitude/latitude degrees which is not always the case or even not usually the case. It seems like it would be better (although not implemented in the AreaDefinition objects already) to get the X/Y coordinates in the area_to_cover CRS then use pyproj to convert to the X/Y of the src_area then convert those to indices. That way lon/lat are never involved.

[djhoese]

I'm realizing now that even the main boundary code assumes lon/lat so maybe this idea of lon/lat being used is less of a problem for the current implementation.

[Zaczero]

As I understand it, the algorithm needs to check every point. This cross-CRS case is now 30% slower. There's some reason why a simple boundary check (not area) would not be correct for every cross-CRS, but I don't remember it now. I'll recheck that and get back to you with concrete information.

[Zaczero]

So yeah, for maximum correctness, all points need to be sampled. GDAL samples only edges by default and provides options to enable sampling of all points on a case-by-case basis. It works most of the time, but not always.

Normally when computing the source raster data to load to generate a particular output area, the warper samples transforms 21 points along each edge of the destination region back onto the source file, and uses this to compute a bounding window on the source image that is sufficient. Depending on the transformation in effect, the source window may be a bit too small, or even missing large areas. Problem situations are those where the transformation is very non-linear or "inside out". Examples are transforming from WGS84 to Polar Stereographic for areas around the pole, or transformations where some of the image is untransformable. The following options provide some additional control to deal with errors in computing the source window:

https://gdal.org/en/stable/doxygen/structGDALWarpOptions.html

So the existing implementation equals to SAMPLE_GRID=YES + SAMPLE_STEPS=ALL in GDAL.

Some other libs also use this 21 number by default.

[Zaczero]

I implemented the 21 sampled approach with 2 new flags:

sample_steps: int | None = 21,
sample_grid: bool = False,

sample_steps=None, sample_grid=True, would work as the before this patch.

In my case it improved the performance 94ms -> 3ms with minimal numerical drift (<=1px).

[djhoese]

This worries me that there isn't a clear "this is why this error happens" definition and is instead "catch any problems". I'd like if it was clear why some of these errors are happening and to try to handle them before they happen. Or at the very least comment on why they would happen.

[Zaczero]

Yes, I tried to follow project patterns with that one. Such a pattern is very common in the codebase. Generic except catches.

[djhoese]

Such a pattern is very common in the codebase

Yes, this project has never had the strictest of coding rules.

[djhoese]

Is this different than get_boundary_lonlats?

pyresample/pyresample/geometry.py

Lines 293 to 294 in 07711c0

	def get_bbox_lonlats(self, vertices_per_side: Optional[int] = None, force_clockwise: bool = True,
	frequency: Optional[int] = None) -> tuple:

[Zaczero]

Yes. _get_edge_lonlat_samples is geostationary-aware, and sampling follows the visibility curve rather than the image rectangular boundary. get_bbox_lonlats could pick samples that are not the true visible boundary. Plus it's more efficient because it's purpose build.

[djhoese]

The pyproj Proj class is less preferred these days (deprecated even?) in favor of Transformer. This is an even stronger argument for maybe looking at not using lons/lats as the intermediate value. By that I mean since you have two define the src and dst CRS in a Transformer you might as well go from source area CRS to area_to_cover CRS. At least that's my gut reaction as I look at this.

[Zaczero]

Transformer.from_proj is deprecated, but Proj is not. Sample points are in lon/lat units, so Proj is the cleanest API to use in this case: we do lon/lat->CRS, not general CRS->CRS. get_projection_coordinates_from_lonlat uses a similar technique.

[djhoese]

Hm I could have sworn there was an explicit deprecation, but can't find it now. However, Proj is limited datum-wise:

https://pyproj4.github.io/pyproj/stable/gotchas.html#proj-not-a-generic-latitude-longitude-to-projection-converter

[djhoese]

I could be wrong but I don't think this logic matches GDAL. I read the GDAL docs as sample_grid=False and sample_steps being ALL (0/None in our case) that that would do all the pixels along the edges. When sample_grid=True then sample_steps=ALL/None is the full dense grid (edges and all internal points).

[Zaczero]

Fixed that

[djhoese]

I've restarted CI as it had a weird failure. Let's see how it does now...

[codecov]

Codecov Report

❌ Patch coverage is 95.28302% with 10 lines in your changes missing coverage. Please review.
✅ Project coverage is 93.71%. Comparing base (4da28b0) to head (7c0cc75).
⚠️ Report is 103 commits behind head on main.

Files with missing lines	Patch %	Lines
pyresample/future/geometry/_subset.py	91.66%	9 Missing ⚠️
pyresample/test/test_geometry/test_area.py	99.00%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #700      +/-   ##
==========================================
+ Coverage   93.67%   93.71%   +0.04%     
==========================================
  Files          89       89              
  Lines       13621    13863     +242     
==========================================
+ Hits        12759    12992     +233     
- Misses        862      871       +9     

Flag	Coverage Δ
unittests	93.71% <95.28%> (+0.04%)	⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[Manny7717]

Verified locally on head 7c0cc75 (checkout of pull/700/head, merge-base d0cecf5 = origin/main):

Bug and fix direction confirmed. For cropped geostationary source areas in cross-projection mode, the old boundary-intersection path under-cropped (existing test_get_area_slices_geos_stereographic expectation moved 56..3659 → 52..3663 — strictly wider, the conservative direction; no narrowed expectations anywhere). The new sampled-coverage path (_get_covered_source_slices) maps sampled destination points into source pixel space and takes the min/max envelope, clamped to the source area, with a clean fallback to the old boundary path when no sample is valid.

Regression-proven: head test file applied to a base worktree → 15/15 new cross-projection tests FAIL on base (sample_steps handling, grid/edge relationships, step=1 rejection, short-circuit, boundary fallback) and PASS on head. Full suites on head: test_geometry/ 560 passed, slice-related tests 89 passed, test_area.py 224 passed.

Design checks:

• _normalize_sample_steps: None/<=0/False → ALL mode; 1 → clear ValueError (no existing caller passes 1; default is 21, so no breakage).
• Dense grid mode is a superset of edge mode (asserted); ALL-edge ⊇ default-edge ⊇ sampled-edge monotonicity asserted in _assert_default_matches_edge_and_is_within_dense.
• Memory-bounded: dense mode chunks at MAX_POINTS_PER_CHUNK=600k; default mode is 84 edge points — cheap.
• Same-CRS fast path preserved; non-geostationary sources unchanged (boundary intersection); _finalize_slices keeps orientation + shape_divisible_by handling identical.
• Fallback safety: sampled path returns None (→ old boundary path) on any projection error or empty coverage.

Consideration (non-blocking, worth documenting): the new default (sample_steps=21 edge sampling) applies to ALL geostationary cross-CRS get_area_slices calls, including downstream satpy usage — it replaces boundary intersection with a 21-point-per-side perimeter approximation. For strongly concave destination areas this is a slight approximation vs the exact polygon intersection; the previous behavior was actively wrong for cropped geostationary sources (this fixes satpy #986), and dense mode is available for exactness. A one-line note in the docstring that edge sampling is a conservative approximation (and when to use sample_steps=None, sample_grid=True) would preempt user confusion.

Also noting the maintainer thread (djhoese restarted CI; AI-assistance disclosed by author and handled in-thread — not my call to weigh in on). APPROVE.

[djhoese]

Thank you for all the work on this PR, but after testing other solutions like #729 and using Claude to migrate it to the area slicing in #757, I don't think we should merge this PR. I've been going back and forth and was concerned by the large number of changes here. I had Claude run some additional test cases and it found some issues which push me even further into closing this PR and not merging it. Here is Claude's summary:

Root cause

Both paths build the source outline with get_geostationary_bounding_box_in_proj_coords. For a partial disk (rapid scan, area-of-interest crops), clipping the disk to the area extent leaves each cut edge as a straight line with only two vertices. In get_area_slices that outline goes into a spherical polygon intersection, where the two-vertex edge becomes one long great-circle arc that bulges toward the pole. Part of the destination falls outside the intersection and the slice comes out too small. #729 fixes this in AreaSlicer by adding vertices along the clipped edges before reprojecting, and the same fix applies to _get_area_boundary here.

Problems with the sampling approach

1. The "edge" samples aren't edge points for non-geostationary destinations. _get_edge_lonlat_samples passes paired row and column index arrays to get_lonlats(data_slice=(rows, cols)), but get_lonlats builds the full row × column grid from them. For a 1024×1024 destination with the default 21 steps, that is an 80×80 grid instead of 80 perimeter points. For example, area.get_lonlats(data_slice=(np.array([0, 4, 4]), np.array([0, 0, 9]))) returns 3×3 arrays.
2. Sampled destination points only bound the footprint when the whole destination lies inside the source's valid region. When the source's extent edge or the Earth's limb passes through the destination's interior, the outermost source pixels fall between samples and the slice is too small.
3. Some cases that were safe are now silently wrong. Where main raises NotImplementedError, Satpy skips data reduction and the output is complete. With this PR, some of those cases return a slice that is too small and data is dropped without any error. test_area_to_cover_all_nan_bounds was changed to accept this.
4. The branch for geostationary destinations has the original bug. It uses get_geostationary_bounding_box_in_lonlats, which has the same two-vertex clipped edges.
5. The new tests miss these cases. They only use a destination that is fully inside the source.

Measured results

Source pixels missing from the get_area_slices result, compared with a dense ground truth (every destination pixel mapped into the source):

Source → destination	main	this PR
MSG RSS strip → euro4 (#728)	3 cols, 274 rows	0
SEVIRI Europe crop → eurol	0	21 cols, 3 rows
RSS strip → lon/lat box crossing the strip's southern edge	NotImplementedError	7 cols, 38 rows
Full disk → global lon/lat	4 cols, 10 rows	92 cols, 2 rows
Full disk → global Mollweide	NotImplementedError	92 cols, 33 rows
RSS strip → global lon/lat	4 cols, 5 rows	149 cols, 102 rows
RSS strip → geostationary strip at lon_0=0	257 cols, 320 rows	1314 cols, 12 rows

Script used for the table

import warnings

import numpy as np
from pyproj import Transformer

from pyresample import AreaDefinition

warnings.filterwarnings("ignore")

GEOS = {"a": 6378169.0, "b": 6356583.8, "h": 35785831.0, "lon_0": 9.5, "proj": "geos", "units": "m"}
FULL = 5568748.2758
LONLAT = {"proj": "longlat", "datum": "WGS84"}

full_disk = AreaDefinition.from_extent("full", GEOS, (3712, 3712), (FULL, FULL, -FULL, -FULL))
# Northern 1392 lines of MSG rapid scan (pyresample#728)
rss = AreaDefinition.from_extent("rss", GEOS, (1392, 3712), (FULL, FULL, -FULL, 1392187.0689))
# A SEVIRI area-of-interest crop over Europe
aoi = AreaDefinition.from_extent("aoi", GEOS, (700, 1500), (2000000.0, 5400000.0, -2500000.0, 3300000.0))
euro4 = AreaDefinition.from_extent(
    "euro4", {"proj": "stere", "ellps": "bessel", "lat_0": 90.0, "lon_0": 14.0, "lat_ts": 60.0}, (1024, 1024),
    (-2717181.7304994687, -5571048.14031214, 1378818.2695005313, -1475048.1403121399))
eurol = AreaDefinition.from_extent(
    "eurol", {"proj": "stere", "ellps": "WGS84", "lat_0": 90.0, "lon_0": 0.0, "lat_ts": 60.0}, (2048, 2560),
    (-3780000.0, -7644000.0, 3900000.0, -1500000.0))
straddle = AreaDefinition.from_extent("straddle", LONLAT, (800, 800), (-10.0, 0.0, 30.0, 40.0))
global_ll = AreaDefinition.from_extent("global", LONLAT, (720, 1440), (-180.0, -90.0, 180.0, 90.0))
moll = AreaDefinition.from_extent("moll", {"proj": "moll"}, (1000, 1000),
                                  (-18000000.0, -9000000.0, 18000000.0, 9000000.0))
geos_0 = AreaDefinition.from_extent("geos_0", dict(GEOS, lon_0=0.0), (400, 1000), (FULL, FULL, -FULL, 1392187.0689))


def needed_slices(src, dst):
    """Source cols/rows hit by any destination pixel center (dense ground truth)."""
    lons, lats = dst.get_lonlats()
    x, y = Transformer.from_crs("EPSG:4326", src.crs, always_xy=True).transform(lons, lats)
    cols, rows = src.get_array_indices_from_projection_coordinates(x, y)
    valid = ~np.ma.getmaskarray(cols) & ~np.ma.getmaskarray(rows)
    cols, rows = np.ma.getdata(cols)[valid], np.ma.getdata(rows)[valid]
    return slice(int(cols.min()), int(cols.max()) + 1), slice(int(rows.min()), int(rows.max()) + 1)


def missing(got, need):
    return max(0, got.start - need.start) + max(0, need.stop - got.stop)


for src, dst in [(rss, euro4), (aoi, eurol), (rss, straddle), (full_disk, global_ll), (full_disk, moll),
                 (rss, global_ll), (rss, geos_0)]:
    need_x, need_y = needed_slices(src, dst)
    try:
        got_x, got_y = src.get_area_slices(dst)
        result = f"missing {missing(got_x, need_x)} cols, {missing(got_y, need_y)} rows"
    except NotImplementedError:
        result = "NotImplementedError"
    print(f"{src.area_id:>4} -> {dst.area_id:<8} {result}")

Smaller concerns

• sample_steps and sample_grid become public arguments on AreaDefinition.get_area_slices (and BaseDefinition) with surprising values: None, <= 0 and False all mean "all points", and 1 raises. They expose an implementation detail that only applies to geostationary sources.
• pyproj.Proj is used instead of a Transformer, so datum differences between the two areas are ignored.
• The broad except (RuntimeError, TypeError, ValueError) silently falls back to the old path, which can hide bugs.
• The change to how the two CRSes are compared isn't related to the fix.
• The existing isinstance(slice_x.start, int) assertions were removed.