;*********************************************** ; radar_NWCT.ncl ; this plots the HI measurements. ; note: there are four levels in the radar reflectivity data in the file I analyzed, ; but only 2 (about 0.5 and 1.5 in the differential reflectivity_HI, etc ; The reflectivity variables without "HI" have more and consistent levels. ; Calculates kdp from differential phase ; plots data from *nc files created by the NOAA Weather & Climate Tool Kit ; There are four levels within each variable (different elevation angles). ; The first level is used (index 0, about 0.5 degrees). ; The next index is the azimuth (angle). ; The last index is the range or along radials (gates). ; the last two maps are based on data that has a simple deattenuation based on ; KDP ; ; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" begin calc_kdp = True plots_dbz = new ( 6, graphic ) plots_zdr = new ( 6, graphic ) plots_cc = new ( 6, graphic ) plots_dp = new ( 6, graphic ) plots_dp_ave = new ( 6, graphic ) plots_kdp = new ( 6, graphic ) plots_dbz_dat = new ( 6, graphic ) plots_zdr_dat = new ( 6, graphic ) ; dir = "./" do i_file = 0,5 ; do i_file = 0,0,1 filenames = ((/"KOKX_V06_20130209_000138.nc","KOKX_V06_20130209_003053.nc", \ "KOKX_V06_20130209_010213.nc", "KOKX_V06_20130209_013125.nc", \ "KOKX_V06_20130209_020046.nc", "KOKX_V06_20130209_023001.nc"/)) ; elev = 0.5 ; i_elev_a = 0 ; i_elev_b = 0 elev = 1.5 i_elev_a = 3 i_elev_b = 1 ; map Reflectivity filename = filenames(i_file) if (i_file.eq.0)then wks = gsn_open_wks("png",filename + "_" + elev) ; send graphics to PNG file end if print("filename = " + filename) ; varname = "Reflectivity_HI" f = addfile(filename,"r") elev := f->elevationR_HI printVarSummary(elev) ; print("elev = " + elev(:,:)) ; end do var3d := f->Reflectivity_HI(:,:,:) var := f->Reflectivity_HI(i_elev_a,:,:) printVarSummary(var3d) dimsvar = dimsizes(var) print("dimsvar = " + dimsvar) scale_factor = var@scale_factor_metadata add_offset = var@add_offset_metadata dbz := byte2flt( (/ var /) ) ;; strip data and don’t apply offset / scale yet. copy_VarAtts(var,dbz) dbz = where(dbz.lt.0, dbz+256,dbz) dbz = (dbz * scale_factor)+add_offset ;; now scale the data angles = f->azimuthR_HI(i_elev_a,:) dims_ang = dimsizes(angles) j_s = 0 j_e = dims_ang(0)-1 dsizes_dbz = dimsizes(dbz) i_s = 0 i_e = 4*dsizes_dbz(0)/5 dbz_sub = dbz(j_s:j_e,i_s:i_e) dbz_o = dbz ;---To see a list of all the variables on the file if (i_file .eq. 0)then print(getfilevarnames(f)) end if ; angles=90-angles angles_sub = angles(j_s:j_e) DEGTORAD = 0.017453292519943 xcenter = 0.0 ycenter = 0.0 ; Create 2D coordinate arrays. x=f->distanceR_HI x_sub = x(i_s:i_e) angles2d = conform(dbz_sub,angles_sub,0) x2d = conform(dbz_sub,x_sub,1)/1000 ; convert to km xarr = xcenter + x2d * cos(DEGTORAD * angles2d) yarr = ycenter + x2d * sin(DEGTORAD * angles2d) ; wks = gsn_open_wks("png",str_join((/filename,varname/),"_")) ; send graphics to PNG file dbres = True dbres@gsnDraw = False dbres@gsnFrame = False dbres@cnInfoLabelOn = False dbres@lbLabelBarOn = False dbres@gsnMaximize = True dbres@sfXArray = xarr dbres@sfYArray = yarr dbres@cnFillOn = True ; dbres@cnFillPalette = "NCV_jet" ; set color map cmap := read_colormap_file("NCV_rainbow2") printVarSummary(cmap) cmap := cmap(50:, :) dbres@cnFillPalette = cmap ; set color map dbres@cnLinesOn = False dbres@cnFillMode = "RasterFill" ; this mode is fastest dbres@trGridType = "TriangularMesh" dbres@tiMainString = filename dbres@tiXAxisString = "Zonal Distance (km)" dbres@tiYAxisString = "Meridional Distance (km)" dbres@lbOrientation = "Vertical" dbres@cnLevelSelectionMode = "ExplicitLevels" ; use explicit levels dbres@cnLevels = ispan(-9,57,3) ; set the contour levels ; dbres@cnLevels = ispan(-16,56,2) ; set the contour levels dbres@cnLineLabelsOn = False ; do not use line labels dbres@cnFillOn = True ; color fill dbres@cnLinesOn = False ; do not draw contour lines contour = gsn_csm_contour(wks,dbz_sub,dbres) plots_dbz(i_file) = contour ; will reuse old dbz variable for differential reflectivity delete(dbz) delete(var) var3dd = f->DifferentialReflectivity_HI(:,:,:) printVarSummary(var3dd) var = f->DifferentialReflectivity_HI(i_elev_b,:,:) ; var_other = f->DifferentialReflectivity(i_elev,:,:) elev_df := f->elevationD_HI ; elev_d := f->elevationD ; print(elev_df) ; print(elev_d) printVarSummary(var) ; printVarSummary(var_other) scale_factor = var@scale_factor_metadata add_offset = var@add_offset_metadata dbz = byte2flt( (/ var /) ) ;; strip data and don’t apply offset / scale yet. copy_VarAtts(var,dbz) dbz = where(dbz.lt.0, dbz+256,dbz) dbz = (dbz * scale_factor)+add_offset ;; now scale the data dbz_sub = dbz(j_s:j_e,i_s:i_e) zdr_o = dbz printVarSummary(zdr_o) angles2d = conform(dbz_sub,angles_sub,0) x2d = conform(dbz_sub,x_sub,1)/1000 ; convert to km ; units = "dB" dfres = True dfres@lbLabelBarOn = False dfres@gsnDraw = False dfres@gsnFrame = False dfres@cnInfoLabelOn = False dfres@gsnMaximize = True dfres@sfXArray = xarr dfres@sfYArray = yarr dfres@cnFillOn = True ; dfres@cnFillPalette = "NCV_jet" ; set color map ; cmap := read_colormap_file("NCV_rainbow2") dfres@cnFillPalette = "NCV_rainbow2" ; printVarSummary(cmap) ; cmap := cmap(75:, :) ; dbres@cnFillPalette := cmap ; set color map dfres@cnLinesOn = False dfres@cnFillMode = "RasterFill" ; this mode is fastest dfres@trGridType = "TriangularMesh" dfres@tiMainString = filename dfres@gsnRightString = "dB" dfres@tiXAxisString = "Zonal Distance (km)" dfres@tiYAxisString = "Meridional Distance (km)" dfres@lbOrientation = "Vertical" ; cnres@tiMainString = str_join((/filename," - ",varname/),"_") dfres@tiMainString = filename dfres@cnLevelSelectionMode = "ExplicitLevels" ; use explicit levels ; dfres@cnLevels = ispan(-2,6,1) ; set the contour levels ; dfres@cnLevels = ispan(-2,5,1) ; set the contour levels dfres@cnLevels = (/-2.0,-1.0,-0.75,-0.5,-0.25,0,0.25,0.5,0.75,1.0,1.25,1.50,1.75,2.0,\ 2.5,3.0,3.5,4.0,4.5,5.0,7.5,10.0/) dfres@cnFillColors = (/97,90,80,70,63,55,108,115,122,129,157,161,169,176,183,190,197,204,211, \ 218,225,232,239,247,252/) dfres@cnLineLabelsOn = False ; do not use line labels dfres@cnFillOn = True ; color fill dfres@cnLinesOn = False ; do not draw contour lines ; cnres@gsnSpreadColorStart= -1 contour = gsn_csm_contour(wks,dbz_sub,dfres) plots_zdr(i_file) = contour ; will reuse old dbz variable for correlation coefficient delete(dbz) delete(var) var_o = f->CorrelationCoefficient_HI printVarSummary(var_o) var = f->CorrelationCoefficient_HI(i_elev_b,:,:) elev_cc = f->elevationC_HI ; print(elev_cc) scale_factor = var@scale_factor_metadata add_offset = var@add_offset_metadata ; dbz = byte2flt( (/ f->Reflectivity_HI(i_elev,:,:) /) ) ;; strip data and don’t apply offset / scale yet. dbz = byte2flt( (/ var /) ) ;; strip data and don’t apply offset / scale yet. copy_VarAtts(var,dbz) dbz = where(dbz.lt.0, dbz+256,dbz) dbz = (dbz * scale_factor)+add_offset ;; now scale the data cc_o = dbz dbz_sub = dbz(j_s:j_e,i_s:i_e) angles2d = conform(dbz_sub,angles_sub,0) x2d = conform(dbz_sub,x_sub,1)/1000 ; convert to km ccres = True ccres@lbLabelBarOn = False ccres@gsnDraw = False ccres@gsnFrame = False ccres@cnInfoLabelOn = False ccres@gsnMaximize = True ccres@sfXArray = xarr ccres@sfYArray = yarr ccres@cnFillOn = True ccres@cnFillPalette = "MPL_gist_ncar" ; set color map ccres@cnLinesOn = False ccres@cnFillMode = "RasterFill" ; this mode is fastest ccres@trGridType = "TriangularMesh" ccres@tiMainString = filename ccres@tiXAxisString = "Zonal Distance (km)" ccres@tiYAxisString = "Meridional Distance (km)" ccres@lbOrientation = "Vertical" ccres@tiMainString = filename ccres@cnLevelSelectionMode = "ExplicitLevels" ; set explicit contour levels ; cmap = read_colormap_file("NCV_jet") ; cmap(0,:) = 1. ; Make first color white ; ccres@cnFillPalette = "NCV_jet" ccres@cnLevels = (/ 0.70,0.71,0.72,0.73,0.74,0.75,0.76,0.77,0.78,0.79,0.80,0.81,\ 0.82,0.83,0.84,0.85,0.86,0.87,0.88,\ 0.89,0.90,0.91,0.92,0.93,0.94,0.95,0.96,\ 0.97,0.98,0.99,1.0/) ccres@cnFillColors = (/ 13,17,21,25,29,33,37,41,\ 45,49,53,57,61,65,69,73,75,77,\ 81 ,83 ,85,\ 87,89,91,93,95,126,122,117,115,113,111,110/) ccres@cnLineLabelsOn = False ; do not use line labels ccres@cnFillOn = True ; color fill ccres@cnLinesOn = False ; do not draw contour lines ; cnres@gsnSpreadColorStart= -1 ; x = ispan (1,10,1) contour = gsn_csm_contour(wks,dbz_sub,ccres) plots_cc(i_file) = contour ; will reuse old dbz variable for Differential Phase delete(dbz) delete(var) var_oo = f->DifferentialPhase_HI printVarSummary(var_oo) var = f->DifferentialPhase_HI(i_elev_b,:,:) elev_dp = f->elevationD_HI ; print(elev_dp) scale_factor = var@scale_factor_metadata add_offset = var@add_offset_metadata ; dbz = byte2flt( (/ f->Reflectivity_HI(i_elev,:,:) /) ) ;; strip data and don’t apply offset / scale yet. dbz = byte2flt( (/ var /) ) ;; strip data and don’t apply offset / scale yet. copy_VarAtts(var,dbz) dbz = where(dbz.lt.0, dbz+256,dbz) dbz = (dbz * scale_factor)+add_offset ;; now scale the data wk_smooth121(dbz) dbz_sub = dbz(j_s:j_e,i_s:i_e) dfp_o = dbz angles2d = conform(dbz_sub,angles_sub,0) x2d = conform(dbz_sub,x_sub,1)/1000 ; convert to km dpres = True dpres@lbLabelBarOn = False dpres@gsnDraw = False dpres@gsnFrame = False dpres@cnInfoLabelOn = False dpres@gsnMaximize = True dpres@sfXArray = xarr dpres@sfYArray = yarr dpres@cnFillOn = True dpres@cnFillPalette = "NCV_jet" ; set color map dpres@cnLinesOn = False dpres@cnFillMode = "RasterFill" ; this mode is fastest dpres@trGridType = "TriangularMesh" dpres@tiMainString = filename dpres@tiXAxisString = "Zonal Distance (km)" dpres@tiYAxisString = "Meridional Distance (km)" dpres@lbOrientation = "Vertical" dpres@tiMainString = filename dpres@cnLevelSelectionMode = "ExplicitLevels" ; set explicit contour levels ; dfres@cnLevels = ispan(15,180,15) ; set the contour levels ; ccres@cnLevels = (/ 0.75,0.76,0.77,0.78,0.79,0.80,0.81,\ ; 0.82,0.83,0.84,0.85,0.86,0.87,0.88,\ ; 0.89,0.90,0.91,0.92,0.93,0.94,0.95,0.96,\ ; 0.97,0.98,0.99,1.0,1.01,1.02,1.03,1.04,1.05/) dpres@cnLevels = (/ 10,20,30,40,50,60,70,80,90,120,180/) ; ccres@cnFillColors = (/ 21,25,31,35,39,43,48,53,\ ; 61,81,85,89,95,99,104,108,115,123,\ ; 130,135,143,148,155,159,167,171,\ ; 180,184,188,192,196,200,205/) dpres@cnLineLabelsOn = False ; do not use line labels dpres@cnFillOn = True ; color fill dpres@cnLinesOn = False ; do not draw contour lines ; cnres@gsnSpreadColorStart= -1 ; x = ispan (1,10,1) contour = gsn_csm_contour(wks,dbz_sub,dpres) plots_dp(i_file) = contour ; filtered dp dfp_rave = dfp_o dfp_rave = 0. dims_o = dimsizes(dfp_o) do j = 0,dims_o(0)-1 dfp_rave(j,:) = runave(dfp_o(j,:),18,0) end do dbz_sub = dfp_rave(j_s:j_e,i_s:i_e) contour = gsn_csm_contour(wks,dbz_sub,dpres) plots_dp_ave(i_file) = contour ; var --> calculate deattenuation ;Bringi et al. (1990): ;Bringi, V.N., V. Chandrasekar, N. Balakrishnan, and D.S. Zrnić, 1990: An Examination of Propagation Effects in Rainfall on Radar Measurements at Microwave Frequencies. J. Atmos. Oceanic Technol., 7, 829–840, https://doi.org/10.1175/1520-0426(1990)007<0829:AEOPEI>2.0.CO;2 ; Calculate KDP (Specific Differential Phase) from Differential Phase. ; the grid spacing is 0.1 km dfp_dif = dfp_o dfp_dif = 0 ; need to relabel kdp with correct units, variable name ; not done here (KDP is specific differential attenuation and has ; units of db/km). ; move away from ground cluster (2 km) do j = 0,dims_o(0)-1 ; do j = 0,0,1 do i = 1,dims_o(1)-1 dfp_dif(j,i) = (dfp_o(j,i) - dfp_o(0,20)) end do end do ; approach from from Snyder et al (2010): http://journals.ametsoc.org/doi/full/10.1175/2010JTECHA1356.1 ; coefficients from Bringi et al, 1991: http://journals.ametsoc.org/doi/pdf/10.1175/1520-0426%281990%29007%3C0829%3AAEOPEI%3E2.0.CO%3B2 ; alpha = 0.313 and alpha_dp 0.0483 alpha = 0.016 gamma_dp = 0.23 ; ratio of alpha_dp/alpha_H ; units = "dB" ; print("calculated PIA") dbz_new = dbz_o zdr_new = zdr_o added = 0 do j = 0,dims_o(0)-1 do i = 1,dims_o(1)-1 check = ismissing(dbz_o(j,i)) ; if (check.eq.False)then if (dfp_dif(j,i)*alpha.gt.0)then dbz_new(j,i) = dbz_o(j,i) + alpha*dfp_dif(j,i) zdr_new(j,i) = zdr_o(j,i) + gamma_dp*alpha*dfp_dif(j,i) added = added + 1 else ; no physical meaning dbz_new(j,i) = dbz_o(j,i) zdr_new(j,i) = zdr_o(j,i) end if ; end if end do end do print("added grid numbers = " + added) dbz_sub = dbz_new(j_s:j_e,i_s:i_e) contour = gsn_csm_contour(wks,dbz_sub,dbres) plots_dbz_dat(i_file) = contour zdr_sub = zdr_new(j_s:j_e,i_s:i_e) dfres@gsnRightString = "dB" contour = gsn_csm_contour(wks,zdr_sub,dfres) plots_zdr_dat(i_file) = contour if (calc_kdp)then ; calculate KDP kdp = dfp_o kdp = 0. zdr_new = zdr_o ; define variable to hold points for determing slope of "line." print("calculating kdp") dfp_x = new(9,float) dfp_y = new(9,float) do j = 0,dims_o(0)-1 do i = 5,dims_o(1)-5 n_dfp = 0 check = ismissing(dfp_rave(j,i)) if (check.eq.False)then ; check = ismissing(cc_o(j,i)) ; if (check.eq.False)then ; if (cc_o(j,i).ge.0.9)then do ii = -4 + i,4+i dfp_x(n_dfp) = ii dfp_y(n_dfp) = dfp_rave(j,ii) n_dfp = n_dfp + 1 end do ; print("dfp_x = " + dfp_x) ; print("dfp_y = " + dfp_y) rc = regline(dfp_x,dfp_y) ; rc = avg(dfp_y) kdp(j,i) = rc ; end if ; end if end if end do end do kdp_sub = kdp(j_s:j_e,i_s:i_e) print("calculated kdp") printMinMax(kdp_sub,False) dkres = True ; units = "deg/km" dkres@gsnLeftString = "Specific Differential Phase" dkres@gsnRightString = "deg/km" dkres@lbLabelBarOn = False dkres@gsnDraw = False dkres@gsnFrame = False dkres@cnInfoLabelOn = False dkres@gsnMaximize = True dkres@sfXArray = xarr dkres@sfYArray = yarr dkres@cnFillOn = True ; dkres@cnFillPalette = "NCV_jet" ; set color map ; dkres@cnFillPalette = "MPL_gist_ncar" ; set color map dkres@cnFillPalette = "NCV_rainbow2" dkres@cnLinesOn = False dkres@cnFillMode = "RasterFill" ; this mode is fastest dkres@trGridType = "TriangularMesh" dkres@tiMainString = filename dkres@tiXAxisString = "Zonal Distance (km)" dkres@tiYAxisString = "Meridional Distance (km)" dkres@lbOrientation = "Vertical" dkres@tiMainString = filename dkres@cnLevelSelectionMode = "ExplicitLevels" ; set explicit contour levels ; dfres@cnLevels = ispan(15,180,15) ; set the contour levels ; ccres@cnLevels = (/ 0.75,0.76,0.77,0.78,0.79,0.80,0.81,\ ; 0.82,0.83,0.84,0.85,0.86,0.87,0.88,\ ; 0.89,0.90,0.91,0.92,0.93,0.94,0.95,0.96,\ ; 0.97,0.98,0.99,1.0,1.01,1.02,1.03,1.04,1.05/) ; dkres@cnLevels = (/ -2,-1,0,1,2,3,4,5,6,7/) dkres@cnLevels = (/-5,-2,-1.0,-0.5,-0.25,0,0.01,0.02,0.06,0.1,0.2,0.3,0.4,0.5,0.6,\ 0.7,0.8,0.9,1.0,2.0,3.0,5.0/) dkres@cnFillColors = (/97,90,80,70,63,55,108,115,122,129,157,161,169,176,183,190,197,204,211, \ 218,225,232,239,247,251,252,253/) ; ccres@cnFillColors = (/ 21,25,31,35,39,43,48,53,\ ; 61,81,85,89,95,99,104,108,115,123,\ ; 130,135,143,148,155,159,167,171,\ ; 180,184,188,192,196,200,205/) dkres@cnLineLabelsOn = False ; do not use line labels dkres@cnFillOn = True ; color fill dkres@cnLinesOn = False ; do not draw contour lines ; cnres@gsnSpreadColorStart= -1 ; x = ispan (1,10,1) contour = gsn_csm_contour(wks,kdp_sub,dkres) ; contour = gsn_csm_contour(wks,kdp_sub,dpres) plots_kdp(i_file) = contour end if end do ; Set Panel Resources pnlres = True ; pnlres@txString = t@description + " (" + t@units + ")" pnlres@txString = "Observed Reflectivity (dbZ)" pnlres@gsnPanelYWhiteSpacePercent = 13 ; Add white space b/w plots. pnlres@gsnPanelLabelBar = True ; Turn on common labelbar pnlres@lbLabelAutoStride = True ; Spacing of lbar labels. pnlres@lbBoxMinorExtentF = 0.13 pnlres@lbLabelFontHeightF = .015 pnlres@lbPerimOn = False pnlres@lbLabelFont = "Helvetica" ; label font pnlres@lbTitleOn = True ; turn on title pnlres@lbTitleString = "dbZ" ; title string ; pnlres@lbTitlePosition = "Top" ; title position pnlres@lbTitlePosition = "Bottom" ; title position pnlres@lbTitleFontHeightF= .015 ; make title smaller pnlres@lbTitleDirection = "Across" ; title direction pnlres@pmLabelBarWidthF = 0.5 pnlres@pmLabelBarHeightF = 0.15 pnlres@lbTopMarginF = 0.001 pnlres@pmLabelBarOrthogonalPosF = 0.02 ; position wrt plot pnlres@lbTitleOffsetF = -0.4 gsn_panel(wks,(/plots_dbz/),(/2,3/),pnlres) pnlres@txString = "Observed ZDR (dB)" pnlres@lbTitleString = "dB" ; title string gsn_panel(wks,(/plots_zdr/),(/2,3/),pnlres) pnlres@txString = "Observed Correlation Coeficient" pnlres@lbTitleString = " " ; title string gsn_panel(wks,(/plots_cc/),(/2,3/),pnlres) pnlres@txString = "Observed Differential Phase (deg)" pnlres@lbTitleString = "deg" ; title string gsn_panel(wks,(/plots_dp/),(/2,3/),pnlres) pnlres@txString = "Observed Corrected Reflectivity (dbZ)" pnlres@lbTitleString = "dbZ" ; title string gsn_panel(wks,(/plots_dbz_dat/),(/2,3/),pnlres) pnlres@txString = "Observed Corrected ZDR (dB)" pnlres@lbTitleString = "dB" ; title string gsn_panel(wks,(/plots_zdr_dat/),(/2,3/),pnlres) pnlres@txString = "Filtered Differential Phase (deg)" pnlres@lbTitleString = "deg" ; title string gsn_panel(wks,(/plots_dp_ave/),(/2,3/),pnlres) pnlres@txString = "Observed KDP (deg/km)" pnlres@lbTitleString = "deg/km" ; title string gsn_panel(wks,(/plots_kdp/),(/2,3/),pnlres) end