CODEWITHSUNDEEP
pythonmath
Samuel
Samuel

Reputation: 119

Python Work out area of a polygon on a spherical surface

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)

Upvotes: 1

Views: 1853

Answers (2)

Yellows
Yellows

Reputation: 723

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.

Upvotes: 2

Samuel
Samuel

Reputation: 119

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)

Upvotes: 0

Related Questions