Python Разработайте область многоугольника на сферической поверхности

Question

У меня есть ряд точек, значений прямого восхождения и склонения. Эти точки соответствуют вершинам многоугольника на поверхности сферы.

Каким будет лучший способ рассчитать площадь, приложенную этими точками? Я бы предположил, что преобразование точек с проекцией на равную площадь, а затем выполнение типичного вычисления многоугольной площади на плоской поверхности было бы подходящим решением.

примечание: я не могу использовать пользовательские библиотеки python. например pyproj или стройная

Пример кода (работает для широты долготы, какие изменения должны будет вступать в силу этого работает с координатами неба?)

def reproject(latitude, longitude): 
    """Returns the x & y coordinates in metres using a sinusoidal projection""" 
    from math import pi, cos, radians 
    earth_radius = 6371009 
    lat_dist = pi * earth_radius/180.0 

    y = [lat * lat_dist for lat in latitude] 
    x = [long * lat_dist * cos(radians(lat)) 
       for lat, long in zip(latitude, longitude)] 
    return x, y 

def area_of_polygon(x, y): 
    """Calculates the area of an arbitrary polygon given its vertices""" 
    area = 0.0 
    for i in xrange(-1, len(x)-1): 
     area += x[i] * (y[i+1] - y[i-1]) 
    return abs(area)/2.0 

dec = [-15.,89.,89.,-15.,-15.]  
ra = [105.,105.,285.,285.,105.] 

x,y = reproject(dec, ra) 
print area_of_polygon(x,y) 

Answer 1

Похоже, я могу относиться к ра и разлу как латы и долго, выработать область на поверхности Земли в m^2 и использовать это значение для преобразования в область в квадратных градусах.

Пожалуйста, дайте мне знать, если решение, которое я предлагаю ниже недостатки:

def reproject(latitude, longitude): 
    """Returns the x & y coordinates in metres using a sinusoidal projection""" 
    from math import pi, cos, radians 
    earth_radius = 6371009 
    lat_dist = pi * earth_radius/180.0 

    y = [lat * lat_dist for lat in latitude] 
    x = [long * lat_dist * cos(radians(lat)) 
       for lat, long in zip(latitude, longitude)] 
    return x, y 

def area_of_polygon(x, y): 
    """Calculates the area of an arbitrary polygon given its vertices""" 
    area = 0.0 
    for i in xrange(-1, len(x)-1): 
     area += x[i] * (y[i+1] - y[i-1]) 
    return ((abs(area)/2.0)/5.10100E14) * 41253 

dec = [-15.,89.,89.,-15.,-15.]  
ra = [105.,105.,285.,285.,105.] 

x,y = reproject(dec, ra) 
print area_of_polygon(x,y)