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

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

Tutorial: Fill between; how to shade area on line graphs

Sometimes when plotting data on simple line graph it is useful to use colors to fill between different curves.  For example, when plotting anomalies it may be helpful to fill the below zero values blue, and the above zero values red (see below figure).  This tutorial will show you how to use the gxout option 'linefill' to fill areas on your x/y graph. 

500mb height anomaly.  Negative values filled blue.

'set gxout linefill' requires two input values, and will fill the space between these value.  It has been used before in previous tutorials here to help you fill in topography for vertical cross sections (see: here for this example).  We will begin by setting up our file and location.  I've been interested in ensemble forecasting recently, so we'll work with the GFS ensemble data.






  'sdfopen http://nomads.ncep.noaa.gov:9090/dods/gens/gens20141109/gep_all_00z'
  'set display color white'
  'clear'
  'set mpdset hires'

  'set t 1 65'
  'set e 1'
  'set lev 500'
  'set lat 44';'set lon -72'

  'set vrange -30 30'







Notice how I only have one varying dimension (time) in this tutorial.  This will allow us to plot a basic line graph, showing the 500mb height anomaly in time.  Now we are ready to move on to plot the data.  We will be plotting the ensemble mean height anomaly.  >0 => Red <0 => Blue

  'set gxout linefill'
  'set lfcols 2 4'
 'd mean(gpa500mb/10.0,e=1,e=21);const(mean(gpa500mb/10.0,e=1,e=21),0,-a)' 

The above code first sets up, then makes the plot.  First use set 'gxout linefill' to get the graphics set up correctly.  The key to getting the different colors is using the command 'set lfcols 2 4'.  What this command does, is set up your plotting environment such that any value greater than the reference value is set to red(2) and anything less than the reference value is set to blue (4).  Once our colors are set, we plot the variable, in this case the ensemble mean 500mb height (divided by 10 to get units: Decameters).

When using the linefill graphics option , you need two inputs separated by a semi-colon.  These values indicate the values that will be filled between.  In this case, the first value is the height anomaly, and the second value uses the const function to set all height anomalies = 0.  So in this case, you are plotting the height anomaly relative to the zero line.  This code will give you the image above.

We can take this one step further.  The 2nd value in the 'linefill' display does not have to be a constant, it can vary.  So to stick with the theme of ensemble prediction, we'll plot the height anomaly again, only now we will fill the + and - 1 standard deviation around the mean.

GFS Ensemble mean Height anomaly and shaded Standard Deviation

The following will produce the above plot.  The method for producing the above plot is very similar to what we just discussed,  only now, instead of using the const function to set the 2nd value in the linefill display, we use two variable values.

   'set lfcols 15 0'

  'd mean(gpa500mb/10.0,e=1,e=21)+sqrt(ave(pow(ave(gpa500mb/10.0,e=1,e=21)-gpa500mb/10.0,2),e=1,e=21));mean(gpa500mb/10.0,e=1,e=21)-sqrt(ave(pow(ave(gpa500mb/10.0,e=1,e=21)-gpa500mb/10.0,2),e=1,e=21))'

Now we fill every value between the first (mean + 1std) and the second (mean - 1std) in gray.  Then to make the plot nicer, we throw on the ensemble mean, and we cap the standard deviation with black lines.



There is one last aspect with linefill, that I feel should be included in this tutorial, and that is how to fill the area under the curve with more than 1 color (see example below).  For this, we use similar code as above, only now we put it in a loop.  Now to do this, it works best if your variable is always greater than zero, due to the way the looping is performed.  So we are going to use total cloud cover as our variable.  Additionally, to save time, I'm only going to take a subset of times from the file, and only a subset of ensembles.

 

GFS Cloud Cover Multi Color

So to start, we set the range of levels that we want to shade in using bars, then we simply loop through all of the levels, use the linefill display, and simply reset the value in the const function to correspond to each level.

      levs='0 10 20 30 40 50 60 70 80 90'
      c=1
      while(subwrd(levs,c)!='')
        'set gxout linefill'
        'set lfcols 'c+1' 0'
        'd mean(tcdcclm,e=1,e=5);const(mean(tcdcclm,e=1,e=5),'subwrd(levs,c)',-a)'
        c=c+1
     endwhile


    'set gxout line'
    'set ccolor 3'
    'set cthick 8'
    'd mean(tcdcclm,e=1,e=5)'

The code above uses the subwrd command to grab the value from the "array" levs.  It breaks the loop if you run into the end of the array.  Then the code sets a new fill color, and voila, you have a multi-colored shaded plot.

So that's about all there is to the linefill function.  Hopefully after this tutorial you have a better understanding of how the linefill function works, and what some of it's capabilities are.  Below is a script with the various examples used here.

Download Example Script

Read More

