Determine if sky polygon winding direction is clockwise or counter clockwise

poly_winding.py

# %matplotlib inline

from __future__ import print_function, division

import logging
import numpy as np
import matplotlib.pyplot as plt


def polygon_to_points(polygon):
    if not polygon.startswith("POLYGON ICRS"):
        logging.error("Polygon string is not in the expected format of 'POLYGON ICRS ra dec [ra dec]+'")
        return []
    
    coords = polygon.split(' ')
    coords = coords[2:]
    points = np.array(coords).reshape(int(len(coords)/2),2).astype(float)
    #print(points)
    return points

    
def verify_polygon(polygon):
    points = polygon_to_points(polygon)
    
    total = 0
    for curr_idx in range(len(points)):
        next_idx = curr_idx+1 if curr_idx < len(points)-1 else 0
        total += (points[next_idx,0] - points[curr_idx,0])*(points[next_idx,1] + points[curr_idx,1])
    
    if (total < 0):
        logging.error('Polygon winding is clockwise, where counter clockwise winding is required for ' + polygon + ' ')
    else:
        print('Polygon winding is counter clockwise for ' + polygon )
    
    return total < 0
  
  
def plot_polygon(polygon):
    points = polygon_to_points(polygon)
    plt.gca().invert_xaxis()
    plt.plot(points[:,0], points[:,1], 'b')
    plt.plot((points[0,0], points[-1,0]), (points[0,1], points[-1,1]), 'b')
    for i in range(points.shape[0]):
        plt.text(points[i,0], points[i,1], str(i+1), fontsize=14, fontweight='bold')
    plt.xlabel('Right Ascension')
    plt.ylabel('Declination')
    plt.show()

    
# Define the polygon to be checked - this can be copied from an obscore query
polygon = "POLYGON ICRS 356.72001015390003 -62.46047089520114 334.15651937444 -63.19263697411308 336.2140028315898 -48.593531680826004 351.68358252472024 -48.090537419574"

verify_polygon(polygon)
plot_polygon(polygon)