Bugfix polygon vertices ordering

Fixes an issue where the edge vertices ordering can be incorrectly flipped.

Change also removes unnecessary arguments from `discretize_all_edge_boundaries`, and propagates this change to methods that call this function.

Author: sedonaprice

regions/_utils/spherical_helpers.py

@@ -316,22 +316,38 @@ def _add_tan_pts_if_in_pa_range(
     return coord_list
 
 
-def _validate_vertices_ordering(verts, gc, centroid):
-    pas_verts = gc.center.position_angle(verts).to(u.deg)
+def _check_edge_lt_pi(pas_verts, wrap_ang):
+    pas_verts_wrap = pas_verts.wrap_at(wrap_ang)
+
+    is_valid_arc_length = ((pas_verts_wrap[1] - pas_verts_wrap[0]).to(u.deg) <= 180 * u.deg)
+
+    return pas_verts_wrap, is_valid_arc_length
+
+
+def _validate_vertices_ordering(verts, gc, gc_center=None):
+    if gc is not None:
+        gc_center = gc.center
+    pas_verts = gc_center.position_angle(verts).to(u.deg)
 
-    # Check ordering of vertices pas:
-    pa_centroid = gc.center.position_angle(centroid).to(u.deg)
     wrap_ang = pas_verts[0]
-    pas_verts_wrap = pas_verts.wrap_at(wrap_ang)
-    pa_centroid_wrap = pa_centroid.wrap_at(wrap_ang)
-    in_range_lon = (pa_centroid_wrap >= pas_verts_wrap[0]) & (
-        pa_centroid_wrap <= pas_verts_wrap[1]
-    )
-    if not in_range_lon:
-        # If not in range, invert the vertices PA ordering:
-        pas_verts = pas_verts[::-1]
-        wrap_ang = pas_verts[0]
-        pas_verts_wrap = pas_verts.wrap_at(wrap_ang)
+
+    # Check ordering of vertices pas:
+
+    # Principle: ALL EDGES must be <= 180 deg of length
+    # So difference of pa[1] - pa[0] <= 180 deg
+    # If not, swap the order.
+
+    pas_verts_wrap, is_valid_arc_length = _check_edge_lt_pi(pas_verts, wrap_ang)
+
+    if not is_valid_arc_length:
+        wrap_ang_opp = pas_verts[-1]
+        pas_verts_wrap_opp, is_valid_arc_length_opp = _check_edge_lt_pi(
+            pas_verts[::-1], wrap_ang_opp
+        )
+        if is_valid_arc_length_opp:
+            return pas_verts_wrap_opp, wrap_ang_opp
+
+        raise ValueError('Invalid arc')  # should never occur
 
     return pas_verts_wrap, wrap_ang
 
@@ -400,7 +416,7 @@ def get_edge_raw_lonlat_bounds_circ_edges(vertices, centroid, gcs):
         # PAs from gc center to vertices:
         verts = SkyCoord(np.concatenate([[vertices[i - 1]], [vertices[i]]]))
 
-        pas_verts_wrap, wrap_ang = _validate_vertices_ordering(verts, gc, centroid)
+        pas_verts_wrap, wrap_ang = _validate_vertices_ordering(verts, gc)
 
         # --------------------------------------------------------
         # Latitude tangent points from bound circle as len 2 SkyCoord:
@@ -409,7 +425,8 @@ def get_edge_raw_lonlat_bounds_circ_edges(vertices, centroid, gcs):
         tan_lat_pts = _get_circle_latitude_tangent_points(gc.center, gc.radius)
 
         lats_list = _add_tan_pts_if_in_pa_range(
-            lats_list, tan_lat_pts, gc, wrap_ang, pas_verts_wrap, coord='lat'
+            lats_list, tan_lat_pts, gc, wrap_ang, pas_verts_wrap,
+            coord='lat'
         )
 
         # --------------------------------------------------------
@@ -419,7 +436,8 @@ def get_edge_raw_lonlat_bounds_circ_edges(vertices, centroid, gcs):
         tan_lon_pts = _get_circle_longitude_tangent_points(gc.center, gc.radius)
         if tan_lon_pts is not None:
             lons_list = _add_tan_pts_if_in_pa_range(
-                lons_list, tan_lon_pts, gc, wrap_ang, pas_verts_wrap, coord='lon'
+                lons_list, tan_lon_pts, gc, wrap_ang, pas_verts_wrap,
+                coord='lon'
             )
 
     lons_arr = [lons_list.min(), lons_list.max()]
@@ -432,7 +450,90 @@ def get_edge_raw_lonlat_bounds_circ_edges(vertices, centroid, gcs):
     return Longitude(lons_arr).to(u.deg), Latitude(lats_arr).to(u.deg)
 
 
-def _discretize_edge_boundary(vertices, circ, centroid, n_points):
+def get_line_edge_raw_lonlat_bounds_circ_edges(coords, center):
+    """
+    Get the raw longitude / latitude bounds from the circle edges of
+    spherical sky region.
+
+    Parameters
+    ----------
+    coords : `~astropy.coordinates.SkyCoord`
+        The start, end of the line as a SkyCoord.
+
+    center : `~astropy.coordinates.SkyCoord`
+        The line center as a SkyCoord.
+
+    Returns
+    -------
+    longitude_limits, latitude_limits: `~astropy.coordinates.Longitude`
+        Length two |Longitude| and |Latitude| with the computed
+        longitude/latitude bounds from the polygon edges.
+    """
+    # Consider lon/lat of vertices: may produce min/max bounds:
+    vrepr = coords.represent_as('spherical')
+
+    # Special handling:
+    # Exclude vertices from longitude bounds if any is on a pole
+    lons_list = []
+    lats_list = []
+    for v in vrepr:
+        if np.abs(v.lat.to(u.deg).deg) < 90:
+            lons_list.append(v.lon)
+            lats_list.append(v.lat)
+    lons_list = Longitude(lons_list, unit=u.radian)
+    lats_list = Latitude(lats_list, unit=u.radian)
+
+    # Need to also check for "out bulging" from edges,
+    # as far as latitude/lon bounds:
+    # eg, 2 vertices at ~60deg: the gc arc goes ~closer to the pole;
+    # a circle centered close to the pole but not extending over it:
+    # WIDE lon bounds
+
+    gc_center = cross_product_skycoord2skycoord(coords[0], coords[1])
+    gc_radius = 90 * u.deg
+
+    pas_verts_wrap, wrap_ang = _validate_vertices_ordering(
+        coords, None, gc_center=gc_center
+    )
+
+    # --------------------------------------------------------
+    # Latitude tangent points from bound circle as len 2 SkyCoord:
+    # Only add to the list if the tangent point is located along this edge
+    tan_lat_pts = _get_circle_latitude_tangent_points(gc_center, gc_radius)
+
+    lats_list = _add_tan_pts_if_in_pa_range(
+        lats_list, tan_lat_pts, None, wrap_ang, pas_verts_wrap, coord='lat',
+        gc_center=gc_center
+    )
+
+    # --------------------------------------------------------
+    # Longitude tangent points from bound circle as len 2 SkyCoord:
+    # Only add to the list if the tangent point is located along this edge
+
+    tan_lon_pts = _get_circle_longitude_tangent_points(gc_center, gc_radius)
+    if tan_lon_pts is not None:
+        lons_list = _add_tan_pts_if_in_pa_range(
+            lons_list, tan_lon_pts, None, wrap_ang, pas_verts_wrap, coord='lon',
+            gc_center=gc_center
+        )
+
+    lons_arr = [lons_list.min(), lons_list.max()]
+    lats_arr = [lats_list.min(), lats_list.max()]
+
+    # --------------------------------------------------------
+    # Invert longitude order if centroid is outside of range:
+    lons_arr = _validate_lon_bounds_ordering(lons_arr, center)
+
+    return Longitude(lons_arr).to(u.deg), Latitude(lats_arr).to(u.deg)
+
+
+def _discretize_edge_boundary(vertices, circ, n_points,
+                              circ_center=None, circ_radius=None):
+
+    if circ is not None:
+        circ_center = circ.center
+        circ_radius = circ.radius
+
     # Discretize an edge boundary defined by a circle, geodetic or not:
     # either great circle arc, or a span of a non-great circle
     # (e.g., constant lat edges of RangeSphericalSkyRegion)
@@ -441,7 +542,9 @@ def _discretize_edge_boundary(vertices, circ, centroid, n_points):
     # connecting circle center to the two vertices bounding that edge:
     pas_verts = circ.center.position_angle(vertices).to(u.deg)
 
-    pas_verts_wrap, wrap_ang = _validate_vertices_ordering(vertices, circ, centroid)
+    pas_verts_wrap, wrap_ang = _validate_vertices_ordering(
+        vertices, circ, gc_center=circ_center,
+    )
 
     # Need to wrap angles to calculate span: wrap at lower value:
     pas_verts_wrap = pas_verts.wrap_at(wrap_ang)
@@ -453,12 +556,47 @@ def _discretize_edge_boundary(vertices, circ, centroid, n_points):
 
     # Calculate directional offsets to get boundary discretization,
     # with vertices as SkyCoords
-    bound_verts = circ.center.directional_offset_by(theta, circ.radius)
+    bound_verts = circ.center.directional_offset_by(theta, circ_radius)
 
     return bound_verts
 
 
-def discretize_all_edge_boundaries(vertices, circs, centroid, n_points):
+def discretize_line_boundary(coords, n_points):
+    """
+    Discretize all edge boundaries for spherical sky regions.
+
+    Parameters
+    ----------
+    coords : `~astropy.coordinates.SkyCoord`
+        The start, end of the line as a SkyCoord.
+
+    n_points : int
+        The number of points for discretization along each edge.
+
+    Returns
+    -------
+    edge_bound_verts : `~astropy.coordinates.SkyCoord`
+        The discretized boundary edge vertices.
+    """
+    # Iterate over full set of vertices & boundary circles for
+    # a region (polygon or range)
+
+    gc_center = cross_product_skycoord2skycoord(coords[0], coords[1])
+    gc_radius = 90 * u.deg
+
+    bound_verts = _discretize_edge_boundary(coords, None, n_points,
+                                            circ_center=gc_center, circ_radius=gc_radius)
+
+    return SkyCoord(
+        SkyCoord(
+            bound_verts,
+            representation_type=UnitSphericalRepresentation,
+        ),
+        representation_type=SphericalRepresentation,
+    )
+
+
+def discretize_all_edge_boundaries(vertices, circs, n_points):
     """
     Discretize all edge boundaries for spherical sky regions.
 
@@ -470,9 +608,6 @@ def discretize_all_edge_boundaries(vertices, circs, centroid, n_points):
     circs : `~regions.CircleSphericalSkyRegion`
         The circle boundaries as CircleSphericalSkyRegion instances.
 
-    centroid : `~astropy.coordinates.SkyCoord`
-        The region centroid as a SkyCoord.
-
     n_points : int
         The number of points for discretization along each edge.
 
@@ -489,7 +624,7 @@ def discretize_all_edge_boundaries(vertices, circs, centroid, n_points):
         # PAs from gc center to vertices:
         verts = SkyCoord(np.concatenate([[vertices[i - 1]], [vertices[i]]]))
 
-        bound_verts = _discretize_edge_boundary(verts, circ, centroid, n_points)
+        bound_verts = _discretize_edge_boundary(verts, circ, n_points)
 
         if all_edge_bound_verts is None:
             all_edge_bound_verts = bound_verts

regions/shapes/lune.py

@@ -184,7 +184,7 @@ def transform_to(self, frame, merge_attributes=True):
 
     def discretize_boundary(self, n_points=100):
         bound_verts = discretize_all_edge_boundaries(
-            self.vertices, self._edge_circs, self.centroid, n_points
+            self.vertices, self._edge_circs, n_points
         )
         return PolygonSphericalSkyRegion(bound_verts)
 

regions/shapes/polygon.py

@@ -597,7 +597,7 @@ def transform_to(self, frame, merge_attributes=True):
 
     def discretize_boundary(self, n_points=10):
         bound_verts = discretize_all_edge_boundaries(
-            self.vertices, self._edge_circs, self.centroid, n_points
+            self.vertices, self._edge_circs, n_points
         )
         return PolygonSphericalSkyRegion(bound_verts)
 

regions/shapes/range.py

@@ -825,7 +825,7 @@ def discretize_boundary(self, n_points=10):
             return self.longitude_bounds.discretize_boundary(n_points=n_points)
         elif (bound_nverts == 4) | (bound_nverts == 3):
             bound_verts = discretize_all_edge_boundaries(
-                self.vertices, self._edge_circs, self.centroid, n_points
+                self.vertices, self._edge_circs, n_points
             )
             return PolygonSphericalSkyRegion(bound_verts)
         elif bound_nverts == 6: