;*********************************************** ; radar_NWCT.ncl ; this plots the non-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 ; Note, need to find slope of kdp line +/- 4 from gate location. Not done yet. ; 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 ;1 KDP ;;;; NOTE, I SHOULD REALLY SET THE DISTANCE USED x_p as a function of angle (see gmail). ; ; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/cd_string.ncl" undef("set_res1_list") function set_res1_list() begin res = True res@gsnDraw = False res@gsnFrame = False res@cnInfoLabelOn = False res@cnFillOn = True res@cnLevelSelectionMode = "ExplicitLevels" ; use explicit levels res@cnLevels = ispan(-9,48,3) ; set the contour levels cmap := read_colormap_file("NCV_rainbow2") printVarSummary(cmap) cmap := cmap(50:, :) res@cnFillPalette = cmap ; set color map ; res@pmLabelBarOrthogonalPosF = -0.1 res@lbOrientation = "vertical" ; vertical label bar res@tiYAxisString = "Height (km)" res@tiXAxisString = "Time" ; opts_fft@lbLabelBarOn = False ; res@gsnDraw = False ; don't draw ; res@gsnFrame = False ; don't advance frame res@cnLinesOn = False ; this gives a warning (once) that it's not a valid setting, but it is. ; I think it happens because I reuse the opts for different variables. res@cnLineLabelsOn = False ; opts_fft@lbLabelBarOn = False res@cnInfoLabelOn = False return(res) end undef("set_res2_list") function set_res2_list() begin res = True ; res@lbLabelBarOn = False res@gsnDraw = False res@gsnFrame = False res@cnInfoLabelOn = False res@cnFillOn = True ; dfres@cnFillPalette = "NCV_jet" ; set color map ; cmap := read_colormap_file("NCV_rainbow2") res@cnLinesOn = False res@cnLevelSelectionMode = "ExplicitLevels" ; use explicit levels cmap = read_colormap_file("NCV_rainbow2") cmap(0,:) = 1. ; Make first color white res@cnFillPalette = cmap res@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/) res@cnFillColors = (/0,90,85,80,70,63,55,108,115,122,129,157,161,169,176,183,190,197,204/) res@cnLineLabelsOn = False ; do not use line labels res@cnFillOn = True ; color fill res@cnLinesOn = False ; do not draw contour lines res@lbOrientation = "vertical" ; vertical label bar res@tiYAxisString = "Height (km)" res@tiXAxisString = "Time" return(res) end undef("set_res3_list") function set_res3_list() begin res = True ; res@lbLabelBarOn = False res@gsnDraw = False res@gsnFrame = False res@cnInfoLabelOn = False res@cnFillOn = True ; dfres@cnFillPalette = "NCV_jet" ; set color map ; cmap := read_colormap_file("NCV_rainbow2") res@cnFillPalette = "MPL_gist_ncar" ; set color map ; printVarSummary(cmap) ; cmap := cmap(75:, :) ; dbres@cnFillPalette := cmap ; set color map res@cnLinesOn = False res@cnLevelSelectionMode = "ExplicitLevels" ; use explicit levels ; dfres@cnLevels = ispan(-2,6,1) ; set the contour levels res@cnLevelSelectionMode = "ExplicitLevels" ; set explicit contour levels res@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/) res@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/) res@cnLineLabelsOn = False ; do not use line labels res@cnFillOn = True ; color fill res@cnLinesOn = False ; do not draw contour lines res@lbOrientation = "vertical" ; vertical label bar res@tiYAxisString = "Height (km)" res@tiXAxisString = "Time" return(res) end undef("set_res4_list") function set_res4_list() begin res = True ; res@lbLabelBarOn = False res@gsnDraw = False res@gsnFrame = False res@cnInfoLabelOn = False res@cnFillOn = True ; dfres@cnFillPalette = "NCV_jet" ; set color map ; cmap := read_colormap_file("NCV_rainbow2") res@cnLinesOn = False res@cnLevelSelectionMode = "ExplicitLevels" ; use explicit levels res@cnFillPalette = "NCV_rainbow2" res@cnLevels = (/-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/) res@cnFillColors = (/0,90,80,70,63,55,108,115,122,129,157,161,169,176,183,190,197,204,211, \ 218,225,232,239,247/) res@cnLineLabelsOn = False ; do not use line labels res@cnFillOn = True ; color fill res@cnLinesOn = False ; do not draw contour lines res@lbOrientation = "vertical" ; vertical label bar res@tiYAxisString = "Height (km)" res@tiXAxisString = "Time" return(res) end begin ;main program plot_rad1 = new ( 4, graphic ) ; type = "png" ;;; ; type = "pdf" type = "png" type@wkWidth = 2500 ; produces a better looking PNG type@wkHeight = 2500 ; type = "pdf" ; type = "ps" ; type = "ncgm" wks = gsn_open_wks(type,"plot_qvp") ; dir = "./" filename = "KOKX_V06_20130209_000138.nc" f = addfile(filename,"r") elev := f->elevationR elev_o = elev ; printVarSummary(elev) n_elev = dimsizes(elev) n_elev2 = n_elev(0)+2 print("n_elev = " + n_elev) print("n_elev2 = " + n_elev2) var := f->Reflectivity(0,:,:) dims_dbz = dimsizes(var) filenames = ((/"qvp_output.201302090001.csv","qvp_output.201302090030.csv", \ "qvp_output.201302090102.csv", "qvp_output.201302090131.csv",\ "qvp_output.201302090200.csv","qvp_output.201302090230.csv","qvp_output.201302090259.csv",\ "qvp_output.201302090328.csv","qvp_output.201302090357.csv"/)) yr = str_get_cols(filenames, 11, 14) year = stringtoint(yr) print("year = " + year) mnth = str_get_cols(filenames, 15, 16) month = stringtoint(mnth) print("month = " + month) dy = str_get_cols(filenames, 17, 18) day = stringtoint(dy) print("day = " + day) hr = str_get_cols(filenames, 19, 20) hour = stringtoint(hr) print("hour = " + hour) mn = str_get_cols(filenames, 21, 22) minute = stringtoint(mn) print("minute = " + minute) n_files = dimsizes(filenames) num_obs = numAsciiRow(filenames(0)) print("num_obs = " + num_obs) num_col = 5 data = asciiread(filenames(0),(/num_obs,num_col/),"float") zo = data(:,0) num_rows = dimsizes(zo)-1 print("num_rows = " + num_rows) z = new((/n_files,num_rows/),float) dbz = new((/n_files,num_rows/),float) zdr = new((/n_files,num_rows/),float) cc = new((/n_files,num_rows/),float) kdp = new((/n_files,num_rows/),float) Times = new((/n_files/),string) ;; opts_rad = set_res1_list() opts_zdr = set_res2_list() opts_crs = set_res3_list() opts_kdp = set_res4_list() opts_rad@gsnLeftString = "Reflectivity" opts_rad@gsnRightString = "dBZ" opts_zdr@gsnLeftString = "Differential Reflectivity" opts_zdr@gsnRightString = "dB" opts_crs@gsnLeftString = "Correlation Coefficient" opts_crs@gsnRightString = " " opts_kdp@gsnLeftString = "Specific Differential Phase" opts_kdp@gsnRightString = "dB/km" do i_file = 0,n_files-1 ; do i_file = 0,0,1 new_time_units = "hours since 1800-01-01 00:00" new_time_array = cd_inv_calendar( year(i_file), month(i_file), day(i_file), \ hour(i_file), minute(i_file), 0, new_time_units,0) print("new_time_array = " + new_time_array) ; format = "" ;; defaults to "%H%M EST %d %c %Y" format = "" ;; defaults to "%H%M UTC %d %c %Y" format = "%H%M UTC %d %c %Y" ; new_time_array = new_time_array + 2 new_time_array = new_time_array + 0 time_var=cd_string(new_time_array,format) print("time_var = " + time_var) Times(i_file) = time_var num_obs = numAsciiRow(filenames(i_file)) data = asciiread(filenames(i_file),(/num_obs,num_col/),"float") print("filenames(i_file) = " + filenames(i_file)) ; printVarSummary(data) ; printVarSummary(z) ; see: https://www.ncl.ucar.edu/Document/Manuals/Getting_Started/basics.shtml#ArraySubscripting ; z(i_file,:) = data(::2,0) ; every other value is written ; print(data) ; zz = data(::2,0); every other value is written ; z = data(0:num_obs-2,0); every other value is written ; print(zz) ; printVarSummary(zz) z(i_file,:) = data(0:num_obs-2,0) dbz(i_file,:) = data(0:num_obs-2,1) print(dbz(i_file,:)) zdr(i_file,:) = data(0:num_obs-2,2) cc(i_file,:) = data(0:num_obs-2,3) kdp(i_file,:) = data(0:num_obs-2,3) ; print(kdp(i_file,:)) ; plot_rad1(0) = gsn_csm_contour(wks,fft_plane(0:zmax_pos,:),opts_fft) ; plot_rad2(0) = gsn_csm_contour(wks,ff1_plane(0:zmax_pos,:),opts_ff1) ; plot_rad3(0) = gsn_csm_contour(wks,fft_plane(0:zmax_pos,:),opts_rad) ; plot_rad4(0) = gsn_csm_contour(wks,qcr_plane(0:zmax_pos,:),opts_mass1) end do zmin = min(z)/1000. zmax = max(z)/1000. nz = 11 dim_z = dimsizes(z) print("dim_z = " + dim_z) zspan = dim_z(1)-1 opts_xy = True opts_xy@tmXTOn = False opts_xy@tmYROn = False opts_xy@tmXBMode = "Explicit" ; opts_xy@tmXBValues := fspan(0,xspan,nx) ; Create tick marks ; opts_xy@tmXBLabels := sprintf("%.1f",fspan(xmin,xmax,nx)) ; Create labels ; print("opts_xy@tmXBLabels" + opts_xy@tmXBLabels) ; opts_xy@tmXBLabelFontHeightF = 0.015 opts_xy@tmYLMode = "Explicit" opts_xy@tmYLValues = fspan(0,zspan,nz) ; Create tick marks opts_xy@tmYLLabels := sprintf("%.1f",fspan(zmin,zmax,nz)) ; Create labels print("opts_xy@tmYLLabels" + opts_xy@tmYLLabels) copy_VarAtts(opts_xy, opts_rad) copy_VarAtts(opts_xy, opts_zdr) copy_VarAtts(opts_xy, opts_crs) copy_VarAtts(opts_xy, opts_kdp) plot_rad1(0) = gsn_csm_contour(wks,dbz,opts_rad) plot_rad1(1) = gsn_csm_contour(wks,zdr,opts_zdr) plot_rad1(2) = gsn_csm_contour(wks,cc,opts_crs) plot_rad1(3) = gsn_csm_contour(wks,kdp,opts_kdp) ; 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@gsnPanelYWhiteSpacePercent = 2 ; Add white space b/w plots. ; pnlres@gsnPanelLabelBar = True ; Turn on common labelbar pnlres@gsnPanelLabelBar = False ; 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,(/plot_rad1/),(/2,2/),pnlres) end