Tutorial: Overlay model resolution onto weather map

There may come a time in your professional or student career when you will have to show a model grid overlaid on a map, to illustrate some sense of your model grid size.  This approximation can be done fairly easily in GrADS using a little bit of spherical geometry, a couple of while loops, and the 'q w2xy' command.

This tutorial will show you how to do this using the topography file from the Grads Data Server.
When you are finished with this tutorial, you will be able to reproduce the plot below.

40km grid overlaid on southwest US.

The first step in this process is to get that data.  One we have the data we set up the domain and do a little formatting.

   file='http://www.monsoondata.org:9090/dods/topo/rose/etopo05'
   'sdfopen 'file
   'set mpdset hires'
   'set display color white'
   'clear'

Now that we have our file open (and if you want, you can do a 'q file' to see the details), since we are plotting over an area bigger than the area we are drawing our model grid on, we need two sets of lat/lon boundaries.  The next bit of code provides the coordinates for the example plot, but you could of course choose your own boundaries.  I like to set these before hand, for easy changes later on.  In addition to your lat/lon boundaries, you need to choose a model resolution.  Once again, I just set this as a variable rather than "hard coding" it in, for quick access to changes.

  lat_bound='32 41.5'
  lon_bound='-115.5 -105'

  lat='32.5 40.25'
  lon='-115.0 -105.5'

   res=40000

Now that we have our initial conditions set, we can move on to plot the topography.  In this example, I use the script checkered.gs, available in the script library, to plot the nice boarder.  So I want to turn the x/y labeling off.  If you have read any of my other tutorials you'll already know this, but I use color.gs to set the color scale.  So if you have these two scripts in your script folder, you can follow along, with the code below.  If not, you can still follow the tutorial, just without the boarder and color scheme.

  'set xlab off';'set ylab off'
  'set gxout shaded'
  'set grads off'
  'set lat 'lat_bound
  'set lon 'lon_bound

  'color 50 5200 100 -kind green->tan->dimgray'
  'd rose'

The above code plots the topography from the topography file.  Remember, when setting the lat/lon boundaries to use your outer boundaries.  So that your inner box will fit inside.

This next part introduces the use of the 'q w2xy' command: World (w) to (2) x/y (xy).  Basically, you are translating lat/lon coordinates to page coordinates.  Before we plot the grid, we are going to first approximate the domain, and I use the term "approximate" because it won't be exact due to the fact that the physical distance between longitudinal coordinates is dependent on latitude, so the simple rectangle isn't an exact representation of the world coordinate system.  Anyway, we simply take the inner lat/lon boundaries and the the lower-left and upper-right corners using 'q w2xy' and then we draw a rectangle using the 'draw rec' command. 'Set line' is used to do a bit of formatting (e.g., the blue color in the example image).  A more detailed explanation of how to do this part is provided in this tutorial

  'q w2xy 'subwrd(lon,1)' 'subwrd(lat,1)
  xpos1=subwrd(result,3);ypos1=subwrd(result,6)

  'q w2xy 'subwrd(lon,2)' 'subwrd(lat,2)
  xpos2=subwrd(result,3);ypos2=subwrd(result,6)

  'set line 4 1 8'
  'draw rec 'xpos1' 'ypos1' 'xpos2' 'ypos2

Now we are ready to draw our model grid.  There are more than a few ways to go about doing this, and  the main thing is you want to get your starting and ending lat/lon points set.  Then we want to set up the constants that we will use in our spherical geometry calculations.  For all intents and purposes, we assume the distance per degree of latitude (mperlat) is constant. While it does have some variation due to the not quite spherical shape of the Earth, it's not a bad assumption that it's constant.

  latstart=subwrd(lat,1);latend=subwrd(lat,2)
  lonstart=subwrd(lon,1);lonend=subwrd(lon,2)

 mperlat=111200
 r_earth=6378000
 pi=3.141592

So, now we have our boundaries set, and our constants we can start to draw our grid.  So the basic idea is that we loop through all longitude and latitude points drawing boxes with sides roughly equal to our resolution.  There are several levels of complexity that can be used to approximate the distance per degree of longitude.  This tutorial will focus on the most simplistic way, however it will touch upon some of the more complex ways to go about it.

The simplest way to approximate your zonal grid resolution, is to just set one value for "delta_lon" for the whole grid, and apply that throughout.  The best way to do this, would be to take a mean latitude for your domain, and calculate delta_lon from that.  So setting delta_lon and delta_lat is reduced to:


   avg_lat=(latstart+latend)/2
   delta_lon=(res/r_earth)*(180/pi)/math_cos(avg_lat*pi/180)
   delta_lat=res/mperlat

