Python Work out area of a polygon on a spherical surface

Question

I have a series of points, of right ascension and declination values. These points correspond to the vertices of a polygon on the surface of a sphere.

What would be the best way to calculate the area enclosed by these points? I would assume that converting the points with an equal-area projection, and then carrying out typical polygonal area calculating on a flat surface would be an appropriate solution.

note: I cannot use custom python libraries. eg pyproj or shapely

Example code (works for latitude longitude, what modifications would be required to enure this works with sky coordinates?)

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)

I think this might belong on [CS.SE](https://cs.stackexchange.com/). You're just asking for an algorithm to do this. — Cyphase, Aug 13 '15 at 09:20

score 2 · Accepted Answer · answered Apr 12 '20 at 19:36

One of the ways is to perform a line integral based on Green's Theorem. See below an implementation, and this question for more details.

def polygon_area(lats, lons, algorithm = 0, radius = 6378137):
    """
    Computes area of spherical polygon, assuming spherical Earth. 
    Returns result in ratio of the sphere's area if the radius is specified.
    Otherwise, in the units of provided radius.
    lats and lons are in degrees.
    """
    from numpy import arctan2, cos, sin, sqrt, pi, power, append, diff, deg2rad
    lats = np.deg2rad(lats)
    lons = np.deg2rad(lons)

    # Line integral based on Green's Theorem, assumes spherical Earth

    #close polygon
    if lats[0]!=lats[-1]:
        lats = append(lats, lats[0])
        lons = append(lons, lons[0])

    #colatitudes relative to (0,0)
    a = sin(lats/2)**2 + cos(lats)* sin(lons/2)**2
    colat = 2*arctan2( sqrt(a), sqrt(1-a) )

    #azimuths relative to (0,0)
    az = arctan2(cos(lats) * sin(lons), sin(lats)) % (2*pi)

    # Calculate diffs
    # daz = diff(az) % (2*pi)
    daz = diff(az)
    daz = (daz + pi) % (2 * pi) - pi

    deltas=diff(colat)/2
    colat=colat[0:-1]+deltas

    # Perform integral
    integrands = (1-cos(colat)) * daz

    # Integrate 
    area = abs(sum(integrands))/(4*pi)

    area = min(area,1-area)
    if radius is not None: #return in units of radius
        return area * 4*pi*radius**2
    else: #return in ratio of sphere total area
        return area

Please find a somewhat more explicit version (and with many more references and TODOs...) here.

There is a flaw in this function. To check the closure of a polygon, you need to consider both the lats and lons, not lats alone. I found this issue because I used it to calculate a rectangle latlon cell and the result is wrong. Other than that, it works great! — Chang, Jun 01 '22 at 19:02
@Chang This is somewhat of an old answer. Could you elaborate? What coordinates precisely do you use, and what do you mean by "check the closure of a polygon"? Any chance that the function computed the complement area of what you intended? There is a well-known ambiguity to the definition of a polygon on a sphere. — Yellows, Jun 03 '22 at 06:53
I made my edit about the polygon closure part: #close polygon if lats[0]==lats[-1] and lons[0]==lons[-1]: pass else: lats = append(lats, lats[0]) lons = append(lons, lons[0]) — Chang, Jun 03 '22 at 23:16

score 0 · Answer 2 · answered Aug 13 '15 at 09:29

0

Looks like I can treat ra and dec like lat and long, work out the area on the Earth's surface in m^2, and use this value to convert into an area in sq degrees.

Please let me know if the solution I propose below is flawed:

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)

answered Aug 13 '15 at 09:29

Samuel

157
6
22

1

You'd better add this to your question and remove this answer. – Eli Korvigo Aug 13 '15 at 09:33
@EliKorvigo Not necessarily: https://stackoverflow.com/help/self-answer – Nulano Mar 31 '19 at 17:30
@Nulano this does not look like a classical self-answer, because the OP is not sure, whether the solution is correct, and asks for further assistance: "*Please let me know if the solution I propose below is flawed*". – Eli Korvigo Mar 31 '19 at 18:48
@EliKorvigo I'm not sure about this one. That's why I wrote "Not necessarily" instead of "You are wrong". But I would argue, that even with a regular self-answer, the OP is always looking for better answers / comments about issues with their answer. – Nulano Apr 01 '19 at 06:55

Python Work out area of a polygon on a spherical surface

2 Answers2