Refactor for boundary lat/lon extraction

docs/source/howtos/spherical_geometry.rst

@@ -72,14 +72,13 @@ satellite passes. See trollschedule_ how to generate a list of satellite overpas
 `area_def` is an :class:`~pyresample.geometry.AreaDefinition` object.
 
  >>> from pyresample.spherical_utils import GetNonOverlapUnions
- >>> from pyresample.boundary import AreaDefBoundary
 
- >>> area_boundary = AreaDefBoundary(area_def, frequency=100) # doctest: +SKIP
- >>> area_boundary = area_boundary.contour_poly # doctest: +SKIP
+ >>> area_boundary = area_def.boundary(vertices_per_side=100) # doctest: +SKIP
+ >>> area_boundary = area_boundary.polygon # doctest: +SKIP
 
  >>> list_of_polygons = []
  >>> for mypass in passes: # doctest: +SKIP
- >>>     list_of_polygons.append(mypass.boundary.contour_poly) # doctest: +SKIP
+ >>>     list_of_polygons.append(mypass.boundary().polygon) # doctest: +SKIP
 
  >>> non_overlaps = GetNonOverlapUnions(list_of_polygons) # doctest: +SKIP
  >>> non_overlaps.merge() # doctest: +SKIP

pyresample/__init__.py

@@ -56,8 +56,13 @@
 
 from .version import get_versions  # noqa
 
-__all__ = ['grid', 'image', 'kd_tree', 'utils', 'plot', 'geo_filter', 'geometry', 'CHUNK_SIZE',
-           'load_area', 'create_area_def', 'get_area_def', 'parse_area_file', 'convert_def_to_yaml']
+_root_path = os.path.dirname(os.path.realpath(__file__))
+
+__all__ = [
+    'grid', 'image', 'kd_tree', 'utils', 'plot', 'geo_filter', 'geometry', 'CHUNK_SIZE',
+    'load_area', 'create_area_def', 'get_area_def', 'parse_area_file', 'convert_def_to_yaml',
+    "_root_path",
+]
 
 __version__ = get_versions()['version']
 del get_versions

pyresample/_config.py

@@ -40,6 +40,7 @@
         "features": {
             "future_geometries": False,
         },
+        "force_boundary_computations": False,
     }],
     paths=_CONFIG_PATHS,
 )

pyresample/bilinear/_base.py

@@ -599,14 +599,14 @@
         target_geo_def, source_lons, source_lats, radius_of_influence):
     """Get valid indices from lonlat boundaries."""
     # Resampling from swath to grid or from grid to grid
-    lonlat_boundary = target_geo_def.get_boundary_lonlats()
-
-    # Combine reduced and legal values
-    return data_reduce.get_valid_index_from_lonlat_boundaries(
-        lonlat_boundary[0],
-        lonlat_boundary[1],
-        source_lons, source_lats,
-        radius_of_influence)
+    try:
+        sides_lons, sides_lats = target_geo_def.boundary().sides
+        valid_indices = data_reduce.get_valid_index_from_lonlat_boundaries(sides_lons, sides_lats,
+                                                                           source_lons, source_lats,
+                                                                           radius_of_influence)
+    except Exception:
+        valid_indices = np.ones(source_lons.size, dtype=bool)
+    return valid_indices
 
 
 def get_slicer(data):

pyresample/boundary/__init__.py

@@ -0,0 +1,32 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# Copyright (c) 2014-2021 Pyresample developers
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+"""The Boundary classes."""
+
+from pyresample.boundary.area_boundary import AreaBoundary, AreaDefBoundary
+from pyresample.boundary.planar_boundary import PlanarBoundary
+from pyresample.boundary.simple_boundary import SimpleBoundary
+from pyresample.boundary.spherical_boundary import SphericalBoundary
+
+__all__ = [
+    "SphericalBoundary",
+    "PlanarBoundary",
+    # Deprecated
+    "SimpleBoundary",
+    "AreaBoundary",
+    "AreaDefBoundary",
+]

pyresample/boundary.py → pyresample/boundary/area_boundary.py

@@ -1,7 +1,7 @@
 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
 #
-# Copyright (c) 2014-2021 Pyresample developers
+# Copyright (c) 2014-2023 Pyresample developers
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
@@ -15,8 +15,7 @@
 #
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
-"""The Boundary classes."""
-
+"""Deprecated Boundary, AreaBoundary, AreaDefBoundary class."""
 import logging
 import warnings
 
