Fix boundary ordering logic for swath and projections !

docs/source/howtos/spherical_geometry.rst

@@ -74,11 +74,11 @@ satellite passes. See trollschedule_ how to generate a list of satellite overpas
  >>> from pyresample.spherical_utils import GetNonOverlapUnions
 
  >>> area_boundary = area_def.geographic_boundary(vertices_per_side=100) # doctest: +SKIP
- >>> area_boundary = area_boundary.contour_poly # 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.geographic_boundary().polygon) # doctest: +SKIP
 
  >>> non_overlaps = GetNonOverlapUnions(list_of_polygons) # doctest: +SKIP
  >>> non_overlaps.merge() # doctest: +SKIP

pyresample/_config.py

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

pyresample/boundary/area_boundary.py

@@ -43,6 +43,9 @@ def contour(self):
     @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))

pyresample/boundary/base_boundary.py

@@ -30,18 +30,25 @@ class BaseBoundary:
     """Base class for boundary objects."""
     __slots__ = ["_sides_x", "_sides_y"]
 
-    def __init__(self, sides_x, sides_y, order=None):
+    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)
+        self.is_clockwise = self._check_is_boundary_clockwise(sides_x, sides_y, area)
         self.is_counterclockwise = not self.is_clockwise
-        self._set_order(order)
+        self._set_order(order=None)  # FIX !
 
-    def _check_is_boundary_clockwise(self, sides_x, sides_y):
+    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:
@@ -112,9 +119,9 @@ def contour(self, closed=False):
             y = np.hstack((y, y[0]))
         return x, y
 
-    def _to_shapely_polygon(self):
+    def to_shapely_polygon(self):
         """Define a Shapely Polygon."""
         from shapely.geometry import Polygon
-        self = self.set_counterclockwise()  # FIXME: add exception for pole wrapping polygons
+        self = self.set_counterclockwise()
         x, y = self.contour(closed=True)
         return Polygon(zip(x, y))

pyresample/boundary/geographic_boundary.py

@@ -24,44 +24,86 @@
 
 from pyresample.boundary.area_boundary import Boundary as OldBoundary
 from pyresample.boundary.base_boundary import BaseBoundary
+from pyresample.spherical import SphPolygon
 
 logger = logging.getLogger(__name__)
 
 
-def _is_corner_is_clockwise(lon1, lat1, corner_lon, corner_lat, lon2, lat2):
-    """Determine if coordinates follow a clockwise path.
+def _get_swath_local_inside_point_from_sides(sides_x, sides_y, start_idx=0):
+    """Retrieve point inside boundary close to the point used to determine order.
+
+    This is required for global areas (spanning more than 180 degrees) and swaths.
+    The start_idx indicates the opposites sides corners (start_idx, start_idx+2) from
+    which to try identify a point inside the polygon.
+    """
+    # NOTE: the name of the object here refers to start_idx=0
+    from pyresample.spherical import Arc, SCoordinate
+    idx1 = start_idx
+    idx2 = (start_idx + 2) % 4
+
+    top_corner = sides_x[idx1][0], sides_y[idx1][0]
+    top_right_point = sides_x[idx1][1], sides_y[idx1][1]
+    bottom_corner = sides_x[idx2][-1], sides_y[idx2][-1]
+    bottom_right_point = sides_x[idx2][-2], sides_y[idx2][-2]
+    point_top_corner = SCoordinate(*np.deg2rad(top_corner))
+    point_top_right_point = SCoordinate(*np.deg2rad(top_right_point))
+    point_bottom_corner = SCoordinate(*np.deg2rad(bottom_corner))
+    point_bottom_right_point = SCoordinate(*np.deg2rad(bottom_right_point))
+    arc1 = Arc(point_top_corner, point_bottom_right_point)
+    arc2 = Arc(point_top_right_point, point_bottom_corner)
+    point_inside = arc1.intersection(arc2)
+    if point_inside is not None:
+        point_inside = point_inside.vertices_in_degrees[0]
+    return point_inside
+
+
+def _try_get_local_inside_point(sides_x, sides_y):
+    """Try to get a local inside point from one of the 4 boundary sides corners."""
+    for start_idx in range(0, 4):
+        point_inside = _get_swath_local_inside_point_from_sides(sides_x, sides_y, start_idx=start_idx)
+        if point_inside is not None:
+            return point_inside, start_idx
+    else:
+        return None, None
+
+
+def _is_clockwise_order(first_point, second_point, point_inside):
+    """Determine if polygon 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.
+    between a polygon arc segment and a point known to be inside the polygon.
 
