[gradsusr] how to find typhoon center via minloc function
sushant puranik
sushantpuranik at gmail.com
Wed May 4 05:47:01 EDT 2011
Hello Charles
Thank you for the script. Script for individual data file runs fine. But I
am not able to run maxwind2.gs using two separate data files uwind05.nc and
vwind05.nc (which are attached here with this mail).
One thing can we check maximum wind at low level. I want to check maximum
wind from 1000 upto 600 or 500hPa only. My study is related with low level
jet for that purpose i need maximum wind at low level only.
Thank you again
On Wed, May 4, 2011 at 1:17 AM, Charles Seman <Charles.Seman at noaa.gov>wrote:
> Sushant,
>
> Glad to hear the script was helpful.
>
> An updated script "maxwind.gs" (attached) for a single wind file and wind
> variable (e.g. ucomp or vcomp) is now able to search over a user-specified
> domain.
>
> Another script "maxwind2.gs" (attached) has been coded based on "
> maxwind.gs" to find the max wind magnitude of 2 wind variables (e.g.
> ucomp, vcomp) over a user-specified domain (although the wind magnitude is
> defined over all (x,y,z) points at each time level).
>
> The scripts were run over a subset of the total (x,y,z,t) domain to verify
> they work for a limited domain (note, you can set these in the scripts
> yourself):
> x1 = 2 ; x2 = nx-1
> y1 = 2 ; y2 = ny-1
> z1 = 2 ; z2 = nz-1
> t1 = 2 ; t2 = nt-1
>
> Run commands:
> cjs: /home/cjs/grads/Sushant_Puranik/ --> set wind_file = (
> ucomp.1981-2000.climatology.nc vcomp.1981-2000.climatology.nc )
> cjs: /home/cjs/grads/Sushant_Puranik/ --> set wind_var = ( ucomp vcomp )
> cjs: /home/cjs/grads/Sushant_Puranik/ --> set wind_units = m/s
> cjs: /home/cjs/grads/Sushant_Puranik/ --> set lev_units = hPa
> cjs: /home/cjs/grads/Sushant_Puranik/ --> grads -lcx "maxwind.gs$wind_nc_file[1] $wind_var[1] $wind_units $lev_units"
> ...
> cjs: /home/cjs/grads/Sushant_Puranik/ --> grads -lcx "maxwind.gs$wind_nc_file[2] $wind_var[2] $wind_units $lev_units"
> ...
> cjs: /home/cjs/grads/Sushant_Puranik/ --> grads -lcx "maxwind2.gs$wind_nc_file[1] $wind_var[1] $wind_nc_file[2] $wind_var[2] $wind_units
> $lev_units"
> ...
>
> Please find attached sample plot files: maxwind_ucomp.pdf,
> maxwind_vcomp.pdf, maxwind2.pdf
>
> The scripts were tested and appear to be working OK. Please let me know if
> you find any problems and/or have any questions.
>
> Thanks,
> Chuck
>
> sushant puranik wrote:
>
>> Thanks again It works fine. Since i am interested in low level wind so can
>> you please tell me how to make it to check for low levels only(i have data
>> at 17 pressure levels). Also i have uwind and vwind data in separate files,
>> how to run script with two files?
>>
>>
>> thank you
>> On Mon, May 2, 2011 at 11:04 PM, Charles Seman <Charles.Seman at noaa.gov<mailto:
>> Charles.Seman at noaa.gov>> wrote:
>>
>> Sushant,
>>
>> Following Kun-Hsuan Chou in the "hurricane_tracking.txt" file (or
>> plot_hurricane_center3.gs <http://plot_hurricane_center3.gs>),
>> please find attached a revised version "maxwind.gs
>> <http://maxwind.gs>" which plots the wind and its max for each time
>> level.
>>
>> To run:
>>
>> % grads -lcx "maxwind.gs <http://maxwind.gs> nc_file wind wind_units"
>>
>> or
>>
>> ga-> run maxwind.gs <http://maxwind.gs> nc_file wind wind_units
>>
>> Where a new argument "wind_units" has been added which is plotted in
>> a title. Be sure there are no blank spaces in the "wind_units" so
>> that GrADS gets the correct string to plot (e.g. m/s or knots,
>> etc.). Note, the script as it is clears the screen before each time
>> level (it plots the wind field itself at the level of maximum wind
>> and an "x" at the location of the max wind). I'm not sure if this
>> is what you want? If not, I would suggest adapting the code from
>> maxwind.gs <http://maxwind.gs> into a modified version of
>> plot_hurricane_center3.gs <http://plot_hurricane_center3.gs> (and
>> rename it to something more appropriate like plot_maxwind.gs
>> <http://plot_maxwind.gs> or whatever you choose). One thing I found
>> when running the script several times on the GrADS command line is
>> that at some point the script exits with a "Memory allocation
>> error!" condition. I'm not sure why this happens. In an attempt to
>> prevent this, 'c events' and 'close 1' statements were added in
>> addition to the previous 'reinit' command, but this did not prevent
>> the problem.
>>
>> I hope this helps,
>> Chuck
>>
>> sushant puranik wrote:
>>
>> Hello Charles
>> I tried maxwind.gs <http://maxwind.gs> <http://maxwind.gs> it
>> works fine as per your suggestion. But i want to know at which
>> lat lon maximum value is instead of x and y location. So can you
>> please tell me how to convert x y in lat lon.
>>
>> Thanks
>>
>> Sushant
>> On Sat, Apr 30, 2011 at 11:43 AM, sushant puranik
>> <sushantpuranik at gmail.com <mailto:sushantpuranik at gmail.com>
>> <mailto:sushantpuranik at gmail.com
>> <mailto:sushantpuranik at gmail.com>>> wrote:
>>
>> Thank you for help i will try it. Please don't mind if error
>> occurs
>> i'll mail you. Thanks again
>>
>>
>> On Sat, Apr 30, 2011 at 1:32 AM, Charles Seman
>> <Charles.Seman at noaa.gov <mailto:Charles.Seman at noaa.gov>
>> <mailto:Charles.Seman at noaa.gov <mailto:Charles.Seman at noaa.gov>>>
>> wrote:
>>
>> Hello Sushant,
>>
>> I think it should be possible. Please find attached an
>> email
>> exchange to this List from several years ago
>> "hurricane_tracking.txt" in which Kun-Hsuan Chou gave a
>> method
>> to track the maximum wind speed using maxloc and max. This
>> technique works on a 2D field. To find the location of
>> the max
>> wind in 4D wind dataset, a short script "maxwind.gs
>> <http://maxwind.gs>
>> <http://maxwind.gs>" was coded that blends the technique
>> used in
>>
>> a local script "plot_hurricane_center3.gs
>> <http://plot_hurricane_center3.gs>
>> <http://plot_hurricane_center3.gs>" and the technique in
>>
>> Kun-Hsuan Chou's "hurricane_tracking.txt".
>> ("plot_hurricane_center3.gs
>> <http://plot_hurricane_center3.gs>
>> <http://plot_hurricane_center3.gs>"
>>
>> itself was developed using Joe Covert's script
>> "plot_tc_shi.gs <http://plot_tc_shi.gs>
>> <http://plot_tc_shi.gs>" and Kun-Hsuan Chou's
>> "hurricane_tracking.txt".) "maxwind.gs
>> <http://maxwind.gs> <http://maxwind.gs>"
>>
>> loops over the time levels and finds the value of the max
>> wind
>> and its (x,y,z) location for each time level. These
>> values are
>> printed out to the screen. Here is a copy of what the
>> "maxwind.gs <http://maxwind.gs> <http://maxwind.gs>"
>> script printed out using GrADS
>>
>> 1.9b4 for a local model "ucomp" climatology dataset:
>> ---
>> ...
>> (nx,ny,nz,nt) = (144,90,23,12)
>> tt = 1
>> max ucomp = 81.833
>> xc = 137
>> yc = 71
>> zc = 23
>> ucomp(x=137,y=71,z=23) = 81.833
>> tt = 2
>> max ucomp = 68.9826
>> xc = 57
>> yc = 61
>> zc = 10
>> ucomp(x=57,y=61,z=10) = 68.9826
>> tt = 3
>> max ucomp = 62.7262
>> xc = 56
>> yc = 60
>> zc = 11
>> ucomp(x=56,y=60,z=11) = 62.7262
>> tt = 4
>> max ucomp = 61.6225
>> xc = 10
>> yc = 15
>> zc = 23
>> ucomp(x=10,y=15,z=23) = 61.6225
>> tt = 5
>> max ucomp = 85.5501
>> xc = 4
>> yc = 18
>> zc = 23
>> ucomp(x=4,y=18,z=23) = 85.5501
>> tt = 6
>> max ucomp = 109.233
>> xc = 4
>> yc = 19
>> zc = 23
>> ucomp(x=4,y=19,z=23) = 109.233
>> tt = 7
>> max ucomp = 119.698
>> xc = 143
>> yc = 17
>> zc = 23
>> ucomp(x=143,y=17,z=23) = 119.698
>> tt = 8
>> max ucomp = 117.018
>> xc = 2
>> yc = 16
>> zc = 22
>> ucomp(x=2,y=16,z=22) = 117.018
>> tt = 9
>> max ucomp = 98.7713
>> xc = 63
>> yc = 9
>> zc = 20
>> ucomp(x=63,y=9,z=20) = 98.7713
>> tt = 10
>> max ucomp = 74.1137
>> xc = 65
>> yc = 7
>> zc = 19
>> ucomp(x=65,y=7,z=19) = 74.1137
>> tt = 11
>> max ucomp = 77.7408
>> xc = 137
>> yc = 67
>> zc = 23
>> ucomp(x=137,y=67,z=23) = 77.7408
>> tt = 12
>> max ucomp = 84.3092
>> xc = 16
>> yc = 70
>> zc = 23
>> ucomp(x=16,y=70,z=23) = 84.3092
>> ...finished.
>> ---
>>
>> It seems the script is working OK. You can give it a try
>> for
>> one of your datasets. To run the script as it is coded,
>> run it
>> in batch mode as follows:
>> % grads -blcx "maxwind.gs <http://maxwind.gs>
>> <http://maxwind.gs> nc_file wind"
>>
>> or on the GrADS command line
>> http://grads.iges.org/grads/gadoc/gradcomdrun.html as (no
>> quotes
>> as with the batch command):
>> ga-> run maxwind.gs <http://maxwind.gs>
>> <http://maxwind.gs> nc_file wind
>>
>> where "nc_file" is the netCDF file containing the "wind"
>> variable name (the "nc_file" and "wind" are generic names;
>> please replace them with the actual netCDF file and wind
>> variable in the netCDF file). Note, please also include the
>> quotes surrounding "maxwind.gs <http://maxwind.gs>
>> <http://maxwind.gs> nc_file wind"
>>
>> for the batch command or else the batch command won't work.
>>
>> I hope this helps. Please let me know if you have any
>> questions.
>>
>> Thanks,
>> Chuck
>>
>> sushant puranik wrote:
>>
>> 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
>> <mailto:Charles.Seman at noaa.gov> <mailto:Charles.Seman at noaa.gov
>> <mailto:Charles.Seman at noaa.gov>>
>> <mailto:Charles.Seman at noaa.gov
>> <mailto:Charles.Seman at noaa.gov>
>> <mailto:Charles.Seman at noaa.gov
>> <mailto: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://plot_tc_shi.gs>
>> <http://plot_tc_shi.gs>
>> <http://plot_tc_shi.gs>":
>> http://caos.iisc.ernet.in/gslib/plot_tc_shi.gs
>>
>>
>> The attached GrADS script
>> plot_hurricane_center3.gs <http://plot_hurricane_center3.gs>
>> <http://plot_hurricane_center3.gs>
>> <http://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 <http://plot_tc_shi.gs>
>> <http://plot_tc_shi.gs>
>> <http://plot_tc_shi.gs>" (Joe.Covert at noaa.gov
>> <mailto:Joe.Covert at noaa.gov>
>> <mailto:Joe.Covert at noaa.gov <mailto:Joe.Covert at noaa.gov
>> >>
>> <mailto:Joe.Covert at noaa.gov
>> <mailto:Joe.Covert at noaa.gov>
>> <mailto:Joe.Covert at noaa.gov
>> <mailto:Joe.Covert at noaa.gov>>>) which was posted to the GrADS
>> Listserv
>>
>> on 9/13/2004 (his script was named plot_tc_shi.gs
>> <http://plot_tc_shi.gs>
>> <http://plot_tc_shi.gs>
>> <http://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
>> <http://plot_hurricane_center3.gs>
>> <http://plot_hurricane_center3.gs>
>> <http://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
>> <http://plot_hurricane_center3.gs>
>> <http://plot_hurricane_center3.gs>
>> <http://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
>> <http://hurricane.nf> <http://hurricane.nf>
>> <http://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>
>> <mailto:totangjie at gmail.com
>> <mailto:totangjie at gmail.com>> <mailto:totangjie at gmail.com
>> <mailto:totangjie at gmail.com>
>> <mailto:totangjie at gmail.com
>> <mailto:totangjie at gmail.com>>>
>> <mailto:totangjie at gmail.com
>> <mailto:totangjie at gmail.com>
>> <mailto:totangjie at gmail.com
>> <mailto:totangjie at gmail.com>> <mailto:totangjie at gmail.com
>> <mailto:totangjie at gmail.com>
>> <mailto:totangjie at gmail.com
>> <mailto:totangjie at gmail.com>>>>
>>
>>
>> Tel: 0086-2154896104
>> Shanghai Typhoon Institute,China
>>
>> ------------------------------------------------------------------------
>>
>> _______________________________________________
>> gradsusr mailing list
>> gradsusr at gradsusr.org
>> <mailto:gradsusr at gradsusr.org> <mailto:gradsusr at gradsusr.org
>> <mailto:gradsusr at gradsusr.org>>
>> <mailto:gradsusr at gradsusr.org
>> <mailto:gradsusr at gradsusr.org> <mailto:gradsusr at gradsusr.org
>> <mailto:gradsusr at gradsusr.org>>>
>>
>>
>> http://gradsusr.org/mailman/listinfo/gradsusr
>> --
>> Please note that Charles.Seman at noaa.gov
>> <mailto:Charles.Seman at noaa.gov>
>> <mailto:Charles.Seman at noaa.gov
>> <mailto:Charles.Seman at noaa.gov>>
>> <mailto:Charles.Seman at noaa.gov
>> <mailto:Charles.Seman at noaa.gov>
>> <mailto:Charles.Seman at noaa.gov
>> <mailto:Charles.Seman at noaa.gov>>> should be considered my NOAA
>> email address, not cjs at gfdl.noaa.gov
>> <mailto:cjs at gfdl.noaa.gov>
>> <mailto:cjs at gfdl.noaa.gov <mailto:cjs at gfdl.noaa.gov>>
>> <mailto:cjs at gfdl.noaa.gov <mailto:cjs at gfdl.noaa.gov>
>>
>> <mailto:cjs at gfdl.noaa.gov <mailto:cjs at gfdl.noaa.gov>>>.
>>
>>
>>
>> ********************************************************************
>> Charles Seman
>> Charles.Seman at noaa.gov <mailto:Charles.Seman at noaa.gov>
>> <mailto:Charles.Seman at noaa.gov <mailto:Charles.Seman at noaa.gov>>
>> <mailto:Charles.Seman at noaa.gov
>> <mailto:Charles.Seman at noaa.gov>
>> <mailto:Charles.Seman at noaa.gov
>> <mailto: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 <http://plot_tc_shi.gs>
>> <http://plot_tc_shi.gs>
>> <http://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
>> <http://plot_hurricane_center2c.gs>
>> <http://plot_hurricane_center2c.gs>
>> <http://plot_hurricane_center2c.gs>, which was
>> adapted from:
>>
>> *
>> * (1) "plot_hurricane_center2[a][b].gs" (from
>> "plot_hurricane_center.gs
>> <http://plot_hurricane_center.gs>
>> <http://plot_hurricane_center.gs>
>> <http://plot_hurricane_center.gs>"),
>>
>> * locally modified versions derived from Joe
>> Covert's script
>> * "plot_tc_shi.gs <http://plot_tc_shi.gs>
>> <http://plot_tc_shi.gs>
>> <http://plot_tc_shi.gs>" (Joe.Covert at noaa.gov
>> <mailto:Joe.Covert at noaa.gov>
>> <mailto:Joe.Covert at noaa.gov <mailto:Joe.Covert at noaa.gov
>> >>
>> <mailto:Joe.Covert at noaa.gov
>> <mailto:Joe.Covert at noaa.gov>
>> <mailto:Joe.Covert at noaa.gov
>> <mailto:Joe.Covert at noaa.gov>>>), which was posted to the
>>
>> * GrADS Listserv on 9/13/2004 (his script was
>> named
>> plot_tc_shi.gs <http://plot_tc_shi.gs>
>> <http://plot_tc_shi.gs>
>> <http://plot_tc_shi.gs>)
>>
>> * on 4/27/06, found via a Google search for
>> "plot_tc_shi.gs <http://plot_tc_shi.gs>
>> <http://plot_tc_shi.gs>
>> <http://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 <http://hurricane.nf>
>> <http://hurricane.nf>
>> <http://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 <http://cbar_line_box.gs>
>> <http://cbar_line_box.gs>
>> <http://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
>> <http://hurricane.nf> <http://hurricane.nf>
>> <http://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
>> <http://hurricane.nf> <http://hurricane.nf>
>> <http://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 <http://marks.nf>
>> <http://marks.nf> <http://marks.nf>
>> lcolor = colors.nf <http://colors.nf>
>> <http://colors.nf> <http://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 <http://hurricane.nf>
>> <http://hurricane.nf>
>> <http://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
>> <http://hurricane.nf> <http://hurricane.nf>
>> <http://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
>> <mailto:gradsusr at gradsusr.org> <mailto:gradsusr at gradsusr.org
>> <mailto:gradsusr at gradsusr.org>>
>> <mailto:gradsusr at gradsusr.org
>> <mailto:gradsusr at gradsusr.org> <mailto:gradsusr at gradsusr.org
>> <mailto:gradsusr at gradsusr.org>>>
>>
>>
>> http://gradsusr.org/mailman/listinfo/gradsusr
>>
>>
>>
>>
>>
>>
>>
>> ------------------------------------------------------------------------
>>
>> _______________________________________________
>> gradsusr mailing list
>> gradsusr at gradsusr.org <mailto:gradsusr at gradsusr.org>
>> <mailto:gradsusr at gradsusr.org <mailto:gradsusr at gradsusr.org>>
>> http://gradsusr.org/mailman/listinfo/gradsusr
>>
>>
>> --
>> Please note that Charles.Seman at noaa.gov
>> <mailto:Charles.Seman at noaa.gov>
>> <mailto:Charles.Seman at noaa.gov
>> <mailto:Charles.Seman at noaa.gov>> should be considered my NOAA
>> email address, not cjs at gfdl.noaa.gov
>> <mailto:cjs at gfdl.noaa.gov> <mailto:cjs at gfdl.noaa.gov
>> <mailto:cjs at gfdl.noaa.gov>>.
>>
>>
>> ********************************************************************
>> Charles Seman
>> Charles.Seman at noaa.gov <mailto:Charles.Seman at noaa.gov>
>> <mailto:Charles.Seman at noaa.gov <mailto: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."
>>
>> *
>> *
>> * Subject: Re: Tracking Hurricanes
>> * From: Kun-Hsuan Chou <cwujou at TYPHOON.AS.NTU.EDU.TW
>> <mailto:cwujou at TYPHOON.AS.NTU.EDU.TW>
>> <mailto:cwujou at TYPHOON.AS.NTU.EDU.TW
>> <mailto:cwujou at TYPHOON.AS.NTU.EDU.TW>>>
>>
>> * Date: Fri, 8 Oct 2004 08:03:20 +0800
>> * To: GRADSUSR at LIST.CINECA.IT
>> <mailto:GRADSUSR at LIST.CINECA.IT> <mailto:GRADSUSR at LIST.CINECA.IT
>> <mailto:GRADSUSR at LIST.CINECA.IT>>
>>
>> *
>> * Hi,
>> *
>> * I write a sample to track the maxminum wind speed, it
>> should
>> be the script that you want.
>>
>> *
>> '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'
>>
>>
>> *============================================================================
>> *Kun-Hsuan Chou
>> *Postdoctoral Researcher
>> *Department of Atmospheric Sciences
>> *National Taiwan University
>> *No. 1, Sec. 4, Roosevelt Rd., Taipei 106, Taiwan.
>> *Phone: 886-2-23625896 ext. 243
>> *Mobile: 0920079247
>> *E-mail: cwujou at typhoon.as.ntu.edu.tw
>> <mailto:cwujou at typhoon.as.ntu.edu.tw>
>> <mailto:cwujou at typhoon.as.ntu.edu.tw
>> <mailto:cwujou at typhoon.as.ntu.edu.tw>>
>>
>> *----- Original Message -----
>> *From: Arturo Caracas Uribe
>> *To: GRADSUSR at LIST.CINECA.IT
>> <mailto:GRADSUSR at LIST.CINECA.IT> <mailto:GRADSUSR at LIST.CINECA.IT
>> <mailto:GRADSUSR at LIST.CINECA.IT>>
>>
>> *Sent: Thursday, October 07, 2004 11:48 PM
>> *Subject: Tracking Hurricanes
>> *
>> *Does anyone know how to get the maximum values of
>> vorticity of
>> a domain? (the value, lat and lon).I need the positions
>> in order
>> to track a hurricane.
>> *Thanks in advance
>> *
>> *Arturo Caracas
>> *Atmospheric Sciences, México.
>>
>>
>>
>> *----------------------------------------------------------------------
>> *
>> * Adapted from plot_hurricane_center2c.gs
>> <http://plot_hurricane_center2c.gs>
>> <http://plot_hurricane_center2c.gs>, which was adapted
>> from:
>> *
>> * (1) "plot_hurricane_center2[a][b].gs" (from
>> "plot_hurricane_center.gs
>> <http://plot_hurricane_center.gs>
>> <http://plot_hurricane_center.gs>"),
>>
>> * locally modified versions derived from Joe
>> Covert's script
>> * "plot_tc_shi.gs <http://plot_tc_shi.gs>
>> <http://plot_tc_shi.gs>"
>> (Joe.Covert at noaa.gov <mailto:Joe.Covert at noaa.gov>
>> <mailto:Joe.Covert at noaa.gov <mailto:Joe.Covert at noaa.gov>>),
>> which was
>> posted to the
>> * GrADS Listserv on 9/13/2004 (his script was named
>> plot_tc_shi.gs <http://plot_tc_shi.gs>
>> <http://plot_tc_shi.gs>)
>>
>> * on 4/27/06, found via a Google search for
>> "plot_tc_shi.gs <http://plot_tc_shi.gs>
>> <http://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<
>>
>> ...
>>
>> [Message clipped] function maxwind(args)
>> 'reinit'
>> nc_file = subwrd(args,1)
>> wind = subwrd(args,2)
>> wind_units = subwrd(args,3)
>> say
>> say nc_file
>> say wind' ('wind_units')'
>> say
>> 'sdfopen 'nc_file
>> 'q file'
>> say result
>> say
>> rec5 = sublin(result,5)
>> nx = subwrd(rec5,3) ; ny = subwrd(rec5,6) ; nz = subwrd(rec5,9) ;
>> nt = subwrd(rec5,12)
>> say
>> say '(nx,ny,nz,nt) = ('nx','ny','nz','nt')'
>> say
>> nxp = nx
>> nyp = ny
>> nzp = nz
>> *
>> * get the x- and y-limits for virtual page
>> *
>> 'q gxinfo'
>> rec3 = sublin(result,3) ; xpl = subwrd(rec3,4) ; xpu = subwrd(rec3,6)
>> rec4 = sublin(result,4) ; ypl = subwrd(rec4,4) ; ypu = subwrd(rec4,6)
>> tt = 1
>> while ( tt <= nt )
>> 'set x 1'
>> 'set y 1'
>> 'set z 1'
>> 'set t 'tt
>> say 'tt = 'tt
>> 'set gxout print'
>> 'd max(max(max('wind',x=1,x='nx'),y=1,y='ny'),z=1,z='nz')'
>> rec = sublin(result,(nyp+1)*nzp+3)
>> max_wind = subwrd(rec,1)
>> say ' max 'wind' = 'max_wind
>> 'd maxloc(max(max('wind',y=1,y='ny'),z=1,z='nz'),x=1,x='nx')'
>> rec = sublin(result,(nzp+1)*nxp+3)
>> xc = subwrd(rec,1)
>> 'd maxloc(max(max('wind',x=1,x='nx'),z=1,z='nz'),y=1,y='ny')'
>> rec = sublin(result,(nzp+1)*nyp+3)
>> yc = subwrd(rec,1)
>> 'd maxloc(max(max('wind',x=1,x='nx'),y=1,y='ny'),z=1,z='nz')'
>> rec = sublin(result,(nyp+1)*nzp+3)
>> zc = subwrd(rec,1)
>> say ' (xc,yc,zc) = ('xc','yc','zc')'
>> 'd 'wind'(x='xc',y='yc',z='zc')'
>> rec = sublin(result,2)
>> wind_xc_yc_zc = subwrd(rec,1)
>> say ' 'wind'(x='xc',y='yc',z='zc') = 'wind_xc_yc_zc
>> *
>> * find "world" coordinates of (xc,yc) and convert "world" coordinates
>> * to "xy" coordinates for plotting track of hurricane center...
>> *
>> 'set x 'xc
>> lonc = subwrd(result,4)
>> 'set y 'yc
>> latc = subwrd(result,4)
>> say ' (lonc,latc,zc) = ('lonc','latc','zc')'
>> 'd 'wind'(lon='lonc',lat='latc',z='zc')'
>> rec = sublin(result,2)
>> wind_lnc_ltc_zc = subwrd(rec,1)
>> say ' 'wind'(lon='lonc',lat='latc',z='zc') = 'wind_lnc_ltc_zc
>> *
>> * plot the wind and place a mark at the max
>> *
>> 'c'
>> 'set vpage 'xpl' 'xpu' 'ypl' 'ypu
>> 'set parea 'xpl+1' 'xpu-1' 'ypl+1' 'ypu-1
>> 'set x 1 'nx
>> 'set y 1 'ny
>> 'set z 'zc
>> 'set gxout shaded'
>> 'set grads off'
>> 'd 'wind
>> 'run cbar 0.7'
>> 'draw title \ Max 'wind' = 'wind_lnc_ltc_zc' ('wind_units')\at
>> (lonc,latc,zc)=('lonc','latc','zc')'
>> 'q w2xy 'lonc' 'latc
>> xpos = subwrd(result,3)
>> ypos = subwrd(result,6)
>> say ' (xpos,ypos) = ('xpos','ypos')'
>> 'set line 1 1 6'
>> 'draw mark 6 'xpos' 'ypos' 0.2'
>> 'c events'
>> pause()
>> tt=tt+1
>> endwhile
>> say '...finished.'
>> 'close 1'
>> return
>>
>> function pause()
>> *
>> * Pauses till ENTER key is pressed
>> *
>> say
>> say '*********************************'
>> say 'Press "ENTER" key to continue...'
>> say '*********************************'
>> pull answer
>> *
>> return
>>
>>
>> _______________________________________________
>> gradsusr mailing list
>> gradsusr at gradsusr.org <mailto:gradsusr at gradsusr.org>
>> http://gradsusr.org/mailman/listinfo/gradsusr
>>
>>
>>
>>
>>
>>
>> ------------------------------------------------------------------------
>>
>> _______________________________________________
>> 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."
>
> function maxwind(args)
> 'reinit'
> wind_nc_file = subwrd(args,1)
> wind_var = subwrd(args,2)
> wind_units = subwrd(args,3)
> lev_units = subwrd(args,4)
> say
> say wind_nc_file
> say wind_var' ('wind_units')'
> say
>
> * open "wind_nc_file" and define variables for (x,y,z,t) dimension limits
>
> 'sdfopen 'wind_nc_file
> 'q file'
> rec5 = sublin(result,5)
> nx = subwrd(rec5,3) ; ny = subwrd(rec5,6) ; nz = subwrd(rec5,9) ; nt =
> subwrd(rec5,12)
> say result
> say '(nx,ny,nz,nt) = ('nx','ny','nz','nt')'
> say
>
> * define (x,y,z,t) lower and upper limits for finding max wind
>
> x1 = 2 ; x2 = nx-1
> y1 = 2 ; y2 = ny-1
> z1 = 2 ; z2 = nz-1
> t1 = 2 ; t2 = nt-1
> say '(x1,x2) = ('x1','x2')'
> say '(y1,y2) = ('y1','y2')'
> say '(z1,z2) = ('z1','z2')'
> say '(t1,t2) = ('t1','t2')'
> say
>
> * get the x- and y-limits for virtual page
>
> 'q gxinfo'
> rec3 = sublin(result,3) ; xpl = subwrd(rec3,4) ; xpu = subwrd(rec3,6)
> rec4 = sublin(result,4) ; ypl = subwrd(rec4,4) ; ypu = subwrd(rec4,6)
> xptime = (xpu-xpl)/2 ; yptime = ypu-0.4
>
> * print plots into an output file?
>
> do_print = true
>
> * plot output
>
> if ( do_print = true )
> gx_file = 'maxwind_'wind_var'.gx'
> ps_file = 'maxwind_'wind_var'.ps'
> 'enable print 'gx_file
> endif
>
> * loop over time
>
> tt = t1
> while ( tt <= t2 )
>
> * find max wind and its location
>
> 'set x 1'
> 'set y 1'
> 'set z 1'
> 'set t 'tt
> 'q time' ; ttime = subwrd(result,3)
> say 'tt = 'tt' ('ttime')'
> 'set gxout print'
> 'd max(max(max('wind_var',x='x1',x='x2'),y='y1',y='y2'),z='z1',z='z2')'
> rec = sublin(result,(y2-y1+2)*(z2-z1+1)+3)
> max_wind = subwrd(rec,1)
> say ' max 'wind_var' = 'max_wind
> 'd
> maxloc(max(max('wind_var',y='y1',y='y2'),z='z1',z='z2'),x='x1',x='x2')'
> rec = sublin(result,(z2-z1+2)*(x2-x1+1)+3)
> xc = subwrd(rec,1)
> 'd
> maxloc(max(max('wind_var',x='x1',x='x2'),z='z1',z='z2'),y='y1',y='y2')'
> rec = sublin(result,(z2-z1+2)*(y2-y1+1)+3)
> yc = subwrd(rec,1)
> 'd
> maxloc(max(max('wind_var',x='x1',x='x2'),y='y1',y='y2'),z='z1',z='z2')'
> rec = sublin(result,(y2-y1+2)*(z2-z1+1)+3)
> zc = subwrd(rec,1)
> say ' (xc,yc,zc) = ('xc','yc','zc')'
> 'd 'wind_var'(x='xc',y='yc',z='zc')'
> rec = sublin(result,2)
> wind_xc_yc_zc = subwrd(rec,1)
> say ' 'wind_var'(x='xc',y='yc',z='zc') = 'wind_xc_yc_zc
>
> * find "world" coordinates of (xc,yc) and convert "world" coordinates
> * to "xy" coordinates for plotting track of hurricane center...
>
> 'set x 'xc
> lonc = subwrd(result,4)
> 'set y 'yc
> latc = subwrd(result,4)
> 'set z 'zc
> levc = subwrd(result,4)
> say ' (lonc,latc,levc) = ('lonc','latc','levc')'
> 'd 'wind_var'(lon='lonc',lat='latc',lev='levc')'
> rec = sublin(result,2)
> wind_lnc_ltc_lvc = subwrd(rec,1)
> say ' 'wind_var'(lon='lonc',lat='latc',lev='levc') = 'wind_lnc_ltc_lvc
>
> * plot the wind at the level of the max and place a mark at the max
>
> 'c'
> 'set vpage 'xpl' 'xpu' 'ypl' 'ypu
> 'set parea 'xpl+1' 'xpu-1' 'ypl+1' 'ypu-1
> 'set x 'x1' 'x2 ; lon1 = subwrd(result,4) ; lon2 = subwrd(result,5)
> 'set y 'y1' 'y2 ; lat1 = subwrd(result,4) ; lat2 = subwrd(result,5)
> 'set z 'zc
> 'set gxout shaded'
> 'set grads off'
> 'd 'wind_var
> 'run cbar 0.7'
> 'set strsiz 0.16'
> 'set string 1 c 6'
> 'draw string 'xptime' 'yptime' 'ttime
> 'draw title 'wind_var' ('wind_units') \ max='wind_lnc_ltc_lvc' at: \
> (x,y,z)=('xc','yc','zc'), (lon,lat,lev)=('lonc'E,'latc'N,'levc''lev_units')'
> 'draw xlab plot: 'lon1'E to 'lon2'E \'
> 'draw ylab \ plot: 'lat1'N to 'lat2'N'
> 'q w2xy 'lonc' 'latc
> xpos = subwrd(result,3)
> ypos = subwrd(result,6)
> say ' (xpos,ypos) = ('xpos','ypos')'
> 'set line 1 1 6'
> 'draw mark 6 'xpos' 'ypos' 0.2'
> if ( do_print = true )
> 'print'
> endif
> 'c events'
> *** pause()
> tt=tt+1
> endwhile
> say '...finished.'
> 'close 1'
> if ( do_print = true )
> say
> 'disable print'
> '!gxps -c -i 'gx_file' -o 'ps_file
> '!ps2pdf 'ps_file
> '!rm 'gx_file' 'ps_file
> endif
> return
>
> function pause()
>
> * Pauses till ENTER key is pressed
>
> say
> say '*********************************'
> say 'Press "ENTER" key to continue...'
> say '*********************************'
> pull answer
> return
>
>
> function maxwind2(args)
> 'reinit'
> wind1_nc_file = subwrd(args,1) ; wind1_var = subwrd(args,2)
> wind2_nc_file = subwrd(args,3) ; wind2_var = subwrd(args,4)
> wind_units = subwrd(args,5)
> lev_units = subwrd(args,6)
> say
> say wind1_nc_file
> say wind1_var
> say wind2_nc_file
> say wind2_var
> say 'wind units: ('wind_units')'
> say ' lev units: ('lev_units')'
> say
>
> * open "wind1_nc_file" and "wind2_nc_file" and
> * define variables for (x,y,z,t) dimension limits
>
> 'sdfopen 'wind1_nc_file
> 'q file 1'
> rec5 = sublin(result,5)
> nx1 = subwrd(rec5,3) ; ny1 = subwrd(rec5,6) ; nz1 = subwrd(rec5,9) ; nt1 =
> subwrd(rec5,12)
> say result
> say '(nx1,ny1,nz1,nt1) = ('nx1','ny1','nz1','nt1')'
> say
> 'sdfopen 'wind2_nc_file
> 'q file 2'
> rec5 = sublin(result,5)
> nx2 = subwrd(rec5,3) ; ny2 = subwrd(rec5,6) ; nz2 = subwrd(rec5,9) ; nt2 =
> subwrd(rec5,12)
> say result
> say '(nx2,ny2,nz2,nt2) = ('nx2','ny2','nz2','nt2')'
> say
>
> * make sure (x,y,z,t) dimensions are the same for "wind1_nc_file" and
> "wind2_nc_file"
> * and, if so, define general variable names for (x,y,z,t) dimension limits
>
> if( nx2 = nx1 )
> nx = nx1
> else
> say
> say 'nx1 = 'nx1
> say 'nx2 = 'nx2
> say '...are not equal'
> say
> exit
> endif
> if( ny2 = ny1 )
> ny = ny1
> else
> say
> say 'ny1 = 'ny1
> say 'ny2 = 'ny2
> say '...are not equal'
> say
> exit
> endif
> if( nz2 = nz1 )
> nz = nz1
> else
> say
> say 'nz1 = 'nz1
> say 'nz2 = 'nz2
> say '...are not equal'
> say
> exit
> endif
> if( nt2 = nt1 )
> nt = nt1
> else
> say
> say 'nt1 = 'nt1
> say 'nt2 = 'nt2
> say '...are not equal'
> say
> exit
> endif
>
> * define (x,y,z,t) lower and upper limits for finding max wind
>
> x1 = 2 ; x2 = nx-1
> y1 = 2 ; y2 = ny-1
> z1 = 2 ; z2 = nz-1
> t1 = 2 ; t2 = nt-1
> say '(x1,x2) = ('x1','x2')'
> say '(y1,y2) = ('y1','y2')'
> say '(z1,z2) = ('z1','z2')'
> say '(t1,t2) = ('t1','t2')'
> say
>
> * get the x- and y-limits for virtual page
>
> 'q gxinfo'
> rec3 = sublin(result,3) ; xpl = subwrd(rec3,4) ; xpu = subwrd(rec3,6)
> rec4 = sublin(result,4) ; ypl = subwrd(rec4,4) ; ypu = subwrd(rec4,6)
> xptime = (xpu-xpl)/2 ; yptime = ypu-0.4
>
> * print plots into an output file?
>
> do_print = true
>
> * plot output
>
> if ( do_print = true )
> gx_file = 'maxwind2.gx'
> ps_file = 'maxwind2.ps'
> 'enable print 'gx_file
> endif
>
> * loop over time
>
> tt = t1
> while ( tt <= t2 )
>
> * use "wind1_nc_file" to define dimensions
>
> 'set dfile 1'
> 'set t 'tt
> 'q time' ; ttime = subwrd(result,3)
> say 'tt = 'tt' ('ttime')'
>
> * define wind magnitude variable for this time level
>
> 'set x 1 'nx
> 'set y 1 'ny
> 'set z 1 'nz
> 'define wind = mag('wind1_var'.1,'wind2_var'.2)'
>
> * find max wind and its location
>
> 'set x 1'
> 'set y 1'
> 'set z 1'
> 'set gxout print'
> 'd max(max(max(wind,x='x1',x='x2'),y='y1',y='y2'),z='z1',z='z2')'
> rec = sublin(result,(y2-y1+2)*(z2-z1+1)+3)
> max_wind = subwrd(rec,1)
> say ' max wind = 'max_wind
> 'd maxloc(max(max(wind,y='y1',y='y2'),z='z1',z='z2'),x='x1',x='x2')'
> rec = sublin(result,(z2-z1+2)*(x2-x1+1)+3)
> xc = subwrd(rec,1)
> 'd maxloc(max(max(wind,x='x1',x='x2'),z='z1',z='z2'),y='y1',y='y2')'
> rec = sublin(result,(z2-z1+2)*(y2-y1+1)+3)
> yc = subwrd(rec,1)
> 'd maxloc(max(max(wind,x='x1',x='x2'),y='y1',y='y2'),z='z1',z='z2')'
> rec = sublin(result,(y2-y1+2)*(z2-z1+1)+3)
> zc = subwrd(rec,1)
> say ' (xc,yc,zc) = ('xc','yc','zc')'
> 'd wind(x='xc',y='yc',z='zc')'
> rec = sublin(result,2)
> wind_xc_yc_zc = subwrd(rec,1)
> say ' wind(x='xc',y='yc',z='zc') = 'wind_xc_yc_zc
>
> * find "world" coordinates of (xc,yc) and convert "world" coordinates
> * to "xy" coordinates for plotting track of hurricane center...
>
> 'set x 'xc
> lonc = subwrd(result,4)
> 'set y 'yc
> latc = subwrd(result,4)
> 'set z 'zc
> levc = subwrd(result,4)
> say ' (lonc,latc,levc) = ('lonc','latc','levc')'
> 'd wind(lon='lonc',lat='latc',lev='levc')'
> rec = sublin(result,2)
> wind_lnc_ltc_lvc = subwrd(rec,1)
> say ' wind(lon='lonc',lat='latc',lev='levc') = 'wind_lnc_ltc_lvc
>
> * plot the wind at the level of the max and place a mark at the max
>
> 'c'
> 'set vpage 'xpl' 'xpu' 'ypl' 'ypu
> 'set parea 'xpl+1' 'xpu-1' 'ypl+1' 'ypu-1
> 'set x 'x1' 'x2 ; lon1 = subwrd(result,4) ; lon2 = subwrd(result,5)
> 'set y 'y1' 'y2 ; lat1 = subwrd(result,4) ; lat2 = subwrd(result,5)
> 'set z 'zc
> 'set gxout shaded'
> 'set grads off'
> 'd 'wind
> 'run cbar 0.7'
> 'set strsiz 0.16'
> 'set string 1 c 6'
> 'draw string 'xptime' 'yptime' 'ttime
> 'draw title Wind Magnitude ('wind_units') \ max='wind_lnc_ltc_lvc' at: \
> (x,y,z)=('xc','yc','zc'), (lon,lat,lev)=('lonc'E,'latc'N,'levc''lev_units')'
> 'draw xlab plot: 'lon1'E to 'lon2'E\'
> 'draw ylab \ plot: 'lat1'N to 'lat2'N'
> 'q w2xy 'lonc' 'latc
> xpos = subwrd(result,3)
> ypos = subwrd(result,6)
> say ' (xpos,ypos) = ('xpos','ypos')'
> 'set line 1 1 6'
> 'draw mark 6 'xpos' 'ypos' 0.2'
> if ( do_print = true )
> 'print'
> endif
> 'c events'
> 'undefine wind'
> *** pause()
> tt=tt+1
> endwhile
> say '...finished.'
> 'close 2' ; 'close 1'
> if ( do_print = true )
> say
> 'disable print'
> '!gxps -c -i 'gx_file' -o 'ps_file
> '!ps2pdf 'ps_file
> '!rm 'gx_file' 'ps_file
> endif
> return
>
> function pause()
>
> * Pauses till ENTER key is pressed
>
> say
> say '*********************************'
> say 'Press "ENTER" key to continue...'
> say '*********************************'
> pull answer
> 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/20110504/5d94e3e7/attachment-0003.html
-------------- next part --------------
A non-text attachment was scrubbed...
Name: uwind05.nc
Type: application/x-netcdf
Size: 1307636 bytes
Desc: not available
Url : http://gradsusr.org/pipermail/gradsusr/attachments/20110504/5d94e3e7/attachment-0006.nc
-------------- next part --------------
A non-text attachment was scrubbed...
Name: vwind05.nc
Type: application/x-netcdf
Size: 1307636 bytes
Desc: not available
Url : http://gradsusr.org/pipermail/gradsusr/attachments/20110504/5d94e3e7/attachment-0007.nc
More information about the gradsusr
mailing list