diff --git a/pyresample/geometry.py b/pyresample/geometry.py index a25839945..6e3de3179 100644 --- a/pyresample/geometry.py +++ b/pyresample/geometry.py @@ -37,7 +37,7 @@ from pyresample import CHUNK_SIZE from pyresample._spatial_mp import Cartesian, Cartesian_MP, Proj_MP from pyresample.area_config import create_area_def -from pyresample.boundary import AreaDefBoundary, Boundary, SimpleBoundary +from pyresample.boundary import Boundary, SimpleBoundary from pyresample.utils import check_slice_orientation, load_cf_area from pyresample.utils.proj4 import ( get_geostationary_height, @@ -2595,14 +2595,10 @@ def get_area_slices(self, area_to_cover, shape_divisible_by=None): y_slice = _ensure_integer_slice(y_slice) return x_slice, y_slice - if not self.is_geostationary: - raise NotImplementedError("Source projection must be 'geos' if " - "source/target projections are not " - "equal.") - - data_boundary = Boundary(*get_geostationary_bounding_box_in_lonlats(self)) - area_boundary = self._get_area_to_cover_boundary(area_to_cover) - intersection = data_boundary.contour_poly.intersection(area_boundary.contour_poly) + data_boundary = _get_area_boundary(self) + area_boundary = _get_area_boundary(area_to_cover) + intersection = data_boundary.contour_poly.intersection( + area_boundary.contour_poly) if intersection is None: logger.debug('Cannot determine appropriate slicing. ' "Data and projection area do not overlap.") @@ -2622,15 +2618,6 @@ def get_area_slices(self, area_to_cover, shape_divisible_by=None): return (check_slice_orientation(x_slice), check_slice_orientation(y_slice)) - @staticmethod - def _get_area_to_cover_boundary(area_to_cover: AreaDefinition) -> Boundary: - try: - if area_to_cover.is_geostationary: - return Boundary(*get_geostationary_bounding_box_in_lonlats(area_to_cover)) - return AreaDefBoundary(area_to_cover, 100) - except ValueError: - raise NotImplementedError("Can't determine boundary of area to cover") - def crop_around(self, other_area): """Crop this area around `other_area`.""" xslice, yslice = self.get_area_slices(other_area) @@ -2730,6 +2717,16 @@ def geocentric_resolution(self, ellps='WGS84', radius=None): return res +def _get_area_boundary(area_to_cover: AreaDefinition) -> Boundary: + try: + if area_to_cover.is_geostationary: + return Boundary(*get_geostationary_bounding_box_in_lonlats(area_to_cover)) + boundary_shape = max(max(*area_to_cover.shape) // 100 + 1, 3) + return area_to_cover.boundary(frequency=boundary_shape, force_clockwise=True) + except ValueError: + raise NotImplementedError("Can't determine boundary of area to cover") + + def _make_slice_divisible(sli, max_size, factor=2): """Make the given slice even in size.""" rem = (sli.stop - sli.start) % factor diff --git a/pyresample/test/test_geometry/test_area.py b/pyresample/test/test_geometry/test_area.py index 0c698a8df..d19b44944 100644 --- a/pyresample/test/test_geometry/test_area.py +++ b/pyresample/test/test_geometry/test_area.py @@ -1821,6 +1821,18 @@ def test_on_flipped_geos_area(self, create_test_area): assert slice_lines == expected_slice_lines assert slice_cols == expected_slice_cols + def test_non_geos_can_be_cropped(self, create_test_area): + """Test that non-geos areas can be cropped also.""" + src_area = create_test_area(dict(proj="utm", zone=33), + 10980, 10980, + (499980.0, 6490200.0, 609780.0, 6600000.0)) + crop_area = create_test_area({'proj': 'latlong'}, + 100, 100, + (15.9689, 58.5284, 16.4346, 58.6995)) + slice_x, slice_y = src_area.get_area_slices(crop_area) + assert slice_x == slice(5630, 8339) + assert slice_y == slice(9261, 10980) + def test_area_to_cover_all_nan_bounds(self, geos_src_area, create_test_area): """Check area slicing when the target doesn't have a valid boundary.""" area_def = geos_src_area