This approximation is pretty good, so long as you are at low to mid-latitudes, and your domain is small. Usually when illustrating your model grid, you want to keep the domain small so that you can see the boxes better (unless your resolution is 200km or something like that, in which case you have bigger problems).  So this approximation isn't bad.  Now that we have delta_lat and delta_lon, we just loop through and draw boxes.  I like to start with my grid-box centered on my lower-left corner, so before entering the loop, I do a quick adjustment to my latstart, and lonstart.  Remember, to reset your starting longitude value within your latitude loop.

   'set line 2 1 1'        ;*Formatting
   latstart=latstart-delta_lat/2

   lat=latstart   
   while(lat < latend)

     lat1=lat
     lat2=lat1+delta_lat
     lon=lonstart-delta_lon/2

     while(lon < lonend)
       lon1=lon
       lon2=lon1+delta_lon


       'q w2xy 'lon1' 'lat1
        xpos1=subwrd(result,3);ypos1=subwrd(result,6)

       'q w2xy 'lon2' 'lat2
        xpos2=subwrd(result,3);ypos2=subwrd(result,6)

       'draw rec 'xpos1' 'ypos1' 'xpos2' 'ypos2
       lon=lon+delta_lon
    endwhile
    lat=lat+delta_lat
  endwhile


That's all there is too it!  You simply use a nested loop set up, to loop through all longitude and latitude values, drawing boxes first in the x-direction, and then in the y-direction.  Now, if you want the nice border around the plot, simply call the 'checkered.gs' script (again more information on this, can be found in the GrADS-aholic! script library).

  'checkered  -s 1.0 -w 0.1 -ss 0.09'

And, you will have your plot.

Now, as mentioned you can get more complex with setting the delta_lon values.  For example, you could re-calculate delta_lon as you loop through latitudes, taking the average latitude value to be the average of each grid-box.  Essentially, you recalculate delta_lon inside of the latitude loop, rather than use the same one for the whole domain.  If you really want to get fancy, you could calculate a delta_lon at the bottom of and the top of your box.  In this case, you are no longer able to use the 'draw rec' command, because you are dealing with 4 unique points, rather than a simple rectangle.  So you would need to use the 'draw line' command to draw your borders.  These methods are better used if you are near the poles, or using large domains, and they are simply trying to refine your approximation of distance between longitude points, and essentially account for the fact that your distance between degrees of longitude decreases as you move to high latitudes.

I hope you enjoyed this tutorial and found it worth while.  As always, a link to a script with this tutorial in it is provided.

Download Example Script

Read More

Tutorial: An introduction to using surface station data in GrADS

In addition to plotting gridded data, GrADS can also be used to plot station data; that is data formatted with information specified at specific latitude/longitude coordinate rather than data gridded in a mesh that spans a region. This tutorial will show you how to use some of the functions in GrADS to work with surface data from both the mesonet and METAR inventories to generate surface maps with and without contoured data. At the end of the tutorial, you will know how to make the images shown below using station data. 

Surface map with METAR data.

Surface data and Contoured Sea level Pressure From METAR data

Surface data and Contoured Sea level Pressure From METAR data

First things first: Getting the station data.  Getting the station data and, more importantly, getting it into the right format is the trickiest part of using station data with GrADS.  To get the data into the right format, you must have knowledge of coding in another language (Fortran, C, Python, etc,.).  The scope of this tutorial is to focus on using the station data in GrADS, and not showing you how to get the data into the right format.  So if you are looking to learn how to format station data such that it can be read into GrADS, this tutorial is not for you. 

Some good information regarding how to format station data can be found here and here.  The basic idea is that you need to take the station data, write out a binary file with a header including information such as station id, latitude and longitude, and then the actual data for each station.  Afterwords you need to write out a control file to tell GrADS how to interpret your binary file.

Since I am not going to show you how to grid station data, I have prepared a few files for you to use with this tutorial.  They can be found here:

Mesonet Data (for temperature map)

• mesotest.ctl
• mesotest.dat 

 METAR Data (For station and SLP maps)

• metartest.ctl
• metartest.dat
  

Once you have downloaded these files and placed them in your GrADS directory (note: Mine is c:/OpenGrADS/, you may need to change the 'dset' variable in the .ctl files to match your GrADS directory).  Now before you can open the file, you need to use the 'stnmap' utility to help map out the data.  This is easy enough to do.  Before you open GrADS, open up a console (windows command prompt) and type the following command to make your station map.

    'stnmap -i mesotest.ctl'     ---> For the mesonet data.

    'stnmap -i metartest.ctl'    ---> For the METAR data.

This makes your station data accessible in GrADS.  For more information on the station map utility you can read more here.

Now you are ready to open up GrADS!  We'll focus on the station data only first.  To do this start by opening up metertest.ctl, and do a 'q file' to see whats in the data.

    ''open metartest.ctl'

    'q file'

