Tutorial: Draw arbitrary cross-section with terrain in pressure coordinates using GrADS

As mentioned in the introductory post, one of the major deficiencies with GrADS, is that drawing arbitrary vertical cross-sections is not an easy task.  There are however work-arounds to do this, but they can be quite complex, especially when including terrain.  This tutorial will guide you step by step through the process of plotting an arbitrary cross section in GrADS with terrain in pressure coordinates.

Essentially, to plot arbitrary vertical cross sections in GrADS, you have to collect your (z-varying) data at each lat/lon point in your cross section, and then re-grid that data before plotting it.  The g2stn(), and the coll2gr() functions in GrADS help you along with doing this.  For this tutorial, I am going to start by borrowing the code from this link and change around a couple of variables.  In this section will provide a basic walk through of what this code is doing.  Once we have gone though this code, I will show you my work around to plotting terrain, which is similar to some of the other work-arounds out there, but I think mine works a little faster.


     'sdfopen http://nomads.ncep.noaa.gov:9090/dods/gfs_hd/gfs_hd20130512/gfs_hd_00z'
     'set x 1'
     'set y 1'
     'set xlog on'
     'set lev 1000 100'
     lon1 = -95.0
     lon2 = -90.0
     lat1 = 55.0
     lat2 = 15.0
     lon = lon1
     'collect 1 free'
     while (lon <= lon2)
       lat = lat1 + (lat2-lat1)*(lon-lon1) / (lon2-lon1)
       'collect 1 gr2stn(rhprs,'lon','lat')'
       lon = lon + 1
     endwhile

     'set lon 'lon1' 'lon2   
     'set clab on'
     'set gxout shaded'
     'color 50 100 2 -kind black->lightgreen->darkgreen->darkblue'

*If you don't have color.gs:
*         'set clevs 40 50 60 70 80 90 100'
*         'set ccols 0 3 4 5 6 7 9 10'
*--------------------------

     'd coll2gr(1,-u)'


So, this code is pretty similar to the example shown in the link provided above, but with a few minor differences, for example we are plotting relative humidity, instead of pv.  Also, I used the world dimension to set the x-axis for the arbitrary cross-section.  This is fairly meaningless, as it will still just plot the values that you have there, regardless of the dimensions.  Lastly, I used color.gs to do the RH color scale.

Lets break down this code a little bit, so that you have an idea as to what you are doing with it.  The first thing you see below the plot scaling is the declaration of your longitude and latitude boundaries.  To make arbitrary cross sections, you must declare your boundaries.  Once the boundaries are declared, you see the 'collect 1 free' command.  This basically, tells GrADS that you are emptying out the grid-collection numbered "1".  So, if there is anything in there, you are freeing this variable.  If you wanted to plot more than one variable, you would put another 'collect' command for a different number, e.g., 'collect 2 free'.  Assuming, you are only plotting relative humidity, the loop will run through (slowly if I may add) collecting the z-varying variable at each lat/lon point.  You can see the latitude is really just a linear interpolation between your latitude boundaries.  Once you have collected all of your data, the dimensions are set on the x-axis (longitude), and the variable is plotted at all vertical levels using the coll2gr() function.  And that is about as simple as this gets, you can see the resulting image below.

Arbitrary cross-section varying in both lat/lon.

As you can see, it's not too much to look at, I did the plot over a small area so to reduce computational time, but this is basically what you get if you use the method described in the above link.

Now, how about adding terrain.  This is tricky business because terrain (or surface pressure in pressure coordinates) is not a z-varying variable.  So, unfortunately, it is not as simple as adding terrain into another collection variable.  Once again, similar to plotting terrain on a log scale, we need to pull some tricks!  Basically, what we are are going to do, is scale the environment, and use the 'q w2xy' command and generate a bunch of vertices to use in a polygon fill.  Then to draw the terrain, we simply draw our polygon.

So, here is the code for drawing a cross section with arbitrary terrain.

    'sdfopen http://nomads.ncep.noaa.gov:9090/dods/gfs_hd/gfs_hd20130512/gfs_hd_00z'

    'clear'
    'set grads off'
    'set zlog on'
    lon1 = -115.0
    lon2 = -95.0
    lat1 = 40.0
    lat2 = 35.0
    lon = lon1

    terrain_arr=''

    'set lat 'lat1
    'set lon 'lon1' 'lon2
    'set lev 1000 100'
    'set gxout shaded'
    'd rhprs'
    'draw title Temporary Scaling Environment'

    'collect 1 free'
    'collect 2 free'
    'set x 1'
    'set y 1'

    say 'collecting station data'

    'q gxinfo'
    xline=sublin(result,3)
    yline=sublin(result,4)
    xs=subwrd(xline,4)' 'subwrd(xline,6)
    ys=subwrd(yline,4)' 'subwrd(yline,6)


    while (lon <= lon2)
       'set lev 1000 100'
       lat = lat1 + (lat2-lat1)*(lon-lon1) / (lon2-lon1)  
       'collect 1 gr2stn(rhprs,'lon','lat')'

       'set gxout print'
       'set lat 'lat
       'set lon 'lon
       'set lev 1000'
       'd pressfc/100'
       press=sublin(result,2);press=subwrd(press,1)
       'q w2xy 'lon' 'press
       x=subwrd(result,3)
       y=subwrd(result,6)

       if(x<=subwrd(xs,1));x=subwrd(xs,1);endif
       if(x>=subwrd(xs,2));x=subwrd(xs,2);endif
       if(y<=subwrd(ys,1));y=subwrd(ys,1);endif
       if(y>=subwrd(ys,2));y=subwrd(ys,2);endif
       terrain_arr=terrain_arr''x' 'y' '

       lon = lon + 1
     endwhile

    say 'plotting data'
    'clear'
    'set lev 1000 100'
    'set lon 'lon1' 'lon2
    'set clab on'
    'set gxout shaded'
    'color 50 100 2 -kind black->lightgreen->darkgreen->darkblue'
    'd coll2gr(1,-u)'

    say 'drawing terrain'

     x1=subwrd(terrain_arr,1)
     y1=subwrd(terrain_arr,2)
    terrain_arr=terrain_arr''subwrd(xs,2)' 'subwrd(ys,1)' 'subwrd(xs,1)' 'subwrd(ys,1)' 'x1' 'y1

    'set line 15 1 1'
    'draw polyf 'terrain_arr


