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