Hello Charles<div>I tried <a href="http://maxwind.gs">maxwind.gs</a> 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.</div>
<div><br></div><div>Thanks </div><div><br></div><div><br></div><div>Sushant <br><br><div class="gmail_quote">On Sat, Apr 30, 2011 at 11:43 AM, sushant puranik <span dir="ltr"><<a href="mailto:sushantpuranik@gmail.com">sushantpuranik@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">Thank you for help i will try it. Please don't mind if error occurs i'll mail you. <div>Thanks again</div><div>
<br><br><div class="gmail_quote"><div><div></div><div class="h5">On Sat, Apr 30, 2011 at 1:32 AM, Charles Seman <span dir="ltr"><<a href="mailto:Charles.Seman@noaa.gov" target="_blank">Charles.Seman@noaa.gov</a>></span> wrote:<br>
</div></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div><div></div><div class="h5">Hello Sushant,<br>
<br>
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 "<a href="http://maxwind.gs" target="_blank">maxwind.gs</a>" was coded that blends the technique used in a local script "<a href="http://plot_hurricane_center3.gs" target="_blank">plot_hurricane_center3.gs</a>" and the technique in Kun-Hsuan Chou's "hurricane_tracking.txt". ("<a href="http://plot_hurricane_center3.gs" target="_blank">plot_hurricane_center3.gs</a>" itself was developed using Joe Covert's script "<a href="http://plot_tc_shi.gs" target="_blank">plot_tc_shi.gs</a>" and Kun-Hsuan Chou's "hurricane_tracking.txt".) "<a href="http://maxwind.gs" target="_blank">maxwind.gs</a>" 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 "<a href="http://maxwind.gs" target="_blank">maxwind.gs</a>" script printed out using GrADS 1.9b4 for a local model "ucomp" climatology dataset:<br>
---<br>
...<br>
(nx,ny,nz,nt) = (144,90,23,12)<br>
tt = 1<br>
max ucomp = 81.833<br>
xc = 137<br>
yc = 71<br>
zc = 23<br>
ucomp(x=137,y=71,z=23) = 81.833<br>
tt = 2<br>
max ucomp = 68.9826<br>
xc = 57<br>
yc = 61<br>
zc = 10<br>
ucomp(x=57,y=61,z=10) = 68.9826<br>
tt = 3<br>
max ucomp = 62.7262<br>
xc = 56<br>
yc = 60<br>
zc = 11<br>
ucomp(x=56,y=60,z=11) = 62.7262<br>
tt = 4<br>
max ucomp = 61.6225<br>
xc = 10<br>
yc = 15<br>
zc = 23<br>
ucomp(x=10,y=15,z=23) = 61.6225<br>
tt = 5<br>
max ucomp = 85.5501<br>
xc = 4<br>
yc = 18<br>
zc = 23<br>
ucomp(x=4,y=18,z=23) = 85.5501<br>
tt = 6<br>
max ucomp = 109.233<br>
xc = 4<br>
yc = 19<br>
zc = 23<br>
ucomp(x=4,y=19,z=23) = 109.233<br>
tt = 7<br>
max ucomp = 119.698<br>
xc = 143<br>
yc = 17<br>
zc = 23<br>
ucomp(x=143,y=17,z=23) = 119.698<br>
tt = 8<br>
max ucomp = 117.018<br>
xc = 2<br>
yc = 16<br>
zc = 22<br>
ucomp(x=2,y=16,z=22) = 117.018<br>
tt = 9<br>
max ucomp = 98.7713<br>
xc = 63<br>
yc = 9<br>
zc = 20<br>
ucomp(x=63,y=9,z=20) = 98.7713<br>
tt = 10<br>
max ucomp = 74.1137<br>
xc = 65<br>
yc = 7<br>
zc = 19<br>
ucomp(x=65,y=7,z=19) = 74.1137<br>
tt = 11<br>
max ucomp = 77.7408<br>
xc = 137<br>
yc = 67<br>
zc = 23<br>
ucomp(x=137,y=67,z=23) = 77.7408<br>
tt = 12<br>
max ucomp = 84.3092<br>
xc = 16<br>
yc = 70<br>
zc = 23<br>
ucomp(x=16,y=70,z=23) = 84.3092<br>
...finished.<br>
---<br>
<br>
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:<br>
% grads -blcx "<a href="http://maxwind.gs" target="_blank">maxwind.gs</a> nc_file wind"<br>
or on the GrADS command line <a href="http://grads.iges.org/grads/gadoc/gradcomdrun.html" target="_blank">http://grads.iges.org/grads/gadoc/gradcomdrun.html</a> as (no quotes as with the batch command):<br>
ga-> run <a href="http://maxwind.gs" target="_blank">maxwind.gs</a> nc_file wind<br>
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 "<a href="http://maxwind.gs" target="_blank">maxwind.gs</a> nc_file wind" for the batch command or else the batch command won't work.<br>
<br>
I hope this helps. Please let me know if you have any questions.<br>
<br>
Thanks,<br>
Chuck<br>
<br>
sushant puranik wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div>
Hello Charles<br>
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.<br>
any suggestion?<br>
<br>
thanks<br>
<br>
Sushant <br></div><div>
On Wed, Jan 5, 2011 at 1:01 AM, Charles Seman <<a href="mailto:Charles.Seman@noaa.gov" target="_blank">Charles.Seman@noaa.gov</a> <mailto:<a href="mailto:Charles.Seman@noaa.gov" target="_blank">Charles.Seman@noaa.gov</a>>> wrote:<br>
<br>
Jie,<br>
<br>
Attached is a file File_Error_1_Grads_1.9-Plot_tc_shi.gs.txt<br>
containing the text from an email exchange to the GrADS Listserv<br>
between Joe Covert and Diane Stokes on 9/13/2004 -- on 4/27/06,<br>
found via a Google search for "<a href="http://plot_tc_shi.gs" target="_blank">plot_tc_shi.gs</a><br></div>
<<a href="http://plot_tc_shi.gs" target="_blank">http://plot_tc_shi.gs</a>>": <a href="http://caos.iisc.ernet.in/gslib/plot_tc_shi.gs" target="_blank">http://caos.iisc.ernet.in/gslib/plot_tc_shi.gs</a><div>
<br>
<br>
The attached GrADS script <a href="http://plot_hurricane_center3.gs" target="_blank">plot_hurricane_center3.gs</a><br></div>
<<a href="http://plot_hurricane_center3.gs" target="_blank">http://plot_hurricane_center3.gs</a>> (along with pause.gsf, aGrADS<div><br>
script function used by the plot script) is a locally modified<br>
version derived from Joe Covert's script "<a href="http://plot_tc_shi.gs" target="_blank">plot_tc_shi.gs</a><br></div>
<<a href="http://plot_tc_shi.gs" target="_blank">http://plot_tc_shi.gs</a>>" (<a href="mailto:Joe.Covert@noaa.gov" target="_blank">Joe.Covert@noaa.gov</a><br>
<mailto:<a href="mailto:Joe.Covert@noaa.gov" target="_blank">Joe.Covert@noaa.gov</a>>) which was posted to the GrADS Listserv<div><br>
on 9/13/2004 (his script was named <a href="http://plot_tc_shi.gs" target="_blank">plot_tc_shi.gs</a><br></div>
<<a href="http://plot_tc_shi.gs" target="_blank">http://plot_tc_shi.gs</a>>) and code from hurricane_tracking.txt (taken<div><br>
from an email exchange posted to the GrADS Listserv by Kun-Hsuan<br>
Chou and Arturo Caracas Uribe in Oct, 2004)...<br>
<br>
Please find below code from hurricane_tracking.txt (taken from an<br>
email exchange posted to the GrADS Listserv by Kun-Hsuan Chou and<br>
Arturo Caracas Uribe in Oct, 2004) illustrating a technique for<br>
finding max wind in a 2D field... and below that some code from<br>
attached script <a href="http://plot_hurricane_center3.gs" target="_blank">plot_hurricane_center3.gs</a><br></div>
<<a href="http://plot_hurricane_center3.gs" target="_blank">http://plot_hurricane_center3.gs</a>> to find the min sea-level<div><br>
pressure in a 2D field...<br>
<br>
code from Kun-Hsuan Chou to find maximum wind speed from<br>
hurricane_tracking.txt:<br>
---<br>
'd maxloc(max(mag(u,v),lon=120,lon=130),lat=15,lat=25)'<br>
line=sublin(result,2)<br>
ygrd=subwrd(line,4)<br>
'd maxloc(max(mag(u,v),lat=15,lat=25),lon=120,lon=130)'<br>
line=sublin(result,2)<br>
xgrd=subwrd(line,4)<br>
'set x 'xgrd<br>
lonval = subwrd(result,4)<br>
'set y 'ygrd<br>
latval = subwrd(result,4)<br>
'q w2xy 'lonval' 'latval<br>
xpos = subwrd(result,3)<br>
ypos = subwrd(result,6)<br>
'draw mark 1 'xpos' 'ypos' .2'<br>
---<br>
<br>
sample code from <a href="http://plot_hurricane_center3.gs" target="_blank">plot_hurricane_center3.gs</a><br></div>
<<a href="http://plot_hurricane_center3.gs" target="_blank">http://plot_hurricane_center3.gs</a>> to find hurricane center (see<div><div></div><div><br>
script for supporting code):<br>
---<br>
*<br>
* find minimum "pressure" within box area (x1,x2), (y1,y2)...<br>
*<br>
'set x 'x1<br>
'set y 'y1<br>
'set z 1'<br>
'set t 'tt<br>
'set gxout print'<br>
nxp = x2-x1+1<br>
* say ' nxp = 'nxp<br>
nyp = y2-y1+1<br>
* say ' nyp = 'nyp<br>
'd minloc(min(psl,y='y1',y='y2'),x='x1',x='x2')'<br>
rec=sublin(result,nxp+3)<br>
xc=subwrd(rec,1)<br>
'd minloc(min(psl,x='x1',x='x2'),y='y1',y='y2')'<br>
rec=sublin(result,nyp+3)<br>
yc=subwrd(rec,1)<br>
say<br>
say 'location of minimum "pressure"...'<br>
say<br>
say ' (xc,yc) = ('xc','yc')'<br>
*<br>
* find "world" coordinates of (xc,yc) and convert "world" coordinates<br>
* to "xy" coordinates for plotting track of hurricane center...<br>
*<br>
'set x 'xc<br>
lonval = subwrd(result,4)<br>
'set y 'yc<br>
latval = subwrd(result,4)<br>
say ' (lonval,latval) = ('lonval','latval')'<br>
'q w2xy 'lonval' 'latval<br>
xpos = subwrd(result,3)<br>
ypos = subwrd(result,6)<br>
say ' (xpos,ypos) = ('xpos','ypos')'<br>
*<br>
* write (xpos,ypos) to output file...<br>
*<br></div></div>
res = write (''<a href="http://hurricane.nf" target="_blank">hurricane.nf</a> <<a href="http://hurricane.nf" target="_blank">http://hurricane.nf</a>>'','t = 'tt'<div><br>
xpos = 'xpos' ypos = 'ypos'')<br>
<br>
pause()<br>
---<br>
<br>
Please let me know if you have any questions.<br>
<br>
Hope this helps,<br>
Chuck<br>
<br>
Jie TANG wrote:<br>
<br>
<br>
hi,grads folks ,<br>
I am using grads v1.9, trying to find the center of typhoon via<br>
minloc function.<br>
the key script is shown as below:<br>
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)<br>
<br>
but grads tell me that "function not found min " and when i<br>
changed my script to be :<br>
tclat=min(minloc(slpt,lon=119,lon=126), lat=19,lat=26)<br>
tclon=min(minloc(slpt, lat=19,lat=26 )lon=119,lon=126)<br>
<br>
<br>
how can i finx my scrpit ? thank you .:)<br>
-- <br>
TANG Jie<br>
Email: <a href="mailto:totangjie@gmail.com" target="_blank">totangjie@gmail.com</a> <mailto:<a href="mailto:totangjie@gmail.com" target="_blank">totangjie@gmail.com</a>><br></div>
<mailto:<a href="mailto:totangjie@gmail.com" target="_blank">totangjie@gmail.com</a> <mailto:<a href="mailto:totangjie@gmail.com" target="_blank">totangjie@gmail.com</a>>><div><br>
<br>
Tel: 0086-2154896104<br>
Shanghai Typhoon Institute,China<br>
------------------------------------------------------------------------<br>
<br>
_______________________________________________<br>
gradsusr mailing list<br></div>
<a href="mailto:gradsusr@gradsusr.org" target="_blank">gradsusr@gradsusr.org</a> <mailto:<a href="mailto:gradsusr@gradsusr.org" target="_blank">gradsusr@gradsusr.org</a>><div><br>
<a href="http://gradsusr.org/mailman/listinfo/gradsusr" target="_blank">http://gradsusr.org/mailman/listinfo/gradsusr</a><br>
<br>
<br>
-- <br>
Please note that <a href="mailto:Charles.Seman@noaa.gov" target="_blank">Charles.Seman@noaa.gov</a><br></div>
<mailto:<a href="mailto:Charles.Seman@noaa.gov" target="_blank">Charles.Seman@noaa.gov</a>> should be considered my NOAA<br>
email address, not <a href="mailto:cjs@gfdl.noaa.gov" target="_blank">cjs@gfdl.noaa.gov</a> <mailto:<a href="mailto:cjs@gfdl.noaa.gov" target="_blank">cjs@gfdl.noaa.gov</a>>.<div><br>
<br>
********************************************************************<br>
Charles Seman <a href="mailto:Charles.Seman@noaa.gov" target="_blank">Charles.Seman@noaa.gov</a><br></div>
<mailto:<a href="mailto:Charles.Seman@noaa.gov" target="_blank">Charles.Seman@noaa.gov</a>><div><br>
U.S. Department of Commerce / NOAA / OAR<br>
Geophysical Fluid Dynamics Laboratory voice: (609) 452-6547<br>
201 Forrestal Road fax: (609) 987-5063<br>
Princeton, NJ 08540-6649 <a href="http://www.gfdl.noaa.gov/~cjs/" target="_blank">http://www.gfdl.noaa.gov/~cjs/</a><br>
********************************************************************<br>
<br>
"The contents of this message are mine personally and do not reflect any<br>
official or unofficial position of the United States Federal Government,<br>
the United States Department of Commerce, or NOAA."<br>
<br>
<br>
<br>
worked great ... thanks. joe<br>
<br>
Diane Stokes wrote:<br>
<br>
> Joe,<br>
><br>
> Read works fine for me in 1.9.<br>
><br>
> I think this is the issue where 'pull' in 1.9 throws in a carriage<br>
> return. The read statement is not getting the exact filename.<br>
><br>
> Try adding:<br>
> file=sublin(file,1)<br>
> after each:<br>
> pull file<br>
><br>
> Diane<br>
><br>
> Joe Covert wrote:<br>
><br>
>> There appears to be a problem with version 1.9 of GrADS "read"<br>
function.<br>
>><br></div>
>> The <a href="http://plot_tc_shi.gs" target="_blank">plot_tc_shi.gs</a> <<a href="http://plot_tc_shi.gs" target="_blank">http://plot_tc_shi.gs</a>> script (see below)<div><div></div><div><br>
uses the read function to open the<br>
>> track data file and the program gives the subject error message<br>
"File<br>
>> Error 1".<br>
>><br>
>> jc<br>
>><br>
>><br>
>> * Script to draw an hurricane-track plot.<br>
>> * Does little error checking on the input file.<br>
>> * Assumes the input file is set up as follows:<br>
>> *<br>
>> * Line 1: Title<br>
>> * Line 2: Drawing primitives for marks: marktype size<br>
>> * Line 3: Drawing primitives for lines: color style thickness<br>
>> * Line 4: Starting hour and the interval of plotting points<br>
>> * e.g., 0 6 means that track starts at 0 hour and mark<br>
>> * will be plotted every 6 hours.<br>
>> * Rest of lines: hour long. lat.<br>
>> * e.g., 0 -70.5 25.0<br>
>> * 6 -71.8 25.2<br>
>> * :<br>
>> * :<br>
>> *<br>
>> * Also assumes that a file has been opened (any file, doesn't<br>
>> * matter -- the set command doesn't work until a file has been<br>
>> * opened).<br>
>> *<br>
>><br>
>> function main()<br>
>><br>
>> * 'clear'<br>
>><br>
>> 'open dummy.ctl'<br>
>> 'set lat 20 50'<br>
>> 'set lon -90 -30'<br>
>> 'set mpdset hires'<br>
>> 'set poli on'<br>
>> 'draw map'<br>
>><br>
>><br>
>> say 'Enter File Name: (type q to stop)'<br>
>> pull fname<br>
>><br>
>> while (fname != 'q')<br>
>><br>
>> * Read the 1st record: Title<br>
>><br>
>> ret = read(fname)<br>
>> rc = sublin(ret,1)<br>
>> if (rc>0)<br>
>> say 'File Error 1'<br>
>> return<br>
>> endif<br>
>> title = sublin(ret,2)<br>
>> say title<br>
>><br>
>> * Read the drawing primitives<br>
>><br>
>> ret = read(fname)<br>
>> rc = sublin(ret,1)<br>
>> if (rc>0)<br>
>> say 'File Error 2'<br>
>> return<br>
>> endif<br>
>> dpline = sublin(ret,2)<br>
>> marktype = subwrd(dpline,1)<br>
>> marksize = subwrd(dpline,2)<br>
>> ret = read(fname)<br>
>> rc = sublin(ret,1)<br>
>> if (rc>0)<br>
>> say 'File Error 3'<br>
>> return<br>
>> endif<br>
>> dpline = sublin(ret,2)<br>
>> lcolor = subwrd(dpline,1)<br>
>> lstyle = subwrd(dpline,2)<br>
>> lthick = subwrd(dpline,3)<br>
>> say ' marktype, marksize, lcolor, lstyle and lthick:'<br>
>> say ' 'marktype ' ' marksize ' ' lcolor ' ' lstyle ' ' lthick<br>
>><br>
>> * Read starting hour and the interval hours of plotting points<br>
>><br>
>> ret = read(fname)<br>
>> rc = sublin(ret,1)<br>
>> if (rc>0)<br>
>> say 'File Error 4'<br>
>> return<br>
>> endif<br>
>> dhour = sublin(ret,2)<br>
>> start = subwrd(dhour,1)<br>
>> jump = subwrd(dhour,2)<br>
>> say ' starting hour and the interval hours of plotting points:'<br>
>> say ' 'start' 'jump<br>
>><br>
>> * Read all data points<br>
>><br>
>> ret = read(fname)<br>
>> rc = sublin(ret,1)<br>
>> while (rc = 0)<br>
>> loc = sublin(ret,2)<br>
>> hour = subwrd(loc,1)<br>
>> dlong.hour = subwrd(loc,2)<br>
>> dlat.hour = subwrd(loc,3)<br>
>> * prompt ' hour ' hour ' are read,'<br>
>> * say ' long=' dlong.hour ' lat=' dlat.hour<br>
>> ret = read(fname)<br>
>> rc = sublin(ret,1)<br>
>> endwhile<br>
>><br>
>> if (rc!=2 & rc!=0)<br>
>> say 'File Error 5, rc=' rc<br>
>> return<br>
>> endif<br>
>><br>
>> endhour = hour<br>
>> say ' endhour=' endhour<br>
>><br>
>> * Plotting<br>
>><br>
>> 'set line 'lcolor' 'lstyle' 'lthick<br>
>> 'query w2xy 'dlong.start' 'dlat.start<br>
>> xprev = subwrd(result,3)<br>
>> yprev = subwrd(result,6)<br>
>> 'draw mark 'marktype' 'xprev' 'yprev' 'marksize<br>
>> next = start+jump<br>
>> while (next <= endhour)<br>
>> * say ' 'dlong.start' 'dlat.start<br>
>> 'query w2xy 'dlong.next' 'dlat.next<br>
>> xnext = subwrd(result,3)<br>
>> ynext = subwrd(result,6)<br>
>> 'draw line 'xprev' 'yprev' 'xnext' 'ynext<br>
>> * say ' 'xprev' 'yprev' 'xnext' 'ynext<br>
>> 'draw mark 'marktype' 'xnext' 'ynext' 'marksize<br>
>> next = next+jump<br>
>> xprev = xnext<br>
>> yprev = ynext<br>
>> endwhile<br>
>><br>
>> say fname ' is working fine.'<br>
>><br>
>> * read in the filename to be plotted next<br>
>><br>
>> say 'Enter File Name: (type q to stop)'<br>
>> pull fname<br>
>><br>
>> endwhile<br>
>><br>
<br>
<br>
*----------------------------------------------------------------------<br>
*<br>
* Adapted from <a href="http://plot_hurricane_center2c.gs" target="_blank">plot_hurricane_center2c.gs</a><br></div></div>
<<a href="http://plot_hurricane_center2c.gs" target="_blank">http://plot_hurricane_center2c.gs</a>>, which was adapted from:<div><br>
*<br>
* (1) "plot_hurricane_center2[a][b].gs" (from<br></div>
"<a href="http://plot_hurricane_center.gs" target="_blank">plot_hurricane_center.gs</a> <<a href="http://plot_hurricane_center.gs" target="_blank">http://plot_hurricane_center.gs</a>>"),<div>
<br>
* locally modified versions derived from Joe Covert's script<br></div>
* "<a href="http://plot_tc_shi.gs" target="_blank">plot_tc_shi.gs</a> <<a href="http://plot_tc_shi.gs" target="_blank">http://plot_tc_shi.gs</a>>" (<a href="mailto:Joe.Covert@noaa.gov" target="_blank">Joe.Covert@noaa.gov</a><br>
<mailto:<a href="mailto:Joe.Covert@noaa.gov" target="_blank">Joe.Covert@noaa.gov</a>>), which was posted to the<div><br>
* GrADS Listserv on 9/13/2004 (his script was named<br></div>
<a href="http://plot_tc_shi.gs" target="_blank">plot_tc_shi.gs</a> <<a href="http://plot_tc_shi.gs" target="_blank">http://plot_tc_shi.gs</a>>)<div><br>
* on 4/27/06, found via a Google search for "<a href="http://plot_tc_shi.gs" target="_blank">plot_tc_shi.gs</a><br></div>
<<a href="http://plot_tc_shi.gs" target="_blank">http://plot_tc_shi.gs</a>>":<div><div></div><div><br>
* <a href="http://caos.iisc.ernet.in/gslib/plot_tc_shi.gs" target="_blank">http://caos.iisc.ernet.in/gslib/plot_tc_shi.gs</a><br>
*<br>
* (2) code in "hurricane_tracking.txt" (taken from an email exchange<br>
* posted to the GrADS Listserv by Kun-Hsuan Chou and<br>
* Arturo Caracas Uribe in Oct, 2004)<br>
*<br>
*<br>
* To use: grads -l<br>
*<br>
*----------------------------------------------------------------------<br>
<br>
function main()<br>
<br>
'reinit'<br>
<br>
************************************************************************<br>
* Allow external GrADS functions<br>
************************************************************************<br>
<br>
rc = gsfallow("on")<br>
<br>
************************************************************************<br>
* Define dataset file information and plotting limits for base map...<br>
************************************************************************<br>
<br>
** nc_file1 = ''<br>
** nc_file2 = ''<br>
** nc_file3 = ''<br>
** nc_file = ''nc_file1' 'nc_file2' 'nc_file3''<br>
nc_file = 'your_netCDF_file'<br>
<br>
lon1 = 298 ; lon2 = 308<br>
lat1 = 20 ; lat2 = 25<br>
<br>
************************************************************************<br>
* Legend captions, colors, and marks for the plots<br>
************************************************************************<br>
<br>
** nmodels = 1<br>
** model.1 = 'isabel_ras'<br>
** model.2 = 'isabel_isotke'<br>
** model.3 = 'isabel_my25'<br>
<br>
** titles.1 = model.1 ; colors.1 = 2 ; marks.1 = 2<br>
** titles.2 = model.2 ; colors.2 = 3 ; marks.2 = 2<br>
** titles.3 = model.3 ; colors.3 = 4 ; marks.3 = 2<br>
***<br>
*** Define the legend plotter input...<br>
***<br>
** titles = '-t' ; colors = '-c' ; marks = '-m' ; lines = '-l'<br>
** n=1<br>
** while ( n <= nmodels )<br>
** titles = ''titles' "'titles.n'"'<br>
** colors = ''colors' 'colors.n<br>
** marks = ''marks' "'marks.n'"'<br>
** lines = ''lines' 1'<br>
** n=n+1<br>
** endwhile<br>
<br>
** legend_info = ''colors' 'lines' 'marks' 'titles' -p'<br>
<br>
nmodels = 1<br>
model.1 = 'title_name_for_your_netCDF_dataset'<br>
<br>
titles.1 = model.1 ; colors.1 = 2 ; marks.1 = 2<br>
*<br>
* Define the legend plotter input...<br>
*<br>
titles = '-t' ; colors = '-c' ; marks = '-m' ; lines = '-l'<br>
n=1<br>
while ( n <= nmodels )<br>
titles = ''titles' "'titles.n'"'<br>
colors = ''colors' 'colors.n<br>
marks = ''marks' "'marks.n'"'<br>
lines = ''lines' 1'<br>
n=n+1<br>
endwhile<br>
<br>
legend_info = ''colors' 'lines' 'marks' 'titles' -p'<br>
<br>
* ...other drawing primitives<br>
<br>
marksize = 0.1<br>
<br>
lstyle = 1 ; lthick = 5<br>
<br>
************************************************************************<br>
* Use "cbar_line" or "cbar_line_box"?<br>
************************************************************************<br>
<br>
legend_plotter = 'cbar_line'<br>
* legend_plotter = 'cbar_line_box'<br>
<br>
************************************************************************<br>
* Plot page limits<br>
************************************************************************<br>
<br>
plot_vpage = '0 11 0 8.5'<br>
plot_area = '1 10 1 7.5'<br>
<br>
xl = subwrd(plot_area,1) ; xr = subwrd(plot_area,2)<br>
yb = subwrd(plot_area,3) ; yt = subwrd(plot_area,4)<br>
<br>
xc = xl + (xr-xl)/2<br>
<br>
************************************************************************<br>
* GrADS metafile output...<br>
************************************************************************<br>
<br>
'enable print hurricane_track.gx'<br>
<br>
*-----------------------------------------------------------------------<br>
* Define and plot track of minimum surface pressure for each dataset...<br>
*-----------------------------------------------------------------------<br>
<br>
nf=1<br>
while ( nf <= nmodels )<br>
<br>
* Open dataset file and define parameters...<br>
<br>
fname = subwrd(nc_file,nf)<br>
'sdfopen 'fname<br>
'set dfile 'nf<br>
'q file'<br>
rec3 = sublin(result,3) ; binary_file = subwrd(rec3,2)<br>
say<br>
say<br>
say '*** 'binary_file' ***'<br>
say<br>
rec5 = sublin(result,5)<br>
nx = subwrd(rec5,3) ; ny = subwrd(rec5,6) ; nz = subwrd(rec5,9) ;<br>
nt = subwrd(rec5,12)<br>
<br>
if( nf = 1 )<br>
<br>
* Plot base map...<br>
<br>
plot_base_map(plot_vpage,plot_area,lon1,lon2,lat1,lat2)<br>
<br>
endif<br>
<br>
* Define locations of "hurricane center" for plotting points, and write<br>
* locations out to ASCII text file for read-in for plotting track<br>
later...<br>
<br></div></div>
<a href="http://hurricane.nf" target="_blank">hurricane.nf</a> <<a href="http://hurricane.nf" target="_blank">http://hurricane.nf</a>> = 'hurricane_'nf'.txt'<div><br>
<br>
ts=2<br>
tt=ts<br>
while ( tt <= nt )<br>
'set t 'tt<br>
say<br>
say '---------'<br>
say ' tt = 'tt<br>
say '---------'<br>
say<br>
'set lon 'lon1' 'lon2<br>
'set lat 'lat1' 'lat2<br>
'set z 1'<br>
*<br>
* here is where you put in the name of your sea-level pressure variable:<br>
*<br>
'define psl = name_of_your_sea_level_pressure_variable'<br>
'set gxout shaded'<br>
'd psl'<br>
*<br>
* following interactive "box location" code adapted from<br></div>
"<a href="http://cbar_line_box.gs" target="_blank">cbar_line_box.gs</a> <<a href="http://cbar_line_box.gs" target="_blank">http://cbar_line_box.gs</a>>"...<div><div></div><div><br>
*<br>
say 'Click where you want the left upper corner of the box'<br>
'query bpos'<br>
xb1o = subwrd(result,3)<br>
yb2o = subwrd(result,4)<br>
say 'Click where you want the right lower corner of the box'<br>
'query bpos'<br>
xb2o = subwrd(result,3)<br>
yb1o = subwrd(result,4)<br>
<br>
say<br>
say '...left upper corner of the box at X Y: 'xb1o' 'yb2o<br>
say '...right lower corner of the box at X Y: 'xb2o' 'yb1o<br>
*<br>
* draw box to show specified area for defining hurricane center...<br>
*<br>
'set line 1 1'<br>
'draw rec 'xb1o' 'yb1o' 'xb2o' 'yb2o<br>
*<br>
* convert box "xy" coordinates to "grid" coordinates<br>
* for finding minimum "pressure" (hurricane center)...<br>
*<br>
'q xy2gr 'xb1o' 'yb2o<br>
x1o = subwrd(result,3) ; x1 = math_int(x1o+0.5)<br>
y2o = subwrd(result,6) ; y2 = math_int(y2o+0.5)<br>
'q xy2gr 'xb2o' 'yb1o<br>
x2o = subwrd(result,3) ; x2 = math_int(x2o+0.5)<br>
y1o = subwrd(result,6) ; y1 = math_int(y1o+0.5)<br>
say<br>
say 'Specified grid coordinates (x1o,x2o), (y1o,y2o) from box,'<br>
say 'rounded to (x1,x2), (y1,y2) for finding minimum "pressure"...'<br>
say<br>
say ' (x1o,x2o) = ('x1o','x2o') --> (x1,x2) = ('x1','x2')'<br>
say ' (y1o,y2o) = ('y1o','y2o') --> (y1,y2) = ('y1','y2')'<br>
*<br>
* find minimum "pressure" within box area (x1,x2), (y1,y2)...<br>
*<br>
'set x 'x1<br>
'set y 'y1<br>
'set z 1'<br>
'set t 'tt<br>
'set gxout print'<br>
nxp = x2-x1+1<br>
* say ' nxp = 'nxp<br>
nyp = y2-y1+1<br>
* say ' nyp = 'nyp<br>
'd minloc(min(psl,y='y1',y='y2'),x='x1',x='x2')'<br>
rec=sublin(result,nxp+3)<br>
xc=subwrd(rec,1)<br>
'd minloc(min(psl,x='x1',x='x2'),y='y1',y='y2')'<br>
rec=sublin(result,nyp+3)<br>
yc=subwrd(rec,1)<br>
say<br>
say 'location of minimum "pressure"...'<br>
say<br>
say ' (xc,yc) = ('xc','yc')'<br>
*<br>
* find "world" coordinates of (xc,yc) and convert "world" coordinates<br>
* to "xy" coordinates for plotting track of hurricane center...<br>
*<br>
'set x 'xc<br>
lonval = subwrd(result,4)<br>
'set y 'yc<br>
latval = subwrd(result,4)<br>
say ' (lonval,latval) = ('lonval','latval')'<br>
'q w2xy 'lonval' 'latval<br>
xpos = subwrd(result,3)<br>
ypos = subwrd(result,6)<br>
say ' (xpos,ypos) = ('xpos','ypos')'<br>
*<br>
* write (xpos,ypos) to output file...<br>
*<br></div></div>
res = write (''<a href="http://hurricane.nf" target="_blank">hurricane.nf</a> <<a href="http://hurricane.nf" target="_blank">http://hurricane.nf</a>>'','t = 'tt'<div>
<br>
xpos = 'xpos' ypos = 'ypos'')<br>
<br>
pause()<br>
*<br>
* erase box showing area for defining hurricane center...<br>
*<br>
* 'set line 0 1'<br>
* 'draw rec 'xb1o' 'yb1o' 'xb2o' 'yb2o<br>
tt=tt+1<br>
endwhile<br>
<br></div>
res = close (''<a href="http://hurricane.nf" target="_blank">hurricane.nf</a> <<a href="http://hurricane.nf" target="_blank">http://hurricane.nf</a>>'')<div><br>
<br>
nf=nf+1<br>
endwhile<br>
<br>
* Clear frame and re-plot base map...<br>
<br>
plot_base_map(plot_vpage,plot_area,lon1,lon2,lat1,lat2)<br>
<br>
* Plot track(s)...<br>
<br>
nf=1<br>
while ( nf <= nmodels )<br>
<br>
* Specify drawing primitives<br>
<br></div>
marktype = <a href="http://marks.nf" target="_blank">marks.nf</a> <<a href="http://marks.nf" target="_blank">http://marks.nf</a>><br>
lcolor = <a href="http://colors.nf" target="_blank">colors.nf</a> <<a href="http://colors.nf" target="_blank">http://colors.nf</a>><div><br>
<br>
say ' marktype, marksize, lcolor, lstyle and lthick:'<br>
say ' 'marktype ' ' marksize ' ' lcolor ' ' lstyle ' ' lthick<br>
<br>
tt=ts<br>
*<br>
* read (xpos,ypos) from ASCII text input file...<br>
*<br></div>
res = read(''<a href="http://hurricane.nf" target="_blank">hurricane.nf</a> <<a href="http://hurricane.nf" target="_blank">http://hurricane.nf</a>>'') ; rc =<div><br>
sublin(res,1)<br>
if( rc != 0 )<br>
say '"read" status = 'rc<br>
exit<br>
endif<br>
rec = sublin(res,2) ; say '"hurricane.'nf'" record: 'rec<br>
xc = subwrd(rec,6) ; yc = subwrd(rec,9)<br>
<br>
'set line 'lcolor' 'lstyle' 'lthick<br>
'draw mark 'marktype' 'xc' 'yc' 'marksize<br>
tt = tt+1<br>
while ( tt <= nt )<br></div>
res = read(''<a href="http://hurricane.nf" target="_blank">hurricane.nf</a> <<a href="http://hurricane.nf" target="_blank">http://hurricane.nf</a>>'') ; rc =<div><div></div><div>
<br>
sublin(res,1)<br>
if( rc != 0 )<br>
say '"read" status = 'rc<br>
exit<br>
endif<br>
rec = sublin(res,2) ; say '"hurricane.'nf'" record: 'rec<br>
xn = subwrd(rec,6) ; yn = subwrd(rec,9)<br>
say ' 'xc' 'yc' 'xn' 'yn<br>
'draw line 'xc' 'yc' 'xn' 'yn<br>
'draw mark 'marktype' 'xn' 'yn' 'marksize<br>
tt=tt+1<br>
xc = xn<br>
yc = yn<br>
endwhile<br>
<br>
nf=nf+1<br>
endwhile<br>
<br>
say<br>
say '...plot legend'<br>
''legend_plotter' -x 3 -y 6 'legend_info<br>
say<br>
<br>
'print'<br>
'disable print'<br>
<br>
say<br>
say '************************************************'<br>
say ' Finished with this script.'<br>
say '************************************************'<br>
say<br>
<br>
return<br>
<br>
function plot_base_map(plot_vpage,plot_area,lon1,lon2,lat1,lat2)<br>
<br>
* Plots base map...<br>
<br>
'c'<br>
'set vpage 'plot_vpage<br>
'set parea 'plot_area<br>
'set lon 'lon1' 'lon2<br>
'set lat 'lat1' 'lat2<br>
'set mproj latlon'<br>
'set mpt * 1 1 5'<br>
'set mpdset hires'<br>
'set poli on'<br>
* 'set grid on'<br>
'set grads off'<br>
'set gxout contour'<br>
'draw map'<br>
<br>
* add lon/lat lines...<br>
<br>
'set ccolor 1'<br>
'set cstyle 5' ; 'd lon'<br>
'set cstyle 5' ; 'd lat'<br>
<br>
* add title, x- and y-labels...<br>
<br>
'draw title \\ Track of Minimum Surface Pressure'<br>
'draw xlab Longitude \\'<br>
'draw ylab \\ Latitude'<br>
<br>
return<br>
<br>
_______________________________________________<br>
gradsusr mailing list<br></div></div>
<a href="mailto:gradsusr@gradsusr.org" target="_blank">gradsusr@gradsusr.org</a> <mailto:<a href="mailto:gradsusr@gradsusr.org" target="_blank">gradsusr@gradsusr.org</a>><div><br>
<a href="http://gradsusr.org/mailman/listinfo/gradsusr" target="_blank">http://gradsusr.org/mailman/listinfo/gradsusr</a><br>
<br>
<br>
<br>
<br>
<br>
<br>
------------------------------------------------------------------------<br>
<br>
_______________________________________________<br>
gradsusr mailing list<br>
<a href="mailto:gradsusr@gradsusr.org" target="_blank">gradsusr@gradsusr.org</a><br>
<a href="http://gradsusr.org/mailman/listinfo/gradsusr" target="_blank">http://gradsusr.org/mailman/listinfo/gradsusr</a><br>
</div></blockquote><div><div></div><div>
<br>
-- <br>
<br>
Please note that <a href="mailto:Charles.Seman@noaa.gov" target="_blank">Charles.Seman@noaa.gov</a> should be considered my NOAA<br>
email address, not <a href="mailto:cjs@gfdl.noaa.gov" target="_blank">cjs@gfdl.noaa.gov</a>.<br>
<br>
********************************************************************<br>
Charles Seman <a href="mailto:Charles.Seman@noaa.gov" target="_blank">Charles.Seman@noaa.gov</a><br>
U.S. Department of Commerce / NOAA / OAR<br>
Geophysical Fluid Dynamics Laboratory voice: (609) 452-6547<br>
201 Forrestal Road fax: (609) 987-5063<br>
Princeton, NJ 08540-6649 <a href="http://www.gfdl.noaa.gov/~cjs/" target="_blank">http://www.gfdl.noaa.gov/~cjs/</a><br>
********************************************************************<br>
<br>
"The contents of this message are mine personally and do not reflect any<br>
official or unofficial position of the United States Federal Government,<br>
the United States Department of Commerce, or NOAA."<br>
</div></div><br></div></div>*<br>
*<br>
* Subject: Re: Tracking Hurricanes<br>
* From: Kun-Hsuan Chou <<a href="mailto:cwujou@TYPHOON.AS.NTU.EDU.TW" target="_blank">cwujou@TYPHOON.AS.NTU.EDU.TW</a>><br>
* Date: Fri, 8 Oct 2004 08:03:20 +0800<br>
* To: <a href="mailto:GRADSUSR@LIST.CINECA.IT" target="_blank">GRADSUSR@LIST.CINECA.IT</a><br>
*<br>
* Hi,<br>
*<br>
* I write a sample to track the maxminum wind speed, it should be the script that you want.<div class="im"><br>
*<br>
'd maxloc(max(mag(u,v),lon=120,lon=130),lat=15,lat=25)'<br>
line=sublin(result,2)<br>
ygrd=subwrd(line,4)<br>
'd maxloc(max(mag(u,v),lat=15,lat=25),lon=120,lon=130)'<br>
line=sublin(result,2)<br>
xgrd=subwrd(line,4)<br>
'set x 'xgrd<br>
lonval = subwrd(result,4)<br>
'set y 'ygrd<br>
latval = subwrd(result,4)<br>
'q w2xy 'lonval' 'latval<br>
xpos = subwrd(result,3)<br>
ypos = subwrd(result,6)<br>
'draw mark 1 'xpos' 'ypos' .2'<br>
<br></div>
*============================================================================<br>
*Kun-Hsuan Chou<br>
*Postdoctoral Researcher<br>
*Department of Atmospheric Sciences<br>
*National Taiwan University<br>
*No. 1, Sec. 4, Roosevelt Rd., Taipei 106, Taiwan.<br>
*Phone: 886-2-23625896 ext. 243<br>
*Mobile: 0920079247<br>
*E-mail: <a href="mailto:cwujou@typhoon.as.ntu.edu.tw" target="_blank">cwujou@typhoon.as.ntu.edu.tw</a><br>
*----- Original Message -----<br>
*From: Arturo Caracas Uribe<br>
*To: <a href="mailto:GRADSUSR@LIST.CINECA.IT" target="_blank">GRADSUSR@LIST.CINECA.IT</a><br>
*Sent: Thursday, October 07, 2004 11:48 PM<br>
*Subject: Tracking Hurricanes<br>
*<br>
*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.<br>
*Thanks in advance<br>
*<br>
*Arturo Caracas<br>
*Atmospheric Sciences, México.<br>
<br>
<br>*----------------------------------------------------------------------<br>
*<br>
* Adapted from <a href="http://plot_hurricane_center2c.gs" target="_blank">plot_hurricane_center2c.gs</a>, which was adapted from:<br>
*<br>
* (1) "plot_hurricane_center2[a][b].gs" (from "<a href="http://plot_hurricane_center.gs" target="_blank">plot_hurricane_center.gs</a>"),<div class="im"><br>
* locally modified versions derived from Joe Covert's script<br></div>
* "<a href="http://plot_tc_shi.gs" target="_blank">plot_tc_shi.gs</a>" (<a href="mailto:Joe.Covert@noaa.gov" target="_blank">Joe.Covert@noaa.gov</a>), which was posted to the<br>
* GrADS Listserv on 9/13/2004 (his script was named <a href="http://plot_tc_shi.gs" target="_blank">plot_tc_shi.gs</a>)<div class="im"><br>
* on 4/27/06, found via a Google search for "<a href="http://plot_tc_shi.gs" target="_blank">plot_tc_shi.gs</a>":<br></div><div><div></div><div class="h5">
* <a href="http://caos.iisc.ernet.in/gslib/plot_tc_shi.gs" target="_blank">http://caos.iisc.ernet.in/gslib/plot_tc_shi.gs</a><br>
*<br>
* (2) code in "hurricane_tracking.txt" (taken from an email exchange<br>
* posted to the GrADS Listserv by Kun-Hsuan Chou and<br>
* Arturo Caracas Uribe in Oct, 2004)<br>
*<br>
*<br>
* To use: grads -l<br>
*<br>
*----------------------------------------------------------------------<br>
<br>
function main()<br>
<br>
'reinit'<br>
<br>
************************************************************************<br>
* Allow external GrADS functions<br>
************************************************************************<br>
<br>
rc = gsfallow("on")<br>
<br>
************************************************************************<br>
* Define dataset file information and plotting limits for base map...<br>
************************************************************************<br>
<br>
** nc_file1 = ''<br>
** nc_file2 = ''<br>
** nc_file3 = ''<br>
** nc_file = ''nc_file1' 'nc_file2' 'nc_file3''<br>
nc_file = 'your_netCDF_file'<br>
<br>
lon1 = 298 ; lon2 = 308<br>
lat1 = 20 ; lat2 = 25<br>
<br>
************************************************************************<br>
* Legend captions, colors, and marks for the plots<br>
************************************************************************<br>
<br>
** nmodels = 1<br>
** model.1 = 'isabel_ras'<br>
** model.2 = 'isabel_isotke'<br>
** model.3 = 'isabel_my25'<br>
<br>
** titles.1 = model.1 ; colors.1 = 2 ; marks.1 = 2<br>
** titles.2 = model.2 ; colors.2 = 3 ; marks.2 = 2<br>
** titles.3 = model.3 ; colors.3 = 4 ; marks.3 = 2<br>
***<br>
*** Define the legend plotter input...<br>
***<br>
** titles = '-t' ; colors = '-c' ; marks = '-m' ; lines = '-l'<br>
** n=1<br>
** while ( n <= nmodels )<br>
** titles = ''titles' "'titles.n'"'<br>
** colors = ''colors' 'colors.n<br>
** marks = ''marks' "'marks.n'"'<br>
** lines = ''lines' 1'<br>
** n=n+1<br>
** endwhile<br>
<br>
** legend_info = ''colors' 'lines' 'marks' 'titles' -p'<br>
<br>
nmodels = 1<br>
model.1 = 'title_name_for_your_netCDF_dataset'<br>
<br>
titles.1 = model.1 ; colors.1 = 2 ; marks.1 = 2<br>
*<br>
* Define the legend plotter input...<br>
*<br>
titles = '-t' ; colors = '-c' ; marks = '-m' ; lines = '-l'<br>
n=1<br>
while ( n <= nmodels )<br>
titles = ''titles' "'titles.n'"'<br>
colors = ''colors' 'colors.n<br>
marks = ''marks' "'marks.n'"'<br>
lines = ''lines' 1'<br>
n=n+1<br>
endwhile<br>
<br>
legend_info = ''colors' 'lines' 'marks' 'titles' -p'<br>
<br>
* ...other drawing primitives<br>
<br>
marksize = 0.1<br>
<br>
lstyle = 1 ; lthick = 5<br>
<br>
************************************************************************<br>
* Use "cbar_line" or "cbar_line_box"?<br>
************************************************************************<br>
<br>
legend_plotter = 'cbar_line'<br>
* legend_plotter = 'cbar_line_box'<br>
<br>
************************************************************************<br>
* Plot page limits<br>
************************************************************************<br>
<br>
plot_vpage = '0 11 0 8.5'<br>
plot_area = '1 10 1 7.5'<br>
<br>
xl = subwrd(plot_area,1) ; xr = subwrd(plot_area,2)<br>
yb = subwrd(plot_area,3) ; yt = subwrd(plot_area,4)<br>
<br>
xc = xl + (xr-xl)/2<br>
<br>
************************************************************************<br>
* GrADS metafile output...<br>
************************************************************************<br>
<br>
'enable print hurricane_track.gx'<br>
<br>
*-----------------------------------------------------------------------<br>
* Define and plot track of minimum surface pressure for each dataset...<br>
*-----------------------------------------------------------------------<br>
<br>
nf=1<br>
while ( nf <= nmodels )<br>
<br>
* Open dataset file and define parameters...<br>
<br>
fname = subwrd(nc_file,nf)<br>
'sdfopen 'fname<br>
'set dfile 'nf<br>
'q file'<br>
rec3 = sublin(result,3) ; binary_file = subwrd(rec3,2)<br>
say<br>
say<br>
say '*** 'binary_file' ***'<br>
say<br>
rec5 = sublin(result,5)<br>
nx = subwrd(rec5,3) ; ny = subwrd(rec5,6) ; nz = subwrd(rec5,9) ; nt = subwrd(rec5,12)<br>
<br>
if( nf = 1 )<br>
<br>
* Plot base map...<br>
<br>
plot_base_map(plot_vpage,plot_area,lon1,lon2,lat1,lat2)<br>
<br>
endif<br>
<br>
* Define locations of "hurricane center" for plotting points, and write<br>
* locations out to ASCII text file for read-in for plotting track later...<br>
<br></div></div>
<a href="http://hurricane.nf" target="_blank">hurricane.nf</a> = 'hurricane_'nf'.txt'<div class="im"><br>
<br>
ts=2<br>
tt=ts<br>
while ( tt <= nt )<br>
'set t 'tt<br>
say<br>
say '---------'<br>
say ' tt = 'tt<br>
say '---------'<br>
say<br>
'set lon 'lon1' 'lon2<br>
'set lat 'lat1' 'lat2<br>
'set z 1'<br>
*<br>
* here is where you put in the name of your sea-level pressure variable:<br>
*<br>
'define psl = name_of_your_sea_level_pressure_variable'<br>
'set gxout shaded'<br>
'd psl'<br>
*<br></div>
* following interactive "box location" code adapted from "<a href="http://cbar_line_box.gs" target="_blank">cbar_line_box.gs</a>"...<div><div></div><div class="h5"><br>
*<br>
say 'Click where you want the left upper corner of the box'<br>
'query bpos'<br>
xb1o = subwrd(result,3)<br>
yb2o = subwrd(result,4)<br>
say 'Click where you want the right lower corner of the box'<br>
'query bpos'<br>
xb2o = subwrd(result,3)<br>
yb1o = subwrd(result,4)<br>
<br>
say<br>
say '...left upper corner of the box at X Y: 'xb1o' 'yb2o<br>
say '...right lower corner of the box at X Y: 'xb2o' 'yb1o<br>
*<br>
* draw box to show specified area for defining hurricane center...<br>
*<br>
'set line 1 1'<br>
'draw rec 'xb1o' 'yb1o' 'xb2o' 'yb2o<br>
*<br>
* convert box "xy" coordinates to "grid" coordinates<br>
* for finding minimum "pressure" (hurricane center)...<br>
*<br>
'q xy2gr 'xb1o' 'yb2o<br>
x1o = subwrd(result,3) ; x1 = math_int(x1o+0.5)<br>
y2o = subwrd(result,6) ; y2 = math_int(y2o+0.5)<br>
'q xy2gr 'xb2o' 'yb1o<br>
x2o = subwrd(result,3) ; x2 = math_int(x2o+0.5)<br>
y1o = subwrd(result,6) ; y1 = math_int(y1o+0.5)<br>
say<br>
say 'Specified grid coordinates (x1o,x2o), (y1o,y2o) from box,'<br>
say 'rounded to (x1,x2), (y1,y2) for finding minimum "pressure"...'<br>
say<br>
say ' (x1o,x2o) = ('x1o','x2o') --> (x1,x2) = ('x1','x2')'<br>
say ' (y1o,y2o) = ('y1o','y2o') --> (y1,y2) = ('y1','y2')'<br>
*<br>
* find minimum "pressure" within box area (x1,x2), (y1,y2)...<br>
*<br>
'set x 'x1<br>
'set y 'y1<br>
'set z 1'<br>
'set t 'tt<br>
'set gxout print'<br>
nxp = x2-x1+1<br>
* say ' nxp = 'nxp<br>
nyp = y2-y1+1<br>
* say ' nyp = 'nyp<br>
'd minloc(min(psl,y='y1',y='y2'),x='x1',x='x2')'<br>
rec=sublin(result,nxp+3)<br>
xc=subwrd(rec,1)<br>
'd minloc(min(psl,x='x1',x='x2'),y='y1',y='y2')'<br>
rec=sublin(result,nyp+3)<br>
yc=subwrd(rec,1)<br>
say<br>
say 'location of minimum "pressure"...'<br>
say<br>
say ' (xc,yc) = ('xc','yc')'<br>
*<br>
* find "world" coordinates of (xc,yc) and convert "world" coordinates<br>
* to "xy" coordinates for plotting track of hurricane center...<br>
*<br>
'set x 'xc<br>
lonval = subwrd(result,4)<br>
'set y 'yc<br>
latval = subwrd(result,4)<br>
say ' (lonval,latval) = ('lonval','latval')'<br>
'q w2xy 'lonval' 'latval<br>
xpos = subwrd(result,3)<br>
ypos = subwrd(result,6)<br>
say ' (xpos,ypos) = ('xpos','ypos')'<br>
*<br>
* write (xpos,ypos) to output file...<br>
*<br></div></div>
res = write (''<a href="http://hurricane.nf" target="_blank">hurricane.nf</a>'','t = 'tt' xpos = 'xpos' ypos = 'ypos'')<div class="im"><br>
<br>
pause()<br>
*<br>
* erase box showing area for defining hurricane center...<br>
*<br>
* 'set line 0 1'<br>
* 'draw rec 'xb1o' 'yb1o' 'xb2o' 'yb2o<br>
tt=tt+1<br>
endwhile<br>
<br></div>
res = close (''<a href="http://hurricane.nf" target="_blank">hurricane.nf</a>'')<div class="im"><br>
<br>
nf=nf+1<br>
endwhile<br>
<br>
* Clear frame and re-plot base map...<br>
<br>
plot_base_map(plot_vpage,plot_area,lon1,lon2,lat1,lat2)<br>
<br>
* Plot track(s)...<br>
<br>
nf=1<br>
while ( nf <= nmodels )<br>
<br>
* Specify drawing primitives<br>
<br>
marktype = <a href="http://marks.nf" target="_blank">marks.nf</a><br></div>
lcolor = <a href="http://colors.nf" target="_blank">colors.nf</a><div class="im"><br>
<br>
say ' marktype, marksize, lcolor, lstyle and lthick:'<br>
say ' 'marktype ' ' marksize ' ' lcolor ' ' lstyle ' ' lthick<br>
<br>
tt=ts<br>
*<br>
* read (xpos,ypos) from ASCII text input file...<br>
*<br></div>
res = read(''<a href="http://hurricane.nf" target="_blank">hurricane.nf</a>'') ; rc = sublin(res,1)<div class="im"><br>
if( rc != 0 )<br>
say '"read" status = 'rc<br>
exit<br>
endif<br>
rec = sublin(res,2) ; say '"hurricane.'nf'" record: 'rec<br>
xc = subwrd(rec,6) ; yc = subwrd(rec,9)<br>
<br>
'set line 'lcolor' 'lstyle' 'lthick<br>
'draw mark 'marktype' 'xc' 'yc' 'marksize<br>
tt = tt+1<br>
while ( tt <= nt )<br></div>
res = read(''<a href="http://hurricane.nf" target="_blank">hurricane.nf</a>'') ; rc = sublin(res,1)<div><div></div><div class="h5"><br>
if( rc != 0 )<br>
say '"read" status = 'rc<br>
exit<br>
endif<br>
rec = sublin(res,2) ; say '"hurricane.'nf'" record: 'rec<br>
xn = subwrd(rec,6) ; yn = subwrd(rec,9)<br>
say ' 'xc' 'yc' 'xn' 'yn<br>
'draw line 'xc' 'yc' 'xn' 'yn<br>
'draw mark 'marktype' 'xn' 'yn' 'marksize<br>
tt=tt+1<br>
xc = xn<br>
yc = yn<br>
endwhile<br>
<br>
nf=nf+1<br>
endwhile<br>
<br>
say<br>
say '...plot legend'<br>
''legend_plotter' -x 3 -y 6 'legend_info<br>
say<br>
<br>
'print'<br>
'disable print'<br>
<br>
say<br>
say '************************************************'<br>
say ' Finished with this script.'<br>
say '************************************************'<br>
say<br>
<br>
return<br>
<br>
function plot_base_map(plot_vpage,plot_area,lon1,lon2,lat1,lat2)<br>
<br>
* Plots base map...<br>
<br>
'c'<br>
'set vpage 'plot_vpage<br>
'set parea 'plot_area<br>
'set lon 'lon1' 'lon2<br>
'set lat 'lat1' 'lat2<br>
'set mproj latlon'<br>
'set mpt * 1 1 5'<br>
'set mpdset hires'<br>
'set poli on'<br>
* 'set grid on'<br>
'set grads off'<br>
'set gxout contour'<br>
'draw map'<br>
<br>
* add lon/lat lines...<br>
<br>
'set ccolor 1'<br>
'set cstyle 5' ; 'd lon'<br>
'set cstyle 5' ; 'd lat'<br>
<br>
* add title, x- and y-labels...<br>
<br>
'draw title \\ Track of Minimum Surface Pressure'<br>
'draw xlab Longitude \\'<br>
'draw ylab \\ Latitude'<br>
<br>
return<br>
<br></div></div>function maxwind(args)<br>
'reinit'<br>
nc_file = subwrd(args,1)<br>
wind = subwrd(args,2)<br>
say nc_file<br>
say wind<br>
'sdfopen 'nc_file<br>
'q file'<br>
say result<div class="im"><br>
rec5 = sublin(result,5)<br>
nx = subwrd(rec5,3) ; ny = subwrd(rec5,6) ; nz = subwrd(rec5,9) ; nt = subwrd(rec5,12)<br></div>
say '(nx,ny,nz,nt) = ('nx','ny','nz','nt')'<br>
nxp = nx<br>
nyp = ny<br>
nzp = nz<br>
'set gxout print'<div class="im"><br>
tt = 1<br>
while ( tt <= nt )<br></div>
'set x 1'<br>
'set y 1'<div class="im"><br>
'set z 1'<br>
'set t 'tt<br></div>
say 'tt = 'tt<br>
'd max(max(max('wind',x=1,x='nx'),y=1,y='ny'),z=1,z='nz')'<br>
rec = sublin(result,(nyp+1)*nzp+3)<br>
max_wind = subwrd(rec,1)<br>
say ' max 'wind' = 'max_wind<br>
'd maxloc(max(max('wind',y=1,y='ny'),z=1,z='nz'),x=1,x='nx')'<br>
rec = sublin(result,(nzp+1)*nxp+3)<br>
xc = subwrd(rec,1)<br>
say ' xc = 'xc<br>
'd maxloc(max(max('wind',x=1,x='nx'),z=1,z='nz'),y=1,y='ny')'<br>
rec = sublin(result,(nzp+1)*nyp+3)<br>
yc = subwrd(rec,1)<br>
say ' yc = 'yc<br>
'd maxloc(max(max('wind',x=1,x='nx'),y=1,y='ny'),z=1,z='nz')'<br>
rec = sublin(result,(nyp+1)*nzp+3)<br>
zc = subwrd(rec,1)<br>
say ' zc = 'zc<br>
'd 'wind'(x='xc',y='yc',z='zc')'<br>
rec = sublin(result,2)<br>
wind_xc_yc_zc = subwrd(rec,1)<br>
say ' 'wind'(x='xc',y='yc',z='zc') = 'wind_xc_yc_zc<br>
tt=tt+1<br>
endwhile<br>
say '...finished.'<div class="im"><br>
return<br>
<br>_______________________________________________<br>
gradsusr mailing list<br>
<a href="mailto:gradsusr@gradsusr.org" target="_blank">gradsusr@gradsusr.org</a><br>
</div><div class="im"><a href="http://gradsusr.org/mailman/listinfo/gradsusr" target="_blank">http://gradsusr.org/mailman/listinfo/gradsusr</a><br>
<br></div></blockquote></div><br><br clear="all"><br>-- <br><br><br><br>
</div>
</blockquote></div><br><br clear="all"><br>-- <br>सुशांत पुराणिक <br>वातावरणशास्त्र व अवकाशशास्त्र विभाग,<br>पुणे विद्यापीठ <br>पुणे-०७, महाराष्ट्र <br>भारत. <br><br><br><br>
</div>