Python | Tutorial: Temperature Advection

Introduction:

Temperature advection is a useful field for meteorologists and weather forecasters since it is an easy way to diagnose the dynamic state of the atmosphere.  Plotted on a map, it is particularly useful for identifying frontal boundaries.  In this blog post, you will learn how to:

• Use netCDF4 to load GFS data into arrays
• Subset data by geographic boundaries
• Use Cartopy transform_points to convert lat/lon to dx and dy
• Use Numpy to compute temperature advection
• Use Cartopy to plot temperature advection
• Use Scipy to smooth data for cleaner figures

Temperature Advection:
I'm not going to demonstrate how to derive the equation for temperature advection, but essentially, it's the -1 multiplied by the dot product between the horizontal wind vector and the temperature gradient:

 Temperature advection equation

To plot this on a map from gridded data, the first, and most important step is to "discretize" the above formula:

Discretized advection

Now, all we've done is changed the partial derivative to a finite change in temperature over a change in x.  So, basically, to plot temperature advection on a map, we just need the u and v components of the wind, and temperature at a specified model level, and Δx/Δy.  The only tricky part is getting Δx/Δy from common lat/lon coordinates associated with most gridded meteorological datasets.

 
Load the Data:
In this example, I'll load data from the GFS model available from the NOMADS data server directly into Python using the netCDF4 package.  The first step in building the script is to import all of the packages needed to plot temperature advection:

import numpy as np
from netCDF4 import Dataset as ds
from datetime import datetime as DT
from cartopy import crs as ccrs
from cartopy import feature as cfeature
from matplotlib import pyplot as plt 

from scipy.ndimage import gaussian_filter

Now that the packages are loaded into Python, we'll use the "datetime" module to make this script dynamic, updating to the current date no matter when it's run.  To do this, we simply copy the url for the GFS data into a string, and replace the date and time attributes with strings formatted to match the URL format.  Note, that there are numerous ways to format strings, I am simply most comfortable with the syntax used in the example below.

hour='00'
now=DT.now().strftime('%Y%m%d')
gfs_url='https://nomads.ncep.noaa.gov:9090/dods/gfs_0p25/gfs%s/gfs_0p25_%sz_anl'%(now,hour)

Basically, we chose the '00' hour model analysis, and put todays date into the yyyymmdd format.  These strings are then passed to the url such that it matches the current date.  Now that the url is set, it's time to load the data using the netCDF4 package.

data=ds(gfs_anal,'r')

For brevity, since this isn't a tutorial on using the netCDF package, I'll just give you the code to load the lat/lon/level variables into arrays.  Basically, we do that by accessing the "variables" dictionary of the dataset.

lon=dataset.variables['lon'][:]
lat=dataset.variables['lat'][:]
levs=dataset.variables['lev'][:]

Now, longitude, latitude and pressure levels are in arrays.  Since this is a large dataset, we'll want to subset the data to match geographic boundaries so we don't calculate temperature advection over the entire globe since we only care about a small region.  Since, in this instance, longitude and latitude are independent 1D vectors this is a simple task, simply choose lat/lon boundaries and use the Numpy "argmin" function which returns the index where a given operation is minimized:

boundbox=[-120+360,-65+360,20,55]
level=850 ##850 mb pressure level
x1,x2=np.argmin(np.abs(lon-boundbox[0])),np.argmin(np.abs(lon-boundbox[1]))
y1,y2=np.argmin(np.abs(lat-boundbox[2])),np.argmin(np.abs(lat-boundbox[3]))
z=np.argmin(np.abs(levs-level))

Basically, the array indices are identifed by taking the value of the absolute difference bewteen the array and the value and seeing where it's closest to zero. Now that the indices are found, load the wind and temperature data:

T=dataset.variables['tmpprs'][0,z,y1:y2+1,x1:x2+1]
U=dataset.variables['ugrdprs'][0,z,y1:y2+1,x1:x2+1]
V=dataset.variables['vgrdprs'][0,z,y1:y2+1,x1:x2+1]

Note that the +1 for the ending index accounts for the indexing structure of Python.

Coordinate Transforms and Computing Advection:
Transforming the lat/lon to physical x/y space to compute dx and dy is the critical challenge in this exercise.  Instead of crudly approximating a constant conversion factor or running a potentially slow and complicated haversine formula, we'll take advantage of cartopys "transform_points" function to transform the cylindrical lat/lon points to a coordinate system that uses x/y distance coordinates.  If you are unfamiliar with Cartopy, I strongly suggest you read through my introduction to cartopy post: Intro to Cartopy.  UTM is a good choice for smaller domains, I like the LambertConformal projection for this domain, simply set the central longitude argument to something close to the center of the domain.

