Tutorial: Use GrADS to draw shapefiles

Here we go with yet another tutorial on using shapefiles in GrADS.  This one is slightly more advanced than the other tutorials I have on this site: here and here.  The reason it is a little more advanced is that, in this tutorial we will now use GrADS to create shapefiles as well as read them.

The ability to create shapefiles can be quite useful, especially if you are would like to output weather data into files that you can easily read into any type of GIS application.  Another possible application (one visited later on here) is the creation of a stippling effect on your GrADS plot.  In this tutorial, we will look at the 500mb geopotential height from the GFS, and make shapefiles from this data.

Now, there are four commands that will be primarily used in setting the shapefile options and drawing our shapefile.

1. 'set shp'
2. 'set shpattr'
3. 'set shpopts'
4. 'set gxout shp'

So, lets start out simple.  The first thing we will do is open up the file, and draw the 500mb geopotential height.

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

    'set lev 500'
    'set lat 20 50'
    'set lon -125 -65'
    'set gxout shaded'
    'd hgtprs'

So, now that our variable is plotted, we have an idea what our shapefile will look like.  Now we get into the shapefile stuff.  The first command we will use is the 'set shp' command.  This command sets up our shapefile type, and our file output.

    filename='shapeout'
    'set shp -ln 'filename

The reason I gave the filename a variable, is because we will call it later on when we open the newly made shapefile, so I made this simple to carry over it by saving it to a variable.  Be sure not to include a file suffix to the filename, as when you draw the shape it makes up four files; .dbf, .shp, .shx, .prj.

The dashed option tells GrADS what kind of shapefile to output, the three available options are "line (-ln), point (-pt) and poly (-poly)."  I do want to point out that I have personally had trouble with using -poly on my windows machine.  You can always check your shape options by doing a 'q shpopts' command.

Once you have these options set, we move onto the 'set shpattr' and the 'set gxout shp' commands.

    'set shpattr Author string GeoHeight500'
    'set gxout shp'
    'd hgtprs'

Now, we have set the shape attributes for the DBF file, set the data output to shapefile, and created the shapefile.  The 'set shpattr' command helps provide a little meta data for the DBF file.  So, now when you do a 'q dbf' on the file, you will see your shapefile labeled with the correct meta data.  Following our example, a 'q dbf' would bring up the following result.

    RECORD#,CREATED_BY,Author,CNTR_VALUE
    0,GrADS-2.0.a9.oga.1,GeoHeight500,5450.000000
    1,GrADS-2.0.a9.oga.1,GeoHeight500,5500.000000
    2,GrADS-2.0.a9.oga.1,GeoHeight500,5550.000000
    3,GrADS-2.0.a9.oga.1,GeoHeight500,5550.000000
    4,GrADS-2.0.a9.oga.1,GeoHeight500,5600.000000
    5,GrADS-2.0.a9.oga.1,GeoHeight500,5600.000000
    6,GrADS-2.0.a9.oga.1,GeoHeight500,5650.000000
    7,GrADS-2.0.a9.oga.1,GeoHeight500,5650.000000
    8,GrADS-2.0.a9.oga.1,GeoHeight500,5700.000000
    9,GrADS-2.0.a9.oga.1,GeoHeight500,5700.000000
    10,GrADS-2.0.a9.oga.1,GeoHeight500,5750.000000
    11,GrADS-2.0.a9.oga.1,GeoHeight500,5800.000000
    12,GrADS-2.0.a9.oga.1,GeoHeight500,5850.000000
    13,GrADS-2.0.a9.oga.1,GeoHeight500,5850.000000

What this shows is the values in which shapes were drawn.  Basically, just lines following the geopotential height values on the far right.  Now, if we want to draw these shapes, we simple use teh 'draw shp' command.

    'draw shp 'filename

And ... voila! Your shapefile is drawn on the map.  The resulting image is below.  Now, this might not look to fancy, and really only look as through you contoured your values over your already shaded values.  But remember this is the most basic application of this process.

500mb Geopotential height from the GFS, shapefile is outlined in white

Now, that you have a basic idea how to draw shapefiles using GrADS, lets move on to a slightly more complex example.

Using GrADS shapefile output to create a stippling effect:

