pytroll / pytroll/pyresample

gradient_search leaves part of the destination empty for a partial geostationary disk

Open
#728 0 comments 0 reactions 0 assignees View on GitHub

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

  1. Read the whole issue, then the project's contributing guide.
  2. Comment on the issue to say you are picking it up — it saves two people doing the same work.
  3. Fork the repository and make your change on a branch.
  4. 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

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.