Initial refactor of AreaBoundary

pyresample/boundary.py

@@ -54,62 +54,175 @@ def draw(self, mapper, options, **more_options):
         self.contour_poly.draw(mapper, options, **more_options)
 
 
+def _is_corner_is_clockwise(lon1, lat1, corner_lon, corner_lat, lon2, lat2):
+    """Determine if coordinates follow a clockwise path.
+
+    This uses :class:`pyresample.spherical.Arc` to determine the angle
+    between the first line segment (Arc) from (lon1, lat1) to
+    (corner_lon, corner_lat) and the second line segment from
+    (corner_lon, corner_lat) to (lon2, lat2). A straight line would
+    produce an angle of 0, a clockwise path would have a negative angle,
+    and a counter-clockwise path would have a positive angle.
+
+    """
+    import math
+
+    from pyresample.spherical import Arc, SCoordinate
+    point1 = SCoordinate(math.radians(lon1), math.radians(lat1))
+    point2 = SCoordinate(math.radians(corner_lon), math.radians(corner_lat))
+    point3 = SCoordinate(math.radians(lon2), math.radians(lat2))
+    arc1 = Arc(point1, point2)
+    arc2 = Arc(point2, point3)
+    angle = arc1.angle(arc2)
+    is_clockwise = -np.pi < angle < 0
+    return is_clockwise
+
+
+def _is_boundary_clockwise(sides_lons, sides_lats):
+    """Determine if the boundary sides are clockwise."""
+    is_clockwise = _is_corner_is_clockwise(
+        lon1=sides_lons[0][-2],
+        lat1=sides_lats[0][-2],
+        corner_lon=sides_lons[0][-1],
+        corner_lat=sides_lats[0][-1],
+        lon2=sides_lons[1][1],
+        lat2=sides_lats[1][1])
+    return is_clockwise
+
+
+def _check_sides_list(sides):
+    if not isinstance(sides, list):
+        raise TypeError("Boundary sides must be a list")
+    if len(sides) != 4:
+        raise ValueError("Boundary sides list must be a list with 4 elements.")
+    # TODO:
+    # - Numpy array elements of at least length 2
+
+
 class AreaBoundary(Boundary):
     """Area boundary objects.
 
     The inputs must be a (lon_coords, lat_coords) tuple for each of the 4 sides.
     """
 
-    def __init__(self, *sides):
-        Boundary.__init__(self)
-        # Check 4 sides are provided
-        if len(sides) != 4:
-            raise ValueError("AreaBoundary expects 4 sides.")
-        # Retrieve sides
-        self.sides_lons, self.sides_lats = zip(*sides)
-        self.sides_lons = list(self.sides_lons)
-        self.sides_lats = list(self.sides_lats)
+    def __init__(self, lon_sides, lat_sides, wished_order=None):
+        _check_sides_list(lon_sides)
+        _check_sides_list(lat_sides)
 
-    @classmethod
-    def from_lonlat_sides(cls, lon_sides, lat_sides):
-        """Define AreaBoundary from list of lon_sides and lat_sides.
+        # Old interface
+        self._contour_poly = None
+        self.sides_lons = lon_sides
+        self.sides_lats = lat_sides
 
-        For an area of shape (m, n), the sides must adhere the format:
+        # New interface
+        # TODO: self.sides (BoundarySide(s))
 
-        sides = [np.array([v00, v01, ..., v0n]),
-                 np.array([v0n, v1n, ..., vmn]),
-                 np.array([vmn, ..., vm1, vm0]),
-                 np.array([vm0, ... ,v10, v00])]
-        """
-        boundary = cls(*zip(lon_sides, lat_sides))
-        return boundary
+        # Check if it is clockwise/counterclockwise
+        self.is_clockwise = _is_boundary_clockwise(sides_lons=lon_sides,
+                                                   sides_lats=lat_sides)
+        self.is_counterclockwise = not self.is_clockwise
+
+        # Define wished order
+        if self.is_clockwise:
+            self._actual_order = "clockwise"
+        else:
+            self._actual_order = "counterclockwise"
+
+        if wished_order is None:
+            self._wished_order = self._actual_order
+        else:
+            if wished_order not in ["clockwise", "counterclockwise"]:
+                raise ValueError("Valid order is 'clockwise' or 'counterclockwise'")
+            self._wished_order = wished_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 lons(self):
+        """Retrieve boundary longitude vertices."""
+        lons = np.concatenate([lns[:-1] for lns in self.sides_lons])
+        if self._wished_order == self._actual_order:
+            return lons
+        else:
+            return lons[::-1]
+
+    @property
+    def lats(self):
+        """Retrieve boundary latitude vertices."""
+        lats = np.concatenate([lts[:-1] for lts in self.sides_lats])
+        if self._wished_order == self._actual_order:
+            return lats
+        else:
+            return lats[::-1]
+
+    @property
+    def vertices(self):
+        """Return boundary vertices 2D array [lon, lat]."""
+        vertices = np.vstack((self.lons, self.lats)).T
+        vertices = vertices.astype(np.float64, copy=False)  # Important for spherical ops.
+        return vertices
 
     def contour(self, closed=False):