Now, as mentioned above, one of the options with the 'set shp' command is the -pt option, which tells GrADS to output your shapes as points, as opposed to lines or polygons.  So if we do that, the shapefile output is much different.  So if we follow the exact same commands, only replacing the -ln with a -pt, we get a completely different looking result.

    'set shp -pt 'filename
    'set shpattr Author string GeoHeight500'
    'set shpopts 15'
    'set gxout shp'
    'd hgtprs'
    'draw shp 'filename

Shapefile output using the -pt option

We can take advantage of this output to create a stippling effect, all we need to do is use the 'q dbf' command and set a range of values we want to stipple.  Now, I won't list the full result here because there are upwards of 7000 values, but here is a small sample.


    RECORD#,CREATED_BY,Author,LONGITUDE,LATITUDE,GRID_VALUE
    0,GrADS-2.0.a9.oga.1,GeoHeight500,-125.000000,20.000000,5856.270020
    1,GrADS-2.0.a9.oga.1,GeoHeight500,-124.500000,20.000000,5854.689941
    2,GrADS-2.0.a9.oga.1,GeoHeight500,-124.000000,20.000000,5853.449707
    3,GrADS-2.0.a9.oga.1,GeoHeight500,-123.500000,20.000000,5852.479980


Now, if you went through this tutorial (also listed at the beginning of this post), you can probably guess what we are going to do.  We need to run a loop through each shape, and grab the value of geopotential height on the far right.  This means we need a block of code to separate the data by comma.  Once you have your value, you basically just run it through a quick conditional statement to see if it falls between some range (5600 and 5750 in this example). Below is a brick of code that does all of this.

    check=1
    div=2
    while(check=1)
    'q dbf 'filename'.dbf'
      line=sublin(result,div)
      com_count=1
      vars=1
      if(line='');check=0;else
         while(vars<=6)
           com_check=1
           len=1
           c=com_count
           while(com_check=1)
             comma=substr(line,c,1)
             if(comma=',');var.vars=substr(line,com_count,len-1);com_check=0;endif
             if(comma='' & vars=6);var.vars=substr(line,com_count,len-1);com_check=0;endif
             len=len+1
             c=c+1
           endwhile
           com_count=com_count+len-1
           vars=vars+1
         endwhile
        val=var.6
        say div' 'val
      endif
        if(val>=5600 & val <=5750)
          'draw shp 'filename' 'div-2
        endif
      div=div+1
    endwhile


Note the value is set to the 6th value on each line, corresponding to the geopotential height value, and also notice that the 'draw shp' command points to div-2, corresponding to the shape ID on the far right, the number before the first comma in the 'q dbf'.  The resulting image is shown here:

Stippling effect using GrADS shapefile output

You should be able to get the idea on this here, you can even get more complicated by using different colors with different ranges using the 'set line' command, but I'll leave that to you.  One last note, with 7000 shapes to loop though, this stippling can take a while.

I do believe that is going to do it for this tutorial.  Hopefully, it provided you with some basic knowledge on how to use the shapefile output abilities of GrADS.  There are certainly a lot more applications to this than the ones I provided here, and those applications I will leave for you to experiment with.

In the example script provided, I put in a brief conditional statement to determine whether or not to use the stippling effect.

Download Example Script

Read More

Script: theta.gs; Defines potential and equlivalent potential temperature given temperature and moisture inputs

This script, is really more of a function, the best way to use this would be to copy the contents of this script into the bottom of the script you wish to use it in.  This function reads in temperature and moisture information and returns a good approximation of both potential and equivalent potential temperature.

This function asks for up to four inputs, though only temperature and moisture inputs are required.
It is capable of calculating  equivalent potential temperature given four different moisture input variables:

1. Relative Humidity (-rh)
2. Mixing Ratio (-r)
3. Vapor Pressure (-e)
4. Specific Humidity (-q)

The default is relative humidity.

Example Usage: theta(tmpprs,rhprs)
                           -This example would define PT and EPT using relative humidity as the input.

In addition to calculating potential and equivalent potential temperature at given pressure levels, this function also allows you to specify surface pressure as the pressure input.  This adjusts for the fact that it would simply calculate everything at fixed pressure levels.

Example Usage: theta(tmp2m,spfhum2m,pressfc,'-q')
                           -This returns the surface PT and EPT calculated using specific humidity as the input.

