geodesic_area.py

import numpy as np
from math import pi, sqrt

def get_area(lat: float, d_lat: float, d_lng: float) -> float:
	"""
	Calculate geodesic area for lat/lng pixels using WGS 1984 as spatial reference.
	Pixel area for a fix pixel width and hight various with latitude.
	Pixels close to the equator have a larger areas. Area decreases towards the poles.
	
	Arguments:
		lat: Latitude of pixel
		d_lat: pixel hight
		d_lng: pixel width	
	Returns:
		The geodesic area of the pixel in m2
	
	"""
	a: float = 6378137.0  # Semi major axis of WGS 1984 ellipsoid
	b: float = 6356752.314245179  # Semi minor axis of WGS 1984 ellipsoid
	
	q: float = d_lng/360
	e: float = sqrt(1 - (b/a)**2)

	area: float = abs(
			(pi * b**2 * ( 
				2*np.arctanh(e*np.sin(np.radians(lat+d_lat))) / 
					(2*e) + 
				np.sin(np.radians(lat+d_lat)) / 
					((1 + e*np.sin(np.radians(lat+d_lat)))*(1 - e*np.sin(np.radians(lat+d_lat)))))) -
			(pi * b**2 * ( 
				2*np.arctanh(e*np.sin(np.radians(lat))) / (2*e) + 
				np.sin(np.radians(lat)) / ((1 + e*np.sin(np.radians(lat)))*(1 - e*np.sin(np.radians(lat))))))) 
			* q
	
	return area