Merge pull request #516 from mraspaud/feature-non-geos-crop

Allow cropping non-geos areas

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

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