Notes on Usage: As always, all input variables must be in MKS units.  If you wish to use a method other than relative humidity, but do not wish to use a specific variable for surface pressure, just type 'lev' or '-' in place of where surface pressure would normally go (e.g., theta(tmp,mixr,'-','-r'))  Be sure to put '' around any variables with a "-" in them, or GrADS will try to subtract.

I have not tested this function with all of the above variables,but this function does work well with both relative humidity and specific humidity.  As always, report any bugs here! 

*Minor bug fixed 09/2013

Download Theta.gs

Read More

Some tips on running GrADS from the Windows Command Line and .bat scripts

If you are using GrADS through a linux operating system, chances are you are familiar with running GrADS from the command line, so this post is mostly aimed at the Windows users.  This post will go through the basic instructions for running GrADS from your command line along with arguments and options.  So lets first begin with a question.

Why would you want to run GrADS from the command line if your mouse is perfectly capable of opening GrADS via the "double click?"

The best answer that I can give you is that running GrADS from the command line vs. double clicking gives you the option to include arguments.  Furthermore, the command syntax used is the same used in .bat scripts.

Now, there are a number of different options available that you can include on the command line when you open the program.  The most important is the -help option.  So if you type:

    'grads -help'

Into the command line, you will get a whole list of options you can include when you run the program.

 The options I have found to have been the most useful are these:

• -l/p: Choosing to open GrADS in landscape or portrait mode
• -b: running GrADS as a batch process (program is hidden from view)
• -c: Allows you to execute a command (likely a script)
• -x: Exits GrADS automatically after executing your command

There are other options that you can see with the -help command, but these are the ones I have found most useful.

Example: grads -lcx 'run script.gs'

This example would open GrADS in landscape mode and run the file "script.gs" and then exit after it was done.

Now, you can actually carry this command line method over to a Windows batch (.bat) script.  Basically, a batch script is similar to a GrADS script, in that it is a file that can execute a bunch of Windows commands at once.  Now, you can include a command to run GrADS with options within your batch script.

Now, this can be very useful as this way you can actually pass arguments from your windows script into your GrADS script.  For example, if you wanted to set your lat/lon boundaries in your windows script, or your display color, you can do that.  All you need to do, is set these variables inside your batch script, and call GrADS with the arguments attached.

Batch Script:
      set display="white"
      set minlat="20"
      set maxlat="50"
      set minlon="-125"
      set maxlon="-65"

     grads -lc 'run script.gs '%display%' '%minlat%' '%maxlat%' '%minlon%' '%maxlon%

GrADS Script:     function script(args)
     display=subwrd(args,1)
     minlat=subwrd(args,2)
     maxlat=subwrd(args,3)
     minlon=subwrd(args,4)
     maxlon=subwrd(args,5)

Now, what you see in the GrADS script, that may be new to you, is the function declaration at the top.  This is required if you want your script to read in the arguments that you laid out in the batch script.  All this line does, is essentially tell GrADS that you are looking for arguments.  You also see in the batch script you have spaces between each argument denoted by the ' '.  This allows you to separate your arguments in GrADS using the subwrd() function.

But the ability to pass arguments from the windows script into the GrADS script can be very helpful, especially if you are running scripts that need variable input (e.g., date or time).  As this method could allow you to run GrADS without having to constantly change the script around.

In any case, I just wanted to provide a few tips on running GrADS from the command line in Windows.  I didn't go too in depth into windows scripting (I leave that to brighter minds), but hopefully you have an idea as to how to run GrADS in a batch script as well as how to pass arguments from one to the other.

Read More

Script: checkered.gs; Draws a black/white checkered border around the outside of your map.

This script helps you create more professional looking area maps by plotting a clean border around your map.  An example is shown here:

World map plotted with checkered border.

Required Arguments: There are no required arguments for checkered.gs

Options: -help(-h)   : Pulls up Help Page
                -size(-s)    : Toggles size/number of bars bordering the plot
                -width(-w): Toggles the width of the border
                -strsiz(-ss) : Toggles the size of the lat/lon coords
                -laboff(-lo): Turns off the drawing of Lat/lon Labels

Example Usage: 'checkered  -s 1.8 -w 0.1'


Note that checkered.gs is currently only compatible with map projections of scaled and latlon.  Also, you must have something plotted in order to have a border drawn.  Please report any bugs so that future versions can be made better.

Download checkered.gs

Read More

