-1

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.

Ocean Scientist
  • 193
  • 1
  • 13

1 Answers1

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()
Ocean Scientist
  • 193
  • 1
  • 13