-        """Get the (lons, lats) tuple of the boundary object.
+        """Return the (lons, lats) 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.
         """
-        lons = np.concatenate([lns[:-1] for lns in self.sides_lons])
-        lats = np.concatenate([lts[:-1] for lts in self.sides_lats])
+        lons = self.lons
+        lats = self.lats
         if closed:
             lons = np.hstack((lons, lons[0]))
             lats = np.hstack((lats, lats[0]))
         return lons, lats
 
-    @property
-    def vertices(self):
-        """Return boundary polygon vertices."""
-        lons, lats = self.contour()
-        vertices = np.vstack((lons, lats)).T
-        vertices = vertices.astype(np.float64, copy=False)  # Important for spherical ops.
-        return vertices
+    def _to_shapely_polygon(self):
+        from shapely.geometry import Polygon
+        self = self.set_counterclockwise()  # TODO: add exception for pole wrapping polygons
+        lons, lats = self.contour(closed=True)
+        return Polygon(zip(lons, lats))
+
+    def _to_spherical_polygon(self):
+        self = self.set_clockwise()  # TODO: add exception for pole wrapping polygons
+        raise NotImplementedError("This will return a SPolygon in pyresample 2.0")
+
+    def polygon(self, shapely=False):
+        """Return the boundary polygon."""
+        if shapely:
+            return self._to_shapely_polygon()
+        else:
+            return self._to_spherical_polygon()
+
+    # For backward compatibility !
+    @classmethod
+    def from_lonlat_sides(cls, lon_sides, lat_sides):
+        """Define AreaBoundary from list of lon_sides and lat_sides.
+
+        For an area of shape (m, n), the sides must adhere the format:
+
+        sides = [np.array([v00, v01, ..., v0n]),
+                 np.array([v0n, v1n, ..., vmn]),
+                 np.array([vmn, ..., vm1, vm0]),
+                 np.array([vm0, ... ,v10, v00])]
+        """
+        warnings.warn("Use `AreaBoundary(lon_sides, lat_sides)` instead of `from_lonlat_sides`",
+                      PendingDeprecationWarning, stacklevel=2)
+        boundary = cls(lon_sides=lon_sides, lat_sides=lat_sides)
+        return boundary
 
     def decimate(self, ratio):
         """Remove some points in the boundaries, but never the corners."""
+        # TODO: to update --> used by AreaDefBoundary
         for i in range(len(self.sides_lons)):
             length = len(self.sides_lons[i])
             start = int((length % ratio) / 2)
@@ -122,33 +235,27 @@ def decimate(self, ratio):
             self.sides_lons[i] = self.sides_lons[i][points]
             self.sides_lats[i] = self.sides_lats[i][points]
 
-    def _to_shapely_polygon(self):
-        from shapely.geometry import Polygon
-        lons, lats = self.contour(closed=True)
-        return Polygon(zip(lons, lats))
-
-    def _to_spherical_polygon(self):
-        raise NotImplementedError("This will return a SPolygon in pyresample 2.0")
+    @property
+    def contour_poly(self):
+        """Return the pyresample SphPolygon."""
+        if self._contour_poly is None:
+            self._contour_poly = SphPolygon(np.deg2rad(self.vertices))
+        return self._contour_poly
 
-    def polygon(self, shapely=False):
-        """Return the boundary polygon."""
-        if shapely:
-            return self._to_shapely_polygon()
-        else:
-            return self._to_spherical_polygon()
+    def draw(self, mapper, options, **more_options):
+        """Draw the current boundary on the *mapper*."""
+        self.contour_poly.draw(mapper, options, **more_options)
 
 
 class AreaDefBoundary(AreaBoundary):
     """Boundaries for area definitions (pyresample)."""
 
     def __init__(self, area, frequency=1):
-        lons, lats = area.get_bbox_lonlats()
+        lon_sides, lat_sides = area.get_bbox_lonlats()
         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, lon_sides=lon_sides, lat_sides=lat_sides)
         if frequency != 1:
             self.decimate(frequency)
 

pyresample/geometry.py

@@ -534,9 +534,15 @@ def boundary(self, *, vertices_per_side=None, force_clockwise=False, frequency=N
             warnings.warn("The `frequency` argument is pending deprecation, use `vertices_per_side` instead",
                           PendingDeprecationWarning, stacklevel=2)
         vertices_per_side = vertices_per_side or frequency
-        lon_sides, lat_sides = self.get_bbox_lonlats(vertices_per_side=vertices_per_side,
-                                                     force_clockwise=force_clockwise)
-        return AreaBoundary.from_lonlat_sides(lon_sides, lat_sides)
+        lon_sides, lat_sides = self._get_boundary_sides(coordinates="geographic",
+                                                        vertices_per_side=vertices_per_side)
+        # TODO: this could be changed but it would breaks backward compatibility
+        # TODO: Implement code to return projection boundary !
+        if force_clockwise:
+            wished_order = "clockwise"
+        else:
+            wished_order = None
+        return AreaBoundary(lon_sides, lat_sides, wished_order=wished_order)
 
     def get_cartesian_coords(self, nprocs=None, data_slice=None, cache=False):
         """Retrieve cartesian coordinates of geometry definition.