Script: us_map.gs; Draws and saves map of the United States with Alaska and Hawaii in the lower left corner

This script gives you the ability to generate the classic US map that has helped convince countless of Americans that Alaska is really an island near Mexico, and not a giant landmass far to the north.


Required Arguments: -var: chosen plot variable

Options:   -help(-h):  Pulls up Help Page
                  -base(-b): Uses Basemap options, you must follow this option with an "O" or "L"
                  -iname   : Sets name of output image, you must follow this with an image name
                  -xsiz      : Toggles the number of pixels in x
                  -ysiz      : Toggles the number of pixels in y
                  -bcol     : Toggles the color of the basemap'

Example Usage: 'us_map rose -b O -iname topomap.png'
                             -Plots the topography with masked out ocean and saves to the image: topomap.png

US Topomap made using us_map.gs

Note, when using us_map.gs, the lat/lon boundaries are set within the script, so you cannot set these before hand.  Also, be sure to set your other plot options before using us_map.gs (e.g., page area, contour levels, etc.).  Lastly, us_map.gs turns off the x and y labels to plot the boxes for Alaska and Hawaii, so you may need to turn them back on after using us_map.gs.  Also, please report any bugs so that this script may be fixed in future versions.

Download us_map.gs

Read More

Tutorial: Shapefile; A look at the DBF file.

This is the 2nd tutorial on using shapefiles in GrADS.  If this is your first time using shapefiles in GrADS, it is recommended you check out this tutorial.  This tutorial will dig a little deeper and look at the 'q dbf' command to look at the different shapefile properties.  For this tutorial, we will take a break from the usual focus on weather data, and instead focus on politics!  What we will do is draw the United States colored by which way each state voted in the 2013 presidential election.

Blue: Democrat, Red: Republican

So, to start you first need to download the required shapefile: Shapefile
Once you have this file put in your shapefile folder, we can begin.  Now we do need to open a file to get the environmental scaling, but that's all we need it for.  I'm going to use a simple GFS file from the GrADS Data Server.

    'sdfopen http://nomads.ncep.noaa.gov:9090/dods/gfs_hd/gfs_hd20130512/gfs_hd_00z'
    'set lat 22 50'
    'set lon -125 -68'
    'set clevs 10000000'
    'd pressfc/100'

Now, we have our file open and our map set up.  Now we are ready to plot our shapefiles.  However, we need to specify the color of each state.  This is where the 'q dfb' command.  The 'q dbf' takes a look a the database (.dbf) file that came in the shapefile zip folder, which is essentially the metadata.  So for this example, we would query the database file of the state shapefile.

    'q dbf Shapefiles/s_06se12.dbf'
   say result

What you will see is a list of information separated by commas.  What this information tells you, is the state abbreviation, shapefile id number, lat/lon, state name.  So, all we need to do, is extract this information and compare it to our previous knowledge of who won the election.  The hardest part of this is dealing with the fact that each value is separated by a comma instead of a space; we can't just use the 'subwrd()' command.  So, since we know that arrays can be set with strings as well as numbers, we can simply set a list of different colors, corresponding to each state abbreviation.

The first few numbers might look like:
    *2=Republican
    *4=Democrat

    col.AL=2
    col.AR=2
    col.AZ=2
    col.CA=4
    ....
    ....


Once, you have your list set, we simply need to loop through each shape within the shapefile, and determine its state abbreviation, then set the shapefile fill color to the color corresponding to the list above.

    i=2
    check=1
    while(check=1)
      'q dbf Shapefiles/s_06se12.dbf'
      line=sublin(result,i)
      if(line='' | line='(NULL)');check=0;else;

        com_count=1
        vars=1
        while(vars<=3)
          com_check=1
          len=1
          c=com_count
          while(com_check=1)
            comma=substr(line,c,1)
            if(comma=',');var.vars=substr(line,com_count,len-1);com_check=0;endif
            if(comma='' & vars=11);var.vars=substr(line,com_count,len-1);com_check=0;endif
            len=len+1
            c=c+1
         endwhile
         com_count=com_count+len-1
         vars=vars+1
       endwhile
       state=subwrd(var.3,1)
       'set shpopts 'col.state
       'draw shp Shapefiles/s_06se12.shp 'i-2
      endif
     i=i+1
   endwhile