proj=ccrs.LambertConformal(central_longitude=-90)
lon,lat=np.meshgrid(lon,lat)
output=proj.transform_points(ccrs.PlateCarree(),lon,lat)
x,y=output[:,:,0],output[:,:,1]

The transform points function is called from a defined projection object (in this example, the Lambert Conformal projection), and takes an input projection (PlateCarree, or cylindrical lat/lon coordinates), and then the x (lon) and y(lat), and (if applicable) z coordinates you are transforming.  It returns an array sized (nx,ny,3) where the 3 cooresponds to the x,y,z coordinates.  Note that prior to input into the function, the "meshgrid" function is used to reset the 1D lat/lon arrays into 2D arrays that have the same shape.  This is required for transform_points to work.

Now that the x,y coordinates are defined, the numpy "gradient" function is used to get the dx and dy for temperature advection (see Intro to Numpy script for more information on "gradient").

gradx=np.gradient(x,axis=1)
grady=np.gradient(y,axis=0)

Now, all the required pieces are required, and all that's left to do is calculate temperature advection:

Tadv=-(U*(np.gradient(T,axis=1)/gradx)+V*(np.gradient(T,axis=0)/grady))

Coordinate Transforms and Computing Advection:
The below code uses Cartopy to plot temperature advection:

skip=7

fig=plt.figure(figsize=(13,12))
ax=plt.subplot(111,projection=ccrs.PlateCarree())
Z=plt.contourf(x,y,Tadv*3600.,cmap='coolwarm',levels=np.linspace(-2.5,2.5,31),transform=proj,extend='both')
C=plt.contour(x,y,T-273.15,cmap='jet',levels=np.linspace(-25,20,10),transform=proj)
plt.barbs(x[::skip,::skip],y[::skip,::skip],U[::skip,::skip],V[::skip,::skip],length=4.5,transform=proj)
plt.colorbar(Z)
ax.set_extent(box,crs=ccrs.PlateCarree())
ax.add_feature(cfeature.COASTLINE.with_scale('50m'))
ax.add_feature(cfeature.STATES.with_scale('50m'))

This produces an image similar to the following:

Temperature Advection ( 00 UTC Nov 24 2019)

 

Note that advection is multiplied by 3600 to convert the units from K/s to a more managable K/hr.  Furthermore, note that the x/y arrays are input into the contourf function with the transform keyword set to the projection used when converting from lat/lon to xy.  It's messy, but you should be able to get an idea as to where the cold and warm air advection is on the map.  One way to clean up this figures is to use the "gaussian_filter" function from the SciPy data-anaylsis package to smooth out the data.  It's pretty straightforward to use, you simply supply the function with the array you want to smooth, and a "sigma" value to determine the degree of smoothing.

Tadv_smooth=gaussian_filter(Tadv, sigma=3.0)
T_smooth=gaussian_filter(T, sigma=3.0)

mslp=dataset.variables['prmslmsl'][0,y1:y2+1,x1:x2+1]
mslp_smooth=gaussian_filter(mslp, sigma=3.0)

In this code the temperature data, and temperature advection data is smoothed, and then the sea-level pressure  variable is loaded and smoothed as well to help visually connect the storm centers to the advection: Tweaking the plotting code to account for these changes:

fig=plt.figure(figsize=(13,12))
ax=plt.subplot(111,projection=ccrs.PlateCarree())
Z=plt.contourf(x,y,Tadv_smooth*3600.,cmap='coolwarm',levels=np.linspace(-2.5,2.5,31),transform=proj,extend='both')
C=plt.contour(x,y,T_smooth-273.15,cmap='jet',levels=np.linspace(-25,20,10),transform=proj)
C=plt.contour(x,y,mslp_smooth/100.,colors='k',levels=np.arange(900,1040,2),transform=proj)
plt.clabel(C,fmt='%i',fontsize=12)
plt.barbs(x[::skip,::skip],y[::skip,::skip],U[::skip,::skip],V[::skip,::skip],length=4.5,transform=proj)
plt.colorbar(Z)
ax.set_extent(box,crs=ccrs.PlateCarree())
ax.add_feature(cfeature.COASTLINE.with_scale('50m'))
ax.add_feature(cfeature.STATES.with_scale('50m'))
plt.show()

And we get the following:

Smoothed Temperature Advection ( 00 UTC Nov 24 2019)

In the above example, a nice cold front is seen extending along the southern Appalachian mountains from a low pressure center over Ohio, with a weak warm front over the Mid-Atlantic states.

That is pretty much it.  There are a few caveats with this calculation that I won't go into great detail over, since they are negligble for the purposes of the synoptic diagnostic tool.  I hope this tutorial was informative, and thanks for reading.