Refactor for boundary lat/lon extraction

pyresample/future/geometry/area.py

@@ -26,9 +26,9 @@
 from pyresample.geometry import AreaDefinition as LegacyAreaDefinition  # noqa
 from pyresample.geometry import (  # noqa
     DynamicAreaDefinition,
+    _get_geostationary_bounding_box_in_lonlats,
     get_full_geostationary_bounding_box_in_proj_coords,
     get_geostationary_angle_extent,
-    get_geostationary_bounding_box_in_lonlats,
     get_geostationary_bounding_box_in_proj_coords,
     ignore_pyproj_proj_warnings,
 )

pyresample/geometry.py

@@ -63,6 +63,7 @@
 HashType = hashlib._hashlib.HASH
 
 
+
 class DimensionError(ValueError):
     """Wrap ValueError."""
 
@@ -279,14 +280,19 @@ def get_boundary_lonlats(self):
         return (SimpleBoundary(s1_lon.squeeze(), s2_lon.squeeze(), s3_lon.squeeze(), s4_lon.squeeze()),
                 SimpleBoundary(s1_lat.squeeze(), s2_lat.squeeze(), s3_lat.squeeze(), s4_lat.squeeze()))
 
+    @property
+    def is_geostationary(self):
+        return False
+
     def get_bbox_lonlats(self, vertices_per_side: Optional[int] = None, force_clockwise: bool = True,
                          frequency: Optional[int] = None) -> tuple:
-        """Return the bounding box lons and lats.
+        """Return the bounding box lons and lats sides.
 
         Args:
             vertices_per_side:
                 The number of points to provide for each side. By default (None)
-                the full width and height will be provided.
+                the full width and height will be provided, except for geostationary
+                projections where by default only 50 points are selected.
             frequency:
                 Deprecated, use vertices_per_side
             force_clockwise:
@@ -318,14 +324,50 @@ def get_bbox_lonlats(self, vertices_per_side: Optional[int] = None, force_clockw
             warnings.warn("The `frequency` argument is pending deprecation, use `vertices_per_side` instead",
                           PendingDeprecationWarning, stacklevel=2)
         vertices_per_side = vertices_per_side or frequency
-        lons, lats = self._get_bbox_elements(self.get_lonlats, vertices_per_side)
+        if self.is_geostationary:
+            lon_sides, lat_sides = self._get_geo_boundary_sides(vertices_per_side=vertices_per_side)
+        else:
+            lon_sides, lat_sides = self._get_boundary_sides(self.get_lonlats, vertices_per_side)
         if force_clockwise and not self._corner_is_clockwise(
-                lons[0][-2], lats[0][-2], lons[0][-1], lats[0][-1], lons[1][1], lats[1][1]):
+                lon_sides[0][-2], lat_sides[0][-2],
+                lon_sides[0][-1], lat_sides[0][-1],
+                lon_sides[1][1], lat_sides[1][1]):
             # going counter-clockwise
-            lons, lats = self._reverse_boundaries(lons, lats)
-        return lons, lats
+            lon_sides, lat_sides = self._reverse_boundaries(lon_sides, lat_sides)
+        return lon_sides, lat_sides
+
+
+    def _get_geo_boundary_sides(self, vertices_per_side=None):
+        """Retrieve the boundary sides list for geostationary projections."""
+        # Define default frequency
+        if vertices_per_side is None:
+            vertices_per_side = 50
+        # Ensure at least 4 points are used
+        if vertices_per_side < 4:
+            vertices_per_side = 4
+        # Ensure an even number of vertices for side creation
+        if (vertices_per_side % 2) != 0:
+            vertices_per_side = vertices_per_side + 1
+        lons, lats = _get_geostationary_bounding_box_in_lonlats(self, nb_points=vertices_per_side)
+
+        # Retrieve dummy sides for GEO (side1 and side3 always of length 2)
+        # - TODO: _get_geostationary_bounding_box_in_lonlats now does not return nb_points !
+        # side02_step = int(vertices_per_side / 2) - 1 # old code
+        side02_step = int(lons.shape[0] / 2) - 1
+        lon_sides = [lons[slice(0, side02_step + 1)],
+                        lons[slice(side02_step, side02_step + 1 + 1)],
+                        lons[slice(side02_step + 1, side02_step * 2 + 1 + 1)],
+                        np.append(lons[side02_step * 2 + 1], lons[0])
+                        ]
+        lat_sides = [lats[slice(0, side02_step + 1)],
+                        lats[slice(side02_step, side02_step + 1 + 1)],
+                        lats[slice(side02_step + 1, side02_step * 2 + 1 + 1)],
+                        np.append(lats[side02_step * 2 + 1], lats[0])
+                        ]
+        return lon_sides, lat_sides
 
