Friday, May 10, 2013

Tutorial:Two methods to plot more than one shaded variable on the same map using GrADS

Plotting more than one shaded variable on the same map can prove frustrating as the new shaded variable continually plots over the old variable, even where the new variable is below the minimum value.  However, there are a couple of ways to plot more than one variable on the same map without this issue.  This tutorial will teach you how to do this.

The first, and perhaps simplest, way to plot more than one shaded variable on the same plot is to use the 'maskout' function.  The 'maskout' function basically scans your variable and only plots values higher than a specified threshold.  Without the maskout function, values below this threshold are plotted at every grid-point, only at the lowest specified color, so the maskout function keeps you from overwriting the previously plotted variable.  A good example use of this method is plotting different precipitation types on the same map (shown below).  This part of the tutorial will take you through the steps of using the maskout function to plot the precipitation types.

Example of precipitation type using 'maskout' function

Most weather models available from the GrADS Data Server include 3-4 variables corresponding to a categorical surface precipitation type: Rain, Ice Pellets, Freezing Rain, and Snow.  The trick to plotting these on the same map is to plot these variables using the maskout function, masking out values less than 0.5.  What this means, is that any value less than 0.5 will not be plotted for each precipitation type.  A short script showing how to do this is below.  This includes changing the colors of each precipitation category.  Just for fun, we'll focus this on Alaska instead of the continuous U.S.

    'sdfopen http://nomads.ncep.noaa.gov:9090/dods/nam/nam20130510/nam_alaskanest_00z'

    'clear'
    'set datawarn off'

    'set clevs  1'
    'set ccols 3 3'
    'd maskout(crainsfc, crainsfc-0.5)'  ;*Rain

    'set clevs  1'   
    'set ccols 2 2' 
    'd maskout(cfrzrsfc, cfrzrsfc-0.5)'  ;*Freezing Rain

    'set clevs  1'
    'set ccols 9 9'   
    'd maskout(cicepsfc, cicepsfc-0.5)'  ;*Ice Pellets


    'set clevs  1' 
    'set ccols 4 4'
    'd maskout(csnowsfc, csnowsfc-0.5)'  ;*Snow


This script will produce the following image, depending on what time step you choose to plot.  The command 'set datawarn off' is there to get rid of the pesky "ENTIRE GRID UNDEFINED" Text that often crops up on your plot.  Basically, often times the model doesn't have any values of the ice variables, and it reads that as undefined when using the maskout function.

Maskout Example: Rain in Green, Snow in Blue

Anyway, you can clearly see two shaded variables plotted on the same map.  So, lets move on to the 2nd method.

The 2nd way to plot two shaded variables on the same plot is to save an image after the first variable is plotted, and then plot the 2nd variable, and save a new image using the first image as a basemap.  To ensure that your newly shaded variable doesn't  overwrite the old one, simply set the minimum threshold color as transparent.  The example used will be plotting model reflectivity overtop clouds.  Note: if you unsure how to save images onto a basemap, check out this tutorial.

Cloud/Background
Reflectivity Over Background


The code below will produce the above images.  Note: to get the colors for reflectivity, you will need to get the scripts colorset.gs and color.gs.  You could also use color.gs to approximate the RADAR colors as well.

    'sdfopen http://nomads.ncep.noaa.gov:9090/dods/nam/nam20130510/nam_alaskanest_00z'

    'set display color white'
    'clear'
    'set gxout shaded'
    'set mpdset hires'

    'clear'
    'color 60 100 2 -kind white->gray'
    'd tcdcclm'
    'printim image1.png'

    'colorset Rad1'
    'd refcclm'
    'printim image2.png -b image1.png -t 16'


You can see from the code above that the cloud variable is plotted first and then saved as "image1.png"  Then reflectivity is plotted using the color scale defined by colorset.gs (which sets the minimum value color to 16).  The last line of this code then saves the newly plotted variable to "image2.png" using "image1.png" as a basemap.

This method provides an alternative to the maskout function and, personally, I find this method superior in a few ways, for one the maskout function causes the variables to have a 'blocky' sort of pixelated look, where as this method keeps the contouring smooth. Also, it removes the need to worry about the 'datawarn' command.

Regardless of which method you choose to use, both work well for plotting more than one shaded variable on the same map.  For the downloadable script, I simply included one script with both methods inside, for your reference.

Download this Script


3 comments:

  1. I apologize for what is probably an extremely simple/stupid question... I noticed on one of the maps here (and a few others on the blog) that you have UTC and local time plotted. I've been trying to figure out how to get local time to plot with data that is in UTC, since not only the time changes, but the date (and day of week) can also be different (e.g. 00 UTC is not only 1800 CST, but is also the previous day of week and date). Any suggestions you might have on working with the display of dates/times would be appreciated. Thank you.

    ReplyDelete
  2. There are two ways to do this: The best way to do this is to simply use the "set t" command to set your time to the offset you want and grab the time (though this can be tricky if your model data isn't timed by 1 hour intervals). For example if you wanted to do an offset of 6 hours, you could simply do a:

    'set t 't-6

    Then you use the 'q dims' command to get the time string information, and you just ignore the UTC in the string. Then you must remember to set the time back to it's original spot. If your time step is 3 hours, then you would do 'set t 't-2, for a six hour offset. If your offset is not divisible by 3, then you have issues with this method, which is why this is tricky if your data isn't in 1 hour intervals.

    The second way to do this is to isolate the different aspects of the time string, (hour, day, month, year) and then do your offset here. This way can be a little more complex since you would need if statements to help you ensure that you change the month and the year around as well, but this isn't that terribly hard, just takes a little time to input the information. So you do this by getting the time dimensions, and then using the substr() command to get each piece of the date. Then its just a matter of applying your offset to the hours, and then deciding to alter the date as needed.

    Hopefully some of that makes sense!

    ReplyDelete
  3. Aha, that is simple, thank you. I stupidly never thought of changing t to get the time I need, then changing it back afterward.

    ReplyDelete