gradient_search leaves part of the destination empty for a partial geostationary disk
Nobody has claimed this yet.
- Dominant language
- Python
- Stars
- 385
- Forks
- 102
- Avg merge
- 4d 2h
- Merged PRs (30d)
- 9
Description
Describe the bug
When the source is a partial geostationary disk (for example MSG SEVIRI rapid scan data, which covers only the northern part of the disk, or a region of interest), resampling to an area in a different projection with gradient_search leaves a large part of the destination empty even where the source clearly has data. The populated region is bounded by a smooth arc.
The cause is in AreaSlicer.get_polygon_to_contain (pyresample/slicer.py). For a geostationary source it builds the bounding polygon of the disk and reprojects it into the destination projection to compute the source slice. When the disk is clipped to a partial area extent, the clip introduces straight chord edges (for instance the southern boundary of a rapid scan strip), each represented by only its two end points. Those chords are not straight in the destination projection, but with only two vertices the reprojected polygon keeps them straight and cuts the corner. The resulting source slice is too small, so resample_blocks fills the uncovered destination pixels with the fill value. A full disk source is unaffected, because its whole boundary is the densely sampled disk arc with no straight chord.
To Reproduce
from pyresample import AreaDefinition
from pyresample.slicer import create_slicer
# Partial geostationary disk: the northern 1392 lines of an MSG rapid scan image
source = AreaDefinition(
"rss", "rss", None,
{"a": 6378169.0, "b": 6356583.8, "h": 35785831.0, "lon_0": 9.5,
"proj": "geos", "units": "m"},
3712, 1392,
(5568748.2758, 5568748.2758, -5568748.2758, 1392187.0689))
# Polar stereographic destination over Europe, fully inside the disk
euro4 = AreaDefinition(
"euro4", "euro4", None,
{"proj": "stere", "ellps": "bessel", "lat_0": 90.0, "lon_0": 14.0,
"lat_ts": 60.0, "units": "m"},
1024, 1024,
(-2717181.7304994687, -5571048.14031214,
1378818.2695005313, -1475048.1403121399))
x_slice, y_slice = create_slicer(source, euro4).get_slices()
print(y_slice)
Expected behavior
The whole euro4 area falls inside the rapid scan strip, so the source slice should cover every source row the destination needs (down to about row 694), and a gradient_search resample should fill the whole destination.
Actual results
slice(1151, 1335, None)
Only the northern rows are kept (the slice should start near row 694). A gradient_search resample of this source to euro4 then fills only about 44 percent of the destination, with the populated area bounded by an arc.
Environment Info:
- pyresample: reproduced with 1.34.2 and 1.35.0 (the relevant code is unchanged on
main) - shapely: 2.x
- pyproj: 3.x
Additional context
Reported downstream in satpy as pytroll/satpy#3403. A fix that densifies the bounding polygon with segmentize before reprojecting it makes the example above fill the whole destination. Happy to open a pull request.
Contributor guide
No contributing guide indexed for this repository
First steps
- Read the whole issue, then the project's contributing guide.
- Comment on the issue to say you are picking it up — it saves two people doing the same work.
- Fork the repository and make your change on a branch.
- Open a pull request that references the issue number.
Research direction
Start by running the reproducer and inspect AreaSlicer.get_polygon_to_contain in pyresample/slicer.py, then follow how create_slicer produces the source slices. Verify that the partial geostationary boundary is represented accurately in the destination projection: the slice should begin near row 694 and a gradient_search resample should populate the whole euro4 destination.
Written by the indexing model from the issue text.
Assessment
- Tech stack
- python
- Domain
- data
- Issue type
- Bug
- Difficulty
- 3/5
- Estimated time
- 1-2 days
- Activity status
- Quiet
- Clarity
- Clearly specified
- Newbie friendliness
- 72/100