-    def _get_bbox_elements(self, coord_fun, vertices_per_side: Optional[int] = None) -> tuple:
+
+    def _get_boundary_sides(self, coord_fun, vertices_per_side: Optional[int] = None) -> tuple:
         s1_slice, s2_slice, s3_slice, s4_slice = self._get_bbox_slices(vertices_per_side)
         s1_dim1, s1_dim2 = coord_fun(data_slice=s1_slice)
         s2_dim1, s2_dim2 = coord_fun(data_slice=s2_slice)
@@ -337,10 +379,10 @@ def _get_bbox_elements(self, coord_fun, vertices_per_side: Optional[int] = None)
                            (s4_dim1.squeeze(), s4_dim2.squeeze())])
         if hasattr(dim1[0], 'compute') and da is not None:
             dim1, dim2 = da.compute(dim1, dim2)
-        clean_dim1, clean_dim2 = self._filter_bbox_nans(dim1, dim2)
+        clean_dim1, clean_dim2 = self._filter_sides_nans(dim1, dim2)
         return clean_dim1, clean_dim2
 
-    def _filter_bbox_nans(
+    def _filter_sides_nans(
             self,
             dim1_sides: list[np.ndarray],
             dim2_sides: list[np.ndarray],
@@ -423,17 +465,18 @@ def get_edge_bbox_in_projection_coordinates(self, vertices_per_side: Optional[in
             warnings.warn("The `frequency` argument is pending deprecation, use `vertices_per_side` instead",
                           PendingDeprecationWarning, stacklevel=2)
         vertices_per_side = vertices_per_side or frequency
-        x, y = self._get_bbox_elements(self.get_proj_coords, vertices_per_side)
+        x, y = self._get_boundary_sides(self.get_proj_coords, vertices_per_side)
         return np.hstack(x), np.hstack(y)
 
-    def boundary(self, vertices_per_side=None, force_clockwise=False, frequency=None):
+    def boundary(self, *, vertices_per_side=None, force_clockwise=False, frequency=None):
         """Retrieve the AreaBoundary object.
 
         Parameters
         ----------
         vertices_per_side:
              (formerly `frequency`) The number of points to provide for each side. By default (None)
-            the full width and height will be provided.
+            the full width and height will be provided, except for geostationary
+            projection where by default only 50 points are selected.
         force_clockwise:
             Perform minimal checks and reordering of coordinates to ensure
             that the returned coordinates follow a clockwise direction.
@@ -1570,65 +1613,6 @@ def is_geostationary(self):
             return False
         return 'geostationary' in coord_operation.method_name.lower()
 
-    def _get_geo_boundary_sides(self, vertices_per_side=None):
-        """Retrieve the boundary sides list for geostationary projections."""
-        # Define default frequency
-        if vertices_per_side is None:
-            vertices_per_side = 50
-        # Ensure at least 4 points are used
-        if vertices_per_side < 4:
-            vertices_per_side = 4
-        # Ensure an even number of vertices for side creation
-        if (vertices_per_side % 2) != 0:
-            vertices_per_side = vertices_per_side + 1
-        lons, lats = get_geostationary_bounding_box_in_lonlats(self, nb_points=vertices_per_side)
-        # Retrieve dummy sides for GEO (side1 and side3 always of length 2)
-        side02_step = int(vertices_per_side / 2) - 1
-        lon_sides = [lons[slice(0, side02_step + 1)],
-                     lons[slice(side02_step, side02_step + 1 + 1)],
-                     lons[slice(side02_step + 1, side02_step * 2 + 1 + 1)],
-                     np.append(lons[side02_step * 2 + 1], lons[0])
-                     ]
-        lat_sides = [lats[slice(0, side02_step + 1)],
-                     lats[slice(side02_step, side02_step + 1 + 1)],
-                     lats[slice(side02_step + 1, side02_step * 2 + 1 + 1)],
-                     np.append(lats[side02_step * 2 + 1], lats[0])
-                     ]
-        return lon_sides, lat_sides
-
-    def boundary(self, *, vertices_per_side=None, force_clockwise=False, frequency=None):
-        """Retrieve the AreaBoundary object.
-
-        Parameters
-        ----------
-        vertices_per_side:
-            The number of points to provide for each side. By default (None)
-            the full width and height will be provided, except for geostationary
-            projection where by default only 50 points are selected.
-        frequency:
-                Deprecated, use vertices_per_side
-        force_clockwise:
-            Perform minimal checks and reordering of coordinates to ensure
-            that the returned coordinates follow a clockwise direction.
-            This is important for compatibility with
-            :class:`pyresample.spherical.SphPolygon` where operations depend
-            on knowing the inside versus the outside of a polygon. These
-            operations assume that coordinates are clockwise.
-            Default is False.
-        """
-        from pyresample.boundary import AreaBoundary
-        if frequency is not None:
-            warnings.warn("The `frequency` argument is pending deprecation, use `vertices_per_side` instead",
-                          PendingDeprecationWarning, stacklevel=2)
-        vertices_per_side = vertices_per_side or frequency
-        if self.is_geostationary:
-            lon_sides, lat_sides = self._get_geo_boundary_sides(vertices_per_side=vertices_per_side)
-        else:
-            lon_sides, lat_sides = self.get_bbox_lonlats(vertices_per_side=vertices_per_side,
-                                                         force_clockwise=force_clockwise)
-        boundary = AreaBoundary.from_lonlat_sides(lon_sides, lat_sides)
-        return boundary
-
     @property
     def area_extent(self):
         """Tuple of this area's extent (xmin, ymin, xmax, ymax)."""
@@ -2743,9 +2727,10 @@ def geocentric_resolution(self, ellps='WGS84', radius=None):
 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)
+            vertices_per_side = None
+        else:
+            vertices_per_side = max(max(*area_to_cover.shape) // 100 + 1, 3)
+        return area_to_cover.boundary(vertices_per_side=vertices_per_side, force_clockwise=True)
     except ValueError:
         raise NotImplementedError("Can't determine boundary of area to cover")
 
@@ -2832,11 +2817,10 @@ def get_full_geostationary_bounding_box_in_proj_coords(geos_area, nb_points=50):
     y = -np.sin(points_around) * (y_max_angle - 0.0001)
     x *= h
     y *= h
-
     return x, y
 
 
-def get_geostationary_bounding_box_in_lonlats(geos_area, nb_points=50):
+def _get_geostationary_bounding_box_in_lonlats(geos_area, nb_points=50):
     """Get the bbox in lon/lats of the valid pixels inside `geos_area`.
 
     Args:
@@ -2847,16 +2831,28 @@ def get_geostationary_bounding_box_in_lonlats(geos_area, nb_points=50):
     return lons, lats
 
 
+def get_geostationary_bounding_box_in_lonlats(geos_area, nb_points=50):
+    """Get the bbox in lon/lats of the valid pixels inside `geos_area`.
+
+    Args:
+      nb_points: Number of points on the polygon
+    """
+    warnings.warn("'get_geostationary_bounding_box_in_lonlats' is deprecated. Please call "
+                  "'area.boundary().contour()' instead.",
+                  DeprecationWarning, stacklevel=2)
+    return _get_geostationary_bounding_box_in_lonlats(geos_area, nb_points)
+
+
 def get_geostationary_bounding_box(geos_area, nb_points=50):
     """Get the bbox in lon/lats of the valid pixels inside `geos_area`.
 
     Args:
       nb_points: Number of points on the polygon
     """
-    warnings.warn("'get_geostationary_bounding_box' is deprecated. Please use "
-                  "'get_geostationary_bounding_box_in_lonlats' instead.",
+    warnings.warn("'get_geostationary_bounding_box' is deprecated. Please call "
+                  "'area.boundary().contour()' instead.",
                   DeprecationWarning, stacklevel=2)
-    return get_geostationary_bounding_box_in_lonlats(geos_area, nb_points)
+    return _get_geostationary_bounding_box_in_lonlats(geos_area, nb_points)
 
 
 def combine_area_extents_vertical(area1, area2):

pyresample/gradient/__init__.py

@@ -39,7 +39,7 @@
 from pyresample.geometry import (
     AreaDefinition,
     SwathDefinition,
-    get_geostationary_bounding_box_in_lonlats,
+    _get_geostationary_bounding_box_in_lonlats,
 )
 from pyresample.gradient._gradient_search import (
     one_step_gradient_indices,
@@ -375,13 +375,13 @@
 
 def get_border_lonlats(geo_def: AreaDefinition):
     """Get the border x- and y-coordinates."""
+    # TODO: we could use geo_def.boundary().contour() here
     if geo_def.is_geostationary:
-        lon_b, lat_b = get_geostationary_bounding_box_in_lonlats(geo_def, 3600)
+        lon_b, lat_b = _get_geostationary_bounding_box_in_lonlats(geo_def, 3600)
     else:
         lons, lats = geo_def.get_boundary_lonlats()
         lon_b = np.concatenate((lons.side1, lons.side2, lons.side3, lons.side4))
         lat_b = np.concatenate((lats.side1, lats.side2, lats.side3, lats.side4))
-
     return lon_b, lat_b
 
 

pyresample/slicer.py

@@ -152,7 +152,7 @@ class AreaSlicer(Slicer):
     def get_polygon_to_contain(self):
         """Get the shapely Polygon corresponding to *area_to_contain* in projection coordinates of *area_to_crop*."""
         from shapely.geometry import Polygon
-        x, y = self.area_to_contain.get_edge_bbox_in_projection_coordinates(frequency=10)
+        x, y = self.area_to_contain.get_edge_bbox_in_projection_coordinates(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._transformer.transform(x_geos, y_geos, direction=TransformDirection.INVERSE)
@@ -180,7 +180,7 @@ def get_slices_from_polygon(self, poly_to_contain):
         except ValueError as err:
             raise InvalidArea(str(err))
         from shapely.geometry import Polygon
-        poly_to_crop = Polygon(zip(*self.area_to_crop.get_edge_bbox_in_projection_coordinates(frequency=10)))
+        poly_to_crop = Polygon(zip(*self.area_to_crop.get_edge_bbox_in_projection_coordinates(vertices_per_side=10)))
         if not poly_to_crop.intersects(buffered_poly):
             raise IncompatibleAreas("Areas not overlapping.")
         bounds = self._sanitize_polygon_bounds(bounds)

pyresample/test/test_geometry/test_area.py

@@ -27,9 +27,9 @@
 from pyresample import geo_filter, parse_area_file
 from pyresample.future.geometry import AreaDefinition, SwathDefinition
 from pyresample.future.geometry.area import (
+    _get_geostationary_bounding_box_in_lonlats,
     get_full_geostationary_bounding_box_in_proj_coords,
     get_geostationary_angle_extent,
-    get_geostationary_bounding_box_in_lonlats,
     get_geostationary_bounding_box_in_proj_coords,
     ignore_pyproj_proj_warnings,
 )
@@ -1493,8 +1493,7 @@ def test_get_full_geostationary_bbox(self, truncated_geos_area):
 
     def test_get_geostationary_bbox_works_with_truncated_area(self, truncated_geos_area):
         """Ensure the geostationary bbox works when truncated."""
-        lon, lat = get_geostationary_bounding_box_in_lonlats(truncated_geos_area, 20)
-
+        lon, lat = _get_geostationary_bounding_box_in_lonlats(truncated_geos_area, 20)
         expected_lon = np.array(
             [-64.24072434653284, -68.69662326361153, -65.92516214783112, -60.726360278290336,
              -47.39851775032484, 9.500000000000018, 66.39851775032487, 79.72636027829033,
@@ -1524,13 +1523,13 @@ def test_get_geostationary_bbox_works_with_truncated_area_proj_coords(self, trun
 
     def test_get_geostationary_bbox_does_not_contain_inf(self, truncated_geos_area):
         """Ensure the geostationary bbox does not contain np.inf."""
-        lon, lat = get_geostationary_bounding_box_in_lonlats(truncated_geos_area, 20)
+        lon, lat = _get_geostationary_bounding_box_in_lonlats(truncated_geos_area, 20)
         assert not any(np.isinf(lon))
         assert not any(np.isinf(lat))
 
     def test_get_geostationary_bbox_returns_empty_lonlats_in_space(self, truncated_geos_area_in_space):
         """Ensure the geostationary bbox is empty when in space."""
-        lon, lat = get_geostationary_bounding_box_in_lonlats(truncated_geos_area_in_space, 20)
+        lon, lat = _get_geostationary_bounding_box_in_lonlats(truncated_geos_area_in_space, 20)
 
         assert len(lon) == 0
         assert len(lat) == 0
@@ -1547,7 +1546,7 @@ def test_get_geostationary_bbox(self):
         geos_area.crs = CRS(proj_dict)
         geos_area.area_extent = [-5500000., -5500000., 5500000., 5500000.]
 
-        lon, lat = get_geostationary_bounding_box_in_lonlats(geos_area, 20)
+        lon, lat = _get_geostationary_bounding_box_in_lonlats(geos_area, 20)
         expected_lon = np.array([-78.19662326, -75.42516215, -70.22636028,
                                  -56.89851775, 0., 56.89851775, 70.22636028,
                                  75.42516215, 78.19662326, 79.23372832, 78.19662326,
@@ -1572,7 +1571,7 @@ def test_get_geostationary_bbox(self):
         geos_area.crs = CRS(proj_dict)
         geos_area.area_extent = [-5500000., -5500000., 5500000., 5500000.]
 
-        lon, lat = get_geostationary_bounding_box_in_lonlats(geos_area, 20)
+        lon, lat = _get_geostationary_bounding_box_in_lonlats(geos_area, 20)
         np.testing.assert_allclose(lon, expected_lon + lon_0)
 
     def test_get_geostationary_angle_extent(self):
@@ -1896,24 +1895,24 @@ def test_geostationary_projection(self, create_test_area):
 
         # Check default boundary shape
         default_n_vertices = 50
-        boundary = areadef.boundary(frequency=None)
+        boundary = areadef.boundary(vertices_per_side=None)
         assert boundary.vertices.shape == (default_n_vertices, 2)
 
         # Check minimum boundary vertices
         n_vertices = 3
         minimum_n_vertices = 4
-        boundary = areadef.boundary(frequency=n_vertices)
+        boundary = areadef.boundary(vertices_per_side=n_vertices)
         assert boundary.vertices.shape == (minimum_n_vertices, 2)
 
-        # Check odd frequency number
+        # Check odd number of vertices per side
         # - Rounded to the sequent even number (to construct the sides)
         n_odd_vertices = 5
-        boundary = areadef.boundary(frequency=n_odd_vertices)
+        boundary = areadef.boundary(vertices_per_side=n_odd_vertices)
         assert boundary.vertices.shape == (n_odd_vertices + 1, 2)
 
         # Check boundary vertices
         n_vertices = 10
-        boundary = areadef.boundary(frequency=n_vertices, force_clockwise=False)
+        boundary = areadef.boundary(vertices_per_side=n_vertices, force_clockwise=False)
 
         # Check boundary vertices is in correct order
         expected_vertices = np.array([[-7.54251621e+01, 3.53432890e+01],