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.

Read More

Python | Tutorial: Intro to Cartopy

Introduction:

There are currently two main Python libraries for plotting geographic data on map: Cartopy and Basemap.  New users should use the Cartopy since support Cartopy will replace Basemap and support for Basemap is expected to wrap up in 2020.  Therefore, this tutorial will focus on Cartopy.  If you want to use Basemap (e.g., for plotting in "3D"), here are a couple of links to existing tutorials and examples to help get you started.  Furthermore, many concepts described here may be useful for understanding Basemap.

• Basemap Examples
• Basemap Tutorials 

Cartopy:

Cartopy is a geospatial plotting library built on top of Numpy and Matplotlib that makes plotting gridded data, shapefiles, and other geographic data on over 30 different map projections.  Furthermore, the Cartopy "transform" functionality makes it straightforward to convert data from one projection to another.  In this tutorial you will learn how to use Cartopy to:

• Plot GFS surface temperature data on a map
• Mask out land surfaces to plot a map of sea surface temperature and ice-cover from the GFS
• Plot a regional map of surface temperature with US states
• Use the transform function to plot GFS data on a different map projection

This tutorial accesses the NOMADS data server using the netCDF4 library, if you are unfamiliar with doing this, I recommend you see my tutorial on reading NetCDF data.

Getting Started - Grabbing GFS Data:
The following code is used to get started with this tutorial, essentially, it imports all of the required modules, pulls the GFS data from the NOMADS data server, and loads it into arrays.  Once the data is loaded, the remainder of the tutorial will focus on Cartopy. 

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

import cartopy.feature as cfeature

import netCDF4 as nc
import datetime as DT 

datetime=(DT.datetime.utcnow()-DT.timedelta(hours=6)).strftime('%Y%m%d')
fhour=20

nomads_path="https://nomads.ncep.noaa.gov:9090/dods/gfs_0p25/gfs{}/gfs_0p25_00z".format(datetime)

datafile=nc.Dataset(nomads_path)
lat=datafile.variables['lat'][:]
lon=datafile.variables['lon'][:] 

sfcT=datafile.variables['tmpsfc'] #Keep meta data for now!
sfcT_data=sfcT[fhour,:]

title=sfcT.long_name
time=datafile.variables['time']
tstart=DT.datetime.strptime(time.minimum,'%Hz%d%b%Y')
timestr=(tstart+DT.timedelta(hours=time.resolution*24.*fhour)).strftime('%c')

Now that your data is loaded, we initialize a figure, and subplot using matplotlibs "projection" keyword to link to the Cartopy library.  Once you've attached the Cartopy package to your subplot, there are a number of additional objects and classes you can add to your figure.  In the first example, I'll demonstrate by calling upon the "coastlines" attribute to plot continents on the map.

fig=plt.figure(figsize=(11,5))
ax=plt.subplot(111,projection=ccrs.PlateCarree())

Now, your figure and subplot are defined.  The Cartopy "PlateCarree()" projection is a basic cylindrical map projection that accepts coordinate information in the form of lat/lon coordinate pairs.  The lat/lon data can be either 1D or 2D, but essentially, if your coordinate information is lat/lon, your projection is PlateCarree().   You can pass 2 keywords into PlateCarree (central_longitude and globe), but in most instances, you won't need to.   Occasionally, I reset the central_longitude to 180, instead of it's default zero.

The code to plot the map is:

Z=plt.pcolormesh(lon,lat,sfcT_data,cmap='jet')

pos=ax.get_position()

plt.title(title,fontweight='bold',loc='left')
plt.title("Valid:{} UTC".format(timestr),loc='right')
ax.coastlines()

cbar_ax=fig.add_axes([pos.x1+0.01,pos.y0,0.015,pos.height])
cbar=plt.colorbar(Z,cax=cbar_ax)
cbar.set_label(title)

datafile.close()

plt.show()

And you'll get an image similar to this one:

GFS surface temperature from Cartopy

Now, lets say you wanted to mask out land areas, and add sea-ice to the map, to focus on the oceans.
Adding the sea-ice is easy, all it requires is for you to pull the icecsfc variable from the GFS data, and plot it overtop of the surface temperature map:

plt.pcolormesh(lon,lat,np.ma.masked_less(ice,0.1),vmin=0,vmax=1,cmap='Greys_r',zorder=3)