pyresample/test/test_boundary.py

@@ -38,6 +38,7 @@ def test_creation_from_lonlat_sides(self):
                      np.array([7.0, 8.0]),
                      np.array([8.0, 8.5, 9.0]),
                      np.array([9.0, 6.0])]
+
         # Define AreaBoundary
         boundary = AreaBoundary.from_lonlat_sides(lon_sides, lat_sides)
 
@@ -50,15 +51,17 @@ def test_creation_from_lonlat_sides(self):
 
     def test_creation(self):
         """Test AreaBoundary creation."""
-        list_sides = [(np.array([1., 1.5, 2.]), np.array([6., 6.5, 7.])),
-                      (np.array([2., 3.]), np.array([7., 8.])),
-                      (np.array([3., 3.5, 4.]), np.array([8., 8.5, 9.])),
-                      (np.array([4., 1.]), np.array([9., 6.]))]
-        lon_sides = [side[0]for side in list_sides]
-        lat_sides = [side[1]for side in list_sides]
+        lon_sides = [np.array([1.0, 1.5, 2.0]),
+                     np.array([2.0, 3.0]),
+                     np.array([3.0, 3.5, 4.0]),
+                     np.array([4.0, 1.0])]
+        lat_sides = [np.array([6.0, 6.5, 7.0]),
+                     np.array([7.0, 8.0]),
+                     np.array([8.0, 8.5, 9.0]),
+                     np.array([9.0, 6.0])]
 
         # Define AreaBoundary
-        boundary = AreaBoundary(*list_sides)
+        boundary = AreaBoundary(lon_sides, lat_sides)
 
         # Assert sides coincides
         for b_lon, src_lon in zip(boundary.sides_lons, lon_sides):
@@ -69,12 +72,14 @@ def test_creation(self):
 
     def test_number_sides_required(self):
         """Test AreaBoundary requires 4 sides ."""