This block of code may look complicated, but mostly it is just looking for commas within the string, and setting variables to the values in-between each comma.  You will also see the option "i-2" added on at the end of the 'drawshp' command.  This is telling GrADS to only plot the shapefile that has the id number i-2.  Since I starts at 2 (the first line of data in the shapefile), and since the first ID number is 0, we need to offset the pointer by 2.  So what this loop does, is loop through all of the shapes in the shapefile, and set the color to the list prescribed at the beginning of the script. 

If you run this script in GrADS, you will see a the shapes start to show up colored on the map.  So, that does it for this tutorial, hopefully it helped you work a little bit with the 'q dbf' command. 

Download Example Script Here

Read More

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

Script: Dynamic: Calculates temperature and vorticity advection, Q-vectors, deformation, geostrophic and ageostrophic winds, and Fronogenesis in GrADS

Usually, more advanced dynamic variables like frontogenesis and deformation are not included in the list of variables in your model output files.  Instead, it is up to the user to calculate them from the available wind and temperature fields.  Calculating these variables is complex and often requires a lot of steps using the various centered difference functions in GrADS.

Dynamic.gs calculates the following variables from the geopotential height, temperature, and horizontal wind.

Variable                           Name             Units                
-Geostrophic Wind     :    ug,vg             [m/s]   
-Ageostrophic Wind   :    ua,va             [m/s] 
-Q-Vectors                :    Q1,Q1           [pa/m2/s]   
-Temp Advection       :    tadv               [K/s]        
-Vort Advection         :    vadv              [-]         
-Frontogenesis         :    F                   [K/m/s]x10^9
-Fn Vector                :    fnx,fny         [K/m/s]x10^9
-Deformation             :    def1,def2     [m]         

The variables are then saved into the name, and you can plot them as normal variables.  However, since there are a lot of calculations to make in this script, the variables are only calculated for the current dimensions.  To calculate these values at all x,y,z,t points, simply set your dimensions to span the entire 4D domain.

Options:   -help: pulls up help page
                  -var : allows the user to point to specified height, temperature, wind variables
                            Default variables are: hgtprs, tmpprs, ugrdprs, vgrdprs


Example Usage: dynamic -var hgtprs tmpprs ugrdprs vgrdprs

This example will use the hgtprs, tmpprs, ugrdprs, vgrdprs variables to calculate the dynamic variables.

Example Image: NARR Shaded: Fronogenesis, Contoured in Green: Fn vector divergence, Vectors: Fn Vectors for 06UTC Feb 15, 2003

850-700mb average: Frontogensis, Fn Vectors, Fn Divergence, plotted for 06UTC Feb 15 2003

I chose this example to compare to the example used in this presentation.  Things seem to match up reasonably well.  Though I do urge some caution here, I am not (nor have I ever been) a dynamics expert, and some of the calculations were taken from either old GrADS scripts, or approximations of textbook equations.  Also, it is possible that the listed units are not 100% correct.  These caveats aside, I have been using this script to successfully diagnose areas of frontogenesis, and ageostrophic circulations from model output for quite some time, so I think overall it serves it's purpose well.  I post it now, with a newly added help page, and slightly cleaner code.

Download Dynamic.gs Here

Download Frontogenesis Example Here

Read More

Tips on using the finite differencing functions in GrADS (e.g., hdivg, and hcurl): Getting rid the undefined strip along the map boundary

The finite differencing functions in GrADS are pretty handy, especially when it comes to plotting variables such as divergence. The finite differencing function used in GrADS is called cdiff(), and what cdiff() does, is simply take the finite difference of a certain variable. For example, if you wanted to take the chance in u-wind in the x-direction, cdiff would be called as such:

    cdiff(u,x)

It is important to note, that this is not du/dx, as it is taking the chance in u from one grid point to the next, and not the physical distance.  If you wanted du/dx, you would have to multiply the denominator by the grid-spacing.  A good example of using the cdiff function can be found here.

GrADS includes a couple of commonly needed functions that use cdiff: hdivg(), and hcurl().  These function use the finite difference method to plot horizontal divergence and vertical vorticity respectively.  A common issue that people may come across while using this function is the appearance of a "blank strip" along the boarder where the data is not plotted.  An example of this strip is below.

GFS 700hPa divergence (hr-1)