We use the numpy mask (ma) functionality to maskout gridcells where the ice cover is less than 10%.
To mask out land areas, we make use of the "Cartopy Feature" module, which handles shapefiles.  The Cartopy Feature module connects seemlessly to the Natural Earth Dataset, and allows for anyone to plot a host of different GIS datasets without needing to download individual shapefiles before hand.  The interface even has many pre-defined functions, including land.  This makes masking out the land areas, a breeze.

ax.add_feature(cfeature.LAND,zorder=4,color='gray')

In the above code, the "add_feature" function is used, and the pre-defined "LAND" feature from the Cartopy Feature package is supplied with a couple of basic keyword arguments.  Putting it all together:

fig=plt.figure(figsize=(11,5))
ax=plt.subplot(111,projection=ccrs.PlateCarree())
Z=plt.pcolormesh(lon,lat,sfcT_data-273.15,cmap='jet',vmin=-1,vmax=30.)
plt.pcolormesh(lon,lat,np.ma.masked_less(ice,0.1),vmin=0,vmax=1,cmap='Greys_r',zorder=3)
ax.add_feature(cfeature.LAND,zorder=4,color='gray')
pos=ax.get_position()

plt.title(title,fontweight='bold',loc='left')
plt.title("Valid:{} UTC".format(timestr),loc='right')
ax.coastlines(color='k',zorder=5)

cbar_ax=fig.add_axes([pos.x1+0.01,pos.y0,0.015,pos.height])
cbar=plt.colorbar(Z,cax=cbar_ax)
cbar.set_label(title)

datafile.close()

plt.show()

Running the above code will give you an image like this:

GFS SST and Sea-ice cover

Now, lets zoom in over the United States, and add states to the map.  To clip the extent, the "set_extent" function is used:

ax.set_extent([-125, -65, 25, 55],crs=ccrs.PlateCarree())

The list corresponds to lon/lat boundaries from [western-most longitude, eastern-most longitude, southern-most latitude, northern-most latitude] and the crs argument indicates which transform is being applied.  To add the states, you need to define the states from Natural Earth database. the "scale" argument

states = cfeature.NaturalEarthFeature(category='cultural', scale='50m', facecolor='none',
                         name='admin_1_states_provinces_shp',edgecolor='k')

 By changing the extent, and adding the states, you'll get a map that looks like this:

GFS SST centered over the United States

You'll notice that the "LAND" feature masks out the Great Lakes, which is unfortunate, and to my knowledge the only way to get both the land mask and the SST over lakes using GFS surface temperature data, is to use a Lakes shapefile with shapefile clipping, a technique that is beyond this simple introduction, but may be the subject of a future tutorial.

Finally, if you want to use a different map projection to plot the data, e.g., a North Polar Sterographic projection, you simply define a different projection when setting up your subplot, and use the "transform" keyword in the plotting function to specify that you are transforming your lon/lat coordinates from the PlateCarree() projection:

ax=plt.subplot(111,projection=ccrs.NorthPolarStereo())
Z=plt.pcolormesh(lon,lat,sfcT_data-273.15,cmap='jet',vmin=-1,vmax=30.,transform=ccrs.PlateCarree())
plt.pcolormesh(lon,lat,np.ma.masked_less(ice,0.1),vmin=0,vmax=1,cmap='Greys_r',zorder=3,transform=ccrs.PlateCarree())

ax.set_extent([-180, 180, 45, 90],crs=ccrs.PlateCarree()) 

Note, that we also reset the extent to focus on the Northern Latitudes.  You should get something that looks like this:

GFS Sea-ice and SST Polar Plot

That is the extent of this introductory tutorial on how to use Cartopy to work with atmospheric data.  You can learn more from the Cartopy Website, and get a list of different map projections here Cartopy Map Projections. Future tutorials on Cartopy will focus on working with shapefiles, transforming individual points, how to work with data without lat/lon coordinates, and how to put fine details on your map.

Read More

Tutorial: Use basemap.gs to make a beautiful map of "SST", Sea Ice and Snow in GrADS

This tutorial will draw upon many of the skills discussed in several other tutorials on this site, e.g., how to handle multiple files at once, or how to use basemap.gs.  There isn't too much "new" information in this tutorial on how to use these different features and functions of GrADS, rather this tutorial will show you how to specifically make a very pretty map of sea-surface temperature, sea ice, and snow.  The final outcome is shown below.