You can see the similarities to the code above, the variable is still collected and regridded using the g2stn() and the coll2gr() functions.  But now there is a lot more code relating to the drawing of the terrain.  Lets break this down.  The first extra command, is the one at the top that sets the variable terrain_arr=''.  This is going to be our polygon array.  We are going to populate this with a bunch of vertices corresponding to the terrain.  But first, we need to clear it.  Once we have cleared that array, we plot a "Dummy" cross-section to scale our environment.  Remember, it doesn't matter what our latitude is, because in the loop, we set the latitude for each variable, and our world coordinates vary in longitude and height.  But we can plot our dummy variable as such.  I also gave it a title, so to specify that it is just a temporary dummy plot.

      'set lat 'lat1
      'set lon 'lon1' 'lon2
      'set lev 1000 100'
      'set gxout shaded'
      'd rhprs'
      'draw title Temporary Scaling Environment'


Now, once we have our environment scaled, we move into our loop.  Here, we set our variables the same way, but then we need to get our terrain value converted from a world coordinate to a page coordinate.  We do that by setting the gxout to be print, and then setting the specific lat/lon coordinate.  Then we simply print out the surface pressure at this point, and use the 'q w2xy' command to convert that into page coordinates.  We then check to make sure the page coordinates for each vertex in our polygon falls within the plot area.  Then we simply build the array from there.

       'set gxout print'
       'set lat 'lat
       'set lon 'lon
       'set lev 1000'
       'd pressfc/100'
       press=sublin(result,2);press=subwrd(press,1)
       'q w2xy 'lon' 'press
       x=subwrd(result,3)
       y=subwrd(result,6)

       if(x<=subwrd(xs,1));x=subwrd(xs,1);endif
       if(x>=subwrd(xs,2));x=subwrd(xs,2);endif
       if(y<=subwrd(ys,1));y=subwrd(ys,1);endif
       if(y>=subwrd(ys,2));y=subwrd(ys,2);endif
       terrain_arr=terrain_arr''x' 'y' '

At the end of the loop, you should have a big long array filled, with x/y coordinates.  But we are not quite done.  We will need to close out our polygon.  So the last thing we do before drawing our polygon is add a couple more vertex values to it, so it runs back down and closes it off.  Then we simply draw the polygon. 


     x1=subwrd(terrain_arr,1)
     y1=subwrd(terrain_arr,2)
    terrain_arr=terrain_arr''subwrd(xs,2)' 'subwrd(ys,1)' 'subwrd(xs,1)' 'subwrd(ys,1)' 'x1' 'y1
    'set line 15 1 1'
    'draw polyf 'terrain_arr

And, voila! Your terrain is drawn on the map.  What's even more fantastic about this, is that you do not even have to pull any special strings to plot this in log coordinates.  This method works both ways.

To prove this works as intended, I am showing you two images below, both are the same E/W cross section with terrain, however one was made using the arbitrary cross-section method, the other, the simple GrADS command.  Also, for simplicity, I plotted it on a non-log scale, so to make the simple GrADS command easier.

Normal Method

Arbitrary Method

So, as you can gather, these plots are identical, and that is good, because if they weren't than the arbitrary method would not be correctly plotting the terrain.

A few final notes before wrapping up: 

•  The main issue with the plotting arbitrary cross-sections, is that it is slow if you are accessing the data online through, for example, the GrADS data server.  The plot above took ~5-10 minutes, this time will increase the more variables that you add to your plot.

• You need to be aware that the lon=lon+1 command can lower the resolution of your plot if your x grid-spacing is less than 1 degree of longitude.  For example, the above plot was made using the 0.5 degree GFS data so to ensure that my plots matchedin up, my longitude was increased by increments of 0.5: lon=lon+0.5.  You could add a simple 'q file' command to your script to get the grid-spacing and always set your longitude to increment that way, but bear in mind this can really increase your computational time, as you are collecting more points.

• There is more than one way to skin a cat: I want to note there are other methods to plot terrain on your arbitrary cross-sections, but I have found them to be slower than this one, as they require you loop through your data points more than once.  This method just requires one loop.  But other methods are perfectly acceptable as well, and will probably give you similar or the same results.
   

So that's it, this was a fairly long tutorial, but hopefully you have a better understanding how  to plot arbitrary cross-sections in GrADS.  Especially when it comes to including terrain in pressure coordinates.  I imagine you could probably carry this method over to height coordinates easily enough though.

Download Example Script