I'm so sorry that I haven't noticed the code has messed up. I have attached my code and data in the attachment. Thanks
--- 2009年12月1日 星期二，Rick Brownrigg <brownrig@ucar.edu> 寫道﹕
寄件人: Rick Brownrigg <brownrig@ucar.edu>
主題: Re: [ncl-talk] Multiple tracks
收件人: "Li Richard" <lcy114@yahoo.com.hk>, ncl-talk@ucar.edu
日期: 2009年12月1日,星期二,上午3:51
Hi,
This script is hard to read, as in my mailer, the indentation reflecting code structure is all messed up. However, I *think* what is going on is that all the polyline drawing code surrounded by an if-clause of the form:
"if (.not.TypeOn) then"...
but prior to that, TypeOn is set to True, so the polyline drawing code would never execute.
Rick
On Tue, 1 Dec 2009 11:14:23 +0800 (HKT)
Li Richard <lcy114@yahoo.com.hk> wrote:
> Hi everyone,
> Below is the code of multiple typhoon best tracks. However, I can only get all the polymarkers plotted, but not the polylines. I wonder if it is related to the delete function at the end of the code. So how can i further modify the code in order to get the multiple polylines drawn on the same plot? The data used can be seen in the attachment.Thanks.
> Regards
> Richard
>
> >>
> 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/"
> fils = systemfunc ("cd " + diri + " ; ls *.txt")
> print(fils)
> print(dimsizes(fils))
> nfil = dimsizes(fils)
> print(dimsizes(nfil))
> ;fili = "bwp182007.txt"
> ;********************************
> ; 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)
> TypeOn = True ; If you just want the legend to have the basic (H, TS, and TD) designators (* and S), set TypeOn to False
> ; NOTE -- this will cause the program to grab a different input file -- one in which all the ; extratropical (E), remnant low (L), and wave cases (W) have been taken out; additionally,
> ; there will not be a distinction between the S and * cases (they will just be plotted according to
> ; intensity
> PlotSubTropicalSegments = True ; If you wish to exclude all subtropical line segments, set to False ; Note: if TypeOn is set to False, the following (W, L, E) cases won't be plotted regardless of their value 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("x11","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 = ssn + " " + cyclone_name + " (" + year(0) + month(0) + ")"
> ;***************************************
> ; plot base map *
> ;***************************************
>
> plot = gsn_csm_map(wks,res) ; draw one of eight map projections ;********************************
> ; get data
> ;********************************
> do nf=0,nfil-1,1 ; loop over each file
> ; do nf=0,2 ; loop over each file
> data = asciiread(diri+fils(nf), -1,"string")
> printVarSummary(data)
> ;data = asciiread(diri+fili, -1,"string") ; Info includes type so cases can be weeded out
>
> ;tChr = stringtochar(data)
> nrows = dimsizes(data)
> 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) don't have a missing value
> 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
> ; create array of dummy graphic variables. This is required, b/c each line
> ; must be associated with a unique dummy variable.
> dum = new(dimsizes(hour),graphic) ic = 0
> do k=0,nrows-2
> if(.not.ismissing(lat(k)))
> res_poly@gsLineDashPattern = 0
> res_poly@gsLineThicknessF = 4 if(.not. TypeOn) then
> ; Plot hurricane segments
> if((vmax(k).gt..64).and.(type(k).eq."*")) then
> res_poly@gsLineColor = "Red" dum(ic) = gsn_add_polyline(wks,plot,lon(k:k+1),lat(k:k+1),res_poly)
> ic = ic + 1
> end if
> ; Plot tropical storm segments if((vmax(k).lt..64).and.(vmax(k).gt.34).and.(type(k).eq."*")) then
> res_poly@gsLineColor = "Yellow" dum(ic) = gsn_add_polyline(wks,plot,lon(k:k+1),lat(k:k+1),res_poly) ic = ic + 1
> end if ; Plot tropical depression segments
> if((vmax(k).lt..34).and.(type(k).eq."*")) then
> res_poly@gsLineColor = "LimeGreen" dum(ic) = gsn_add_polyline(wks,plot,lon(k:k+1),lat(k:k+1),res_poly)
> ic = ic + 1
> end if ; Now plot the subtropical segments with dashed lines
> ; Plot hurricane segments of subtropical systems with dashed red lines
> if((vmax(k).gt..64).and.(type(k).eq."S")) then
> res_poly@gsLineColor = "Red" res_poly@gsLineDashPattern = 2
> res_poly@gsLineDashSegLenF = 0.1
> dum(ic) = gsn_add_polyline(wks,plot,lon(k:k+1),lat(k:k+1),res_poly)
> ic = ic + 1
> end if
> ; Plot tropical storm segments of subtropical systems with dashed yellow lines if((vmax(k).lt..64).and.(vmax(k).gt.34).and.(type(k).eq."S")) then
> res_poly@gsLineColor = "Yellow"
> res_poly@gsLineDashPattern = 2
> res_poly@gsLineDashSegLenF = 0.1 dum(ic) = gsn_add_polyline(wks,plot,lon(k:k+1),lat(k:k+1),res_poly) ic = ic + 1
> end if ; Plot tropical depression segments of subtropical systems with dashed dark LimeGreen lines
> if((vmax(k).lt..34).and.(type(k).eq."S")) then
> res_poly@gsLineColor = "LimeGreen" res_poly@gsLineDashPattern = 2
> res_poly@gsLineDashSegLenF = 0.1 dum(ic) = gsn_add_polyline(wks,plot,lon(k:k+1),lat(k:k+1),res_poly)
> ic = ic + 1
> end if else ; if(TypeOn) -- plot pure tropical segments as solid lines (red for hurricane, etc.). ; Plot super typhoon segments
> if((vmax(k).ge..64).and.(type(k).eq."ST")) then
> res_poly@gsLineColor = "Purple" dum(ic) = gsn_add_polyline(wks,plot,lon(k:k+1),lat(k:k+1),res_poly)
> ic = ic + 1
> end if
> ; Plot hurricane segments
> if((vmax(k).ge..64).and.(type(k).eq."TY")) then
> res_poly@gsLineColor = "Red" dum(ic) = gsn_add_polyline(wks,plot,lon(k:k+1),lat(k:k+1),res_poly)
> ic = ic + 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 = "Orange" dum(ic) = gsn_add_polyline(wks,plot,lon(k:k+1),lat(k:k+1),res_poly) ic = ic + 1
> end if ; Plot tropical depression segments
> if((vmax(k).lt..34).and.(type(k).eq."TD")) then
> res_poly@gsLineColor = "Yellow" dum(ic) = gsn_add_polyline(wks,plot,lon(k:k+1),lat(k:k+1),res_poly)
> ic = ic + 1
> end if ; Plot tropical disturbance segments
> if((type(k).eq.."DB").and.(vmax(k).lt.34)) then
> res_poly@gsLineColor = "Green" dum(ic) = gsn_add_polyline(wks,plot,lon(k:k+1),lat(k:k+1),res_poly)
> ic = ic + 1 end if if(PlotWaveSegments) then ; Plot tropical wave segments
> if(type(k).eq."WV") then
> res_poly@gsLineColor = "LimeGreen" res_poly@gsLineDashPattern = 2
> res_poly@gsLineDashSegLenF = 0.1
> dum(ic) = gsn_add_polyline(wks,plot,lon(k:k+1),lat(k:k+1),res_poly)
> ic = ic + 1
> end if end if
> if(PlotLowSegments) then ; Plot remanant low segments
> if(type(k).eq."LO") then
> res_poly@gsLineColor = "SlateBlue" res_poly@gsLineDashPattern = 2
> res_poly@gsLineDashSegLenF = 0.1
> dum(ic) = gsn_add_polyline(wks,plot,lon(k:k+1),lat(k:k+1),res_poly)
> ic = ic + 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" dum(ic) = gsn_add_polyline(wks,plot,lon(k:k+1),lat(k:k+1),res_poly) ic = ic + 1
> end if ; Plot subtropical depression segments
> if((vmax(k).lt..34).and.(type(k).eq."SD")) then
> res_poly@gsLineColor = "Blue1" dum(ic) = gsn_add_polyline(wks,plot,lon(k:k+1),lat(k:k+1),res_poly)
> ic = ic + 1
> end if end if
>
> if(PlotExtraTropicalSegments) then ; Plot extratropical segments
> if(type(k).eq."EX") then
> res_poly@gsLineColor = "Black" res_poly@gsLineDashPattern = 2
> res_poly@gsLineDashSegLenF = 0.1
> dum(ic) = gsn_add_polyline(wks,plot,lon(k:k+1),lat(k:k+1),res_poly)
> ic = ic + 1
> end if end if
> end if
> end if
> end do
>
> ;***************************************
> ; Draw black dots at lat/lon itersections
> ;***************************************
> ;***************************************
> ; 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.02 ; polymarker size - was 0.012
> res_mark@gsMarkerColor = "Black" ; change marker color
> duma = new(1,graphic)
> duma = gsn_add_polymarker(wks,plot,lon0,lat0,res_mark)
> ; Now draw white polymarkers at 12z locations
> res_mark@gsMarkerIndex = 1
> res_mark@gsMarkerSizeF = 0.02 ; polymarker size - was 0.012
> res_mark@gsMarkerColor = "White" ; change marker color
> dumb = new(1,graphic)
> dumb = gsn_add_polymarker(wks,plot,lon12,lat12,res_mark)
>
> ;***************************
> ; Plot some text labels *
> ;***************************
> txres = True
> ; Now draw the day labels at each 00Z marker
> txres@txFontHeightF = 0.015 ; 0.008
> txres@txJust = "TopCenter"
> dumz = new(nrows,graphic)
> n12 = dimsizes(day12)
> ;;do k=0,nrows-1
> do k=0,n12-1
> dumz(k) = gsn_add_text(wks,plot,day12(k),lon12(k),lat12(k)-0.25,txres)
> end do
>
>
> ; Draw storm number (actually these are strings) at beginning of track txres@txFontHeightF = 0.02 ; 0.0044
> txres@txJust = "CenterLeft"
> txres@txPerimOn = True
> txres@txPerimColor = "Black"
> txres@txPerimThicknessF = 1.0
> txres@txPerimSpaceF = 0.4
> txres@txBackgroundFillColor = "White"
> ; gsn_text(wks,plot,ssn,lon(0)+0.6,lat(0),txres)
>
> ; Now draw them at the end
> txres@txJust = "BottomRight"
> ; gsn_text(wks,plot,ssn,lon(nrows-1),lat(nrows-1),txres)
>
> if (.not. TypeOn) then
> txres@txFontHeightF = 0.02 ; 0.006
> txres@txJust = "CenterCenter"
> txres@txPerimSpaceF = 0.7
> txres@txPerimThicknessF = 1.5
> txres@txConstantSpacingF = 1.0
> txres@txFontAspectF = 1.5
> txres@txFontThicknessF = 0.8
> gsn_text_ndc(wks,"~F02~ Remnant low, tropical wave, and ~C~ extratropical stages are not shown. ~C~~C~ Subtropical stages are dashed.",0.715,0.680,txres)
> end if
> txres@txFontHeightF = 0.005
> txres@txPerimThicknessF = 1.5
> txres@txPerimSpaceF = 1.0
> txres@txConstantSpacingF = 1.0
> txres@txFontAspectF = 1.5
> txres@txFontThicknessF = 0.8
>
> ;******************
> ; Plot a legend *
> ;******************
> lgres = False
> lgres@lgAutoManage = False
> lgres@vpWidthF = 0.2 ; was 0.1 ; width of legend (NDC)
> lgres@vpHeightF = 0..2 ; was 0.1 ; 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.66 ; 0.820 ylegend = 0.35 ; 0.290
>
> if(TypeOn) then
> lgres@lgLabelFont = 0
> lgres@lgLabelFontHeightF = 0..012 ; 0.004
> 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, 4.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",\
> "SlateBlue1", "LimeGreen", "Black", "Blue1", "Dark Orange", "LimeGreen", "Yellow", "Red" /) ; colors for legend lines
> lgres@lgDashIndexes = (/ 0, 0, 0,\
> 2, 2, 2, 0, 0, 0, 0, 0 /) ; dash indexes
> ; gsn_legend_ndc(wks,10, (/"Tropical Cyclone No.","1200 UTC Position/Date","0000 UTC Position",\
> ; "Remnant Low","Tropical Wave","Extratropical","Subtropical Dep.","Subtropical Storm","Tropical Dep.","Tropical Storm (T)","Hurricane (H)"/),xlegend,ylegend,lgres)
> yoffset_day = -0.0858
> yoffset_storm = -0.0910
> else
> lgres@lgLabelFont = 0
> lgres@lgLabelFontHeightF = 0..004
> lgres@lgLabelFontAspectF = 1..2
> lgres@lgLabelConstantSpacingF = 0.0
> lgres@lgItemCount = 6
> lgres@lgItemTypes = (/"Markers","Markers","Markers","Lines" ,"Lines" ,"Lines"/)
> lgres@lgLineColors = (/"White", "Black" ,"Black" ,"LimeGreen","Yellow","Red" /) ; colors for legend lines
> lgres@lgLineThicknesses = (/ 0.1, 0.1, 0.1, 4.0, 4.0, 4.0/)
> lgres@lgDashIndexes = (/ 0, 0, 0, 0, 0, 0/) lgres@lgMarkerIndexes = (/ 0, 4, 16, 0, 0, 0/) lgres@lgMarkerColors = (/"White", "Black" , "Black"/)
> lgres@lgMarkerSizes = (/ 0.0001, 0.004, 0.004/)
> ; gsn_legend_ndc(wks,5,(/"Tropical Cyclone No.","1200 UTC Position/Date","0000 UTC Position",\
> ; "Tropical Dep.","Tropical Storm (T)","Hurricane (H)"/),xlegend,ylegend,lgres) yoffset_day = -0.0750 ; offsets for the case where we don't plot all the extra cases
> yoffset_storm = -0.0875
> end if
> ; Now draw a day label on the legend
> txres@txPerimOn = False
> txres@txFontHeightF = 0.01
> txres@txPerimSpaceF = 0.0
> txres@txJust = "CenterCenter"
> ; gsn_text_ndc(wks,"21",xlegend + 0.0183, ylegend + yoffset_day,txres)
>
> ; Draw storm number on the legend txres@txFontHeightF = 0.004
> txres@txPerimOn = True
> txres@txPerimColor = "Black"
> txres@txPerimThicknessF = 1.0
> txres@txPerimSpaceF = 0.4
> txres@txBackgroundFillColor = "White"
> ; gsn_text_ndc(wks,"3",xlegend + 0.0142, ylegend + yoffset_storm,txres) ; drawNDCGrid(wks) ; Use to put on a grid for fine-tuning placement
> delete(region)
> delete(data) ; delete all arrays that may change size
> 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)
> delete(dum)
> delete(dumz)
> end do
> draw(wks)
> frame(wks)
> end
>
>
>
> Yahoo!香港提供網上安全攻略，教你如何防範黑客! 請前往 http://hk.promo.yahoo.com/security/ 了解更多!
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
This archive was generated by hypermail 2.1.8 : Thu Dec 03 2009 - 09:53:22 MST