It may be hard to notice here, but you can clearly see a strip of white along the boundary.  The reason this appears is simple enough, basically you are plotting a difference, and at the edge, there is only one point, so there is nothing to plot on the edge points.  The way around this is very simple.  All you need to do, is set the domain one grid point larger than you intend to plot, define your variable, then fix the domain to your plot size, and plot it.  That way, you define the finite difference at the edge of your domain before plotting.  Of course, this method relies on having the data outside of your desired domain.  Anyway, the code to do this is simple enough:

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

    'set lat 20 52'
    'set lon -127 -66'
    'set gxout shaded'
    'set mpdset hires'
    'set display color white'
    'clear'

    'set gxout shaded'
    'set lev 700'
    'define var = hdivg(ugrdprs,vgrdprs)*3600'

    'set lat 21 51'
    'set lon -126 -67'

    'd var'
    'cbar'

As you can see the lat/lon boundaries extend out just a little bit farther when defining the variable, then are shrunk when making the actual plot.  The resulting image is:

GFS 700hPA divergence (hr-1)

This image is very similar to the image above, but you now clearly see that there is no "blank" strip along the borders of the image.

The finite difference functions in GrADS can be useful for a whole number of different purposes.  Hopefully these quick tips will help you better use these functions in the future.

Read More

Tutorial: Use the tcorr function to create a map of regression/correlation coefficients in GrADS

Making a correlation map in GrADS is actually quite simple.  A regression/correlation map, is basically a map of regression and correlation coefficients plotted against a time-series of some variable.  The classic example is the monthly average time-series of El-Nino vs. surface temperature on the globe.  The example used in this tutorial will be 500mb height correlated with the 500mb height over the northeast U.S.  For this tutorial we will use the NCEP reanalysis data.

GrADS makes correlation maps very simple.  All you need is a data set that varies in both space and time.  The first step is to define your time-series variable.  In the interest of time and simplicity, I will just take the 500mb height at one point rather than average over a small spatial area, as is often done when doing these maps.  So lets open the data file and set our time-series variable!

    'sdfopen http://monsoondata.org:9090/dods/rean3d'
    'set t 1 636'
    'set lev 500'
    'set lat 44'
    'set lon 282'
    'ts_var=z'

The variable "ts_var" is our time-series variable, and basically we saved the time-series of geopotential height at 500mb at the lat/lon point of 44/282.  Now, we are going to use two intrinsic GrADS functions to plot the regression coefficients (shaded) and the correlation coefficients (contoured above 0.5).

So, to plot our coefficients, we use the 'tregr()' and the 'tcorr()' functions.  Both of these functions do essentially the same thing; you feed them a time-series variable, and a spatial variable, as well as time-constraints, and the function plots your coefficients.

Example:  tregr(x, y, t=1, t=50)

This variable takes the regression coefficients of spatial variable y and time-series x through the first 50 time-steps.  In our example of 500mb height our code looks like:
  
    'set gxout shaded'
    'set t 1'
    'set lat 0 90'
    'set lon 180 360'
    'd tregr(ts_var,z,t=1,t=636)'

    'set gxout contour'
    'set cstyle 2'
    'set cthick 6'
    'set ccolor 1'
    'set clevs 0.5 0.6 0.7 0.8 0.9 1.0'
    'd tcorr(ts_var,z,t=1,t=636)'

Most of the above code is for display of the contour variable; style, thickness, levels, but you can see the tregr() and the tcorr() functions used correctly.  If you are following along at home, the resultant map is plotted.

Shaded:Regression Coefficients, Contoured: Correlation Coefficients

So, what you can gather from this map is that the 500mb geopotential height is generally well correlated above 20 degrees N (r>0.8).  Furthermore, you can pick out the Rossby wave pattern in the regression coefficients.  In anycase, this example show give you an idea of how to make regression/correlation maps in GrADS.

Download Example Script

Read More

A look at the different map projections in GrADS

There are a number of different possible map projections you can use in GrADS, and you may be familiar with some, or all of them, or perhaps you are only familiar with the standard default lat/lon projection used.  This entry will be less of a tutorial and more of just a detailed look at the different map projections you can use in GrADS.  However, before we get started, we do need a little tutorial section, just so you know the syntax for setting the map projection.  Setting the map project is extremely simple: Use the 'set mproj' command and provide it with a map projection from the list below:

