diff --git a/pyresample/slicer.py b/pyresample/slicer.py index 6a6ed9a4..0d32b07a 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 6efef48e..9c486a6d 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,