This will tell you what variables are in the file.  In the data file provided, cld and wxsym are dummy variables and contain no relevant information;  they are just place holders.  What we are going to do is plot out station plots with temperature and dew point.  This is very simple, all that is needed is to first 'set gxout model' and then display the data.  To commands below will reproduce the first image seen (without the title).

    'set display color white'
    'set mpdset hires'
    'clear'
    'set lat 30.52 40.24'
    'set lon -100.39 -90.11'
    'set gxout model'
    'd u;v;ts;td'

The 'd u;v;ts;td' tells GrADS to plot a station model with wind components u and v, along with the surface temperature and dew point.  Note: without u;v components nothing will plot when the graphics is set to 'model'.  To add more, you simply increase the number of variables in the chain (in proper order).  For example to add SLP the command would be: 'd u;v;ts;td;slp'.  Missing variables are represented by a 0.0.  The expected order of variables is well described here.

Now lets say we want to do a similar station plot, only now with the contoured sea-level pressure.  To do this, you need to fit your station data to a preexisting grid.  In this tutorial, we will use the RAP model as our base grid.  It doesn't matter which data set you use, but it is recommended you pick something that is of similar grid spacing to your station plot (e.g., it wouldn't be good to use grid METAR data to 2km grid spacing).  Since we are using METAR data, the commands will be almost the same as the ones above with the only difference being a domain change.

    ''open metartest.ctl'

    'set display color white'
    'set mpdset hires'
    'clear'
    'set lat 36 49'
    'set lon -84 -65'
    'set gxout model'
    'd u;v;ts;td'

Once this is done, we need to open up the RAP model data using the 'sdfopen' command.

   'sdfopen http://nomads.ncep.noaa.gov:9090/dods/rap/rap20140216/rap32_12z'

This file is now open as file 2.  From here it is simple, we just use the oacres function to fit the data to the grid.  Now here, it does not matter which variable from the RAP you use, since it is just a dummy variable to get the grid to which the station data will be interpolated too using the Cressman objective analysis on. 

First some formatting...

    'set cthick 6'

    'set clevs 990 994 998 1002 1006 1010 1014 1018 1022 1026 1030'

    

Then we can use the oacres function to interpolate the station data.

    'd oacres(rh2m.2,slp.1,45,30,20,10,5)*33.86'

And voila, you have your station SLP contoured on your map.  To break down what exactly is going on in the above command:

The first argument (rh2m.2) is the variable from the RAP model, essentially the grid you are projecting your data onto.  The suffix .2 is simply there to specify that you want the variable from file number 2 (if you are unfamiliar with how to handle multiple files in GrADS, check out this tutorial).  The second argument (slp.1) is the station variable you want interpolated to the grid, in this case Sea Level Pressure.  The arguments after that specify different radii of influence that you are using to interpolate the station data with.  This will take a little experimentation on your part to get right, but the selection used above seemed to work fairly well with this data set.  I am personally not 100% familiar with the exact math behind the Cressman OA, so if you want to know more on that, you are on your own.  The factor of 33.86 is to convert the station SLP reported in inHg to mb.

The final lesson in this tutorial is very similar to the 2nd lesson except now, we are using mesonet data, which has a higher density than METAR data.  So again we will start by opening the file (in this case the mesotest.ctl).

    'open mesotest.ctl'

    'set mpdset hires'
    'clear'
    'set lat 36 49'
    'set lon -84 -65'

Instead of displaying the data as a station model, we are doing a shaded contour underneath the station data.  So we need to do the Cressman OA first, which means we need to open the RAP file and set up our formatting.  In this example, I used color.gs and xcbar.gs to set the color scale and color bar, but you can do that however you would like.

  'sdfopen http://nomads.ncep.noaa.gov:9090/dods/rap/rap20140216/rap32_12z'
  'color -20 100 2 -kind fuchsia->darkblue->lime->yellow->red'
  'set gxout shaded'
  'd oacres(rh2m.2,ts.1,45,30,20,10,5)'
  'xcbar 10.10 10.20 0.5 8.0 -fs 5 -line on'

This will draw your shaded station temperature data fit to the RAP grid.  Once again, you may have to play with the radii numbers to get it exactly to your liking.  Same goes for your color location and your color scale.  Once this is done, simply use the 'set gxout value' command and display your temperature data.

  'set gxout value'
  'd ts'

And that about does it for your temperature plot, and for that matter, this tutorial!  I hope you enjoyed it and found it interesting.  The script below has the examples described here and is set for you to use and experiment with.  

As a final note:  The station ids were not plotted with the data in the examples in this tutorial.  However, this can be toggled in GrADS by using the 'set stnid' command.  For example, to turn on the plotting of the station ids: 'set stnid on' .   

Thanks again for reading!

Download Example Script

Read More