Hi Li (and anyone else interested in plotting typhoon tracks),
Sorry for not responding sooner (I was in Egypt last week), but it looks
like Dennis has given you some helpful advice in the meantime. Like Rick
pointed out, the TypeOn boolean was used in the original script to
control the behavior for setting the polyline styles (colors and dashed
vs. solid) for the track lines depending on storm intensity and/or type
(e.g. tropical, subtropical, extratropical). Also, one reason your
script wasn't working is because each polyline segment must be
associated with a unique identifier, so the dummy variables need to be
created outside of the loop over each storm and be incremented
appropriately. Your script was creating them inside the loop so things
were getting "lost".
I've gone ahead and revised the original script to fix this and a few
other issues. It is somewhat similar to the unique_1.ncl example on the
NCL web site, but specific to Western Pacific storms (the unique_1.ncl
script plots a season's worth of Atlantic tropical cyclones and had some
fancier annotations, but used an older data format). The new script gets
rid of the TypeOn option and sets some some good default line styles for
Western Pacific systems. Note that in the JTWC best tracks, they tend to
use "disturbance" (DB) rather than "remnant low" or "tropical wave"
stages - I've set the line styles accordingly and updated the legend.
I've also added a Supertyphoon stage for tropical systems with winds at
or above 130 kt.
Also, a heads up on the ATCF b-deck files. The read method that Dennis
showed you (by columns) should work in your case, but be aware that the
b-decks aren't always formatted this nicely. Sometimes empty fields
aren't padded with empty spaces so the column alignment isn't consistent
from line to line. I checked your files and this won't be a problem in
this specific instance, but it could become an issue when trying to plot
from previous years or from agencies which don't make their b-decks look
pretty. Since the abr-decks are a comma-delimited format, it'd be best
to extract the information based on field number rather than by
character columns. I have some more sophisticated reader routines that I
could make available if there is interest - but I won't have time to do
this in the next month as I'm in the final stage of my PhD program. But
for right now, you should be fine.
The example plot was too big to send, but you can see it at:
http://euler.atmos.colostate.edu/~vigh/outgoing/WP2007.png
Hopefully this is helpful. You can tweak this further to suit your
purposes, but the basic structure is in place now.
Best regards,
Jonathan
;********************************************************
; typhoon_tracks.ncl
;
; The purpose of this script is to plot the best tracks
; for a given season of storms for the Western Pacific
; basin. All stages of the storm are plotted (subtropical
; storms, depressions, extratropical lows, etc.). The input
; data are individual b-decks files (ATCF format).
;
; This script is a somewhat simplified version of
; the unique_1.ncl script which can be found at:
; http://www.ncl.ucar.edu/Applications/Scripts/unique_1.ncl
; (that script plotted a season of Atlantic storms and did
; some fancier annotations than here, but did not read the newer
; ATCF format).
;
; This script was originally written by Jonathan Vigh with
; assistance from Mary Haley. Dennis Shea and Li Richard
; have also contributed.
;********************************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
;********************************************************
begin
diri = "C:/typhoon-best-track-data/2007-wp-full/"
; diri = "/home/vigh/SCRIPTS/NCL/TRACKS/TYPHOONS/"
fils = systemfunc ("cd " + diri + " ; ls *.txt")
nstorms = dimsizes(fils)
;********************************
; select sub-regions
;********************************
minlat = 0.0 ;min(lat) ; 0.0 ; deg N ; was 0
maxlat = 50.0 ;max(lat) ; 57.0 ; deg N ; was 57
minlon = 100.0 ; min(lon) -.;-100.0 ; - deg E = deg W ; was -100
maxlon = 165.0 ; max(lon) +10.; -10.0 ; - deg E = deg W ; was -20, then -15
;***** Define some NICE color maps to use.
colors_8 = (/"White","Black","Black","MediumPurple1","MediumPurple3","Blue1","CadetBlue3","Aquamarine2","Gold","Tan1","Sienna1","Tomato","VioletRed1", \
"Yellow","LimeGreen","Grey37","Red","Orange","GoldenRod1","DarkOrange","SteelBlue1","SlateBlue1","SlateGray1","LightSlateBlue","Magenta", \
"DodgerBlue","LightSteelBlue1","Moccasin","LightYellow","LemonChiffon1","CornSilk","LightGoldenrodYellow","Tan","PaleTurquoise3"/)
; Next, set the plotting parameters (which cases to plot)
PlotSubTropicalSegments = True ; If you wish to exclude all subtropical stages, set to False
PlotWaveSegments = True ; If you wish to exclude all wave segments, set to False
PlotLowSegments = True ; If you wish to exclude all remnant low segments, set to False
PlotExtraTropicalSegments = True ; If you wish to exclude all extratropical line sements, set to False
;********************************
; create plot
;********************************
wks = gsn_open_wks("ps","WP2007")
gsn_define_colormap(wks,colors_8)
res = True
res@gsnMaximize = True
res@gsnDraw = False ; so we can add poly stuff
res@gsnFrame = False ; do not advance frame
; res@gsnPaperOrientation = "portrait"
res@mpDataBaseVersion = "Ncarg4_1" ; Alias 'MediumRes'
res@mpDataSetName = "Earth..1"
; res@mpProjection = "LambertConformal"
; res@mpLambertParallel1F = 20.0
; res@mpLambertParallel2F = 40.0
; res@mpLambertMeridianF = -60.0
res@mpLimitMode = "LatLon"
res@mpMinLatF = minlat
res@mpMaxLatF = maxlat
res@mpMinLonF = minlon
res@mpMaxLonF = maxlon
res@mpFillOn = True ; False to turn off gray continents
res@mpOutlineOn = True ; turn on continental outline
res@mpOutlineBoundarySets = "AllBoundaries"
res@mpLandFillColor = "Tan" ; was "GoldenRod1"
res@mpInlandWaterFillColor = "PaleTurquoise3" ; was "LightBlue1" ; was "PaleTurquoise3"
res@mpOceanFillColor = "PaleTurquoise3" ; was "LightBlue1"
res@mpGeophysicalLineColor = "Grey37"
res@mpGeophysicalLineThicknessF = 0.5
res@mpUSStateLineColor = "Grey37"
res@mpUSStateLineThicknessF = 0.5
res@mpNationalLineColor = "Grey37"
res@mpNationalLineThicknessF = 0.5
res@mpGridAndLimbOn = "True"
res@mpGridAndLimbDrawOrder = "Draw"
res@mpGridMaskMode = "MaskLand"
res@mpGridSpacingF = 5.0
res@mpGridLineColor = "Grey37"
res@tmXBLabelFontHeightF = 0.012 ; 0.005
res@tmXBMajorLengthF = -0.001
res@pmTickMarkDisplayMode = "Always"
res@tiMainString = "2007 Season"
;***************************************
; plot base map *
;***************************************
plot = gsn_csm_map(wks,res) ; draw one of eight map projections
;***************************************
; get data and plot data for each storm
;***************************************
; create arrays of dummy graphic variables. This is required, b/c each line
; must be associated with a unique dummy variable.
dumline = new(nstorms*40*4,"graphic") ; allow for at least 40 days worth of data
dumblack = new(nstorms*40,graphic)
dumwhite = new(nstorms*40,graphic)
dumdaylabel = new(nstorms*40,graphic)
dumstormlabel = new(nstorms*2,graphic)
iline = 0
iblack = 0
iwhite = 0
idaylabel = 0
istormlabel = 0
do istorm = 0, nstorms-1 ; loop over each file
print(" istorm = " + istorm)
data = asciiread(diri+fils(istorm),-1,"string")
; printVarSummary(data)
nrows = dimsizes(data)
filename_chars = stringtochar(fils(istorm))
ssn = chartostring(filename_chars(3:4))
region = str_get_cols( data, 0, 1)
; print(region)
year = str_get_cols( data, 8, 11) ; year,...hour are type "string"
month = str_get_cols( data, 12, 13) ; Returns an array of substrings, given a start and end index
day = str_get_cols( data, 14, 15)
hour = str_get_cols( data, 16, 17)
vmax = stringtoint( str_get_cols( data, 48, 50) )
; print (vmax)
lat = stringtofloat(str_get_cols( data, 35, 37))*0.1
;if (str_get_cols( data, 38, 38).eq."S") then
; lat = -lat
;end if
; print(lat)
lon = stringtofloat(str_get_cols( data, 41, 44))*0.1
;if (str_get_cols( data, 45, 45).eq."W") then
; lon = -lon
;end if
; print(lon)
type = str_get_cols( data, 59, 60) ; "DB", "TD"
name = str_get_cols( data,149,158)
; print(type+" "+name)
; Now create arrays that only hold the points for the 00z locations
ind0 = ind(hour.eq."00") ; returns the integer indice if hour = "00"
; print(ind0)
; print(ind0(0))
if (.not.ismissing(ind0(0))) then ; if ind0(0) isn't missing, then use these values
lat0 = lat(ind0)
lon0 = lon(ind0)
end if
; print("lat0="+lat0+" lon0="+lon0)
ind12 = ind(hour.eq."12")
if (.not.ismissing(ind12(0))) then
lat12 = lat(ind12)
lon12 = lon(ind12)
day12 = day(ind12)
end if
; print("lat12="+lat12+" lon12="+lon12+" day12="+day12)
;***************************************
; Draw best track history as polylines *
;***************************************
res_poly = True ; polyline mods desired
do k=0,nrows-2
if(.not.ismissing(lat(k)))
res_poly@gsLineDashPattern = 0 ; default to solid line unless otherwise set below
res_poly@gsLineThicknessF = 4 ; default to a thick line unless otherwise set below
; Plot super typhoon segments
if((vmax(k).ge.130).and.(type(k).eq."ST")) then
res_poly@gsLineColor = "Purple"
dumline(iline) = gsn_add_polyline(wks,plot,lon(k:k+1),lat(k:k+1),res_poly)
iline = iline + 1
end if
; Plot hurricane segments
if((vmax(k).ge.64).and.(type(k).eq."TY")) then
res_poly@gsLineColor = "Red"
dumline(iline) = gsn_add_polyline(wks,plot,lon(k:k+1),lat(k:k+1),res_poly)
iline = iline + 1
end if
; Plot tropical storm segments
if((vmax(k).lt.64).and.(type(k).eq."TS").and.(vmax(k).ge.34)) then
res_poly@gsLineColor = "Yellow"
dumline(iline) = gsn_add_polyline(wks,plot,lon(k:k+1),lat(k:k+1),res_poly)
iline = iline + 1
end if
; Plot tropical depression segments
if((vmax(k).lt.34).and.(type(k).eq."TD")) then
res_poly@gsLineColor = "Green"
dumline(iline) = gsn_add_polyline(wks,plot,lon(k:k+1),lat(k:k+1),res_poly)
iline = iline + 1
end if
; Plot tropical disturbance segments
if((type(k).eq."DB").and.(vmax(k).lt.34)) then
res_poly@gsLineColor = "Brown"
res_poly@gsLineThicknessF = 2
res_poly@gsLineDashPattern = 2
res_poly@gsLineDashSegLenF = 0.1
dumline(iline) = gsn_add_polyline(wks,plot,lon(k:k+1),lat(k:k+1),res_poly)
iline = iline + 1
end if
; Plot tropical wave segments
if(PlotWaveSegments) then
if(type(k).eq."WV") then
res_poly@gsLineColor = "Brown"
res_poly@gsLineThicknessF = 2
res_poly@gsLineDashPattern = 2
res_poly@gsLineDashSegLenF = 0.1
dumline(iline) = gsn_add_polyline(wks,plot,lon(k:k+1),lat(k:k+1),res_poly)
iline = iline + 1
end if
end if
; Plot remanant low segments
if(PlotLowSegments) then
if(type(k).eq."LO") then
res_poly@gsLineColor = "Brown"
res_poly@gsLineThicknessF = 2
res_poly@gsLineDashPattern = 2
res_poly@gsLineDashSegLenF = 0.1
dumline(iline) = gsn_add_polyline(wks,plot,lon(k:k+1),lat(k:k+1),res_poly)
iline = iline + 1
end if
end if
if(PlotSubTropicalSegments) then
; Plot subtropical storm segments
if((vmax(k).gt.34).and.(type(k).eq."SS")) then
res_poly@gsLineColor = "DarkOrange"
dumline(iline) = gsn_add_polyline(wks,plot,lon(k:k+1),lat(k:k+1),res_poly)
iline = iline + 1
end if
; Plot subtropical depression segments
if((vmax(k).lt.34).and.(type(k).eq."SD")) then
res_poly@gsLineColor = "Blue1"
dumline(iline) = gsn_add_polyline(wks,plot,lon(k:k+1),lat(k:k+1),res_poly)
iline = iline + 1
end if
end if
; Plot extratropical segments
if(PlotExtraTropicalSegments) then
if(type(k).eq."EX") then
res_poly@gsLineColor = "Black"
res_poly@gsLineDashPattern = 2
res_poly@gsLineDashSegLenF = 0.1
dumline(iline) = gsn_add_polyline(wks,plot,lon(k:k+1),lat(k:k+1),res_poly)
iline = iline + 1
end if
end if
end if ; end check on whether the data are present
end do ; end loop over the times for each individual storm
;***************************************
; Draw black polymarkers at 00Z locations and white polymarkers at 12z locations
;***************************************
res_mark = True
res_mark@gsMarkerIndex = 1 ; polymarker style
res_mark@gsMarkerSizeF = 0.012 ; polymarker size - was 0.012
res_mark@gsMarkerColor = "Black" ; change marker color
do iday = 0, dimsizes(lon0)-1
dumblack(iblack) = gsn_add_polymarker(wks,plot,lon0(iday),lat0(iday),res_mark)
iblack = iblack + 1
end do
; Now draw white polymarkers at 12z locations
res_mark@gsMarkerIndex = 1
res_mark@gsMarkerSizeF = 0.012 ; polymarker size - was 0.012
res_mark@gsMarkerColor = "White" ; change marker color
do iday = 0, dimsizes(lon12)-1
dumwhite(iwhite) = gsn_add_polymarker(wks,plot,lon12(iday),lat12(iday),res_mark)
iwhite = iwhite + 1
end do
;***************************
; Plot some text labels *
;***************************
txres_day_labels = True
; Now draw the day labels at each 00Z marker
txres_day_labels@txFontHeightF = 0.006 ; 0.008
txres_day_labels@txJust = "TopCenter"
do iday=0,dimsizes(day12)-1
dumdaylabel(idaylabel) = gsn_add_text(wks,plot,day12(iday),lon12(iday),lat12(iday)-0.25,txres_day_labels)
idaylabel = idaylabel + 1
end do
; Draw storm number (actually these are strings) at beginning of track
txres_storm_labels = True
txres_storm_labels@txFontHeightF = 0.0054 ; 0.0044
txres_storm_labels@txFontThicknessF = 0.8
txres_storm_labels@txConstantSpacingF = 1.0
txres_storm_labels@txFontAspectF = 1.5
txres_storm_labels@txPerimOn = True
txres_storm_labels@txPerimColor = "Black"
txres_storm_labels@txPerimThicknessF = 1.0
txres_storm_labels@txPerimSpaceF = 0.4
txres_storm_labels@txBackgroundFillColor = "White"
txres_storm_labels@txJust = "CenterLeft"
dumstormlabel(istormlabel) = gsn_add_text(wks,plot,ssn,lon(0)+0.6,lat(0),txres_storm_labels)
istormlabel = istormlabel + 1
; Now draw them at the end
txres_storm_labels@txJust = "BottomRight"
dumstormlabel(istormlabel) = gsn_add_text(wks,plot,ssn,lon(nrows-1),lat(nrows-1),txres_storm_labels)
istormlabel = istormlabel + 1
; delete all arrays that may change size from storm to storm
delete(region)
delete(data)
delete(year)
delete(month)
delete(day)
delete(hour)
delete(vmax)
delete(lat)
delete(lon)
delete(type)
delete(name)
delete(ind0)
delete(lat0)
delete(lon0)
delete(lat12)
delete(lon12)
delete(day12)
delete(ind12)
end do ; end loop over each storm
;************************************************************************
; Create a legend that will be later added to the plot as an annotation *
;************************************************************************
lgres = True
lgres@lgAutoManage = False
lgres@vpWidthF = 0.1 ; was 0.08 ; width of legend (NDC)
lgres@vpHeightF = 0.1 ; was 0.08 ; height of legend (NDC)
; lgres@lgBottomMarginF = 0.17 ; was 0.25
lgres@lgPerimFill = 0 ; Use solid fill (0) instead of the default hollow fill
lgres@lgPerimFillColor = "Background"
lgres@lgBoxMajorExtentF = 0.4
lgres@lgBoxMinorExtentF = 0.2 ; controls how wide the box holding the legend items (lines and markers) can be in relation to legend
; lgres@lgBoxBackground = "PaleTurquoise3"
lgres@lgMonoItemType = False ; indicates that we wish to set the item types individually
lgres@lgMonoMarkerIndex = False
lgres@lgMonoLineThickness = False
lgres@lgMonoMarkerThickness = False
lgres@lgMonoMarkerSize = False
; Position fine elements of legend relative to some positional values:
xlegend = 0.680
ylegend = 0.390
lgres@lgLabelFont = 0
lgres@lgLabelFontHeightF = 0.04
lgres@lgLabelFontAspectF = 1.2
lgres@lgLabelConstantSpacingF = 0.0
lgres@lgItemCount = 11
lgres@lgLineDashSegLenF = 0.1
lgres@lgItemTypes = (/"Markers","Markers","Markers","Lines","Lines","Lines","Lines","Lines","Lines","Lines","Lines"/)
lgres@lgMarkerIndexes = (/ 1, 4, 16, 0, 0, 0, 0, 0, 0, 0, 0/)
lgres@lgLineThicknesses = (/ 0.1, 0.1, 0.1, 2.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0/)
lgres@lgMarkerColors = (/ "White", "Black", "Black"/)
lgres@lgMarkerSizes = (/ 0.0001, 0.004, 0.004/)
lgres@lgLineColors = (/ "White", "Black", "Black",\
"Brown", "Black", "Blue1", "DarkOrange", "Green", "Yellow", "Red" , "Purple" /) ; colors for legend lines
lgres@lgDashIndexes = (/ 0, 0, 0,\
2, 2, 0, 0, 0, 0, 0 , 0 /) ; dash indexes
legend_labels = (/"Tropical Cyclone No.","1200 UTC Position/Date","0000 UTC Position", \
"Disturbance (DB)","Extratropical (EX)","Subtropical Dep. (SD)","Subtropical Storm (SS)","Tropical Dep. (TD)","Tropical Storm (TS)","Typhoon (TY)","Supertyphoon (ST)"/)
legend = gsn_create_legend(wks,10,legend_labels,lgres)
yoffset_day = -0.0858
yoffset_storm = -0.0910
; Now draw a day label on the legend
txres = True
txres@txPerimOn = False
txres@txFontHeightF = 0.0022
txres@txPerimSpaceF = 0.0
txres@txJust = "CenterCenter"
txres@gsnDraw = False
;
; Note: the X and Y positions in this case don't matter, because the
; text will get repositioned in the "gsn_add_anotation" function.
;
number_text = gsn_create_text(wks,"21",txres)
; Draw storm number on the legend
txres@txFontHeightF = 0.0028
txres@txPerimOn = True
txres@txPerimColor = "Black"
txres@txPerimThicknessF = 1.0
txres@txPerimSpaceF = 0.4
txres@txBackgroundFillColor = "White"
three_text = gsn_create_text(wks,"3",txres)
;
; Add the various text and legend annotations to the plot.
;
amres = True
amres@amParallelPosF = 0.383
amres@amOrthogonalPosF = 0.470
amres@amJust = "BottomRight"
txid5 = gsn_add_annotation(plot,number_text,amres)
amres@amParallelPosF = 0.3815
amres@amOrthogonalPosF = 0.479
amres@amJust = "BottomRight"
txid6 = gsn_add_annotation(plot,three_text,amres)
amres@amParallelPosF = 0.49
amres@amOrthogonalPosF = 0.49
amres@amJust = "BottomRight"
annoid1 = gsn_add_annotation(plot,legend,amres)
; finally draw the plot and frame the workstation
draw(plot)
frame(wks)
end
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Thu Dec 3 14:14:58 2009
This archive was generated by hypermail 2.1.8 : Mon Nov 15 2010 - 12:48:08 MST
