; Example script to produce plots for a WRF real-data run, ; with the ARW coordinate dynamics option. load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/NCL_GETools.ncl" begin year = 2012 month = 6 day = 22 dayst = day dirname = "/glade/p/mmm/barthm/IDL/Flight-Tracks/" ; convert single digit months, days, and hours to double digits if(month .lt. 10) then cmon = "0"+month else cmon = month end if if(day .lt. 10) then cday = "0"+day else cday = day end if if(dayst .lt. 10) then cdayst = "0"+dayst else cdayst = dayst end if ; create output file name fout = "DBZ_nexrad_inflow2_"+cmon+"-"+cday type = "ps" wks = gsn_open_wks(type,fout) setvalues NhlGetWorkspaceObjectId() "wsMaximumSize" : 120000000 end setvalues ; Set some Basic Plot Information ; We would like to switch off the Initial Time and Footer info on the plot ; We also do not want titles on the color label bar latmin = 38.8 latmax = 42.5 lonmin = 360.-105.34 lonmax = 360.-99.36 centerlon = (lonmin + lonmax)/2. llonmin = -105.34 llonmax = -99.36 txres = True txres@txFontHeightF = 0.015 txres@txFont = 21 txres@txJust = "CenterLeft" txres@txPerimOn = False res = True res@gsnFrame = False res@gsnDraw = False ; res@gsnAddCyclic = False res@mpMinLatF = latmin ; min(lat2d) res@mpMaxLatF = latmax ; max(lat2d) res@mpMinLonF = lonmin ; min(lon2d) res@mpMaxLonF = lonmax ; max(lon2d) res@mpCenterLonF = centerlon pltres = True pltres@PanelPlot = True pltres@CommonTitle = True pltres1 = True pltres1@FramePlot = False mpres = res mpres@mpFillOn = False mpres@mpGridAndLimbOn = False mpres@mpOutlineOn = True mpres@mpDataBaseVersion = "MediumRes" mpres@mpOutlineBoundarySets = "GeophysicalAndUSStates" mpres@mpGeophysicalLineColor = "Black" mpres@mpNationalLineColor = "Black" mpres@mpUSStateLineColor = "Black" mpres@mpGridLineColor = "Blue" mpres@mpLimbLineColor = "Green" mpres@mpPerimLineColor = "Black" mpres@mpUSStateLineThicknessF = 2.0 mpres@mpGridLatSpacingF = 1. mpres@mpGridLonSpacingF = 1. mpres@tmXBLabelFontHeightF = 0.015 mpres@tmYLLabelFontHeightF = 0.015 mpres@tmXBMode = "Explicit" mpres@tmXBValues = (/255,256,257,258,259,260/) mpres@tmXBLabels = (/"105W","104W","103W","102W","101W","100W"/) res@tmXBMinorValues = ispan(255,260,1) mpres@tmYLMode = "Explicit" mpres@tmYLValues = (/39,40,41,42/) mpres@tmYLLabels = (/"39N","40N","41N","42N"/) res@tmYLMinorValues = ispan(39,42,1) plots = new ( 2, graphic ) ; Change color map to something that has a grey scale gsn_define_colormap(wks,"radar_1") opts = mpres opts@lbTitleOn = True opts@lbTitleString = "" opts@lbTitleFontHeightF = 0.015 opts@cnFillOn = True ; turn on color opts@cnLinesOn = False ; no contour lines opts@lbLabelFontHeightF = 0.015 opts@lbLabelBarOn = True ;False opts@cnLevelSelectionMode = "ManualLevels" ; set manual contour levels opts@cnMinLevelValF = 5 ; set min contour level opts@cnMaxLevelValF = 75 ; 80 ; set max contour level opts@cnLevelSpacingF = 5 ; contour spacing ; opts@gsnSpreadColorStart = 0 ; with colormap "radar" opts@gsnSpreadColorStart = 9 ; with colormap "radar1" ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; ; The NEXRAD input files after decluttered by Cameron Homeyer's IDL routine ifs = 23 ife = 23 do ifl=ifs,ife hr = ifl wrfhr = ifl minutes = 15 if(minutes .lt. 10) then minute = "0"+minutes else minute = minutes end if if(hr .lt. 10) then hour = "0"+hr else hour = hr end if if(hr .gt. 23) then day = dayst + 1 hr = hr-24 if(hr .lt. 10) then delete(hour) hour = "0"+hr else hour = hr end if end if if(month .lt. 10) then cmon = "0"+month else cmon = month end if if(day .lt. 10) then cday = "0"+day else cday = day end if sdate = year+"-"+cmon+"-"+cday+" "+hour+":"+minute+" UTC" s_start = "0000" s_stop = "0000" if(minute .eq. "00") then ; 10 minutes before / after time s_start = "25:50" s_stop = "25:10" else minute_10 = minutes-10 minutep10 = minutes+10 s_start = hour+":"+minute_10 s_stop = hour+":"+minutep10 end if ; nexrad_radar_20120623_0120.nc filename = "nexrad_radar_"+year+cmon+cday+"_"+hour+minute+".nc" fname1 = dirname+filename print(" reading from file "+fname1) a = addfile(fname1, "r") time = 0 lon1d = a->Longitude(0,:) ; (time,nx) lat1d = a->Latitude(0,:) ; (time,ny) alt1d = a->Altitude(0,:) ; (time,nz) refl = a->Reflectivity(0,:,:,:) ; (time, nz,ny,nx) diml = dimsizes(refl) ; print(" dimsizes refl "+diml) nz = diml(0) ny = diml(1) nx = diml(2) lon2d = new( (/ny,nx/), float) lat2d = new( (/ny,nx/), float) do j=0,ny-1 lon2d(j,:) = lon1d end do do i=0,nx-1 lat2d(:,i) = lat1d end do lat2d@units = "degrees_N" lon2d@units = "degrees_E" refl@lon2d = lon2d refl@lat2d = lat2d mrefl = new( (/ny,nx/), float) mrefl = 0. do i = 0,nx-1 do j = 0,ny-1 mrefl(j,i) = max(refl(:,j,i)) end do end do mrefl@description = "f) Maximum Reflectivity at "+sdate mrefl@units = "" ; "(dBZ)" mrefl@lon2d = lon2d mrefl@lat2d = lat2d print(" max refl of mrefl "+max(mrefl)) plot = gsn_csm_contour_map_ce(wks, mrefl, opts) draw(plot) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Change color map for aircraft data gsn_define_colormap(wks,"BlAqGrYeOrReVi200") ; set colors for overlaying aircraft data onto radar map colors = new( (/10/), integer) do imkr = 0,9 colors(imkr) = 11 + imkr*20 ; 11, 31, 51, 71, 91, 111, 131, 151, 171, 191 end do colors(0) = 0 colors(1) = 25 ; put red colors in middle of colormap, blue colors at end, green at beginning to be able to see aircraft data on top of radar data colors(0:7) = colors(2:9) ; colors(5) = 11 + 4*20 ; colors(6) = 11 + 3*20 ; colors(7) = 11 + 2*20 colors(8) = 25 colors(9) = 11 + 1*20 ; print(" colors "+colors) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Get DC-8 Aircraft Data ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; vars = (/ "co", "o3", "h2o2", "hcho", "hno3", "ch3ooh", "no", "kappa" /) spcs = (/ "CO", "O3", "H2O2", "HCHO", "HNO3", "CH3OOH", "NO", "Kappa" /) idc8 = (/ 62, 71, 257, 67, 260, 254, 70, 236 /) ; NO2 = 72 TDLIF, 68 ESRL dc8scl=(/ 1., 1., .001, .001, .001, .001, 1., 1. /) ; scale to ppbv from data file units ilvls =(/ 1, 1, 0, 0, 0, 0, 0, 2 /) ivar = 7 csp = spcs(ivar) var = vars(ivar) print(" variable: "+var+" "+csp+" "+idc8(ivar)+" "+dc8scl(ivar)+" "+ilvls(ivar) ) obsdirname = "/glade/p/mmm/barthm/DC3/Aircraft/DC8_NASA_Merge/" fnamedc8 = "dc3-mrg01-dc8_merge_2012"+cmon+cdayst+"_R7.ict" print("DC8 file: "+obsdirname+fnamedc8) ictVar_DC8 = read_ict(obsdirname+fnamedc8) dc8_time_s = doubletofloat(ictVar_DC8(:,0)) ; Time, s dc8_u_s = doubletofloat(ictVar_DC8(:,43)) ; U_MMS dc8_v_s = doubletofloat(ictVar_DC8(:,44)) ; V_MMS dc8_lat_s = doubletofloat(ictVar_DC8(:,48)) ; GPS_LAT_MMS dc8_lon_s = doubletofloat(ictVar_DC8(:,49)) ; GPS_LON_MMS dc8_galt_s = doubletofloat(ictVar_DC8(:,50)) ; GPS_ALT_MMS iv = idc8(ivar) dc8_co_s = doubletofloat(ictVar_DC8(:,iv)) if(var .eq. "no") then if(month .eq. 6 .and. day .ge. 11) then dc8_no2_s = doubletofloat(ictVar_DC8(:,68)) ; ESRL data starting 6/11/2012, ppbv else dc8_no2_s = doubletofloat(ictVar_DC8(:,72)) ; TDLIF data dc8_no2_s = dc8_no2_s * 0.001 ; ppbv end if dc8_co_s = dc8_co_s + dc8_no2_s ; ppbv end if dc8_var_s = dc8_co_s *dc8scl(ivar) dc8_var_s@_FillValue = -999. dc8_galt_s@_FillValue = -999. dc8_lon_s@_FillValue = -999. dc8_lat_s@_FillValue = -999. print(" max DC8 "+var+" "+max(dc8_var_s) ) avgtime = 30. iavgt = 30 dimda = dimsizes(dc8_galt_s) dimdc = dimsizes(dc8_var_s) ndc8_min = floattointeger(dimdc(0)/avgtime) dc8_time = new( (/ndc8_min/), float) dc8_galt = new( (/ndc8_min/), float) dc8_var = new( (/ndc8_min/), float) dc8_lat = new( (/ndc8_min/), float) dc8_lon = new( (/ndc8_min/), float) dc8_u = new( (/ndc8_min/), float) dc8_v = new( (/ndc8_min/), float) dc8_galt@_FillValue = -999. dc8_var@_FillValue = -999. dc8_lat@_FillValue = -999. dc8_lon@_FillValue = -999. dc8_u@_FillValue = -999. dc8_v@_FillValue = -999. do im = 0,ndc8_min-1 i1 = im*iavgt i2 = i1+iavgt-1 i3 = i1+iavgt/2 if(all(ismissing(dc8_var_s(i1:i2)))) then dc8_time(im) = avg(dc8_time_s(i1:i2)) dc8_galt(im) = avg(dc8_galt_s(i1:i2)) dc8_var(im) = dc8_var_s@_FillValue dc8_lat(im) = avg(dc8_lat_s(i1:i2)) dc8_lon(im) = avg(dc8_lon_s(i1:i2)) dc8_u(im) = avg(dc8_u_s(i1:i2)) dc8_v(im) = avg(dc8_v_s(i1:i2)) else dc8_time(im) = avg(dc8_time_s(i1:i2)) dc8_galt(im) = avg(dc8_galt_s(i1:i2)) dc8_var(im) = avg(dc8_var_s(i1:i2)) dc8_lat(im) = avg(dc8_lat_s(i1:i2)) dc8_lon(im) = avg(dc8_lon_s(i1:i2)) dc8_u(im) = avg(dc8_u_s(i1:i2)) dc8_v(im) = avg(dc8_v_s(i1:i2)) ; print(" "+im+" "+dc8_time_s(i3)+" "+dc8_galt(im)+" "+dc8_var(im) +" "+dc8_lat(im) +" "+dc8_lon(im) ) end if end do ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Add vectors to the plots ; Do this before aircraft data so that they are below the aircraft data on the plot u_stm = 0. ; It may be important to give a storm-relative wind vector to show air is going into storm. However, v_stm = 0. ; we would need to determine storm speed ahead of time. I don't think it is an issue with 6/22. dc8_u = dc8_u - u_stm dc8_v = dc8_v - v_stm startsec = 23.*3600. + 10.*60. ; from inflow times stopsec = 23.*3600. + 20.*60. ; from inflow times c = ind(dc8_time .ge. startsec) im3 = c(0) delete(c) c = ind(dc8_time .ge. stopsec) im4 = c(0) delete(c) wmsetp("vrs - reference vector size",1.) wmsetp("vrn - NDC size corresponding to vrs",0.002) wmsetp("vcw - vector line width scale",3.) ; wmsetp("vcc - vector color", 200) wmsetp("vch - arrow head size",0.015) ; Either you can plot all the vectors (at 30-s intervals) for inflow time period, or every other vector. Here, I do the second method and comment out the first method. vind = ispan(im3,im4,2) wmvectmap(wks, dc8_lat(vind), dc8_lon(vind), dc8_u(vind), dc8_v(vind)) ; Plot vectors. ;wmvectmap(wks, dc8_lat(im3:im4), dc8_lon(im3:im4), dc8_u(im3:im4), dc8_v(im3:im4)) ; Plot vectors. ; wspd = sqrt(dc8_u(im3:im4)^2 + dc8_v(im3:im4)^2) ; print(" wind "+dc8_time(im3:im4)+" "+dc8_lat(im3:im4) +" "+dc8_lon(im3:im4)+" "+dc8_u(im3:im4) +" "+dc8_v(im3:im4) +" "+wspd ) wspd = sqrt(dc8_u(vind)^2 + dc8_v(vind)^2) print(" wind "+dc8_time(vind)+" "+dc8_lat(vind) +" "+dc8_lon(vind)+" "+dc8_u(vind) +" "+dc8_v(vind) +" "+wspd ) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Plot var data gsres = True ; find start, stop times for hour of interest wrfsec = wrfhr*3600. + minutes*60. ; This is set above to 23:15 for radar image startsec = wrfsec - 600. ; inflow was at 23:10-23:20 (see wind vector times). stopsec = wrfsec + 600. ; --> will keep these times (23:05-23:25) for flight overlay c = ind(dc8_time .ge. startsec) im1 = c(0) delete(c) c = ind(dc8_time .ge. stopsec) im2 = c(0) delete(c) if(ismissing(im2)) then im2 = ndc8_min-1 end if print(" wrfsec,start,stop, im1,im2 "+wrfsec+" "+startsec+" -- "+stopsec+" "+im1+" , "+im2 ) print(" min/max var during inflow "+min(dc8_var(im1:im2))+" "+max(dc8_var(im1:im2)) ) do im = im1,im2 ; 0,ndc8_min-1 if(ismissing(dc8_var(im))) then print(" var "+im+" "+dc8_time(im)+" "+dc8_galt(im)+" "+dc8_lat(im) +" "+dc8_lon(im)+" "+dc8_var(im) +" " ) else gsres@gsMarkerIndex = 4 ; 8=upward triangle; 4=open circle; 16=filled dots for markers. gsres@gsMarkerThicknessF = 4. gsres@gsMarkerColor = "black" gsn_polymarker(wks,plot,dc8_lon(im),dc8_lat(im),gsres) gsres@gsMarkerIndex = 16 ; 7=downward triangle; 4=open circle; 16=filled dots for markers. gsres@gsMarkerThicknessF = 2. ; contour levels indices are defined at the top (in "table") if(ilvls(ivar) .eq. 0) then contourlvls = (/ .05, .07, .10, .12, .15, .17, .2, .25, .3 /) else if(ilvls(ivar) .eq. 1) then contourlvls = (/ 50., 70., 80., 90., 100., 110., 120., 130., 150./) else contourlvls = (/ .1, .2, .5, .7, .9, 1.0, 1.2, 1.3, 1.5/) end if end if gsres@gsMarkerColor = colors(0) if(dc8_var(im) .ge. contourlvls(0)) then gsres@gsMarkerColor = colors(1) end if if(dc8_var(im) .ge. contourlvls(1)) then gsres@gsMarkerColor = colors(2) end if if(dc8_var(im) .ge. contourlvls(2)) then gsres@gsMarkerColor = colors(3) end if if(dc8_var(im) .ge. contourlvls(3)) then gsres@gsMarkerColor = colors(4) end if if(dc8_var(im) .ge. contourlvls(4)) then gsres@gsMarkerColor = colors(5) end if if(dc8_var(im) .ge. contourlvls(5)) then gsres@gsMarkerColor = colors(6) end if if(dc8_var(im) .ge. contourlvls(6)) then gsres@gsMarkerColor = colors(7) end if if(dc8_var(im) .ge. contourlvls(7)) then gsres@gsMarkerColor = colors(8) end if if(dc8_var(im) .ge. contourlvls(8)) then gsres@gsMarkerColor = colors(9) end if gsn_polymarker(wks,plot,dc8_lon(im),dc8_lat(im),gsres) print(" var "+im+" "+dc8_time(im)+" "+dc8_galt(im)+" "+dc8_lat(im) +" "+dc8_lon(im)+" "+dc8_var(im) +" "+gsres@gsMarkerColor ) end if end do txres@txFontColor = "black" xpt = 0.80 ypt = 0.148 gsn_text_ndc(wks, "dBZ", xpt, ypt, txres) print(" im1 lon/lat "+im1+" "+dc8_lon(im1)+" "+dc8_lat(im1)+" "+xpt+" "+ypt+" s_start="+s_start ) ; This next section keeps crashing on me (although it worked fine for my paper). It somehow lost the values of lat/lon. ; If you want to have end times of flight track printed on plot, then you will need to not skip this section. iskip = 1 if(iskip .eq. 0) then xpt = (dc8_lon(im1) - llonmin)/(llonmax-llonmin) - 0.08 ; locating the position of the text can be an iterative process ypt = (dc8_lat(im1) - latmin)/(latmax-latmin) - 0.04 print(" im1 lon/lat "+im1+" "+dc8_lon(im1)+" "+dc8_lat(im1)+" "+xpt+" "+ypt+" s_start="+s_start ) gsn_text_ndc(wks, s_start, xpt, ypt, txres) xpt = (dc8_lon(im2) - llonmin)/(llonmax-llonmin) - 0.05 ypt = (dc8_lat(im2) - latmin)/(latmax-latmin) - 0.181 print(" im2 lon/lat "+dc8_lon(im2)+" "+dc8_lat(im2)+" "+xpt+" "+ypt+" s_stop="+s_stop ) gsn_text_ndc(wks, s_stop, xpt, ypt, txres) end if ; Plot legend for flight track colors. clvls = 1.*contourlvls ltxt = flt2string(clvls) do i = 0,9 xpt = 360.-100.15 ypt = 39.05 + 0.235*(i) gsres@gsMarkerIndex = 4 ; 8=upward triangle; 4=open circle; 16=filled dots for markers. gsres@gsMarkerThicknessF = 4. gsres@gsMarkerColor = "black" gsn_polymarker(wks,plot,xpt,ypt,gsres) gsres@gsMarkerIndex = 16 ; 7=downward triangle; 4=open circle; 16=filled dots for markers. gsres@gsMarkerThicknessF = 2. gsres@gsMarkerColor = colors(i) gsn_polymarker(wks,plot,xpt,ypt,gsres) end do do i = 0,8 xpt = (-100.07 - llonmin)/(llonmax-llonmin) - 0.11 ypt = (0.234*i )*0.12000 + 0.325 print(" i lon/lat "+dc8_lon(im2)+" "+dc8_lat(im2)+" "+xpt+" "+ypt+" "+ltxt(i) ) gsn_text_ndc(wks, ltxt(i), xpt, ypt, txres) if(i.eq.8) then gsn_text_ndc(wks, "", xpt+0.04, ypt-0.003, txres) end if end do ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Close the Plot frame(wks) end do ; ifl loop ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; end