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

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: Shapefile Basics

The capability of GrADS to read shapefiles brings a whole new dimension to the plotting software.  It provides you with a whole new range of possibilities when it comes to plotting data.  With the use of shapefiles, you are no longer limited to generic maps with data plotted.  Shapefiles provide you the possibility of plotting everything from county boundaries, to roads and rivers, to congressional districts, weather watches and warnings, and much more.  This tutorial will give you a very basic introduction into using shapefiles in GrADS.

It is recommended (not required): that you have both color.gs and cbar.gs in your script library to plot the RADAR color scale and associated colorbar.

For this tutorial, we will simply plot county boundaries on a map of the US.  Before we get started, we will first need to download the proper shapefiles.  The county shapefiles are available from NOAA here.  This site should provide you with a table showing recent versions of the shapefile, download and unzip it.  The list of files included within the zipped folder shows four files, these are the files from the April 2013 version:

• c_02ap13.prj
• c_02ap13.shp
• c_02ap13.shx
• c_02ap13.dbf 

Most shapefiles I have encountered have some combination of these files, and because each shape tends to come with multiple files, I have found it useful to create a separate folder to store my shapefiles in, aptly named: Shapefiles.

Now that you have the county shapefiles in your proper folder we can begin with the tutorial.  The first, and possibly most important thing to know about shapefiles is that you need to plot something from your GrADS data before you can draw your shapefiles.  If you don't you will get an error that says: "No Scaling Environment."   

So, in order to plot our counties, we need to plot a variable first, which means we will need a file.  For this example we will use the RAP model using the 'sdfopen' command.

     'sdfopen http://nomads.ncep.noaa.gov:9090/dods/rap/rap20130507/rap_00z'

Now that our data file is open, the rest of the script is very short and most of it is just setting up the domain and display options.  If you have read through any previous tutorials, then next block of code will look familiar.


    'clear'
    'set lat 22 51'
    'set lon -128 -65'
    'set gxout shaded'
    'set mpdset hires'
    'set map 0'     ;*Sets map outline color to Black

With the domain set up, we just need to plot the variable.  For this example, we will plot the model predicted RADAR reflectivity.  I included in this block of code with a good approximation of the NWS color scale for reflectivity using color.gs.

 *If you have color.gs*
     'color 5 75 5 -kind white->dodgerblue->blue->lime->green->darkgreen->yellow->goldenrod->orange->firebrick->red->darkred->fuchsia->indigo'
**********************
     'd refcclm'
     'cbar'
  
Now that the RADAR is plotted, its time to draw the shapefile.  This is perhaps the easiest part of the tutorial.  You simply use the 'draw shp' command and point to the file path.  In this example:

    'draw shp Shapefiles/c_03ap12.shp'

This command is all you need to draw the counties on the map.  If you wish to change the properties of the shapefile you can do so using a few additional commands.  The outline color can be set using the 'set line command.  The fill properties can be set using the 'set shpopts' command.  By default, shapefiles are not filled.

    'set line 0'       ;*Where the outline color is set to 0
     'set shpopts 15'    ;*Fills shapefile, where 15 is the number of the desired fill color. A value of -1 turns off fill.

Depending on your exact output options, your final map should look similar to this.

RAP Reflectivity map with counties

Using shapefiles in GrADS can get much more complex.  However, for the purposes of this tutorial, I wanted to provide a real basic introduction to shapefiles in GrADS.  I may, at a later time, write up more advanced tutorials for shapefiles.  In the mean time, if you are interested in learning more about shapefiles in GrADS, I encourage you to check out this resource to get a more comprehensive look.

Again, you can download shapefiles for pretty much everything you could imagine, here is just a small list of shapefile resources to get you started.

• http://www.nws.noaa.gov/geodata/ 
• http://www.nws.noaa.gov/gis/shapepage.htm 
• http://www.statsilk.com/maps/download-free-shapefile-maps 

 
Download this Script

Read More