Model Grid

In this notebook, we will explore how to create the grid of a model.

First, import packages to re-create and visualize the model grid here:

import numpy as np
import matplotlib.pyplot as plt
import cartopy
import cartopy.crs as ccrs

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[1], line 3
      1 import numpy as np
      2 import matplotlib.pyplot as plt
----> 3 import cartopy
      4 import cartopy.crs as ccrs

ModuleNotFoundError: No module named 'cartopy'

The California Coast

The grid for this model will be located on the west coast of California covering 135-115\(^{\circ}\)W in longitude and 29-52\(^{\circ}\)N in latitude. The grid spacing will be \(1/12^{\circ}\) in the zonal (east-west) direction and \(1/16^{\circ}\) in the meridional (north-south) direction, covering a grid of 360 rows and 240 columns.

In the data file for my model, I will specifiy the following parameters in the PARM04 namelist in the data file:

usingSphericalPolarGrid=.TRUE.,
delX=240*0.083333,
delY=360*0.0625,
xgOrigin=-135,
ygOrigin=29,

This grid can be recreated in Python as follows:

# define the parameters that will be used in the data file
delX = 1/12
delY = 1/16
xgOrigin = -135
ygOrigin = 29
n_rows = 360
n_cols = 240

# recreate the grids that will be used in the model
xc = np.arange(xgOrigin+delX/2, xgOrigin+n_cols*delX+delX/2, delX)
yc = np.arange(ygOrigin+delY/2, ygOrigin+n_rows*delY+delY/2, delY)
XC, YC = np.meshgrid(xc, yc)

Visualizing the Grid

The grids above can be visualized as follows:

# make a plot of XC and YC
plt.figure(figsize=(10,7))

plt.subplot(1,2,1)
C = plt.pcolormesh(XC)
plt.colorbar(C, orientation = 'horizontal')
plt.title('XC (Longitude)')

plt.subplot(1,2,2)
C = plt.pcolormesh(YC)
plt.colorbar(C, orientation = 'horizontal')
plt.title('YC (Latitude)')

plt.show()

Visualizing the Grid on a Map with Cartopy

To get a sense of where the model is located on the globe, cartopy can be be used to plot the domain on the globe:

plt.figure(figsize=(5,5))
ax = plt.axes(projection=ccrs.Orthographic(-130,10))
ax.plot(XC[:,0], YC[:,0], 'g-', transform=ccrs.PlateCarree())
ax.plot(XC[:,-1], YC[:,-1], 'g-', transform=ccrs.PlateCarree())
ax.plot(XC[0,:], YC[0,:], 'g-', transform=ccrs.PlateCarree())
ax.plot(XC[-1,:], YC[-1,:], 'g-', transform=ccrs.PlateCarree())
ax.coastlines()
ax.set_global()
plt.show()

Visualizing the Grid Spacing

The model grid is defined in terms of units in longitude and latitude although it is useful to quantify the grid spacing in terms of more familiar units, such as meters. The following great_circle_distance function can be used to quantify this distance:

def great_circle_distance(lon_ref, lat_ref, Lon, Lat):
    earth_radius = 6371000
    lon_ref_radians = np.radians(lon_ref)
    lat_ref_radians = np.radians(lat_ref)
    lons_radians = np.radians(Lon)
    lats_radians = np.radians(Lat)
    lat_diff = lats_radians - lat_ref_radians
    lon_diff = lons_radians - lon_ref_radians
    d = np.sin(lat_diff * 0.5) ** 2 + np.cos(lat_ref_radians) * np.cos(lats_radians) * np.sin(lon_diff * 0.5) ** 2
    h = 2 * earth_radius * np.arcsin(np.sqrt(d))
    return(h)

Then, loop through the points to generate inter-point distances in the horizontal (dXC) and vertical (dYC) directions:

dXC = np.zeros((np.shape(XC)[0], np.shape(XC)[1]-1))
for row in range(np.shape(XC)[0]):
    for col in range(np.shape(XC)[1]-1):
        dXC[row,col] = great_circle_distance(XC[row,col], YC[row,col], XC[row,col+1], YC[row,col+1])

dYC = np.zeros((np.shape(YC)[0]-1, np.shape(YC)[1]))
for row in range(np.shape(XC)[0]-1):
    for col in range(np.shape(XC)[1]):
        dYC[row,col] = great_circle_distance(XC[row,col], YC[row,col], XC[row+1,col], YC[row+1,col])

Finally, make a plot of the inter-point distances:

# make a plot of XC and YC
plt.figure(figsize=(10,5))

plt.subplot(1,2,1)
C = plt.pcolormesh(dXC)
plt.colorbar(C, orientation = 'horizontal')
plt.title('dXC')

plt.subplot(1,2,2)
C = plt.pcolormesh(dYC.round(3))
plt.colorbar(C, orientation = 'horizontal')
plt.title('YC')

plt.show()

As we can see the grid has a resolution of about 7 km, although there is a north-south gradient in horizontal distances (in other words, points further north are closer together).