For this tutorial you will need:

• color.gs 
• xcbar.gs
• colormaps.gs
• basemap.gs (plus the accompanying .asc files)

Snow Cover and Sea Ice on the Robinson Map Projection

So before we start, I want to give a full disclaimer: The plot is not actually showing any observed SST.  In this tutorial I'm going to use the 0.25 degree GFS surface temperature.  In practice the surface temperature should be roughly equal to the SST over the Oceans.

The first thing we need to do is open both the GFS data and the sea ice data.  This data will come from NOMADS and will be opened using the 'sdfopen' command.  Note, the files below are for February 2015, so keep in mind, that you will need to change the date on files if you wish to copy paste the example code below.


   'reinit'
   gfsfile='http://nomads.ncep.noaa.gov:9090/dods/gfs_0p25/gfs20150218/gfs_0p25_18z'
   icefile='http://nomads.ncep.noaa.gov:9090/dods/ice/ice20150218/ice.00z'
  'sdfopen 'gfsfile
  'sdfopen 'icefile


Now that the files are open we will simply set up the map.  Since we are using the Robinson Projection, we need to set the longitude to range from -180 to 180 and the latitude -90 to 90.  

   'set gxout shaded'
   'set mpdset hires'
   'set lon -180 180'
   'set lat -90 90'
   'set mproj robinson'
   'colormaps -l 272 307 0.5 -map jet' ;*Note the use of colormaps.gs


Then we display the variable:

   'd tmpsfc'
   'xcbar -fs 4'
  
After a moment, the surface temperature will be displayed following the "jet" color map.
Now, we are going to mask out the land using the land basemap from basemap.gs.

Note, in my version of basemap.gs, I did not want to include the frame around the map projection.  So I opened up basemap.gs and commented out the lines at the bottom that draw the rectangle:

   * Draw a new square frame around the plot
   * If you have 'set frame off' or 'set frame circle' before running basemap,
   * you may want to comment out the next 2 lines
   *'set line 1 1 6'
   *'draw rec 'x1' 'y1' 'x2' 'y2

We are going to use the "medium (M)" resolution option for the coastlines, and we are going to define a custom dark gray color for the land.

   'set rgb 73 80 80 80'
   'basemap L 73 0 M'         ;*L = Land, 73 = Fill Color, 0 = outline Color, M=Medium resolution

SST only!

  This will fill in the land areas with a gray map and you will have an image like this:

Now that we are done with the "SST" and the land mask, we are going to open up the sea-ice file and plot the sea ice.  It is important that you use the maskout function when plotting the sea-ice, otherwise you will overwrite all of the SST data with zeros.

  'set gxout shaded'
  'set dfile 2'                   ;*Sets Default file to file 2 (Sea ice)
  'set z 1'                        ;*Set z to 1
  'set t 1'                        ;*Set t to 1
  'set map 0 1 6'   
  'color 0 0.6 0.1 -kind dimgray->seashell->white'           ;*Note use of color.gs
  'd maskout(icecmsl,icecmsl-0.1)'                                     ;*Need to mask out ice.


After a few moments, the sea ice will plot and you will have an image like this:

Sea Ice and SST only!

Now this is all well and good, but our knowledge is limited to the ice over the oceans, we don't really get any sense of how ice and snow covers the land masses.  So the final touch that we'll add to the map is to plot the GFS snow water equivalent (SWE) using a similar color scale as ice cover.  To do that, we need to reset the file to our first file and then reset our time our vertical level.  Lastly, we plot SWE (again using maskout to avoid plotting over everything).

  'set dfile 1'
  'color 10 100 5 -kind dimgray->seashell->white'
  'set z 1'
  'set t 1'
  'set map 0 1 6'
  'd maskout(weasdsfc,weasdsfc-10)'

And after a few moments, you have a map like the one shown at the beginning of the tutorial.

You can play around with your minimum values for ice and snow, I used 0.1 for sea-ice and 10 mm for SWE.  That seemed to produce a good looking map.  Experiment as you like!  This plot is really basic, no special formatting outside of the actual plot (titles, etc).  Hopefully you enjoyed this tutorial, I know it was oddly specific, and didn't really present anything new, hopefully you learned something anyway!

Download Example Script

Read More