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