这里有一个近似解

def surf_area_on_ellipsoid(lat0, lat1, delta_lon, a, b):
    """
    Function to be used to get approximate grid cell areas
    
    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
    """
    assert(b <= a)
    d2r = np.pi / 180.
    # latitude width in radians
    delta_lat = d2r * (lat1 - lat0)
    # central latitude in radians
    lat = .5 * d2r * (lat1 + lat0)
    # arc length in latitude direction
    ds_lat = delta_lat * np.hypot(a * np.sin(lat), b * np.cos(lat))
    # arc length in longitude direction
    ds_lon = a * np.cos(lat) * d2r * delta_lon
    return ds_lon * ds_lat

下面是一个完整的椭圆体的精确解来检查它（从https://math.stackexchange.com/questions/2050158/deriving-formula-for-surface-area-of-an-ellipsoid）：

def surf_area_ellipsoid(a, b):
    """
    Exact formula for total area of an ellipsoid (for testing)
    Oblate spheroid case: https://math.stackexchange.com/questions/2050158/deriving-formula-for-surface-area-of-an-ellipsoid
    
    Parameters
    ----------
    a : float
        major semi-axis
    b : float
        minor semi-axis
    
    Returns
    -------
    area : float
    """
    if a == b:
        return 4 * np.pi * a ** 2
    t = np.sqrt(a ** 2 - b ** 2)
    s = b ** 2 * np.arcsinh(t / b) / t
    return 2 * np.pi * a * (a +  s)

在半度分辨率下，数值方法与精确的总面积相当一致：

a = 1
for b in [1, .99, .9, .5, .1]:
    exact = surf_area_ellipsoid(a=a, b=b)
    lat = np.linspace(-90., 90., 360)
    sa = np.sum(surf_area_on_ellipsoid(lat[:-1], lat[1:], 360, a=a, b=b))
    print(sa, exact)

12.566410711226226 12.566370614359172
12.482719164468133 12.482679067595733
11.737538672449274 11.737498575531172
8.671922800413613 8.671882703345052
6.472242602987875 6.472202505854831

它可能已经足够好了，但它可能值得尝试用测试的方法来做。

这段代码现在给出了精确的面积，而不会太复杂（只需要对球体进行不同的处理）

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))