[gradsusr] how to find typhoon center via minloc function

sushant puranik sushantpuranik at gmail.com
Fri Apr 29 08:32:14 EDT 2011


Hello Charles
I want to make a jet stream plot. Can i make it with the help of hurricane
center script by changing minloc and min to maxloc and max.
any suggestion?

thanks

Sushant

On Wed, Jan 5, 2011 at 1:01 AM, Charles Seman <Charles.Seman at noaa.gov>wrote:

> Jie,
>
> Attached is a file File_Error_1_Grads_1.9-Plot_tc_shi.gs.txt containing the
> text from an email exchange to the GrADS Listserv between Joe Covert and
> Diane Stokes on 9/13/2004 -- on 4/27/06, found via a Google search for "
> plot_tc_shi.gs": http://caos.iisc.ernet.in/gslib/plot_tc_shi.gs
>
> The attached GrADS script plot_hurricane_center3.gs (along with pause.gsf,
> aGrADS script function used by the plot script) is a locally modified
> version derived from Joe Covert's script "plot_tc_shi.gs" (
> Joe.Covert at noaa.gov) which was posted to the GrADS Listserv on 9/13/2004
> (his script was named plot_tc_shi.gs) and code from hurricane_tracking.txt
> (taken from an email exchange posted to the GrADS Listserv by Kun-Hsuan Chou
> and Arturo Caracas Uribe in Oct, 2004)...
>
> Please find below code from hurricane_tracking.txt (taken from an email
> exchange posted to the GrADS Listserv by Kun-Hsuan Chou and Arturo Caracas
> Uribe in Oct, 2004) illustrating a technique for finding max wind in a 2D
> field... and below that some code from attached script
> plot_hurricane_center3.gs to find the min sea-level pressure in a 2D
> field...
>
> code from Kun-Hsuan Chou to find maximum wind speed from
> hurricane_tracking.txt:
> ---
> 'd maxloc(max(mag(u,v),lon=120,lon=130),lat=15,lat=25)'
>  line=sublin(result,2)
>  ygrd=subwrd(line,4)
> 'd maxloc(max(mag(u,v),lat=15,lat=25),lon=120,lon=130)'
>  line=sublin(result,2)
>  xgrd=subwrd(line,4)
> 'set x 'xgrd
> lonval = subwrd(result,4)
> 'set y 'ygrd
> latval = subwrd(result,4)
> 'q w2xy 'lonval' 'latval
> xpos = subwrd(result,3)
> ypos = subwrd(result,6)
> 'draw mark 1 'xpos' 'ypos' .2'
> ---
>
> sample code from plot_hurricane_center3.gs to find hurricane center (see
> script for supporting code):
> ---
> *
> * find minimum "pressure" within box area (x1,x2), (y1,y2)...
> *
>     'set x 'x1
>     'set y 'y1
>     'set z 1'
>     'set t 'tt
>     'set gxout print'
>     nxp = x2-x1+1
> *      say '    nxp = 'nxp
>     nyp = y2-y1+1
> *      say '    nyp = 'nyp
>     'd minloc(min(psl,y='y1',y='y2'),x='x1',x='x2')'
>       rec=sublin(result,nxp+3)
>       xc=subwrd(rec,1)
>     'd minloc(min(psl,x='x1',x='x2'),y='y1',y='y2')'
>       rec=sublin(result,nyp+3)
>       yc=subwrd(rec,1)
>     say
>     say 'location of minimum "pressure"...'
>     say
>     say ' (xc,yc) = ('xc','yc')'
> *
> * find "world" coordinates of (xc,yc) and convert "world" coordinates
> * to "xy" coordinates for plotting track of hurricane center...
> *
>     'set x 'xc
>     lonval = subwrd(result,4)
>     'set y 'yc
>     latval = subwrd(result,4)
>     say ' (lonval,latval) = ('lonval','latval')'
>     'q w2xy 'lonval' 'latval
>     xpos = subwrd(result,3)
>     ypos = subwrd(result,6)
>     say ' (xpos,ypos) = ('xpos','ypos')'
> *
> * write (xpos,ypos) to output file...
> *
>     res = write (''hurricane.nf'','t = 'tt' xpos = 'xpos' ypos = 'ypos'')
>
>     pause()
> ---
>
> Please let me know if you have any questions.
>
> Hope this helps,
> Chuck
>
> Jie TANG wrote:
>
>>
>> hi,grads folks ,
>>  I am using grads v1.9, trying to find the center of typhoon via minloc
>> function.
>> the key script is shown as below:
>> slplat=min(slp,lat=19,lat=26)
>> tclon=minloc(slplat,lon=119,lon=126)
>> slplon=min(slp,lon=119,lon=126)
>>  tclat=minloc(slplon,lat=19,lat=26)
>>
>> but grads tell me that "function not found min " and when i changed my
>> script to be :
>> tclat=min(minloc(slpt,lon=119,lon=126), lat=19,lat=26)
>> tclon=min(minloc(slpt, lat=19,lat=26 )lon=119,lon=126)
>>
>>
>> how can i finx my scrpit ? thank you .:)
>> --
>>
>> TANG Jie
>> Email: totangjie at gmail.com <mailto:totangjie at gmail.com>
>>
>> Tel: 0086-2154896104
>> Shanghai Typhoon Institute,China
>> ------------------------------------------------------------------------
>>
>> _______________________________________________
>> gradsusr mailing list
>> gradsusr at gradsusr.org
>> http://gradsusr.org/mailman/listinfo/gradsusr
>>
>>
>
> --
>
> Please note that Charles.Seman at noaa.gov should be considered my NOAA
> email address, not cjs at gfdl.noaa.gov.
>
> ********************************************************************
> Charles Seman                                Charles.Seman at noaa.gov
> U.S. Department of Commerce / NOAA / OAR
> Geophysical Fluid Dynamics Laboratory         voice: (609) 452-6547
> 201 Forrestal Road                              fax: (609) 987-5063
> Princeton, NJ  08540-6649            http://www.gfdl.noaa.gov/~cjs/
> ********************************************************************
>
> "The contents of this message are mine personally and do not reflect any
> official or unofficial position of the United States Federal Government,
> the United States Department of Commerce, or NOAA."
>
>
>
> worked great ... thanks.  joe
>
> Diane Stokes wrote:
>
> > Joe,
> >
> > Read works fine for me in 1.9.
> >
> > I think this is the issue where 'pull' in 1.9 throws in a carriage
> > return.  The read statement is not getting the exact filename.
> >
> > Try adding:
> >   file=sublin(file,1)
> > after each:
> >   pull file
> >
> > Diane
> >
> > Joe Covert wrote:
> >
> >> There appears to be a problem with version 1.9 of GrADS "read" function.
> >>
> >> The plot_tc_shi.gs script (see below) uses the read function to open
> the
> >> track data file and the program gives the subject error message "File
> >> Error 1".
> >>
> >> jc
> >>
> >>
> >> *  Script to draw an hurricane-track plot.
> >> *  Does little error checking on the input file.
> >> *  Assumes the input file is set up as follows:
> >> *
> >> *    Line 1:  Title
> >> *    Line 2:  Drawing primitives for marks: marktype size
> >> *    Line 3:  Drawing primitives for lines: color style thickness
> >> *    Line 4:  Starting hour and the interval of plotting points
> >> *             e.g., 0 6 means that track starts at 0 hour and mark
> >> *                   will be plotted every 6 hours.
> >> *    Rest of lines:  hour  long.  lat.
> >> *             e.g.,   0    -70.5  25.0
> >> *                     6    -71.8  25.2
> >> *                            :
> >> *                            :
> >> *
> >> *  Also assumes that a file has been opened (any file, doesn't
> >> *  matter -- the set command doesn't work until a file has been
> >> *  opened).
> >> *
> >>
> >> function main()
> >>
> >> *  'clear'
> >>
> >>   'open dummy.ctl'
> >>   'set lat 20 50'
> >>   'set lon -90 -30'
> >>   'set mpdset hires'
> >>   'set poli on'
> >>   'draw map'
> >>
> >>
> >> say 'Enter File Name: (type q to stop)'
> >> pull fname
> >>
> >> while (fname != 'q')
> >>
> >> *  Read the 1st record: Title
> >>
> >>   ret = read(fname)
> >>   rc = sublin(ret,1)
> >>   if (rc>0)
> >>       say 'File Error 1'
> >>       return
> >>   endif
> >>   title = sublin(ret,2)
> >>   say title
> >>
> >> *  Read the drawing primitives
> >>
> >>   ret = read(fname)
> >>   rc = sublin(ret,1)
> >>   if (rc>0)
> >>      say 'File Error 2'
> >>      return
> >>   endif
> >>   dpline = sublin(ret,2)
> >>   marktype = subwrd(dpline,1)
> >>   marksize = subwrd(dpline,2)
> >>   ret = read(fname)
> >>   rc = sublin(ret,1)
> >>   if (rc>0)
> >>      say 'File Error 3'
> >>      return
> >>   endif
> >>   dpline = sublin(ret,2)
> >>   lcolor = subwrd(dpline,1)
> >>   lstyle = subwrd(dpline,2)
> >>   lthick = subwrd(dpline,3)
> >>   say ' marktype, marksize, lcolor, lstyle and lthick:'
> >>   say ' 'marktype ' ' marksize ' ' lcolor ' ' lstyle ' ' lthick
> >>
> >> * Read starting hour and the interval hours of plotting points
> >>
> >>   ret = read(fname)
> >>   rc = sublin(ret,1)
> >>   if (rc>0)
> >>      say 'File Error 4'
> >>      return
> >>   endif
> >>   dhour = sublin(ret,2)
> >>   start = subwrd(dhour,1)
> >>   jump = subwrd(dhour,2)
> >>   say ' starting hour and the interval hours of plotting points:'
> >>   say '  'start' 'jump
> >>
> >> *  Read all data points
> >>
> >>   ret = read(fname)
> >>   rc = sublin(ret,1)
> >>   while (rc = 0)
> >>       loc = sublin(ret,2)
> >>       hour = subwrd(loc,1)
> >>       dlong.hour = subwrd(loc,2)
> >>       dlat.hour = subwrd(loc,3)
> >> *      prompt ' hour ' hour ' are read,'
> >> *      say ' long=' dlong.hour '    lat=' dlat.hour
> >>       ret = read(fname)
> >>       rc = sublin(ret,1)
> >>   endwhile
> >>
> >>   if (rc!=2 & rc!=0)
> >>          say 'File Error 5, rc=' rc
> >>          return
> >>   endif
> >>
> >>   endhour = hour
> >>   say ' endhour=' endhour
> >>
> >> * Plotting
> >>
> >>   'set line 'lcolor' 'lstyle' 'lthick
> >>   'query w2xy 'dlong.start' 'dlat.start
> >>   xprev = subwrd(result,3)
> >>   yprev = subwrd(result,6)
> >>   'draw mark 'marktype' 'xprev' 'yprev' 'marksize
> >>   next = start+jump
> >>   while (next <= endhour)
> >> *      say ' 'dlong.start' 'dlat.start
> >>       'query w2xy 'dlong.next' 'dlat.next
> >>       xnext = subwrd(result,3)
> >>       ynext = subwrd(result,6)
> >>       'draw line 'xprev' 'yprev' 'xnext' 'ynext
> >> *      say ' 'xprev' 'yprev' 'xnext' 'ynext
> >>       'draw mark 'marktype' 'xnext' 'ynext' 'marksize
> >>       next = next+jump
> >>       xprev = xnext
> >>       yprev = ynext
> >>   endwhile
> >>
> >>   say fname ' is working fine.'
> >>
> >> * read in the filename to be plotted next
> >>
> >>   say 'Enter File Name: (type q to stop)'
> >>   pull fname
> >>
> >> endwhile
> >>
>
>
> *----------------------------------------------------------------------
> *
> *  Adapted from plot_hurricane_center2c.gs, which was adapted from:
> *
> *  (1) "plot_hurricane_center2[a][b].gs" (from "plot_hurricane_center.gs
> "),
> *      locally modified versions derived from Joe Covert's script
> *      "plot_tc_shi.gs" (Joe.Covert at noaa.gov), which was posted to the
> *      GrADS Listserv on 9/13/2004 (his script was named plot_tc_shi.gs)
> *      on 4/27/06, found via a Google search for "plot_tc_shi.gs":
> *      http://caos.iisc.ernet.in/gslib/plot_tc_shi.gs
> *
> *  (2) code in "hurricane_tracking.txt" (taken from an email exchange
> *      posted to the GrADS Listserv by Kun-Hsuan Chou and
> *      Arturo Caracas Uribe in Oct, 2004)
> *
> *
> *  To use: grads -l
> *
> *----------------------------------------------------------------------
>
> function main()
>
>  'reinit'
>
> ************************************************************************
> *  Allow external GrADS functions
> ************************************************************************
>
>  rc = gsfallow("on")
>
> ************************************************************************
> *  Define dataset file information and plotting limits for base map...
> ************************************************************************
>
> **  nc_file1 = ''
> **  nc_file2 = ''
> **  nc_file3 = ''
> **  nc_file = ''nc_file1' 'nc_file2' 'nc_file3''
>  nc_file = 'your_netCDF_file'
>
>  lon1 = 298 ; lon2 = 308
>  lat1 =  20 ; lat2 =  25
>
> ************************************************************************
> *  Legend captions, colors, and marks for the plots
> ************************************************************************
>
> **  nmodels = 1
> **  model.1 = 'isabel_ras'
> **  model.2 = 'isabel_isotke'
> **  model.3 = 'isabel_my25'
>
> **  titles.1 = model.1 ; colors.1 = 2 ; marks.1 = 2
> **  titles.2 = model.2 ; colors.2 = 3 ; marks.2 = 2
> **  titles.3 = model.3 ; colors.3 = 4 ; marks.3 = 2
> ***
> ***  Define the legend plotter input...
> ***
> **  titles = '-t' ; colors = '-c' ; marks = '-m' ; lines = '-l'
> **  n=1
> **  while ( n <= nmodels )
> **    titles = ''titles' "'titles.n'"'
> **    colors = ''colors' 'colors.n
> **    marks = ''marks' "'marks.n'"'
> **    lines = ''lines' 1'
> **    n=n+1
> **  endwhile
>
> **  legend_info = ''colors' 'lines' 'marks' 'titles' -p'
>
>  nmodels = 1
>  model.1 = 'title_name_for_your_netCDF_dataset'
>
>  titles.1 = model.1 ; colors.1 = 2 ; marks.1 = 2
> *
> *  Define the legend plotter input...
> *
>  titles = '-t' ; colors = '-c' ; marks = '-m' ; lines = '-l'
>  n=1
>  while ( n <= nmodels )
>    titles = ''titles' "'titles.n'"'
>    colors = ''colors' 'colors.n
>    marks = ''marks' "'marks.n'"'
>    lines = ''lines' 1'
>    n=n+1
>  endwhile
>
>  legend_info = ''colors' 'lines' 'marks' 'titles' -p'
>
> *  ...other drawing primitives
>
>  marksize = 0.1
>
>  lstyle = 1 ; lthick = 5
>
> ************************************************************************
> *  Use "cbar_line" or "cbar_line_box"?
> ************************************************************************
>
>  legend_plotter = 'cbar_line'
> *  legend_plotter = 'cbar_line_box'
>
> ************************************************************************
> *  Plot page limits
> ************************************************************************
>
>  plot_vpage = '0 11 0 8.5'
>  plot_area  = '1 10 1 7.5'
>
>  xl = subwrd(plot_area,1) ; xr = subwrd(plot_area,2)
>  yb = subwrd(plot_area,3) ; yt = subwrd(plot_area,4)
>
>  xc = xl + (xr-xl)/2
>
> ************************************************************************
> *  GrADS metafile output...
> ************************************************************************
>
>  'enable print hurricane_track.gx'
>
> *-----------------------------------------------------------------------
> *  Define and plot track of minimum surface pressure for each dataset...
> *-----------------------------------------------------------------------
>
>  nf=1
>  while ( nf <= nmodels )
>
> * Open dataset file and define parameters...
>
>    fname = subwrd(nc_file,nf)
>    'sdfopen 'fname
>    'set dfile 'nf
>    'q file'
>    rec3 = sublin(result,3) ; binary_file = subwrd(rec3,2)
>    say
>    say
>    say '*** 'binary_file' ***'
>    say
>    rec5 = sublin(result,5)
>    nx = subwrd(rec5,3) ; ny = subwrd(rec5,6) ; nz = subwrd(rec5,9) ; nt =
> subwrd(rec5,12)
>
>    if( nf = 1 )
>
> * Plot base map...
>
>      plot_base_map(plot_vpage,plot_area,lon1,lon2,lat1,lat2)
>
>    endif
>
> * Define locations of "hurricane center" for plotting points, and write
> * locations out to ASCII text file for read-in for plotting track later...
>
>    hurricane.nf = 'hurricane_'nf'.txt'
>
>    ts=2
>    tt=ts
>    while ( tt <= nt )
>      'set t 'tt
>      say
>      say '---------'
>      say ' tt = 'tt
>      say '---------'
>      say
>      'set lon 'lon1' 'lon2
>      'set lat 'lat1' 'lat2
>      'set z 1'
> *
> * here is where you put in the name of your sea-level pressure variable:
> *
>      'define psl = name_of_your_sea_level_pressure_variable'
>      'set gxout shaded'
>      'd psl'
> *
> * following interactive "box location" code adapted from "cbar_line_box.gs
> "...
> *
>      say 'Click where you want the left upper corner of the box'
>      'query bpos'
>      xb1o = subwrd(result,3)
>      yb2o = subwrd(result,4)
>      say 'Click where you want the right lower corner of the box'
>      'query bpos'
>      xb2o = subwrd(result,3)
>      yb1o = subwrd(result,4)
>
>      say
>      say '...left upper corner of the box at X Y: 'xb1o' 'yb2o
>      say '...right lower corner of the box at X Y: 'xb2o' 'yb1o
> *
> * draw box to show specified area for defining hurricane center...
> *
>      'set line 1 1'
>      'draw rec 'xb1o' 'yb1o' 'xb2o' 'yb2o
> *
> * convert box "xy" coordinates to "grid" coordinates
> * for finding minimum "pressure" (hurricane center)...
> *
>      'q xy2gr 'xb1o' 'yb2o
>      x1o = subwrd(result,3) ; x1 = math_int(x1o+0.5)
>      y2o = subwrd(result,6) ; y2 = math_int(y2o+0.5)
>      'q xy2gr 'xb2o' 'yb1o
>      x2o = subwrd(result,3) ; x2 = math_int(x2o+0.5)
>      y1o = subwrd(result,6) ; y1 = math_int(y1o+0.5)
>      say
>      say 'Specified grid coordinates (x1o,x2o), (y1o,y2o) from box,'
>      say 'rounded to (x1,x2), (y1,y2) for finding minimum "pressure"...'
>      say
>      say ' (x1o,x2o) = ('x1o','x2o')  -->  (x1,x2) = ('x1','x2')'
>      say ' (y1o,y2o) = ('y1o','y2o')  -->  (y1,y2) = ('y1','y2')'
> *
> * find minimum "pressure" within box area (x1,x2), (y1,y2)...
> *
>      'set x 'x1
>      'set y 'y1
>      'set z 1'
>      'set t 'tt
>      'set gxout print'
>      nxp = x2-x1+1
> *      say '    nxp = 'nxp
>      nyp = y2-y1+1
> *      say '    nyp = 'nyp
>      'd minloc(min(psl,y='y1',y='y2'),x='x1',x='x2')'
>        rec=sublin(result,nxp+3)
>        xc=subwrd(rec,1)
>      'd minloc(min(psl,x='x1',x='x2'),y='y1',y='y2')'
>        rec=sublin(result,nyp+3)
>        yc=subwrd(rec,1)
>      say
>      say 'location of minimum "pressure"...'
>      say
>      say ' (xc,yc) = ('xc','yc')'
> *
> * find "world" coordinates of (xc,yc) and convert "world" coordinates
> * to "xy" coordinates for plotting track of hurricane center...
> *
>      'set x 'xc
>      lonval = subwrd(result,4)
>      'set y 'yc
>      latval = subwrd(result,4)
>      say ' (lonval,latval) = ('lonval','latval')'
>      'q w2xy 'lonval' 'latval
>      xpos = subwrd(result,3)
>      ypos = subwrd(result,6)
>      say ' (xpos,ypos) = ('xpos','ypos')'
> *
> * write (xpos,ypos) to output file...
> *
>      res = write (''hurricane.nf'','t = 'tt' xpos = 'xpos' ypos = 'ypos'')
>
>      pause()
> *
> * erase box showing area for defining hurricane center...
> *
> *      'set line 0 1'
> *      'draw rec 'xb1o' 'yb1o' 'xb2o' 'yb2o
>      tt=tt+1
>    endwhile
>
>    res = close (''hurricane.nf'')
>
>    nf=nf+1
>  endwhile
>
> *  Clear frame and re-plot base map...
>
>  plot_base_map(plot_vpage,plot_area,lon1,lon2,lat1,lat2)
>
> *  Plot track(s)...
>
>  nf=1
>  while ( nf <= nmodels )
>
> *  Specify drawing primitives
>
>    marktype = marks.nf
>    lcolor = colors.nf
>
>    say ' marktype, marksize, lcolor, lstyle and lthick:'
>    say ' 'marktype ' ' marksize ' ' lcolor ' ' lstyle ' ' lthick
>
>    tt=ts
> *
> *  read (xpos,ypos) from ASCII text input file...
> *
>    res = read(''hurricane.nf'') ; rc = sublin(res,1)
>    if( rc != 0 )
>      say '"read" status = 'rc
>      exit
>    endif
>    rec = sublin(res,2) ; say '"hurricane.'nf'" record: 'rec
>    xc = subwrd(rec,6) ; yc = subwrd(rec,9)
>
>    'set line 'lcolor' 'lstyle' 'lthick
>    'draw mark 'marktype' 'xc' 'yc' 'marksize
>    tt = tt+1
>    while ( tt <= nt )
>      res = read(''hurricane.nf'') ; rc = sublin(res,1)
>      if( rc != 0 )
>        say '"read" status = 'rc
>        exit
>      endif
>      rec = sublin(res,2) ; say '"hurricane.'nf'" record: 'rec
>      xn = subwrd(rec,6) ; yn = subwrd(rec,9)
>      say ' 'xc' 'yc' 'xn' 'yn
>      'draw line 'xc' 'yc' 'xn' 'yn
>      'draw mark 'marktype' 'xn' 'yn' 'marksize
>      tt=tt+1
>      xc = xn
>      yc = yn
>    endwhile
>
>    nf=nf+1
>  endwhile
>
>  say
>  say '...plot legend'
>  ''legend_plotter' -x 3 -y 6 'legend_info
>  say
>
>  'print'
>  'disable print'
>
>  say
>  say '************************************************'
>  say ' Finished with this script.'
>  say '************************************************'
>  say
>
> return
>
> function plot_base_map(plot_vpage,plot_area,lon1,lon2,lat1,lat2)
>
> *  Plots base map...
>
>  'c'
>  'set vpage 'plot_vpage
>  'set parea 'plot_area
>  'set lon 'lon1' 'lon2
>  'set lat 'lat1' 'lat2
>  'set mproj latlon'
>  'set mpt * 1 1 5'
>  'set mpdset hires'
>  'set poli on'
> *  'set grid on'
>  'set grads off'
>  'set gxout contour'
>  'draw map'
>
> *  add lon/lat lines...
>
>  'set ccolor 1'
>  'set cstyle 5' ; 'd lon'
>  'set cstyle 5' ; 'd lat'
>
> *  add title, x- and y-labels...
>
>  'draw title \\ Track of Minimum Surface Pressure'
>  'draw xlab Longitude \\'
>  'draw ylab \\ Latitude'
>
> return
>
> _______________________________________________
> gradsusr mailing list
> gradsusr at gradsusr.org
> http://gradsusr.org/mailman/listinfo/gradsusr
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://gradsusr.org/pipermail/gradsusr/attachments/20110429/e2936a21/attachment-0003.html 


More information about the gradsusr mailing list