Calculates the area enclosed by an arbitrary polygon on the sphere.

Parameters
----------

lat : iterable
    The latitudes, in degrees, of the vertex locations of the polygon, in
    clockwise order.
lon : iterable
    The longitudes, in degrees, of the vertex locations of the polygon, in
    clockwise order.

Returns
-------

area : float
    The desired spherical area in the same units as r_sphere.

Notes
-----

This function assumes the vertices form a valid polygon (edges do not
intersect each other).

**References**

Computing the Area of a Spherical Polygon of Arbitrary Shape
Bevis and Cambareri (1987)
Mathematical Geology, vol.19, Issue 4, pp 335-346

def area_poly_sphere(lat, lon, r_sphere):
    '''
Calculates the area enclosed by an arbitrary polygon on the sphere.

Parameters
----------

lat : iterable
    The latitudes, in degrees, of the vertex locations of the polygon, in
    clockwise order.
lon : iterable
    The longitudes, in degrees, of the vertex locations of the polygon, in
    clockwise order.

Returns
-------

area : float
    The desired spherical area in the same units as r_sphere.

Notes
-----

This function assumes the vertices form a valid polygon (edges do not
intersect each other).

**References**

Computing the Area of a Spherical Polygon of Arbitrary Shape
Bevis and Cambareri (1987)
Mathematical Geology, vol.19, Issue 4, pp 335-346
    '''
    dtr = np.pi/180.

    def _tranlon(plat, plon, qlat, qlon):
        t = np.sin((qlon-plon)*dtr)*np.cos(qlat*dtr)
        b = (np.sin(qlat*dtr)*np.cos(plat*dtr) -
             np.cos(qlat*dtr)*np.sin(plat*dtr)*np.cos((qlon-plon)*dtr))
        return np.arctan2(t, b)

    if len(lat) < 3:
        raise ValueError('lat must have at least 3 vertices')
    if len(lat) != len(lon):
        raise ValueError('lat and lon must have the same length')
    total = 0.
    for i in range(-1, len(lat)-1):
        fang = _tranlon(lat[i], lon[i], lat[i+1], lon[i+1])
        bang = _tranlon(lat[i], lon[i], lat[i-1], lon[i-1])
        fvb = bang - fang
        if fvb < 0:
            fvb += 2.*np.pi
        total += fvb
    return (total - np.pi*float(len(lat)-2))*r_sphere**2