Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
13 changes: 12 additions & 1 deletion pyresample/slicer.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
19 changes: 19 additions & 0 deletions pyresample/test/test_slicer.py
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
Loading