@@ -44,6 +43,9 @@
     @property
     def contour_poly(self):
         """Get the Spherical polygon corresponding to the Boundary."""
+        warnings.warn("'contour_poly' is deprecated." +
+                      "Use the 'boundary().polygon' property instead!.",
+                      PendingDeprecationWarning, stacklevel=2)
         if self._contour_poly is None:
             self._contour_poly = SphPolygon(
                 np.deg2rad(np.vstack(self.contour()).T))
@@ -62,6 +64,9 @@
 
     def __init__(self, *sides):
         Boundary.__init__(self)
+        warnings.warn("'AreaBoundary' will be removed in the future. " +
+                      "Use the Swath/AreaDefinition 'boundary' method instead!.",
+                      PendingDeprecationWarning, stacklevel=2)
         # Check 4 sides are provided
         if len(sides) != 4:
             raise ValueError("AreaBoundary expects 4 sides.")
@@ -117,28 +122,13 @@
 
 
 class AreaDefBoundary(AreaBoundary):
-    """Boundaries for area definitions (pyresample)."""
+    """Boundaries for a pyresample AreaDefinition."""
 
     def __init__(self, area, frequency=1):
-        lons, lats = area.get_bbox_lonlats()
+        sides_lons, sides_lats = area.boundary().sides
         warnings.warn("'AreaDefBoundary' will be removed in the future. " +
                       "Use the Swath/AreaDefinition 'boundary' method instead!.",
                       PendingDeprecationWarning, stacklevel=2)
-        AreaBoundary.__init__(self,
-                              *zip(lons, lats))
-
+        AreaBoundary.__init__(self, *zip(sides_lons, sides_lats))
         if frequency != 1:
             self.decimate(frequency)
-
-
-class SimpleBoundary(object):
-    """Container for geometry boundary.
-
-    Labelling starts in upper left corner and proceeds clockwise
-    """
-
-    def __init__(self, side1, side2, side3, side4):
-        self.side1 = side1
-        self.side2 = side2
-        self.side3 = side3
-        self.side4 = side4

pyresample/boundary/base_boundary.py

@@ -0,0 +1,127 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# Copyright (c) 2014-2023 Pyresample developers
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+"""Define the BaseBoundary class."""
+
+import logging
+
+import numpy as np
+
+from pyresample.boundary.sides import BoundarySides
+
+logger = logging.getLogger(__name__)
+
+
+class BaseBoundary:
+    """Base class for boundary objects."""
+    __slots__ = ["_sides_x", "_sides_y"]
+
+    def __init__(self, area, vertices_per_side=None):
+
+        sides_x, sides_y = self._compute_boundary_sides(area, vertices_per_side)
+        self._sides_x = BoundarySides(sides_x)
+        self._sides_y = BoundarySides(sides_y)
+        self._area = area
+
+        self.is_clockwise = self._check_is_boundary_clockwise(sides_x, sides_y, area)
+        self.is_counterclockwise = not self.is_clockwise
+        self._set_order(order=None)  # FIX !
+
+    def _check_is_boundary_clockwise(self, sides_x, sides_y, area):
+        """Check if the boundary is clockwise or counterclockwise."""
+        raise NotImplementedError()
+
+    def _compute_boundary_sides(self, area, vertices_per_side):
+        """Compute boundary sides."""
+        raise NotImplementedError()
+
+    def _set_order(self, order):
+        """Set the order of the boundary vertices."""
+        if self.is_clockwise:
+            self._actual_order = "clockwise"
+        else:
+            self._actual_order = "counterclockwise"
+
+        if order is None:
+            self._wished_order = self._actual_order
+        else:
+            if order not in ["clockwise", "counterclockwise"]:
+                raise ValueError("Valid 'order' is 'clockwise' or 'counterclockwise'")
+            self._wished_order = order
+
+    def set_clockwise(self):
+        """Set clockwise order for vertices retrieval."""
+        self._wished_order = "clockwise"
+        return self
+
+    def set_counterclockwise(self):
+        """Set counterclockwise order for vertices retrieval."""
+        self._wished_order = "counterclockwise"
+        return self
+
+    @property
+    def sides(self):
+        """Return the boundary sides as a tuple of (sides_x, sides_y) arrays."""
+        return self._sides_x, self._sides_y
+
+    @property
+    def _x(self):
+        """Retrieve boundary x vertices."""
+        xs = self._sides_x.vertices
+        if self._wished_order == self._actual_order:
+            return xs
+        else:
+            return xs[::-1]
+
+    @property
+    def _y(self):
+        """Retrieve boundary y vertices."""
+        ys = self._sides_y.vertices
+        if self._wished_order == self._actual_order:
+            return ys
+        else:
+            return ys[::-1]
+
+    @property
+    def vertices(self):
+        """Return boundary vertices 2D array [x, y]."""
+        vertices = np.vstack((self._x, self._y)).T
+        vertices = vertices.astype(np.float64, copy=False)  # Important for spherical ops.
+        return vertices
+
+    def contour(self, closed=False):
+        """Return the (x, y) tuple of the boundary object.
+
+        If excludes the last element of each side because it's included in the next side.
+        If closed=False (the default), the last vertex is not equal to the first vertex
+        If closed=True, the last vertex is set to be equal to the first
+        closed=True is required for shapely Polygon creation.
+        closed=False is required for pyresample SPolygon creation.
+        """
+        x = self._x
+        y = self._y
+        if closed:
+            x = np.hstack((x, x[0]))
+            y = np.hstack((y, y[0]))
+        return x, y
+
+    def to_shapely_polygon(self):
+        """Define a Shapely Polygon."""
+        from shapely.geometry import Polygon
+        self = self.set_counterclockwise()
+        x, y = self.contour(closed=True)
+        return Polygon(zip(x, y))

