; Plots storm tracks (best tracks) from tcvitals info from G5NR ; Storm exist over Pacific and Atlantic basins ; ; ; Andrew Kren ; June 16, 2017 begin ; begin user's settings ;--------------------------------------------------------- path_g5nr = "/scratch4/BMC/qosap/Andrew.Kren/ncl_scripts/G5NR/" ; list files in directory to determine how many dropsondes to read fils_g5nr = systemfunc("csh -c 'cd " + path_g5nr + " ; ls *.tc'") output_file = "storm_tracks_G5NR" ; output figure file output_type = "png" ; "x11", "ps", etc input_file_g5nr = path_g5nr + fils_g5nr scale_factor = 0.1 scale_factor_west = -0.1 plot_basemap = 0 ; if zero, first iteration of loop dumt2 = new(50000,graphic); used for plotting dumt3 = new(50000,graphic) dumt1 = new(50000,graphic) dumt4 = new(50000,graphic) dum = new(50000,graphic) duma = new(50000,graphic) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; end of user's settings ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; loop over each of the 25 storms, plotting the storm track at each iteration do i=0,dimsizes(fils_g5nr)-1 ; read in data for this storm data1 = readAsciiTable(input_file_g5nr(i),19,"string",1) ; 19 column array data_1d = ndtooned(data1) ; loop over data_1d and read data using str_get_cols count = num(.not.ismissing(data_1d)) lat = tofloat(str_get_cols(data_1d,116,118)) * scale_factor ; latitudes lon = tofloat(str_get_cols(data_1d,121,124)) ; longitudes deg_we = str_get_cols(data_1d,125,125) do j=0,count-1 if (deg_we(j).eq."E") then lon(j) = lon(j) * scale_factor ; deg W end if if (deg_we(j).eq."W") then lon(j) = lon(j) * scale_factor_west ; deg E end if end do pres = tofloat(str_get_cols(data_1d,135,138)) ; pressure in hPa at time and location storm_name = str_get_cols(data_1d,92,95) ; Storm name in TC vitals file storm_cat = tofloat(str_get_cols(data_1d,176,177)) ; Storm category, e.g. cat 1 hurricane, 2, 3, 0 is TD/TC ymd = str_get_cols(data_1d,102,109) ; YMD hr = str_get_cols(data_1d,111,112) ; HH date = ymd ; initialize month = ymd day = ymd do j=0,count-1 temp_concat = (/ymd(j),hr(j)/) date(j) = str_concat(temp_concat) month(j) = str_get_cols(date(j),4,5) day(j) = str_get_cols(date(j),6,7) end do if (plot_basemap.eq.0) then ;******************************** ; create plot ;******************************** wks = gsn_open_wks("png","storm_tracks_G5NR") ; send graphics to PNG file res = True res@gsnMaximize = True res@gsnDraw = False ; so we can add poly stuff res@gsnFrame = False ; do not advance frame res@mpDataBaseVersion = "MediumRes" ; 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 = 0.0 res@mpMaxLatF = 57.0 res@mpMinLonF = -100.0 res@mpMaxLonF = -20.0 res@mpFillOn = True ; False to turn off gray continents res@mpOutlineOn = True ; turn on continental outline res@mpOutlineBoundarySets = "AllBoundaries" res@mpLandFillColor = "GoldenRod1" res@mpInlandWaterFillColor = "PaleTurquoise3" res@mpOceanFillColor = "PaleTurquoise3" 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.005 res@tmXBMajorLengthF = -0.001 res@pmTickMarkDisplayMode = "Always" ;*************************************** ; plot base map * ;*************************************** plot = gsn_csm_map(wks,res) ; draw one of eight map projections end if res_poly = True counter_1 = 0 counter_2 = 0 do j=0,count-2 ; draw trajectories if (storm_cat(j).eq.0) then ; tropical storm res_poly@gsLineColor = "DarkGreen" dum(counter_1) = gsn_add_polyline(wks,plot,lon(j:j+1),lat(j:j+1),res_poly) end if if (storm_cat(j).eq.1) then ; Category 1 res_poly@gsLineColor = "Yellow" dum(counter_1) = gsn_add_polyline(wks,plot,lon(j:j+1),lat(j:j+1),res_poly) end if if (storm_cat(j).eq.2) then ; Category 2 res_poly@gsLineColor = "Orange" dum(counter_1) = gsn_add_polyline(wks,plot,lon(j:j+1),lat(j:j+1),res_poly) end if if (storm_cat(j).eq.3) then ; Category 3 res_poly@gsLineColor = "Red" dum(counter_1) = gsn_add_polyline(wks,plot,lon(j:j+1),lat(j:j+1),res_poly) end if if (storm_cat(j).eq.4 .or. storm_cat(j).eq.5) then ; Category 4 and 5 res_poly@gsLineColor = "DarkOrchid4" dum(counter_1) = gsn_add_polyline(wks,plot,lon(j:j+1),lat(j:j+1),res_poly) end if counter_1 = counter_1 + 1 end do ; draw polymarkers at every other time step res_mark = True res_mark@gsMarkerIndex = 1 ; polymarker style res_mark@gsMarkerSizeF = 0.012 ; polymarker size res_mark@gsMarkerColor = "Black" ; change marker color do j=0,count-2,2 ; every other plot duma(counter_1) = gsn_add_polymarker(wks,plot,lon(j),lat(j),res_mark) counter_1 = counter_1 + 1 end do ; add starting and ending date at beginning of track txres = True txres@txFontHeightF = 0.0035 txres@txJust = "TopCenter" txres@txPerimOn = True txres@txPerimColor = "Black" txres@txPerimThicknessF = 1.0 txres@txPerimSpaceF = 0.4 txres@txBackgroundFillColor = "White" dumt1(counter_2) = gsn_add_text(wks,plot,month(0)+"/"+day(0)+" "+hr(0)+"z",lon(0),lat(0)-0.45,txres) dumt4(counter_2) = gsn_add_text(wks,plot,month(0)+"/"+day(0)+" "+hr(0)+"z",lon(count-1),lat(count-1)-0.45,txres) ; draw storm name at beginning and end txres@txFontHeightF = 0.0044 txres@txJust = "CenterLeft" txres@txPerimOn = True txres@txPerimColor = "Black" txres@txPerimThicknessF = 1.0 txres@txPerimSpaceF = 0.4 txres@txBackgroundFillColor = "White" dumt2(counter_2) = gsn_add_text(wks,plot,storm_name(2),lon(2),lat(2)+0.6,txres) dumt3(counter_2) = gsn_add_text(wks,plot,storm_name(count-3),lon(count-3),lat(count-2)+0.6,txres) ; delete arrays before going to next file since TC vital file length will vary delete([/lat,lon,pres,storm_name,storm_cat,temp_concat,date,ymd,hr,data1,data_1d,count,deg_we,month,day/]) plot_basemap = plot_basemap + 1 counter_2 = counter_2 + 1 end do draw(plot) frame(wks) print("GREAT SUCCESS!!") end