From 99f9c93c589442f92f096f5694640234ab49dc11 Mon Sep 17 00:00:00 2001 From: Lorenzo Cingano Date: Mon, 29 Jun 2026 17:51:59 +0200 Subject: [PATCH] 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,