From 99f9c93c589442f92f096f5694640234ab49dc11 Mon Sep 17 00:00:00 2001 From: Lorenzo Cingano Date: Mon, 29 Jun 2026 17:51:59 +0200 Subject: [PATCH 1/4] Densify the geostationary bounding polygon before reprojecting in AreaSlicer AreaSlicer.get_polygon_to_contain builds the bounding polygon of a geostationary source disk and reprojects it into the destination projection to compute the source slice. When the source is a partial disk (rapid scan or region of interest data), clipping the disk to its area extent introduces straight chord edges with only two vertices. Those chords are not straight in the destination projection, so the reprojected polygon cuts the corner and the source slice comes out too small, leaving part of the destination without data. Densify the polygon with segmentize before reprojecting it, so the chord edges keep enough vertices to follow their true shape in the destination projection. --- pyresample/slicer.py | 13 ++++++++++++- pyresample/test/test_slicer.py | 19 +++++++++++++++++++ 2 files changed, 31 insertions(+), 1 deletion(-) diff --git a/pyresample/slicer.py b/pyresample/slicer.py index 6a6ed9a4e..0d32b07a4 100644 --- a/pyresample/slicer.py +++ b/pyresample/slicer.py @@ -175,7 +175,18 @@ def get_polygon_to_contain(self): x, y = self.area_to_contain.get_edge_lonlats(vertices_per_side=10) if self.area_to_crop.is_geostationary: x_geos, y_geos = get_geostationary_bounding_box_in_proj_coords(self.area_to_crop, 360) - x_geos, y_geos = self._source_transformer.transform(x_geos, y_geos, direction=TransformDirection.INVERSE) + geos_poly = Polygon(zip(x_geos, y_geos, strict=True)) + # Densify the polygon before reprojecting it to the destination crs. + # Clipping a partial disk (for example rapid scan or region of interest + # data) to its area extent introduces straight chord edges. Those chords + # are not straight in the destination crs, so without intermediate + # vertices the reprojected polygon cuts the corner and the computed slice + # is too small, leaving part of the destination without data. + segment_length = np.max(np.abs(self.area_to_crop.area_extent)) / 100 + geos_poly = geos_poly.segmentize(segment_length) + x_geos, y_geos = geos_poly.exterior.coords.xy + x_geos, y_geos = self._source_transformer.transform( + np.asarray(x_geos), np.asarray(y_geos), direction=TransformDirection.INVERSE) geos_poly = Polygon(zip(x_geos, y_geos, strict=True)) poly = Polygon(zip(x, y, strict=True)) poly = poly.intersection(geos_poly) diff --git a/pyresample/test/test_slicer.py b/pyresample/test/test_slicer.py index 6efef48e7..9c486a6db 100644 --- a/pyresample/test/test_slicer.py +++ b/pyresample/test/test_slicer.py @@ -64,6 +64,25 @@ def test_source_area_does_not_cover_dest_area_entirely(self): assert x_slice.start > 0 and x_slice.stop < 100 assert y_slice.start > 0 and y_slice.stop >= 100 + def test_partial_geostationary_disk_covering_dest_area_is_not_truncated(self): + """Test that a partial geostationary disk covering the destination is not truncated. + + The southern boundary of a partial disk is a straight chord in the source + projection. It has to be densified before being reprojected to the + destination crs, otherwise the computed source slice stops short and part + of the destination is left without data. + """ + src_area = AreaDefinition('rss', 'rss area', None, + {'ellps': 'WGS84', 'h': '35785831', 'proj': 'geos', 'lon_0': 9.5}, + 3712, 1392, + (5550000.0, 5550000.0, -5550000.0, 2000000.0)) + slicer = create_slicer(src_area, self.dst_area) + x_slice, y_slice = slicer.get_slices() + # The whole destination falls inside the strip; its southern edge maps to + # source row 584, so the slice must reach well below the strip top rows. + assert y_slice.start < 700 + assert y_slice.stop >= 1328 + def test_source_area_does_not_cover_dest_area_at_all(self): """Test source area does not cover dest area at all.""" src_area = AreaDefinition('dst', 'dst area', None, From b4172cd804c76f18c70c49a188bd018d77c54f5c Mon Sep 17 00:00:00 2001 From: David Hoese Date: Wed, 7 Oct 2026 16:24:11 -0500 Subject: [PATCH 2/4] Require shapely 2.0+ --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 61c2cd072..4d553e6db 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -19,7 +19,7 @@ dependencies = [ "pykdtree>=1.3.1", "pyyaml", "numpy>=1.25.0", - "shapely", + "shapely>=2.0", "donfig", "platformdirs", ] From 29176d02cb309f0c8d836db61b3b5a5f87b2b7b5 Mon Sep 17 00:00:00 2001 From: David Hoese Date: Wed, 7 Oct 2026 16:24:47 -0500 Subject: [PATCH 3/4] Apply suggestion from @djhoese --- pyresample/slicer.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pyresample/slicer.py b/pyresample/slicer.py index 0d32b07a4..bede7988b 100644 --- a/pyresample/slicer.py +++ b/pyresample/slicer.py @@ -182,7 +182,8 @@ def get_polygon_to_contain(self): # are not straight in the destination crs, so without intermediate # vertices the reprojected polygon cuts the corner and the computed slice # is too small, leaving part of the destination without data. - segment_length = np.max(np.abs(self.area_to_crop.area_extent)) / 100 + atc_extent = self.area_to_crop.area_extent + segment_length = max(abs(atc_extent[2] - atc_extent[0]), abs(atc_extent[3] - atc_extent[1])) / 100 geos_poly = geos_poly.segmentize(segment_length) x_geos, y_geos = geos_poly.exterior.coords.xy x_geos, y_geos = self._source_transformer.transform( From 106529172e2c5e266db4a912131a1c4c853a40fd Mon Sep 17 00:00:00 2001 From: David Hoese Date: Wed, 7 Oct 2026 21:03:08 -0500 Subject: [PATCH 4/4] Add stricter slice check --- pyresample/test/test_slicer.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/pyresample/test/test_slicer.py b/pyresample/test/test_slicer.py index 9c486a6db..e891e2704 100644 --- a/pyresample/test/test_slicer.py +++ b/pyresample/test/test_slicer.py @@ -80,8 +80,10 @@ def test_partial_geostationary_disk_covering_dest_area_is_not_truncated(self): x_slice, y_slice = slicer.get_slices() # The whole destination falls inside the strip; its southern edge maps to # source row 584, so the slice must reach well below the strip top rows. - assert y_slice.start < 700 - assert y_slice.stop >= 1328 + assert y_slice.start <= 584 + assert y_slice.stop >= 1329 + assert x_slice.start <= 1357 + assert x_slice.stop >= 2620 def test_source_area_does_not_cover_dest_area_at_all(self): """Test source area does not cover dest area at all."""