-        list_sides = [(np.array([1., 1.5, 2.]), np.array([6., 6.5, 7.])),
-                      (np.array([2., 3.]), np.array([7., 8.])),
-                      (np.array([3., 3.5, 4.]), np.array([8., 8.5, 9.])),
-                      (np.array([4., 1.]), np.array([9., 6.]))]
+        lon_sides = [np.array([1.0, 1.5, 2.0]),
+                     np.array([2.0, 3.0]),
+                     np.array([4.0, 1.0])]
+        lat_sides = [np.array([6.0, 6.5, 7.0]),
+                     np.array([7.0, 8.0]),
+                     np.array([9.0, 6.0])]
         with pytest.raises(ValueError):
-            AreaBoundary(*list_sides[0:3])
+            AreaBoundary(lon_sides, lat_sides)
 
     def test_vertices_property(self):
         """Test AreaBoundary vertices property."""
@@ -87,7 +92,7 @@ def test_vertices_property(self):
                      np.array([8.0, 8.5, 9.0]),
                      np.array([9.0, 6.0])]
         # Define AreaBoundary
-        boundary = AreaBoundary.from_lonlat_sides(lon_sides, lat_sides)
+        boundary = AreaBoundary(lon_sides, lat_sides)
 
         # Assert vertices
         expected_vertices = np.array([[1., 6.],
@@ -100,22 +105,32 @@ def test_vertices_property(self):
 
     def test_contour(self):
         """Test that AreaBoundary.contour(closed=False) returns the correct (lon,lat) tuple."""
-        list_sides = [(np.array([1., 1.5, 2.]), np.array([6., 6.5, 7.])),
-                      (np.array([2., 3.]), np.array([7., 8.])),
-                      (np.array([3., 3.5, 4.]), np.array([8., 8.5, 9.])),
-                      (np.array([4., 1.]), np.array([9., 6.]))]
-        boundary = AreaBoundary(*list_sides)
+        lon_sides = [np.array([1.0, 1.5, 2.0]),
+                     np.array([2.0, 3.0]),
+                     np.array([3.0, 3.5, 4.0]),
+                     np.array([4.0, 1.0])]
+        lat_sides = [np.array([6.0, 6.5, 7.0]),
+                     np.array([7.0, 8.0]),
+                     np.array([8.0, 8.5, 9.0]),
+                     np.array([9.0, 6.0])]
+        # Define AreaBoundary
+        boundary = AreaBoundary(lon_sides, lat_sides)
         lons, lats = boundary.contour()
         assert np.allclose(lons, np.array([1., 1.5, 2., 3., 3.5, 4.]))
         assert np.allclose(lats, np.array([6., 6.5, 7., 8., 8.5, 9.]))
 
     def test_contour_closed(self):
         """Test that AreaBoundary.contour(closed=True) returns the correct (lon,lat) tuple."""
-        list_sides = [(np.array([1., 1.5, 2.]), np.array([6., 6.5, 7.])),
-                      (np.array([2., 3.]), np.array([7., 8.])),
-                      (np.array([3., 3.5, 4.]), np.array([8., 8.5, 9.])),
-                      (np.array([4., 1.]), np.array([9., 6.]))]
-        boundary = AreaBoundary(*list_sides)
+        lon_sides = [np.array([1.0, 1.5, 2.0]),
+                     np.array([2.0, 3.0]),
+                     np.array([3.0, 3.5, 4.0]),
+                     np.array([4.0, 1.0])]
+        lat_sides = [np.array([6.0, 6.5, 7.0]),
+                     np.array([7.0, 8.0]),
+                     np.array([8.0, 8.5, 9.0]),
+                     np.array([9.0, 6.0])]
+        # Define AreaBoundary
+        boundary = AreaBoundary(lon_sides, lat_sides)
         lons, lats = boundary.contour(closed=True)
         assert np.allclose(lons, np.array([1., 1.5, 2., 3., 3.5, 4., 1.]))
         assert np.allclose(lats, np.array([6., 6.5, 7., 8., 8.5, 9., 6.]))