+    Note: pyresample.spherical assumes angles are positive if counterclockwise.
+    Note: if the longitude distance between the first_point/second_point and point_inside is
+    larger than 180°, the function likely return a wrong unexpected result !
     """
-    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)
+    point1 = SCoordinate(*np.deg2rad(first_point))
+    point2 = SCoordinate(*np.deg2rad(second_point))
+    point3 = SCoordinate(*np.deg2rad(point_inside))
+    arc12 = Arc(point1, point2)
+    arc23 = Arc(point2, point3)
+    angle = arc12.angle(arc23)
     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_is_clockwise(area, sides_x, sides_y):
+    from pyresample import SwathDefinition
+
+    if isinstance(area, SwathDefinition):
+        point_inside, start_idx = _try_get_local_inside_point(sides_x, sides_y)
+        first_point = sides_x[start_idx][0], sides_y[start_idx][0]
+        second_point = sides_x[start_idx][1], sides_y[start_idx][1]
+        return _is_clockwise_order(first_point, second_point, point_inside)
+    else:
+        if area.is_geostationary:
+            point_inside = area.get_lonlat(row=int(area.shape[0] / 2), col=int(area.shape[1] / 2))
+            first_point = sides_x[0][0], sides_y[0][0]
+            second_point = sides_x[0][1], sides_y[0][1]
+            return _is_clockwise_order(first_point, second_point, point_inside)
+        else:
+            return True
 
 
 class GeographicBoundary(BaseBoundary, OldBoundary):
@@ -74,17 +116,35 @@ class GeographicBoundary(BaseBoundary, OldBoundary):
     #   from the old interface for compatibility to AreaBoundary
 
     @classmethod
-    def _check_is_boundary_clockwise(cls, sides_x, sides_y):
+    def _check_is_boundary_clockwise(cls, sides_x, sides_y, area):
         """GeographicBoundary specific implementation."""
-        return _is_boundary_clockwise(sides_lons=sides_x, sides_lats=sides_y)
+        return _check_is_clockwise(area, sides_x, sides_y)
+
+    @classmethod
+    def _compute_boundary_sides(cls, area, vertices_per_side):
+        sides_lons, sides_lats = area._get_geographic_sides(vertices_per_side=vertices_per_side)
+        return sides_lons, sides_lats
 
-    def __init__(self, sides_lons, sides_lats, order=None, crs=None):
-        super().__init__(sides_x=sides_lons, sides_y=sides_lats, order=order)
+    def __init__(self, area, vertices_per_side=None):
+        super().__init__(area=area, vertices_per_side=vertices_per_side)
 
         self.sides_lons = self._sides_x
         self.sides_lats = self._sides_y
-        self.crs = crs or CRS(proj="longlat", ellps="WGS84")
-        self._contour_poly = None  # Backcompatibility with old AreaBoundary
+
+        # Define CRS
+        if self.is_swath:
+            crs = self._area.crs
+        else:
+            crs = CRS(proj="longlat", ellps="WGS84")  # FIXME: AreaDefinition.get_lonlat for geographic projections?
+        self.crs = crs
+
+        # Backcompatibility with old AreaBoundary
+        self._contour_poly = None
+
+    @property
+    def is_swath(self):
+        """Determine if is the boundary of a swath."""
+        return self._area.__class__.__name__ == "SwathDefinition"
 
     @property
     def lons(self):
@@ -96,28 +156,22 @@ def lats(self):
         """Retrieve boundary latitude vertices."""
         return self._y
 
-    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."""
-        # ALTERNATIVE:
-        # - shapely: to_shapely_polygon(), to_shapely_line(),
-        # - pyresample spherical: to_polygon(), to_line(), polygon, line
-        if shapely:
-            return self._to_shapely_polygon()
-        else:
-            return self._to_spherical_polygon()
-
-    def plot(self, ax=None, subplot_kw=None, **kwargs):
+    @property
+    def polygon(self):
+        """Return the boundary spherical polygon."""
+        self = self.set_clockwise()
+        if self._contour_poly is None:
+            self._contour_poly = SphPolygon(np.deg2rad(self.vertices))
+        return self._contour_poly
+
+    def plot(self, ax=None, subplot_kw=None, alpha=0.6, **kwargs):
         """Plot the the boundary."""
         import cartopy.crs as ccrs
 
         from pyresample.visualization.geometries import plot_geometries
 
-        geom = self.polygon(shapely=True)
+        geom = self.to_shapely_polygon()
         crs = ccrs.Geodetic()
         p = plot_geometries(geometries=[geom], crs=crs,
-                            ax=ax, subplot_kw=subplot_kw, **kwargs)
+                            ax=ax, subplot_kw=subplot_kw, alpha=alpha, **kwargs)
         return p

pyresample/boundary/projection_boundary.py

@@ -45,17 +45,27 @@ class ProjectionBoundary(BaseBoundary):
     """
 
     @classmethod
-    def _check_is_boundary_clockwise(cls, sides_x, sides_y):
+    def _check_is_boundary_clockwise(cls, sides_x, sides_y, area=None):
         """GeographicBoundary specific implementation."""
         return _is_projection_boundary_clockwise(sides_x=sides_x, sides_y=sides_y)
 
-    def __init__(self, sides_x, sides_y, crs, order=None, cartopy_crs=None):
-        super().__init__(sides_x=sides_x, sides_y=sides_y, order=order)
+    @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 = crs
-        self.cartopy_crs = cartopy_crs
+        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):
@@ -67,23 +77,18 @@ def y(self):
         """Retrieve boundary y vertices."""
         return self._y
 
-    def polygon(self, shapely=True):
-        """Return the boundary polygon."""
-        if shapely:
-            return self._to_shapely_polygon()
-        else:
-            raise NotImplementedError("Only shapely polygon available.")
-
-    def plot(self, ax=None, subplot_kw=None, crs=None, **kwargs):
+    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.polygon(shapely=True)
+        geom = self.to_shapely_polygon()
         p = plot_geometries(geometries=[geom], crs=crs,
-                            ax=ax, subplot_kw=subplot_kw, **kwargs)
+                            ax=ax, subplot_kw=subplot_kw, alpha=alpha, **kwargs)
         return p