from pyproj import Geod
from shapely.geometry import LineString, Point, Polygon
import pandas as pd

# function takes lat+lon gridcell centres as well as the resolution as input
def get_latlon_areas(latcentres, loncentres, latres, lonres):
    
    """Calculate the area of gridcells on a regular grid
       using polygons

     Parameters:
     latcentres (array): The latitude centres of the gridcells
     loncentres (array): The longitude centres of the gridcells
     latres (float): The latitude resolution
     lonres (float): The longitude resolution
    
     Returns:
     xarray.DataArray: DataArray of areas in metres squared
    
    """

    geod = Geod(ellps="WGS84")
    
    areas = pd.Series(dtype=float, index=latcentres)
    for clat in latcentres:
        # get bounds from coords
        lat_south, lat_north = clat - latres/2, clat + latres/2
        # (only need to do it for one longitude)
        lon_west, lon_east = 0 - lonres/2, 0 + lonres/2

        # create polygon
        poly = Polygon(

               LineString([
                   Point(lon_east, lat_south),
                   Point(lon_east, lat_north),
                   Point(lon_west, lat_north),
                   Point(lon_west, lat_south)
                   ])

            )
        
        # this computes the area in metres ^2
        area, _ = geod.geometry_area_perimeter(poly)
        areas.loc[clat] = area
    
    # make a dataarray
    da = xr.DataArray(areas, dims={"lat":latcentres})
    # add longitudes in
    da = da.expand_dims({"lon":loncentres})
    
    return da

import numpy as np
import xarray as xr

res = 1
global_grid = xr.DataArray(coords=[("lon", np.arange(-180+res/2, 180, res)),
                                   ("lat", np.arange(-90+res/2, 90, res))])

print(global_grid)

<xarray.DataArray (lon: 360, lat: 180)>
array([[nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       ...,
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan]])
Coordinates:
  * lon      (lon) float64 -179.5 -178.5 -177.5 -176.5 ... 177.5 178.5 179.5
  * lat      (lat) float64 -89.5 -88.5 -87.5 -86.5 -85.5 ... 86.5 87.5 88.5 89.5

# calculate area of earth

areas = get_latlon_areas(latcentres=global_grid.lat.values, 
                         loncentres=global_grid.lon.values,
                         latres=res, lonres=res)
earth_area_km2 = areas.sum()/1e6

# according to wikipedia...
wiki_earth_area = 510072000

print(f'The area of the Earth is {earth_area_km2} km²')
print(f'This is only {(wiki_earth_area-earth_area_km2)/wiki_earth_area*100:.5f}%',
     'different to the value on wikipedia')

The area of the Earth is <xarray.DataArray ()>
array(5.10065622e+08) km²
This is only 0.00125% different to the value on wikipedia

Calculating the area of gridcells

Code