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