1: latlon      Lat/lon projection with aspect ratio maintained (default)
2: scaled      Lat/lon aspect ratio is not maintained; plot fills entire plotting area
3: nps         North polar stereographic
4: sps         South polar stereographic
5: lambert     Lambert conformal conic projection
6: mollweide   Mollweide projection
7: orthogr     Orthographic projection
8: robinson    Robinson projection, requires set lon -180 180, set lat -90 90
9: off         No map is drawn; axis labels are not interpreted as lat/lon

Example:
    'set mproj nps'

That's all there is to it.  Once you have set your projection, you set the lat/lon boundaries as you normally would using the default settings, although certain map projections have intrinsic limits associated with them, which we will explore here.  So, in order to make the examples shown below, we will very simply plot vorticity from the GFS on a bunch of different map projections.  So starting out with the following code we will get this basic image, on the default cylindrical map projection,

    'sdfopen http://nomads.ncep.noaa.gov:9090/dods/gfs_hd/gfs_hd20130513/gfs_hd_00z'
    'set gxout shaded'
    'set lev 600'
    'set mpdset hires'
    'd abs(absvprs)'

Vorticity on a Defualt latlon map projection

So, now that we have that out of the way, lets look at the other map projections.

Scaled:

The scaled map projection does not differ much from the default cylindrical projection, except that it gets rid of the aspect ratio.  Basically, using this option will squish, or elongate the map to fit your specified page area.  The image below shows the transition.  To show the different more dramatically, I set the lat/lon coordinates to -40 to 40 instead of -90 to 90.

latlon vs. scaled map projections

Personally, I don't have much use for the "scaled" map projection, but it is always an option if you need to fit stuff together.

North/South Polar Stereographic NPS/SPS:

The NPS map projection can be useful if you are looking down at the north pole.  Pretend basically that you have the northern hemisphere, and you flatten it.  This would be the north polar stereographic projection.  For the southern hemisphere, simply take the same concept and apply it.  Now, while this map projection does not set latitude boundaries, my rule of thumb has been to set the minimum(maximum, depending on your hemisphere) latitude to be the equator.  I have trouble getting information if I look at southern hemisphere values using an nps projection.  The example below shows each full hemisphere using each projection.  Latitude set from 0 to 90 and -90 to 0.

NPS vs. SPS

Now, you don't need to have longitude set to cover the entire globe either, or the entire hemisphere.

Lambert:

The best way to think about the Lambert map projection is as if the Earth was projected onto a cone.  In fact, often time this projection is referred to as a "conic" projection as well.  Anyway, imagine the Earth is projected onto a cone with the north pole at it's point.  The Lambert projection is basically what the map would look like if the projection was unwrapped and then flattened.  The Lambert projection requires, that your latitude boundaries exist in the same hemisphere, otherwise the projection won't work.  The below example, shows the Lambert projection applied to the northern hemisphere for the longitude coordinates from -180 to 0 (western hemisphere)

Lambert Projection

Mollweide:

The Mollweide projection is often used to show global maps, but in GrADS there are no specific requirements for the latitude/longitude boundaries.  This projection can be neat, as it sort or combines a spherical look with a cylindrical look, with the meridians converging at the pole.  This gives a neat almost 3D appearance.

Mollweide projection

Orthographic Projection:

The orthographic projection is another neat projection that gives a 3D spherical appearance.  In GrADS there are strict requirements for the orthographic projection.  Your latitude boundaries must span -90 to 90.  Additionally, your longitude coordinates must span 180 degrees, no more, no less.  You have the option of choosing your longitude boundaries, so long as they cover 180 degrees of territory.  The result is an image similar to this one.

Orthographic projection

Robinson Projection:

The last map projection I am going to discuss is the "Robinson" Projection (The "off" projection is not really interesting, so I'm not going to talk about it) The Robinson projection is similar to the Mollweide projection, but with the poles kind of flattened out.  This is another commonly used map projection.  In GrADS this projection requires that lat/lon boundaries cover the entire globe.  Lat=-90 to 90, and Lon=-180 to 180 (0 to 360).

Robinson projection

So that's it for the different map projections in GrADS.  As always, each projection is useful depending on the situation, the map, the variable plotted, etc., etc.  I hope this post has given you a little more information regarding each different type of projection.  Just remember that some of them require specific lat/lon boundaries.  Aside from that, picking the right projection is just a matter of trial and error.

Read More