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

 

Read More

Tutorial: Vertical Cross section with topography in pressure coordinates

Vertical cross sections are very easy to plot in GrADS, requiring only that you fix either the latitude or longitude value, and vary the height levels.  For example, an west to east cross section is easily set up using the following commands.

     'set lat 40.0'
     'set lon -125 -68'
     'set lev 1000 100'

A lot of cross sections however like to add in the topography, giving the cross section greater structure.  The method to do this in GrADS is not terribly hard or difficult, it is just a little tricky.

The first step is to fix your page area using the 'set parea' command.  This is because the in order to overlay topography onto your cross section, we need to change around the domain, but keep domains plotting in the same spot.  Setting the page area assures this.  Once your page are is set, just plot your desired cross section variable as usual.

Now it's time to overlay the topography.  In pressure coordinates, the topography can be represented by the surface pressure (surface pressure decreases with increased topographical height).  Since the surface pressure is a 2D and not a 3D variable, you have to fix the height in GrADS to the lowest model level (usually done by 'set lev 1000', or 'set z 1').  Now that you have done this, you need to fix the y-axis of the plot to match the y-axis of your cross section.  This is done with the 'set vrange' command.  The best way to do this is to set it exactly as you would have set the levs when plotting your cross section variable.

For example, if you typed in 'set lev 1000 100' for your cross section then you would type in 'set vrange 1000 100' for your topography.

This way, your topography is set to plot on the same scale as your cross section.  Now, we are ready to plot.  To do this, we make use of the 'set gxout linefill' option.  Then we simply fill the gap between the surface pressure and the baseline pressure of 1000.  The block of code to do this looks is shown below.

    'set lev 1000'
    'set vrange 1000 100'
    'set gxout linefill'
    'set lfcols 1 1'
    'd pressfc/100;const(pressfc/100,1000)'

The 'set lfcols' command is the syntax for setting the color of the topography.  The end result for an east west cross-section through the United States will look something like this:

Example of Cross section with topography

An additional way to plot topography is by using the 'set gxout bar' command, instead of the 'set gxout linefill' command.  The syntax is similar, but with minor differences.  You set the vrange the same, but instead of having to fill between two values, you now only need to plot the variable.  Visually, this gives your topography a more discrete, blocky look to it.

Syntax for plotting topography with 'set gxout bar':

     'set gxout bar'
     'set baropts filled'
     'set ccolor 1'
     'set lev 1000'
     'set vrange 1000 100'
     'd pressfc/100'

Resultant Image:

The real difficulty comes with rying to plot a vertical cross section on a log scale.  Since surface pressure is only a 2D variable, it is unaffected when setting the y-axis to be on a log scale.  As a result, neither of the above methods are capable of correctly plotting the topography on a log scale.  So in order to do this, we need to essentially trick GrADS into plotting the topography using a different method.

This is a little tricky to explain, but the basic idea is you simply plot the geopotential height at contour intervals so close together that they overlap and appear as one entity.  Since we are focused on topography, we subtract the surface geopotential height (elevation) from the geopotential height.  Also, since we are only interested in plotting the topography, any value >0 is masked out.  Hopefully, this short paragraph gives you some indication of what this next block of code is doing.

    'set gxout contour'
    'set clab off'
    size=750
    int=10
    i=0
   'hgt=hgtprs-hgtsfc'
  
    while(i*size>=-5000)
       j=i-1
      'set ccolor 1'
      'set cthick 10'
      'set cstyle 1'
      'set cmax 'i*size
      'set cmin 'j*size
      'set cint 'int
      'd hgt'
      i=i-1
   endwhile

Now that you can look at the code, I can explain it a little more thoroughly.  Basically, we need a while loop because GrADS cannot set enough contour levels at one time, so we have to continuously reset the contour levels until we cover the full range of values.  Within the loop we reset the line thickness and the line style to ensure continuity through the plotting.

The result of this code is a correctly plotted vertical cross section on a log scale with topography.

Vertical cross section (with terrain) plotted on a log scale

Hopefully this tutorial has provided you with a greater understanding of how to plot vertical cross sections in GrADS (at least in pressure coordinates).


Download this script

Read More