How to calculate size in m2 of each lat/long grid square

Question

I have estimates for data in units m2 at gridsquare resolution. I need to calculate the number of m2 in each latitude / longitude grid cell?

Cell sizes are much smaller near the poles than at the equator so this is important.

I would like a netcdf file or array of the number of square meters in each grid square.

Answer 1

In case anyone would like a netcdf of the number of square meters in each lat long grid cell. This is probably not the cleanest solution, but will create a netcdf (earth_m2.nc) of m2 in each grid using xarray.

The gridsize() function is adapted from another stack overflow question. We can then make a dummy array and create a earth field of m2s using the longitude distances at each location.

"""
This will create a global grid of the approximate size of each grid square.
"""
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt

def gridsize(lat1):
   #https://en.wikipedia.org/wiki/Haversine_formula
   #https://stackoverflow.com/questions/639695/how-to-convert-latitude-or-longitude-to-meters/11172685#11172685
   lon1=200
   import math
   lat2=lat1
   lon2=lon1+1

   R = 6378.137 # // Radius of earth in km
   dLat = lat2 * np.pi / 180 - lat1 * np.pi / 180
   dLon = lon2 * np.pi / 180 - lon1 * np.pi / 180
   a = np.sin(dLat/2) * np.sin(dLat/2) + np.cos(lat1 * np.pi / 180) * np.cos(lat2 * np.pi / 180) * np.sin(dLon/2) * np.sin(dLon/2)
   c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1-a))
   d = R * c
   return d * 1000 #; // meters


boxlo,boxla=np.array(np.meshgrid(np.arange(-179.5,179.5,1),np.arange(-89.5,89.5,1)))
sizes=np.ones(boxlo.shape)
grid=gridsize(boxla)

grid_nc=xr.DataArray(grid,coords={'lat':boxla[:,1],'lon':boxlo[1,:]},dims=['lat','lon'])
lat_size=110567 #in m
grid_nc['m2']=grid_nc*lat_size
grid_nc=grid_nc['m2']
grid_nc.to_netcdf('earth_m2.nc')

plt.pcolormesh(boxlo[1,:],boxla[:,1],grid_nc)
plt.colorbar()
plt.show()

Answer 2

By spherical trigonometry, the surface area of a triangle on a sphere of radius R with 1 vertex at the North pole (latitude: π/2) and 2 vertices at the same latitude -π/2 < x < π/2 separated (longitudinally) by d radians is

S(x, d) = (cos⁻¹((cos(b) - cos²(a))/sin²(a)) + 2cos⁻¹((cos(a) - cos(a)cos(b))/(sin(a)sin(b))) - π)R² where a = R(π/2 - x) and b = Rd

So, the surface area of a grid rectangle on a sphere of radius R between lines of longitude separated by d radians and latitudes x₁ > x₂ is

S(x₂, d) - S(x₁, d)

Answer 3

One option is to transform your cells into a coordinate reference system (CRS) that has units in, say, meters rather than degrees. Then the area calculation is simple.

I assume your coordinates are in WGS84.

For the target CRS there are choices especially if you know the locality of the points, but a common collection of global CRSs like this are Universal Transverse Mercator (UTM), or near the poles Universal Polar Stereographic

For example, for UTM, assuming a list of points of the form [lon, lat] where the last point is equal to the first:

import pyproj
from shapely.geometry import Polygon
from shapely.ops import transform

def utm_epsg(lon: float, lat: float) -> int:
  """
  Return the UTM EPSG code for the given lon-lat.
  """
  offset = int(round((183 + lon) / 6.0))
  return 32600 + offset if lat > 0 else 32700 + offset

for lat in range(-79, 83):
    for lon in range(-179, 179):
        polygon = Polygon([
            [lon, lat],
            [lon+1, lat],
            [lon+1, lat+1],
            [lon, lat+1],
            [lon, lat],
        ])
        
        src_crs = pyproj.CRS.from_epsg(4326)
        tgt_crs = pyproj.CRS.from_epsg(utm_epsg(polygon.centroid.x, polygon.centroid.y))

        project = pyproj.Transformer.from_crs(src_crs, tgt_crs, always_xy=True).transform

        utm_polygon = transform(project, polygon)

        # aggregate into some result.  Here just printed to stdout.
        print(polygon.centroid, utm_polygon.area)

It's worth noting that UTM isn't defined south of 80°S and north of 84°N.