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