Thursday, July 25, 2013

Script: cone.gs; Draws "Cone of Uncertainty" for projected storm path using the MFHILO function

As we start to get into the meat of the northern hemisphere Hurricane season, I figured I would release a couple of entries related to hurricanes and tropical storms.  The first of which is a script I recently put together called "cone.gs"  This script helps plot one of those "uncertainty" cones showing the predicted storm track with potential deviations in it.

Now, I should mention that if you want to just plot the official Hurricane uncertainty cone, you can always read in the shapefile provided here. But, if you want to plot your own uncertainty track this script does a nice job!

Example: Tropical Storm Dorian, 6 day storm track

This script basically uses the MFHILO function to find the storm center, and then draws a cone around the storm track based on your input options for error and everything else.

The example above also uses the 'basemap' script to plot the basemap.  All that is required in the example call is the variable you want to use (in this case,sea level pressure).

Example:
    'cone prmslmsl -fill -fcol 2 -circ -error 50 -dunit km -unit day -end 6 -tint 8'

This example uses many of the built in options:
   -fill: Fills cone
   -fcol: Fill color
   -circ: Draws semi-circle on the end
   -error: Radial error in units: dunit/tunit - in this example 50 km/day
   -dunit: Distance unit for error
   -tunit: Time unit for error
   -end: Total number of time-steps to plot
   -tint: Time interval to plot in timesteps.

The example above uses the GFS data (which runs on a 3 hour time step) to plot the cone.  This block of code will produce the above image.

    file='http://nomads.ncep.noaa.gov:9090/dods/gfs_hd/gfs_hd20130725/gfs_hd_00z'

    'sdfopen 'file
    'set display color white'
    'clear'
    'set mpdset hires'
    'set lon -80 1'
    'set lat 5 30'
    'd hgtprs'
    'rgbset'

    'basemap L 74 1 M'
    'basemap O 45 1 M'

    'printim hurricanebase.png'

    'clear'
    'set t 3'

    'set lat 5 30'
    'set lon -80 1'
    'set clevs 10000000'
    'set ccolor 0'
    'd hgtprs'

    'set map 1'
    'draw map'

    'cone prmslmsl -mnlt 10 -maxlat 25 -mnln -70 -mxln -20 -tunit day -fcol 2'
    'printim ucone.png -b hurricanebase.png -t 0'


The Important thing to note is the use of the maxlat/mnlat/etc... options.  These options define an invisible boundary within your domain for which to search for the storm center, without these you can easily find min values that are not part of your storm, especially with larger domains.  It is recommended that you use these options to help draw your cone.  Also you will need the script basemap.gs and the associated .asc files.

That's about it, here are links to cone.gs and the example script above!  Hope you all find this useful!  Be sure to report any bugs here!  Stay tuned for a tutorial on plotting numerous model predicted storm tracks on the same plot!

Download cone.gs

Download example script

