Given an ellipsoid with major semi-axis a and minor semi-axis b, and a regular lon-lat grid, what is the best way to get the areas of each grid cell? I would like to use python and preferably cartopy.
This code now gives the exact areas without being too much more complicated (just need to treat the sphere differently)
def surf_area_on_sphere(lat0, lat1, delta_lon, r):
"""
Function to get the area of a portion of a sphere
Parameters
----------
lat0 : float or array
lower limit(s) of latitude interval(s) in degrees
lat1 : float or array
upper limit(s) of latitude interval(s) in degrees
delta_lon : float or array
width of longitude interval(s) in degrees
r : float
radius
Returns
-------
area : float or array
"""
d2r = np.pi / 180.
return d2r * delta_lon * r ** 2 * (
np.sin(d2r * lat1) - np.sin(d2r * lat0))
def anti_deriv(lat, a, b):
"""
Parameters
----------
lat : float or numpy array
Latitude in radians
a : float
major semi-axis
b : float
minor semi-axis
Returns
-------
ad : float or numpy array
value of the antiderivative with respect to z for
1 / (2 * b) * sqrt(1 + ((a / b) ** 2 - 1) * (z / b) ** 2)
"""
r = np.sin(lat) # z /b
g = np.sqrt(a ** 2 - b ** 2) / b
return np.arcsinh(g * r) / g + r * np.sqrt(1 + (g * r) ** 2)
def surf_area_on_ellipsoid(lat0, lat1, delta_lon, a, b):
"""
Function to get the area of a portion of an ellipsoid.
Based on https://math.stackexchange.com/questions/2050158/deriving-formula-for-surface-area-of-an-ellipsoid
Parameters
----------
lat0 : float or array
lower limit(s) of latitude interval(s) in degrees
lat1 : float or array
upper limit(s) of latitude interval(s) in degrees
delta_lon : float or array
width of longitude interval(s) in degrees
a : float
major semi-axis
b : float
minor semi-axis
Returns
-------
area : float or array
"""
if a == b:
return surf_area_on_sphere(lat0, lat1, delta_lon, a)
assert(b < a)
d2r = np.pi / 180.
return .5 * d2r * delta_lon * a * b * (
anti_deriv(d2r * lat1, a , b) - anti_deriv(d2r * lat0, a, b))