pyresample/boundary/planar_boundary.py

@@ -0,0 +1,94 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# Copyright (c) 2014-2023 Pyresample developers
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+"""Define the PlanarBoundary class."""
+
+import logging
+
+import numpy as np
+
+from pyresample.boundary.base_boundary import BaseBoundary
+
+logger = logging.getLogger(__name__)
+
+
+def _is_projection_boundary_clockwise(sides_x, sides_y):
+    """Determine if the boundary is clockwise-defined in planar coordinates."""
+    from shapely.geometry import Polygon
+    x = np.concatenate([xs[:-1] for xs in sides_x])
+    y = np.concatenate([ys[:-1] for ys in sides_y])
+    x = np.hstack((x, x[0]))
+    y = np.hstack((y, y[0]))
+    polygon = Polygon(zip(x, y))
+    return not polygon.exterior.is_ccw
+
+
+class PlanarBoundary(BaseBoundary):
+    """Projection Boundary object.
+
+    The inputs must be the x and y sides of the projection.
+    It expects the projection coordinates to be planar (i.e. metric, radians).
+    """
+
+    @classmethod
+    def _check_is_boundary_clockwise(cls, sides_x, sides_y, area=None):
+        """SphericalBoundary specific implementation."""
+        return _is_projection_boundary_clockwise(sides_x=sides_x, sides_y=sides_y)
+
+    @classmethod
+    def _compute_boundary_sides(cls, area, vertices_per_side):
+        sides_x, sides_y = area._get_projection_sides(vertices_per_side=vertices_per_side)
+        return sides_x, sides_y
+
+    def __init__(self, area, vertices_per_side=None):
+        super().__init__(area=area, vertices_per_side=vertices_per_side)
+
+        self.sides_x = self._sides_x
+        self.sides_y = self._sides_y
+        self.crs = self._area.crs
+        self.cartopy_crs = self._area.to_cartopy_crs()
+
+    @property
+    def _point_inside(self):
+        x, y = self._area.get_proj_coords(data_slice=(int(self.shape[0] / 2), int(self.shape[1] / 2)))
+        return x.squeeze(), y.squeeze()
+
+    @property
+    def x(self):
+        """Retrieve boundary x vertices."""
+        return self._x
+
+    @property
+    def y(self):
+        """Retrieve boundary y vertices."""
+        return self._y
+
+    def plot(self, ax=None, subplot_kw=None, crs=None, alpha=0.6, **kwargs):
+        """Plot the the boundary. crs must be a Cartopy CRS !"""
+        from pyresample.visualization.geometries import plot_geometries
+
+        if self.cartopy_crs is None and crs is None:
+            raise ValueError("Projection Cartopy 'crs' is required to display projection boundary.")
+        if crs is None:
+            crs = self.cartopy_crs
+        if subplot_kw is None:
+            subplot_kw = {"projection": crs}
+
+        geom = self.to_shapely_polygon()
+        p = plot_geometries(geometries=[geom], crs=crs,
+                            ax=ax, subplot_kw=subplot_kw, alpha=alpha, **kwargs)
+        return p