;_______________________________________________________________________________________ ;To run the script type: ; ncl CreateTracks.ncl {options} ; ; e.g. ; ncl CreateTracks.ncl 'inputFILE="track_2005082800_12k"' 'stormNAME="Katrina"' ; ; ncl CreateTracks.ncl 'inputFILE="files"' 'stormNAME="Katrina"' 'wksTYPE="x11"' \ ; plotTYPE=3 bestTRACK=True 'btFILE="../BestTrack.Katrina"' \ ; 'btCOLOR="black"' lonMIN=-92.0 ; ;_______________________________________________________________________________________ ;The following options are available ; ; inputFILE ; The input file to be used ; For summary plots, the file containing a list of track data files, and ; the dates associated with the data ; No filename convensions required ; MUST be supplied ; ; stormNAME ; Name of the storm ; Optional argument ; The word "UnKnown" will be use for filenames without storm names ; ; dataSOURCE ; RIP4 (tracks generated by RIP4 code) ; ARW (tracks generated by WRF model) ; Default is "ARW" ; ; startTIME ; For RIP4 data, a start time needs to be supplied, as this is not available ; in the output data ; Required for RIP4 ; ; wksTYPE ; X11 / ps / pdf/ ncgm ; Optional argument. Default is pdf ; ; dataINT ; Real value ; Used only for RIP4 data ; Default is data input interval of 3 hours, set this if input interval is different ; ; plotTYPE ; Default (0) is single run plots ; Summary plots: ; 1 - similar to daily plots ; 2 - similar to 1, but without info on trackline. Tracks are numbered ; 3 - colors are use for hurricane CAT. No hurricane symbols ; 4 - different colors for each forecast. No hurricane symbols ; CONVENSION ; For RIP4 data, file names and start dates for each file needs to be supplied ; in the following manner ; filename : YYYYMMDDHH ; ; outputFILE ; Name of output file to be used ; Optional argument. Default is "hur_track" ; ; rotPNG ; Rotate the PNG file created from pdf file ; Default is =90 ; ; bestTRACK ; Should a Best Track be plotted. ; Default is False ; ; btFILE ; Best Track input file name ; Default is "BestTrack" ; ; btCOLOR ; Best Track line color ; Default is Red ; ; latMIN / latMAX / lonMIN / lonMAX ; The min window size to use ; Optional argument ; Defaults are: 25N / 40N / 90W / 75W ; If overwritten by user, the plot may be outside the plot. ; ; disEDGE ; Distance in degrees between track and edge of plot ; Optional argument ; Default is 2 degrees ; ; textANGLE ; To rotate angle of text on track line ; Default is no rotation ; Set to 45.0 for SW-NE orientation ; ; StartLineColor ; Used for summary plots 4 ; Changes the start color for the track lines. ; Default is Blue (5) - see src/TrackColors.ncl for a list of colors ; or to add colors ; ;_______________________________________________________________________________________ load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" ;_______________________________________________________________________________________ ; White and Black are background colors ; LightBlue and SandyBrown are used for the map background ; Azure4 (darkgrey) is used for the map borders and frame ; Blue, DarkGreen, Yellow, Red and Purple3 are used for the hurricane symbols (5,6,7,8,9) ; Add colors below "Purple3" if more are needed undef("TrackColors") function TrackColors() begin track_colors = (/"White","Black", \ "Azure4", \ "LightBlue", \ "SandyBrown", \ "Blue", \ "DarkGreen", \ "Yellow", \ "Red", \ "Purple3", \ "CadetBlue4", \ "Olivedrab1", \ "Green", \ "Orchid1", \ "Cyan", \ "CornFlowerBlue"/) return(track_colors) end ;_______________________________________________________________________________________ ; Make sure the final window is bigger than the data spread undef("MapWindow") procedure MapWindow(lat,lon,wnd,tkres) begin do i=0,dimsizes(wnd)-1 if ( wnd(i) .eq. -1.0 ) then lat(i) = 30. lon(i) = -80. end if end do min_lat = min(lat) max_lat = max(lat) min_lon = min(lon) max_lon = max(lon) if ( min_lat .lt. tkres@latMIN .and. .not. tkres@latMINset ) then min_lat = min_lat - tkres@disEDGE else min_lat = tkres@latMIN end if tkres@min_lat = min_lat if ( max_lat .gt. tkres@latMAX .and. .not. tkres@latMAXset ) then max_lat = max_lat + tkres@disEDGE else max_lat = tkres@latMAX end if tkres@max_lat = max_lat if ( min_lon .lt. tkres@lonMIN .and. .not. tkres@lonMINset ) then min_lon = min_lon - tkres@disEDGE else min_lon = tkres@lonMIN end if tkres@min_lon = min_lon if ( max_lon .gt. tkres@lonMAX .and. .not. tkres@lonMAXset ) then max_lon = max_lon + tkres@disEDGE else max_lon = tkres@lonMAX end if tkres@max_lon = max_lon end ;_______________________________________________________________________________________ ; Set up the basic map resources undef("MapResources") procedure MapResources(res,tkres) begin res = True res@gsnMaximize = True ; Maximize plot size in the frame. res@gsnDraw = False ; Don't draw the plot just yet. res@gsnFrame = False ; Don't advance the frame just yet. res@pmTickMarkDisplayMode = "Always" res@mpMinLatF = tkres@min_lat res@mpMaxLatF = tkres@max_lat res@mpMinLonF = tkres@min_lon res@mpMaxLonF = tkres@max_lon res@mpOutlineBoundarySets = "GeophysicalAndUSStates" res@mpDataBaseVersion = "MediumRes" res@mpFillColors = (/"background","LightBlue","SandyBrown","LightBlue", \ "transparent"/) res@mpGeophysicalLineColor = "Azure4" res@mpGeophysicalLineThicknessF = 0.5 res@mpGridLineColor = "Azure4" res@mpGridLineThicknessF = 0.5 res@mpGridMaskMode = 3 res@mpGridSpacingF = 5 res@mpLimbLineColor = "Azure4" res@mpLimbLineThicknessF = 0.5 res@mpNationalLineColor = "Azure4" res@mpNationalLineThicknessF = 0.5 res@mpUSStateLineColor = "Azure4" res@mpUSStateLineThicknessF = 0.5 res@mpGridAndLimbOn = True res@mpOutlineOn = True ; ; Tickmark stuff. ; res@tmXBMajorOutwardLengthF = 0.0 res@tmXBMajorLengthF = 0.0 res@tmXBLabelFontHeightF = 0.007 end ;_______________________________________________________________________________________ ; Add a title header topleft undef("PlotHeaderL") procedure PlotHeaderL(wks,map,tkres) begin ;Add header at top txt0 = create "MainPlotTille" textItemClass wks "txString" : tkres@header "txFontHeightF" : 0.018 "txBackgroundFillColor" : "White" "txPerimOn" : True "txPerimColor" : "Black" "txFont" : "helvetica-bold" end create anno = NhlAddAnnotation(map,txt0) setvalues anno "amZone" : 0 "amSide" : "Top" "amJust" : "TopLeft" "amParallelPosF" : -0.5 "amOrthogonalPosF" : 0.5 "amResizeNotify" : False end setvalues end ;_______________________________________________________________________________________ ; Add time information topright undef("PlotHeaderR") procedure PlotHeaderR(wks,map,header,tkres) begin ; Set return character cr = "~C~" ; Placement of the time info slon = tkres@max_lon - 0.2 slat = tkres@max_lat - 0.2 txres = True txres@txFont = "helvetica-bold" txres@txFontColor = "Black" txres@txFontHeightF = 0.01 txres@txPerimOn = True txres@txPerimColor = "Black" txres@txJust = "TopRight" txres@txBackgroundFillColor = "White" if (tkres@plotTYPE .eq. 0) then timestring = "TIME:" do i = 0, dimsizes(header)-1 timestring = timestring + cr + header(i) end do else timestring = "Forecasts:" do i = 0, dimsizes(header)-1 ii = i+1 timestring = timestring + cr + ii + ". " + header(i) end do end if gsn_text(wks,map,timestring,slon,slat,txres) end ;_______________________________________________________________________________________ ; Functions to plot different symbols undef("create_hurricane_symbol") function create_hurricane_symbol(wks) begin mstring = "p" fontnum = 37 xoffset = 0.0 yoffset = 0.0 ratio = 1.0 size = 1.0 angle = 0.0 new_index = NhlNewMarker(wks, mstring, fontnum, xoffset, yoffset, \ ratio, size, angle) return(new_index) end undef("create_diamond_symbol") function create_diamond_symbol(wks) begin mstring = "q" fontnum = 35 xoffset = 0.0 yoffset = 0.0 ratio = 1.0 size = 1.0 angle = 0.0 new_index = NhlNewMarker(wks, mstring, fontnum, xoffset, yoffset, \ ratio, size, angle) return(new_index) end ;_______________________________________________________________________________________ ; get lat / lon / lev / wnd from input data ; input filenames are hardcorded - may need to change this later ; return an array (4,*), where the second dim is the number of ; input data to be plotted. ; The first dim contain the fields lat-lon-lev-wnd undef("GetTrackData") function GetTrackData(tkres) begin if ( tkres@dataSOURCE .eq. "RIP4" ) then filename = tkres@filename ; Read lat/lon speed and pressure tmp = systemfunc("cut -c4-10 " + filename) lat = stringtofloat(tmp(tkres@StartGoodData:tkres@EndGoodData)) tmp = systemfunc("cut -c13-20 " + filename) lon = stringtofloat(tmp(tkres@StartGoodData:tkres@EndGoodData)) tmp = systemfunc("cut -c53-60 " + filename) lev = stringtofloat(tmp(tkres@StartGoodData:tkres@EndGoodData)) tmp = systemfunc("cut -c65-71 " + filename) wnd = stringtofloat(tmp(tkres@StartGoodData:tkres@EndGoodData)) ; need to know how many points we have fno = lat tkres@numGridPoints = dimsizes(lat) do it = 0,tkres@numGridPoints-1 fno(it) = it end do end if if ( tkres@dataSOURCE .eq. "ARW" ) then filename = tkres@filename ; Read lat/lon speed and pressure lat = stringtofloat(systemfunc("cut -c29-34 " + filename)) lon = stringtofloat(systemfunc("cut -c35-43 " + filename)) lev = stringtofloat(systemfunc("cut -c44-50 " + filename)) wnd = stringtofloat(systemfunc("cut -c51-57 " + filename)) dummy = dimsizes(lat) fno = new( dummy , float ) fno = 0.0 end if npoints = dimsizes(lat) ll_lev_wnd = new( (/5,npoints/) , float ) ll_lev_wnd(0,:) = lat ll_lev_wnd(1,:) = lon ll_lev_wnd(2,:) = lev ll_lev_wnd(3,:) = wnd ll_lev_wnd(4,:) = fno return (ll_lev_wnd) end ;_______________________________________________________________________________________ ; get date and time information for lables on plot undef("Get_dd_hh") function Get_dd_hh(fno,time_start_end,tkres) begin if ( tkres@dataSOURCE .eq. "RIP4" ) then tm = new (dimsizes(fno), string) MonthEnd = new ( 12 , integer ) MonthEnd = (/ 32,29,32,31,32,31,32,32,31,32,31,32 /) yy = (systemfunc("echo " + time_start_end(0) + " | cut -c1-4" ) ) mm = stringtointeger(systemfunc("echo " + time_start_end(0) + " | cut -c5-6" ) ) dd = stringtointeger(systemfunc("echo " + time_start_end(0) + " | cut -c7-8" ) ) hh = stringtointeger(systemfunc("echo " + time_start_end(0) + " | cut -c9-10" ) ) fno = fno*tkres@interval fnoi = floattoint(fno) do i = 0,dimsizes(fnoi)-1 if ( i .eq. 0 ) then hh = hh + fnoi(0) else hh = hh + fnoi(i) - fnoi(i-1) end if if ( hh .eq. 24 ) then dd = dd + 1 hh = 0 if ( dd .eq. MonthEnd(mm) ) then mm = mm + 1 dd = 1 end if end if if ( mm .lt. 10 ) then time_start_end(1) = yy + "0" + mm else time_start_end(1) = yy + mm end if if ( dd .lt. 10 ) then tm(i) = "0" + dd time_start_end(1) = time_start_end(1) + "0" + dd else tm(i) = dd time_start_end(1) = time_start_end(1) + dd end if if ( hh .lt. 10 ) then tm(i) = tm(i) + "/0" + hh time_start_end(1) = time_start_end(1) + "0" + hh else tm(i) = tm(i) + "/" + hh time_start_end(1) = time_start_end(1) + hh end if end do end if if ( tkres@dataSOURCE .eq. "ARW" ) then ; Read Time stamps dd = (systemfunc("cut -c15-16 " + tkres@filename)) hh = (systemfunc("cut -c18-19 " + tkres@filename)) tm = dd + "/" + hh end if return (tm) end ;_______________________________________________________________________________________ ; get interval of the input data undef("GetInterval") function GetInterval(tkres) begin if ( tkres@dataSOURCE .eq. "RIP4" ) then interval = tkres@dataINT end if if ( tkres@dataSOURCE .eq. "ARW" ) then ; This dataset has a full date, so we can calculate the interval of ; incoming data filename = tkres@filename mins = stringtofloat(systemfunc("cut -c21-22 " + filename)) ; As the first time may be wrong, we use the 2nd and 3rd times for calculations interval = mins(3) - mins(2) ; We need to know the answer as a fraction of an hour interval = interval / 60. end if print("Data interval is " + interval + " hour(s)" ) return(interval) end ;_______________________________________________________________________________________ ; get start and end times to place on the map as information ; return an array (2), where the the first string is the start ; time and the second one the end time undef("GetTimeInfo") function GetTimeInfo(tkres) begin if ( tkres@dataSOURCE .eq. "RIP4" ) then tstart = tkres@startTIME tend = tkres@startTIME end if if ( tkres@dataSOURCE .eq. "ARW" ) then tstart= systemfunc("head -n 1 " + tkres@filename + " | cut -c7-19 ") tend = systemfunc("tail -n 1 " + tkres@filename + " | cut -c7-19 ") end if time_start_end = new(2,string) time_start_end(0) = tstart time_start_end(1) = tend return(time_start_end) end ;_______________________________________________________________________________________ ; get DataType and Hurricane Categories ; return an array (6), where the first dimension is the DataType ; and the other 5 is cat1 to cat5 undef("GetCatInfo") function GetCatInfo(tkres) begin if ( tkres@dataSOURCE .eq. "RIP4" ) then ; we know this data is in m/s DataType = "m/s" cat1 = 33.0 cat2 = 42.0 cat3 = 51.0 cat4 = 59.0 cat5 = 69.0 end if if ( tkres@dataSOURCE .eq. "ARW" ) then ; we know this data is in kts DataType = "kts" cat1 = 64.0 cat2 = 84.0 cat3 = 97.0 cat4 = 114.0 cat5 = 135.0 end if cat = new(6,string) cat(0) = DataType cat(1) = cat1 cat(2) = cat2 cat(3) = cat3 cat(4) = cat4 cat(5) = cat5 return(cat) end ;_______________________________________________________________________________________ ; Plot best track if asked for ; currently only a standard plot, but this can be changed undef("PlotBestTrack") procedure PlotBestTrack(wks,map,tkres) begin print("Now attempting to plot bestTRACK from data file " + tkres@btFILE) lat = stringtofloat(systemfunc("cut -c18-21 " + tkres@btFILE)) latsign = systemfunc("cut -c22-22 " + tkres@btFILE) lon = stringtofloat(systemfunc("cut -c24-28 " + tkres@btFILE)) lonsign = systemfunc("cut -c29-29 " + tkres@btFILE) bthh = systemfunc("cut -c1-2 " + tkres@btFILE) btdd = systemfunc("cut -c11-12 " + tkres@btFILE) btdate = btdd + "/" + bthh btcat = stringtoint(systemfunc("cut -c33-35 " + tkres@btFILE)) ; Ckect to see if we have S/N ; W/E data checkDATA = 0 do isign=0,dimsizes(lat)-1 if (latsign(isign) .eq. "S") then lat(isign) = -1.0 * lat(isign) else if (latsign(isign) .ne. "N") then checkDATA = checkDATA + 1 end if end if if (lonsign(isign) .eq. "W") then lon(isign) = -1.0 * lon(isign) else if (lonsign(isign) .ne. "E") then checkDATA = checkDATA + 1 end if end if end do if ( checkDATA .gt. 0 ) then print("###########################################################") print("Looks like some of your best TRACK data is not in columns.") print("Please ensure all data is in columns and re-run the script.") print("###########################################################") end if lnres = True lnres@gsLineThicknessF = 2.6 do i=0,dimsizes(lat)-2 lnres@gsLineColor = tkres@btCOLOR ;if ( lon(i) .lt. -75 ) then gsn_polyline(wks,map,(/lon(i),lon(i+1)/),(/lat(i),lat(i+1)/),lnres) ;end if end do print(" We have " + dimsizes(lat) + " data points") print(" Track start point (lon/lat) = "+ lon(0) +" "+ lat(0)) print(" Track end point (lon/lat) = "+ lon(dimsizes(lat)-1) +" "+ lat(dimsizes(lat)-1)) txres2 = True txres2@txFont = "helvetica-bold" txres2@txFontHeightF = 0.008 txres2@txJust = "CenterLeft" txres2@txAngleF = 45.00 lnres2 = True lnres2@gsLineColor = "Black" lnres2@gsLineThicknessF = 0.7 pmres = True pmres@gsMarkerIndex = create_hurricane_symbol(wks) do i=0,dimsizes(lat)-1 pmres@gsMarkerColor = "Black" pmres@gsMarkerIndex = 16 pmres@gsMarkerSizeF = 0.005 if (btcat(i) .ge. 74) then pmres@gsMarkerColor = 5 end if if (btcat(i) .ge. 96) then pmres@gsMarkerColor = 6 end if if (btcat(i) .ge. 111) then pmres@gsMarkerColor = 7 end if if (btcat(i) .ge. 131) then pmres@gsMarkerColor = 8 end if if (btcat(i) .ge. 155) then pmres@gsMarkerColor = 9 end if if (btcat(i) .ge. 74) then pmres@gsMarkerIndex = create_hurricane_symbol(wks) pmres@gsMarkerSizeF = 0.01 end if gsn_polyline(wks,map,(/lon(i),lon(i)+1.0/),(/lat(i),lat(i)+1.0/),lnres2) gsn_polymarker(wks,map,lon(i)+1.0,lat(i)+1.0,pmres) if ( bthh(i) .eq. "03" )then gsn_text(wks,map,btdate(i),lon(i)+1.2,lat(i)+1.2,txres2) end if end do end ;_______________________________________________________________________________________ ; Get some line number if we are dealing with RIP4 data ; This is to strip out all junk from the input data ; ; line_breaks will hold 3 pieces of information ; start of good data ; end of good data ; number of lines in file ; undef("GetLineNums") function GetLineNums(tkres) begin if ( tkres@dataSOURCE .eq. "RIP4" ) then filename = tkres@filename s_tester = (systemfunc("cut -c4-10 " + filename)) tester = new(dimsizes(s_tester),float) line_breaks = new(3,integer) line_breaks = 0 line_breaks(1) = dimsizes(s_tester) - 1 line_breaks(2) = dimsizes(s_tester) tkres@bad_lines = False found_start = False do i =0,dimsizes(s_tester)-1 tester(i) = stringtofloat(s_tester(i)) if ( .not. tkres@bad_lines .and. ismissing(tester(i)) ) then tkres@bad_lines = True end if if ( .not. found_start .and. .not. ismissing(tester(i)) ) then line_breaks(0) = i found_start = True end if if ( found_start .and. ismissing(tester(i)) ) then line_breaks(1) = i - 1 return (line_breaks) end if end do end if return (line_breaks) end ;_______________________________________________________________________________________ ; Draw the track line on the map that we have created previously undef("PlotTrackLine") procedure PlotTrackLine(wks,map,lat,lon,wnd,cat,tkres) begin ; We are going to draw the line every 3 hour, or every timestep, ; Depending on frequency of input data line_int = floattoint(3/tkres@interval) if ( line_int .le. 1 ) line_int = 1 end if ; Line resources lnres = True lnres@gsLineColor = "Black" lnres@gsLineThicknessF = 2.7 if (tkres@plotTYPE .eq. 4) lnres@gsLineColor = tkres@StartLineColor + tkres@iNumFil lnres@gsLineThicknessF = 3.0 end if ; Draw the line do i=0,dimsizes(wnd)-line_int-1,line_int if ( wnd(i) .ne. -1.0 .and. wnd(i+line_int) .ne. -1.0 ) then if (tkres@plotTYPE .eq. 3) lnres@gsLineColor = "Azure4" if (wnd(i) .gt. stringtofloat(cat(1))) lnres@gsLineColor = 5 end if if (wnd(i) .gt. stringtofloat(cat(2))) lnres@gsLineColor = 6 end if if (wnd(i) .gt. stringtofloat(cat(3))) lnres@gsLineColor = 7 end if if (wnd(i) .gt. stringtofloat(cat(4))) lnres@gsLineColor = 8 end if if (wnd(i) .gt. stringtofloat(cat(5))) lnres@gsLineColor =9 end if end if gsn_polyline(wks,map,(/lon(i),lon(i+line_int)/),(/lat(i),lat(i+line_int)/),lnres) end if end do end ;_______________________________________________________________________________________ ; Draw pressure and wind information to left of track line ; Plot date information to the right of the track line undef("PlotTrackInfo") procedure PlotTrackInfo(wks,map,lat,lon,wnd,wndlev,dd_hh,tkres) begin ; We are going to plot the information every 12 hours, or every timestep, ; Depending on frequency of input data label_int= floattoint(12/tkres@interval) if ( label_int .le. 1 ) label_int = 1 end if if (tkres@plotTYPE .le. 1) ; Plot the Pressure/Wind Info to left of track line txres = True txres@txFont = "helvetica-bold" txres@txFontHeightF = 0.009 txres@txJust = "CenterRight" txres@txAngleF = tkres@textANGLE do i=0,dimsizes(wnd)-1,label_int if ( wnd(i) .ne. -1.0) then gsn_text(wks,map,wndlev(i),lon(i),lat(i),txres) end if end do delete(txres) end if ; Plot date information to right of track line txres = True txres@txFont = "helvetica-bold" txres@txJust = "CenterLeft" txres@txAngleF = tkres@textANGLE dis_marker = 0.0 if (tkres@plotTYPE .le. 1) txres@txBackgroundFillColor = "White" txres@txPerimOn = True txres@txPerimColor = "Black" txres@txFontHeightF = 0.008 dis_marker = 0.3 ; how far must time stamp be drawn from marker if (tkres@textANGLE .gt. 25.) dis_marker = 0.1 ; how far must time stamp be drawn from marker end if end if if (tkres@plotTYPE .ge. 2) txres@txFontHeightF = 0.007 dis_marker = 0.2 ; how far must time stamp be drawn from marker label_int= floattoint(24/tkres@interval) ; we will only plot these every 24 hours end if lon_dd_hh = lon + dis_marker lat_dd_hh = lat + dis_marker do i=0,dimsizes(wnd)-1,label_int if ( wnd(i) .ne. -1.0) then gsn_text(wks,map,dd_hh(i),lon_dd_hh(i),lat_dd_hh(i),txres) end if end do delete(txres) ; Add forecast numbers if (tkres@plotTYPE .ge. 2) txres = True txres@txFont = "helvetica-bold" txres@txJust = "BottomLeft" txres@txBackgroundFillColor = "White" txres@txPerimOn = True txres@txPerimColor = "Black" txres@txFontHeightF = 0.008 forecastnumber = 0 do i=0,dimsizes(wnd)-1,floattoint(3/tkres@interval) ; this is to find end of track line itest = dimsizes(wnd)-1-i if ( forecastnumber .eq. 0 .and. wnd(itest) .ne. -1.0 ) then gsn_text(wks,map,tkres@iNumFil+1,lon(itest),lat(itest),txres) forecastnumber = 1 end if end do end if end ;_______________________________________________________________________________________ ; Draw Sysmbols over the track line undef("PlotTrackSymbols") procedure PlotTrackSymbols(wks,map,lat,lon,wnd,cat,tkres) begin ; We are going to draw the symbols every 6 hours, or every timestep, ; Depending on frequency of input data symbol_int= floattoint(6/tkres@interval) if ( symbol_int .le. 1 ) symbol_int = 1 end if pmres = True if (tkres@plotTYPE .ge. 3 ) then do i=0,dimsizes(wnd)-1,symbol_int if ( wnd(i) .ne. -1.0) then pmres@gsMarkerColor = "White" pmres@gsMarkerIndex = 16 pmres@gsMarkerSizeF = 0.003 gsn_polymarker(wks,map,lon(i),lat(i),pmres) end if end do return end if do i=0,dimsizes(wnd)-1,symbol_int if ( wnd(i) .ne. -1.0) then pmres@gsMarkerColor = "Black" if ( wnd(i) .gt. stringtofloat(cat(1))) then pmres@gsMarkerIndex = create_hurricane_symbol(wks) pmres@gsMarkerSizeF = 0.01 else ;pmres@gsMarkerIndex = 4 pmres@gsMarkerIndex = 16 pmres@gsMarkerSizeF = 0.006 end if if (wnd(i) .gt. stringtofloat(cat(1))) then pmres@gsMarkerColor = "Blue" end if if (wnd(i) .gt. stringtofloat(cat(2))) then pmres@gsMarkerColor = "DarkGreen" end if if (wnd(i) .gt. stringtofloat(cat(3))) then pmres@gsMarkerColor = "Yellow" end if if (wnd(i) .gt. stringtofloat(cat(4))) then pmres@gsMarkerColor = "Red" end if if (wnd(i) .gt. stringtofloat(cat(5))) then pmres@gsMarkerColor = "Purple3" end if gsn_polymarker(wks,map,lon(i),lat(i),pmres) end if end do end ;_______________________________________________________________________________________ ; Add Hurricane sysmbols (color) and cat's bottomleft undef("HurricaneLegend") procedure HurricaneLegend(wks,map,tkres) begin if (tkres@plotTYPE .eq. 4) return end if ; Placement of the legdends slon = tkres@min_lon + 0.5 slat = tkres@min_lat + 0.5 pmres = True pmres@gsMarkerSizeF = 0.013 pmres@gsMarkerIndex = create_hurricane_symbol(wks) txres = True txres@txFont = "helvetica-bold" txres@txFontColor = "White" txres@txFontHeightF = 0.008 start_color = 5 ; our color table has Blue as color number 5 do j=0,4 ; plot 5 color symbols pmres@gsMarkerColor = start_color + j gsn_polymarker(wks,map,slon,slat,pmres) gsn_text(wks,map,1+j,slon,slat,txres) slon = slon + 0.5 end do ; Draw a box around the legend lnres = True lnres@gsLineColor = "Black" lnres@gsLineThicknessF = 1.7 slat = slat + 0.4 gsn_polyline(wks,map,(/tkres@min_lon,slon/),(/slat+0.4,slat+0.4/),lnres) gsn_polyline(wks,map,(/slon,slon/),(/slat+0.4,tkres@min_lat/),lnres) slon = slon - 0.5 - 0.5 - 0.5 ; Add some text txres@txFontColor = "Black" txres@txFontHeightF = 0.009 txres@txJust = "BottomCenter" gsn_text(wks,map,"hur cat:",slon,slat,txres) end ;_______________________________________________________________________________________ ; Add pressure/wind and dd/hh legend to bottomright undef("SymbolLegend") procedure SymbolLegend(wks,map,tkres) begin if ( tkres@plotTYPE .ge. 2) return end if ; Set return character cr = "~C~" ; Placement of the ledgends slon = tkres@max_lon - 2.0 slat = tkres@min_lat + 0.5 pmres = True pmres@gsMarkerSizeF = 0.008 pmres@gsMarkerIndex = create_hurricane_symbol(wks) pmres@gsMarkerColor = "Black" gsn_polymarker(wks,map,slon,slat,pmres) txres = True txres@txFont = "helvetica-bold" txres@txFontHeightF = 0.008 txres@txJust = "CenterRight" sslon = slon - 0.5 left = "pressure" + cr + "max wnd ("+tkres@DataType+")" gsn_text(wks,map,left,sslon,slat,txres) txres@txBackgroundFillColor = "White" txres@txPerimOn = True txres@txJust = "CenterLeft" right = "dd/hh" slon = slon + 0.5 gsn_text(wks,map,right,slon,slat,txres) end ;_______________________________________________________________________________________ ;_______________________________________________________________________________________ ;_______________________________________________________________________________________ ; Main script begin ;---------------------------------------------------------------------- ; Get dataSOURCE and file info from command line ;---------------------------------------------------------------------- tkres = True ; Make sure we have a datafile to work with if (.not. isvar("inputFILE") ) then print(" ") print(" ### MUST SUPPLY a inputFILE ### ") print(" ") print(" Something like: ") print(" ncl CreateTracks.ncl inputFILE=track_2005082800_12k stormNAME=Katrina dataSOURCE=RIP4 wksTYPE=X11 ") print(" REMEMBER TO ADD QUOTES" ) print(" Refer to the information at the top od this file for more info and syntax" ) exit end if tkres@inputFILE = inputFILE ; dataSOURCE = ARW (data from WRF-ARW model) ; dataSOURCE = RIP4 (data from RIP) if (.not. isvar("dataSOURCE") ) then dataSOURCE = "ARW" end if if (dataSOURCE .ne. "RIP4" .and. dataSOURCE .ne. "ARW") then print("'" + dataSOURCE + "'" + " is not a valid dataSOURCE") print(" Current options are :") print(" RIP4 for data generated by rip4") print(" ARW for data generated by the WRF-ARW model") exit end if tkres@dataSOURCE = dataSOURCE print("Input is " + tkres@dataSOURCE + " type data") ; Do we have any special data that comes in at a different interval from default ; As we only care about this for RIP4 data (ARW has all the info in the file) ; We are going to make the default 3h if (.not. isvar("dataINT") .or. dataINT .eq. 0.0) then dataINT = 3.0 if (dataSOURCE .eq. "RIP4" ) then print(" Using a default of " +dataINT+ " hours between data points" ) print(" If this is not correct, use dataINT on the command line, to override" ) end if end if tkres@dataINT = dataINT ;---------------------------------------------------------------------- ; Are we dealing with SUMMARIES ;---------------------------------------------------------------------- if (.not. isvar("plotTYPE") ) then plotTYPE = 0 end if if (plotTYPE .eq. 0) then filename = inputFILE tkres@plotTYPE = 0 print("Doing a single time plot") else filename = (systemfunc("cut -d':' -f1 " + inputFILE)) filedates = (systemfunc("cut -d':' -f2 " + inputFILE)) TRHeader = new(dimsizes(filename),string) tkres@plotTYPE = plotTYPE print("Doing summary plot type " + tkres@plotTYPE) if ( dimsizes(filename) .gt. 10 ) then print("########################################################") print("Number of input files for summary is: " + dimsizes(filename) ) print("If code is hanging or this looks wrong,") print("make sure you are not pointing to a track data file,") print("but to a file containing a list if track data files.") print("This will be true even for one track plot.") print("########################################################") end if end if ; Do we have any special data that does not contain date info ; As we only care about this for RIP4 data (ARW has all the info in the file) ; No default - we need this set if (dataSOURCE .eq. "RIP4" .and. plotTYPE .eq. 0 .and. .not. isvar("startTIME") ) then print(" ") print(" ### MUST SUPPLY a startTIME for RIP4 data ### ") print(" ") exit end if if (isvar("startTIME") ) then tkres@startTIME = startTIME end if ;---------------------------------------------------------------------- ; Basic Setup ;---------------------------------------------------------------------- ; open the workstation if (.not. isvar("wksTYPE") ) then wksTYPE = "pdf" end if tkres@wksTYPE = wksTYPE if (.not. isvar("outputFILE") ) then outputFILE = "hur_track" end if tkres@outputFILE = outputFILE wks = gsn_open_wks(wksTYPE,outputFILE) ; set up header if (.not. isvar("stormNAME") ) then header = "ARW Forecast" stormNAME = "UnKnown" else header = "ARW Forecast: " + stormNAME end if tkres@header = header tkres@stormNAME = stormNAME ; are we doing the best track if (.not. isvar("bestTRACK") ) then bestTRACK = False end if tkres@bestTRACK = bestTRACK if (.not. isvar("btFILE") ) then btFILE = "BestTrack" end if tkres@btFILE = btFILE if (.not. isvar("btCOLOR") ) then btCOLOR = "Red" end if tkres@btCOLOR = btCOLOR ; less common settings ; rotate the text angle on the track line if (.not. isvar("textANGLE") ) then textANGLE=0.0 end if tkres@textANGLE = textANGLE ; rotate the output .pdf file if (.not. isvar("rot") ) then rot=-90 end if tkres@rotPDF = rot ; change the start number for summary 4 plots if (.not. isvar("StartLineColor") ) then StartLineColor=5 end if tkres@StartLineColor = StartLineColor ; Set return character cr = "~C~" ;---------------------------------------------------------------------- ; Read storm track and info ;---------------------------------------------------------------------- do iNumFil=0,dimsizes(filename)-1 tkres@iNumFil = iNumFil tkres@filename = filename(iNumFil) if ( plotTYPE .ne. 0 ) then tkres@startTIME = filedates(iNumFil) end if if ( dataSOURCE .eq. "RIP4" ) then ; need to clean up the datafile ; First figure out some line numbers line_breaks = GetLineNums(tkres) tkres@StartGoodData = line_breaks(0) tkres@EndGoodData = line_breaks(1) if ( tkres@bad_lines ) then print(" ") print("Warning messages above this point is fine, and can be ignored") print("=============================================================") print(" ") end if end if ll_lev_wnd = GetTrackData(tkres) lat = ll_lev_wnd(0,:) lon = ll_lev_wnd(1,:) lev = ll_lev_wnd(2,:) wnd = ll_lev_wnd(3,:) fno = ll_lev_wnd(4,:) interval = GetInterval(tkres) tkres@interval = interval time_start_end = GetTimeInfo(tkres) dd_hh = Get_dd_hh(fno,time_start_end,tkres) cat = GetCatInfo(tkres) tkres@DataType = cat(0) ; set up the data the way we want to plot it slev = sprintf(" %4.0f",lev) swnd = sprintf(" %5.1f",wnd) wndlev = slev + " " + cr + swnd + " " ; Create top right header if (tkres@plotTYPE .eq. 0) then TRHeader = time_start_end else TRHeader(iNumFil) = time_start_end(0) end if ;---------------------------------------------------------------------- ; Set up the window ;---------------------------------------------------------------------- ; Set user over ride of plot domain to true. If user does not use ; the override commands, this will be set to False again tkres@latMINset = True tkres@latMAXset = True tkres@lonMINset = True tkres@lonMAXset = True if (iNumFil .eq. 0) then ; defaults if not overwritten on the command line if (.not. isvar("latMIN") ) then latMIN = 25. tkres@latMINset = False end if tkres@latMIN = latMIN if (.not. isvar("latMAX") ) then latMAX = 40. tkres@latMAXset = False end if tkres@latMAX = latMAX if (.not. isvar("lonMIN") ) then lonMIN = -90. tkres@lonMINset = False end if tkres@lonMIN = lonMIN if (.not. isvar("lonMAX") ) then lonMAX = -75. tkres@lonMAXset = False end if tkres@lonMAX = lonMAX if (.not. isvar("disEDGE") ) then disEDGE = 2. end if tkres@disEDGE= disEDGE ; Make sure the final window is bigger than the data spread MapWindow(lat,lon,wnd,tkres) print("The Track Plot Window is set to") print(" Latitute from " + tkres@min_lat + " to " + tkres@max_lat) print(" Longitute from " + tkres@min_lon + " to " + tkres@max_lon) print("The buffer between the data and the plot edge is " + tkres@disEDGE + " degrees") ;---------------------------------------------------------------------- ; Set up resources and create map ;---------------------------------------------------------------------- ; Get colors to use track_colors = TrackColors() gsn_define_colormap(wks,track_colors) ;gsn_draw_colormap(wks) ; Set some map resources and zoom in res = False MapResources(res,tkres) ; Create the map. map = gsn_csm_map(wks,res) ; Add Header top left PlotHeaderL(wks,map,tkres) draw(map) end if ;---------------------------------------------------------------------- ; The plotting section. ;---------------------------------------------------------------------- ; Draw the track line PlotTrackLine(wks,map,lat,lon,wnd,cat,tkres) ; Add Symbols over line PlotTrackSymbols(wks,map,lat,lon,wnd,cat,tkres) ; Add pressure/wind information on track line ; Add date information on track line ; Only when needed - determine by plotTYPE PlotTrackInfo(wks,map,lat,lon,wnd,wndlev,dd_hh,tkres) delete(ll_lev_wnd) delete(interval) delete(fno) delete(cat) delete(lat) delete(lon) delete(lev) delete(wnd) delete(slev) delete(swnd) delete(wndlev) delete(dd_hh) end do ; Add time information top right PlotHeaderR(wks,map,TRHeader,tkres) ; Add hurricane legend in bottom left HurricaneLegend(wks,map,tkres) ; Add pressure/wind and dd/hh legend to bottom right SymbolLegend(wks,map,tkres) ; Finally - are we doing a best track if (tkres@bestTRACK) then PlotBestTrack(wks,map,tkres) end if ; We are done frame(wks) delete(wks) end