16 comments:

  1. Good day. I've been using GrADS for quite some time now, yet I'm still new to the program. I'm having difficulties on maxlt/mnlat etc. variables. I've been tracking LPAs and tropical cyclones based on the Philippines using GrADS. I don't know the values for the maxlt/mnlat etc. to be used for my location. Hope you would help me with that. Thank you for your further response. God Bless!

    ReplyDelete
    Replies
    1. Are you talking about the minlat/maxlat values for this script? Or are you talking about finding the lat/lon coordinates of your low pressure center.

      If you are talking about the values in this script, the values are your values to choose. What they do is put a "bounding" box around the region you are interested in, so that you avoid jumping from one low pressure area to another. If you used this script over an larger domain, there is a greater chance that you will get crazy results from jumping around to different low pressure areas. But if you restrict the region you are focused on, you will simply track the same LPA. So choose these values as you need them.

      If you are looking for a way to find the coordinates of you minimum pressure, you can just dig around through this script. The function used is the MFHILO function, and that finds your minimum value in a given domain. I believe you could also nest a "minloc" function in x and y to find the coordinates as well, though as I recall this was less intuitive than using the mfhillo method.

      Hope this answers your question!
      Good luck!

      Delete
  2. Hi there ,
    Great GrADS script to plot Cone Graphic. I try to use that script to plot cone graphic for my imaginary typhoon using my storm surge model. Below are the best track data i'd created :

    66666 1501 16 0002 1501 0 6 INDERA 20150101
    15051500 002 2 65 1131 1000 030 00000 0000 60070 0070

    15051512 002 3 47 1102 998 055 00000 0000 60070 0070

    15051600 002 4 51 1075 990 085 00000 0000 60100 0080 15051612 002 4 46 1039 970 091 00000 0000 60100 0080

    15051700 002 5 43 1015 948 128 00000 0000 60120 0090
    15051712 002 5 54 0973 930 141 00000 0000 60130 0100
    15051800 002 5 59 0954 950 125 00000 0000 60150 0120

    15051812 002 4 58 0951 998 078 00000 0000 60150 0120

    It uses the JMA RSMC Best Track text data format.

    So i use my control file :
    dset ^jma_ssm.grb
    index ^jma_ssm.grb.idx
    undef 9.999E+20
    title jma_ssm.grb
    * produced by grib2ctl v0.9.13
    * command line options: -verf jma_ssm.grb
    dtype grib 255
    ydef 901 linear -6.000000 0.017
    xdef 1681 linear 92.000000 0.017000
    tdef 15 linear 00Z15may2015 6hr
    zdef 1 linear 1 1
    vars 4
    DSLMsfc 0 82,1,0 ** surface Deviation of sea level from mean [m]
    PRMSLmsl 0 2,102,0 ** mean-sea level Pressure reduced to MSL [Pa]
    UGRD10m 0 33,105,10 ** 10 m above ground u wind [m/s]
    VGRD10m 0 34,105,10 ** 10 m above ground v wind [m/s]
    ENDVARS

    Here are the result when i run the "q file" in grads :
    ga-> q file
    File 1 : jma_ssm.grb
    Descriptor: jma_ssm.ctl
    Binary: jma_ssm.grb
    Type = Gridded
    Xsize = 1681 Ysize = 901 Zsize = 1 Tsize = 15 Esize = 1
    Number of Variables = 4
    dslmsfc 0 82 ** surface Deviation of sea level from mean [m]
    prmslmsl 0 2 ** mean-sea level Pressure reduced to MSL [Pa]
    ugrd10m 0 33 ** 10 m above ground u wind [m/s]
    vgrd10m 0 34 ** 10 m above ground v wind [m/s]

    I use script file cone_indera.gs as below :
    *'sdfopen 'file
    'open jma_ssm.ctl'
    'set display color white'
    'clear'
    'set gxout contour'
    'set mpdset hires'
    *'set lon 92 120'
    *'set lat -6 9'
    'd dslmsfc'
    *'rgbset'

    'basemap L 74 1 M'
    'basemap O 45 1 M'

    'printim indera.png'

    'clear'
    'set t 3'
    *'set lon 92 120'
    *'set lat -6 9'
    *'set lat 5 30'
    *'set lon -80 1'
    *'set lev 850'
    'set gxout contour'
    *'set clevs 10000000'
    'set ccolor 0'
    'd PRMSLmsl/100'
    'set map 1'
    'draw map'


    *'cone PRMSLmsl -e 50 -tunit day -dunit km -e_reduct 8 -end 3 -tint 4 -fcol 2'
    'cone prmslmsl/100 -e 50 -tunit day -dunit km ucone.png-fcol 2'
    *'printim ucone.png -b indera.png -t 0'

    Do you have any idea?. Your prompt comment on these is highly appreciated.

    Thanks,
    KChang

    ReplyDelete
  3. By the way, below are the error messages i got :
    ******nlo exceeds nmin=2000
    ******nlo exceeds nmin=2000
    ******nlo exceeds nmin=2000
    ******nlo exceeds nmin=2000
    ******nlo exceeds nmin=2000
    ******nlo exceeds nmin=2000
    ******nlo exceeds nmin=2000
    ******nlo exceeds nmin=2000
    ******nlo exceeds nmin=2000
    ******nlo exceeds nmin=2000
    ******nlo exceeds nmin=2000
    ******nlo exceeds nmin=2000
    ******nlo exceeds nmin=2000
    ******nlo exceeds nmin=2000
    ******nlo exceeds nmin=2000
    ******nlo exceeds nmin=2000
    ******nlo exceeds nmin=2000
    ******nlo exceeds nmin=2000
    ******nlo exceeds nmin=2000
    ******nlo exceeds nmin=2000
    ******nlo exceeds nmin=2000
    ******nlo exceeds nmin=2000
    ******nlo exceeds nmin=2000
    ******nlo exceeds nmin=2000
    ******nlo exceeds nmin=2000
    ******nlo exceeds nmin=2000
    mfhilo CL method N: 1 maxmin: -1 radinf: 1000.0 [nm] cintinf: 300.0 latc: 1.65 lonc: 106.28
    L -0.59 96.96 M: 1010.51 D: 937.54
    Data Request Warning: Request is completely outside file limits
    mfhilo CL method N: 0 maxmin: -1 radinf: 1000.0 [nm] cintinf: 300.0 latc: 1.65 lonc: 106.28
    Query Error: Syntax is QUERY W2XY Lon Lat
    Non-numeric args to numeric operation
    Error occurred on line 337
    In file /opt/opengrads/Contents/Resources/Scripts/cone.gs


    Do you have any idea?. Your prompt comment on these is highly appreciated.

    Thanks,
    KChang

    ReplyDelete
    Replies
    1. I'll have to look more intently at the script to give you a definitive answer, but at first glance is appears that there is an issue when the script calls the "mflhilo" function, which is the function use to find the location of the minimum SLP. The outside of file limits request makes me think you are setting your boundaries wrong, I most commonly do this when I try to set a negative longitude when the data is gridded from 0-360, not -180-180.

      Also it looks like there is a weird call in this :
      'cone prmslmsl/100 -e 50 -tunit day -dunit km ucone.png-fcol 2'

      I don't think the ucone.png should be in there (though honestly, it's been a while since I dug around this script so its entirely possible I put something looking for an image in there...).

      I will glance at the script later and maybe I can give you a better answer. But at first look, those are my thoughts.

      Delete
    2. ok..then let me first look at your comments and try to resolve the issues...

      Delete
  4. I modified the script but still received the following error message :
    ga-> cone_Indera.gs
    mfhilo CL method N: 1 maxmin: -1 radinf: 1000.0 [nm] cintinf: 300.0 latc: 5.50 lonc: 104.00
    L 4.49 101.67 M: 948 D: 151.86
    Tambahan KC :
    4.49
    101.67
    6
    Data Request Warning: Request is completely outside file limits
    mfhilo CL method N: 0 maxmin: -1 radinf: 1000.0 [nm] cintinf: 300.0 latc: 5.50 lonc: 104.00
    Tambahan KC :
    4.49 grid
    101.67 contents
    6
    Query Error: Syntax is QUERY W2XY Lon Lat
    Non-numeric args to numeric operation
    Error occurred on line 341
    In file /opt/opengrads/Contents/Resources/Scripts/cone.gs

    ReplyDelete
  5. Here is my script :

    *'sdfopen 'file
    'open jma_ssm.ctl'
    'set display color white'
    'clear'
    'set gxout contour'
    'set mpdset hires'
    'set t 1'
    'set lon 92 120'
    'set lat -6 9'
    'd PRMSLmsl/100'
    *'rgbset'

    'basemap L 74 1 M'
    'basemap O 45 1 M'

    'printim indera.png'

    'clear'
    'set t 5'
    'set lon 92 120'
    'set lat -6 9'
    *'set lat 5 30'
    *'set lon -80 1'
    *'set lev 850'
    'set gxout contour'
    'set clevs 10000000'
    'set ccolor 0'
    'd PRMSLmsl/100'
    'set map 1'
    'draw map'
    *'cone PRMSLmsl -e 50 -tunit day -dunit km -e_reduct 8 -end 3 -tint 4 -fcol 2'
    'cone prmslmsl/100 -e 50 -mnlt 4 -maxlat 7 -mnln 94 -mxln 114 -tunit50 -tunit day -dunit km -fcol 2'
    *-mnlt 10 -maxlat 25 -mnln -70 -mxln -20 -tunit
    'printim ucone.png -b indera.png -t 0'

    ReplyDelete
  6. The output :
    Tambahan KC :
    4.49
    101.67
    6

    Is just additional lines that i print to check values..

    ReplyDelete
  7. Hi there,

    Did you manage to resolve the errors in my scripts?. Your opinion to solve this issue is highly appreciated.

    ReplyDelete
  8. Hi there,

    I managed to solve the problems and had successfully plot the cone graphic. Just the -tint and the -maxtime options. But do you know how can i control the transperacy of filled colour?

    ReplyDelete
    Replies
    1. Sorry I didn't get back to your earlier, I am glad you solved the problem and that it works okay! Unfortunately, one of the deficiencies of GrADS is that it does not do transparency unless you have version 2.1 or newer. If you have 2.1 or newer try following the last example on this page. I know it's for shapefiles, but it may work for the polygons as well.

      http://www.iges.org/grads/gadoc/shapefiles.html

      Delete
    2. Hi KC, may I get your code? I got the same error

      Delete
  9. Ok.thanks,i will use the set rgb function in grads 2.1 above. Maybe i can modify your cone.gs to show model date in every point the cone plot

    ReplyDelete
  10. Ok.thanks,i will use the set rgb function in grads 2.1 above. Maybe i can modify your cone.gs to show model date in every point the cone plot

    ReplyDelete
  11. I am attempting to use the mflhilo functions but cannot get it to work. I realize I need to enable the user defined extensions but am failing at getting this setup. I have grads v2.1.a3. Any help would be appreciated

    ReplyDelete