load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" 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/cd_string.ncl" external AREA_PROB "./area_prob_value.so" begin ; define number of panel plots plot = new(16,graphic) a_asc = (/"wrf_03_13_2017_times.file"/) ; note plotted time is first, but I plot only last time step ;---Read the data in as a small 1D array values = asciiread(a_asc,-1,"float") nvals = dimsizes(values) print("nvals = " + nvals) ; nx=dims(1) ; ny=dims(0) ;---Calculate ntime based on the number of values we read in. n_l_times =121 print("n_times = " + n_l_times) ; ntxy = nblk+1 nvars = 5 nxy = 1 nblk = nxy * nvars ; create variables to hold the ascii file variables. ; read times asc_time_file = a_asc num_col = 5 ; num_times = n_l_times num_times = numAsciiRow(asc_time_file) year = new((/num_times/),integer) month = new((/num_times/),integer) day = new((/num_times/),integer) hour = new((/num_times/),integer) minute = new((/num_times/),integer) dims = dimsizes(year) print(dims) z1 = asciiread(asc_time_file,(/num_times,num_col/),"integer") year = z1(:,0) month = z1(:,1) day = z1(:,2) hour = z1(:,3) minute = z1(:,4) ; set bounding box ; bounding box for January 23 2016 lats = (/ 39.00, 43.50 /) lons = (/ -77.5, -69.50 /) ; bounding box for January 26 2015 ; lats = (/ 37.75, 44.00 /) ; lons = (/ -77.0, -68.00 /) ; bounding box for February 08 2013 ; lats = (/ 38.00, 42.50 /) ; lons = (/ -77.0, -68.75 /) ; set it_beg, it_inc ; 4 km it_beg = 0 it_inc = 12 ; my resolution ; my_date = "02_08_2013" ; my_date = "01_22_2016" my_date = "03_13_2017" ; my_date = "01_26_2016" ; my_date = "" ; outfile="wrf_snowfall_4.0km_" + my_date outfile="wrf_obs_snowfall_" + my_date ; outfile="wrf_snowfall_12.0km_" + my_date wtype = "png" ; wtype = "pdf" ; for png --> gives "error" for pdf plot wtype@wkWidth = 7500 ; produces a better looking PNG wtype@wkHeight = 7500 wks = gsn_open_wks(wtype,outfile) DATADir = "/data4/wrf/" + my_date + "/GEFS/" ; DATADir = "/barry3/barry/" + my_date + "/GEFS/" print("DATADir = " + DATADir) ; for when there is no date ; find directories with this date, etc cmd = "find " + DATADir + "/GEFS* -type d -print" dirWRF = systemfunc (cmd) print(dirWRF) ; this is the total number of directories within the listed directory header ; it can be just the number of, for example, "WRF_XX" directories or also all subdirectories within n_dirWRF = dimsizes(dirWRF) ; print("n_dirWR" + n_dirWRF) ; this is the total number of directories within the listed directory header ; it can be just the number of, for example, "WRF_XX" directories or also all subdirectories within n_dirWRF = dimsizes(dirWRF) ; ; this is the name of the file ntimes = 10 n_loops = 09 n_loopss = 10 it_old = 0 wrfname2 ="wrfout_d02_2017-03-13_00:00:00_2dom" wrfname3 ="wrfout_d03_2017-03-13_00:00:00_3dom" wrfname4 ="wrfout_d04_2017-03-13_00:00:00" do i_file = 0,2 print("dirWRF(0) = " + dirWRF(0)) if (i_file.eq.0)then filename = dirWRF(0)+"/"+wrfname2 resol = "12.0km" end if if (i_file.eq.1)then filename = dirWRF(0)+"/"+wrfname3 resol = "4.0km" end if if (i_file.eq.2)then filename = dirWRF(0)+"/"+wrfname4 resol = "1.3km" end if print("filename at top = " + filename) a = addfile(filename,"r") times = wrf_user_list_times(a) ; get times in the file ntimes = dimsizes(times) ; number of times in the file print ("ntimes = " + ntimes) rain_old := wrf_user_getvar(a,"SNOWH",it_old) it_val = ntimes-1 ; do it = it_beg,ntimes,it_inc do it = ntimes-1,ntimes-1,1 ; do it = it_beg,it_beg,it_inc print("it = " + it) print("it_old = " + it_old) i_loop = 0 do n=0,n_dirWRF -1,15 ; do n=0,171,19 ; do n=0,0,1 if (i_file.eq.0)then filename = dirWRF(n)+"/"+wrfname2 end if if (i_file.eq.1)then filename = dirWRF(n)+"/"+wrfname3 end if if (i_file.eq.2)then filename = dirWRF(n)+"/"+wrfname4 end if print("i_loop " + i_loop) print("dirWRF(n) = " + dirWRF(n)) print("filename = " + filename) print("here 4") a = addfile(filename,"r") ntimes = dimsizes(times) ; number of times in the file print ("ntimes = " + ntimes) print("here 5") ; print("it = " + it) ; print("it_old = " + it_old) rainnc := wrf_user_getvar(a,"SNOWH",it) if (i_file.eq.2)then rainnc4 = wrf_user_getvar(a,"SNOWH",it) printMinMax(rainnc4,False) end if ; rain_acc := rainnc - rain_old ; rain_acc := rain_acc*100./2.54 rain_acc := rainnc*100./2.54 print("snowh") printMinMax(rain_acc,False) print("rainnc") printMinMax(rainnc,False) print("rainold") printMinMax(rain_old,False) if (i_loop.eq.0)then liq_equiv_ave:=rain_acc/(n_loops+1) end if if (i_loop.gt.0)then liq_equiv_ave:=liq_equiv_ave+rain_acc/(n_loops+1) end if dims2d=dimsizes(rainnc) print("dims2d = " + dims2d) if (i_loop.eq.0)then prob_100:=rainnc prob_50:=rainnc prob_25:=rainnc prob_100 = 0 prob_50 = 0 prob_25 = 0 do j=1,dims2d(0)-1 do i=1,dims2d(1)-1 if (rain_acc(j,i).gt.24)then prob_100(j,i) = 100./n_loopss end if if (rain_acc(j,i).gt.18)then prob_50(j,i) = 100./n_loopss end if if (rain_acc(j,i).gt.12)then prob_25(j,i) = 100./n_loopss end if end do end do end if if (i_loop.gt.0)then do j=1,dims2d(0)-1 do i=1,dims2d(1)-1 if (rain_acc(j,i).gt.24)then prob_100(j,i) = prob_100(j,i) + 100./n_loopss end if if (rain_acc(j,i).gt.18)then prob_50(j,i) = prob_50(j,i) + 100./n_loopss end if if (rain_acc(j,i).gt.12)then prob_25(j,i) = prob_25(j,i) + 100./n_loopss end if end do end do end if i_loop = i_loop + 1 if (i_loop-1.eq.n_loops)then xlat2d:=wrf_user_getvar(a,"XLAT",it) xlon2d:=wrf_user_getvar(a,"XLONG",it) prob_100@lat2d = xlat2d ; this is for plotting prob_100@lon2d = xlon2d ; using gsn_csm routines prob_50@lat2d = xlat2d ; this is for plotting prob_50@lon2d = xlon2d ; using gsn_csm routines prob_25@lat2d = xlat2d ; this is for plotting prob_25@lon2d = xlon2d ; using gsn_csm routines liq_equiv_ave@lat2d = xlat2d ; this is for plotting liq_equiv_ave@lon2d = xlon2d ; using gsn_csm routines ; loc(0,;) is west-east (x) ; loc(1,:) is south-north (y) ; subtract one since we want to use it as an index in NCL loc = wrf_user_ll_to_ij(a, lons, lats, True) x_start = loc(0,0) - 1 x_end = loc(0,1) - 1 y_start = loc(1,0) - 1 y_end = loc(1,1) - 1 ;;;; ;;; create plot (pdf) ;;;; ; copy coordinate variables ; copy_VarCoords(rainnc,liq_equiv_ave) ; copy_VarCoords(rainnc,rain_tot) gsn_define_colormap(wks,"BkBlAqGrYeOrReViWh200") ; load color table ;;; Font Style res1 = True ; plots modification on res1@txFont = "Helvetica" res1@txFontQuality = "High" ;;; Title ; first calculate time new_time_units = "hours since 1800-01-01 00:00" new_time_units = new_time_units new_time_array = cd_inv_calendar( year(it_val), month(it_val), day(it_val), \ hour(it_val), minute(it_val), 0, new_time_units,0) it_val = it_val + it_inc print("new_time_array = " + new_time_array) ; printVarSummary(new_time_array) format = "" ;; defaults to "%H%M EST %d %c %Y" format = "%H%M UTC %d %c %Y" ; 2.16667 = + 2 hours ;+ 10 minutes (for first time on map in EST) ; new_time_array = new_time_array + 2.16667 ; new_time_array = new_time_array -5 new_time_array = new_time_array time_var=cd_string(new_time_array,format) print("time_var = " + time_var) res1@gsnDraw = False ; do not draw the plot res1@gsnFrame = False ; do not advance the frame ; res1@tiMainString = time_var res1@tiMainString = " " res1@gsnLeftString = "WRF " + resol + " Mean Snowfall" res1@gsnRightString = "in" ;;; Map res1 = True res1@mpDataBaseVersion = "Ncarg4_1" ; choose more recent database res1@mpDataSetName = "Earth..4" ; high resolution res1@mpDataBaseVersion = "MediumRes" ; res1@mpOutlineBoundarySets = "National" ; National borders ; res1@mpOutlineBoundarySets = "AllBoundaries" ; National borders res1@mpOutlineBoundarySets = "USStates" ; National borders res1@mpGeophysicalLineColor = "Black" res1@mpNationalLineColor = "Black" res1@mpUSStateLineColor = "Black" res1@mpNationalLineColor = "Black" res1@mpGridLineColor = 0 res1@mpLimbLineColor = "Black" res1@mpPerimLineColor = "Black" res1@mpPerimLineThicknessF = 25 res1@mpUSStateLineThicknessF = 12.0 res1@mpGeophysicalLineThicknessF = 12.0 res1@mpGridAndLimbOn = False res1@pmTickMarkDisplayMode = "Always" ;;; Prec res1@cnFillOn = True ; turn on color fill res1@cnLinesOn = False ; turn of contour lines ;;; color fill tresholds and colors res1@cnLevelSelectionMode = "ExplicitLevels" ; set explicit contour levels res1@cnLevels = (/ 0, 1,4., 8, \ 12,18,24,30,36,50/) res1@cnFillColors = (/"White","White","AntiqueWhite1",\ "AntiqueWhite3",\ "SpringGreen1","SpringGreen4",\ "Yellow","Orange","Red","HotPink3",\ "Violet"/) ;;; label bars style ;;; if not set standard values are used. ;;; so you normally can delete these lines without any problems res1@cnInfoLabelOn = False ; turn off contour info label res1@lbAutoManage = False ; control label bar res1@pmLabelBarDisplayMode = "Always" ; turns on label bar res1@lbOrientation = "Vertical" res1@pmLabelBarSide = "Right" res1@lbLabelAutoStride = True res1@pmLabelBarWidthF = 0.20 res1@pmLabelBarHeightF = 0.7 res1@lbLabelFontHeightF = .020 res1@lbPerimOn = False res1@lbLabelFont = "Helvetica" ; label font res1@lbTitleOn = True ; turn on title res1@lbTitleString = "in" ; title string res1@lbTitlePosition = "Bottom" ; title position res1@lbTitleFontHeightF= .030 ; make title smaller res1@lbTitleDirection = "Across" ; title direction ;;; Map Projection (min and max Lon and Lat) res1@gsnAddCyclic = True ; data already has cyclic point ; this must also be set for any zoom min_lat = 35.5 max_lat = 44 min_lon = -77 max_lon = -69 ; res1@mpProjection = "LambertConformal" res1@mpMinLatF = min_lat res1@mpMaxLatF = max_lat res1@mpMinLonF = min_lon res1@mpMaxLonF = max_lon res1@mpLimitMode="LatLon" ; to have a common label bar, both plots should be set to the same interval ; b/c the label bar is drawn by default from the interval of the first plot. ; res@cnLevelSelectionMode = "ManualLevels" ; res@cnMinLevelValF = -10. ; res@cnMaxLevelValF = 45. ; res@cnLevelSpacingF = 5. res2 = True res2@gsnDraw = False ; do not draw the plot res2@gsnFrame = False ; do not advance the frame res2@gsnLeftString = "WRF " + resol + "Prob > 12 in" res2@gsnRightString = "%" res2@cnFillOn = True ; turn on color fill res2@cnLinesOn = False ; turn of contour lines ;;; Map res2@mpDataBaseVersion = "Ncarg4_1" ; choose more recent database res2@mpDataSetName = "Earth..4" ; high resolution res2@mpDataBaseVersion = "MediumRes" ; res1@mpOutlineBoundarySets = "National" ; National borders ; res1@mpOutlineBoundarySets = "AllBoundaries" ; National borders res2@mpOutlineBoundarySets = "USStates" ; National borders res2@mpNationalLineColor = "Black" ; National borders color res2@mpGeophysicalLineColor = "Black" res2@mpNationalLineColor = "Black" res2@mpUSStateLineColor = "Black" res2@mpNationalLineColor = "Black" res2@mpGridLineColor = 0 res2@mpLimbLineColor = "Black" res2@mpPerimLineColor = "Black" res2@mpUSStateLineThicknessF = 3.0 res2@mpGeophysicalLineThicknessF = 1.5 res2@pmTickMarkDisplayMode = "Always" ;;; Prec res2@cnFillOn = True ; turn on color fill res2@cnLinesOn = False ; turn of contour lines ;;; color fill tresholds and colors res2@cnLevelSelectionMode = "ExplicitLevels" ; set explicit contour levels res2@cnLevels = (/ 1,10., \ 20, 30,40,50,60,70,80,90/) res2@cnFillColors = (/"White","AntiqueWhite","AntiqueWhite3", \ "chartreuse", \ "chartreuse3","ForestGreen", \ "Yellow","Orange","Red","HotPink3","Violet"/) ;;; label bars style ;;; if not set standard values are used. ;;; so you normally can delete these lines without any problems res2@cnInfoLabelOn = False ; turn off contour info label res2@lbAutoManage = False ; control label bar res2@pmLabelBarDisplayMode = "Always" ; turns on label bar res2@lbOrientation = "Vertical" res2@pmLabelBarSide = "Right" res2@lbLabelAutoStride = True res2@pmLabelBarWidthF = 0.20 res2@pmLabelBarHeightF = 0.7 res2@lbLabelFontHeightF = .020 res2@lbPerimOn = False res2@lbLabelFont = "Helvetica" ; label font res2@lbTitleOn = True ; turn on title res2@lbTitleString = "(%)" ; title string res2@lbTitlePosition = "Bottom" ; title position res2@lbTitleFontHeightF= .030 ; make title smaller res2@lbTitleDirection = "Across" ; title direction ;;; Map Projection (min and max Lon and Lat) res2@gsnAddCyclic = True ; data already has cyclic point ; this must also be set for any zoom ; res2@mpProjection = "LambertConformal" res2@mpMinLatF = min_lat res2@mpMaxLatF = max_lat res2@mpMinLonF = min_lon res2@mpMaxLonF = max_lon res2@mpLimitMode="LatLon" ; res3 = res2 ; res3@gsnLeftString = "GEFS Prob > 18 in" res3@gsnLeftString = "WRF " + resol + "Prob > 18 in" res3@lbTitleString = "(%)" ; title string res4 = res2 ; res4@gsnLeftString = "GEFS Prob > 1.00 in" res4@gsnLeftString = "WRF " + resol + "Prob > 24 in" res4@lbTitleString = "(%)" ; title string ; plot(0) = gsn_csm_contour_map(wks,u,res) ; plot(1) = gsn_csm_contour_map(wks,v,res) ;************************************************ ; create panel ;************************************************ resP = True ; modify the panel plot ; resP@gsnPanelMainString = "A plot with a common label bar" ; use this for NCL V6.3.0 and earlier ; resP@txString = "A plot with a common label bar" ; resP@gsnPanelLabelBar = True ; add common colorbar ; resP@lbLabelFontHeightF = 0.007 ; make labels smaller end if end do ; end of directory loop end do ; end of file read plot(i_file) = gsn_csm_contour_map(wks,liq_equiv_ave,res1) ; create plot plot(i_file+4) = gsn_csm_contour_map(wks,prob_25,res2) ; create plot plot(i_file+8) = gsn_csm_contour_map(wks,prob_50,res3) ; create plot plot(i_file+12) = gsn_csm_contour_map(wks,prob_100,res4) ; create plot ; end of i_file end do ; ; Example of using panels with WRF data ; dx =1.3 ; dx =12 dx =4 filename = dirWRF(0)+"/"+wrfname3 a = addfile(filename,"r") xlat2d:=wrf_user_getvar(a,"XLAT",it) xlon2d:=wrf_user_getvar(a,"XLONG",it) dims2d=dimsizes(xlat2d) xlat2d_m=xlat2d xlon2d_m=xlon2d xland:=wrf_user_getvar(a,"XLAND",0) do itime = 0,3 ; t = wrf_user_getvar(a,"T2",itime) ; T2 in Kelvin ; slp = wrf_user_getvar(a,"slp",itime) ; slp asc_file = "../OBS/March13-16Snowfall_nolbl.csv" if (itime.eq.1)then date = "> 12 in" interval = 12 end if if (itime.eq.2)then date = "> 18 in" interval = 18 end if if (itime.eq.3)then date = "> 24 in" interval = 24 end if ; if (itime.ne.3)then num_col = 4 num_obs = numAsciiRow(asc_file) ; end if ; if (itime.eq.3)then ; num_col = 4 ; num_obs = numAsciiRow(asc_file) ; end if z1 := asciiread(asc_file,(/num_obs,num_col/),"float") lat_o := z1(:,0) lon_o := z1(:,1) obs_o := z1(:,3) ; you needs these units so you don't get plot warnings obs_o@long_name = "observations of some sort" obs_o@units = "not sure" lat_o@long_name = "observed latitudes" lat_o@units = "degrees_north" lon_o@long_name = "observed longitudes" lon_o@units = "degrees_east" ;---Get slice of lat/lon data at halfway points. nlat = dims2d(0) nlon = dims2d(1) xlat_m = xlat2d_m(:,nlon/2) xlon_m = xlon2d_m(nlat/2,:) rscan = (/0.8,0.4,0.20/) grid = obj_anal_ic(lon_o,lat_o,obs_o,xlon_m,xlat_m,rscan,False) grid@lat2d = xlat2d ; this is for plotting grid@lon2d = xlon2d ; this is for plotting printVarSummary(grid) f_prob = grid var2d = grid if (itime.ne.0)then do j=0,dims2d(0)-1 do i=0,dims2d(1)-1 check = ismissing(grid(j,i)) if (check.eq.False)then if (grid(j,i).gt.interval)then var2d(j,i) = 1 else var2d(j,i) = 0 end if end if end do end do AREA_PROB::area_prob_value(var2d,dx,f_prob,dims2d(0),dims2d(1)) do j=0,dims2d(0)-1 do i=0,dims2d(1)-1 if (xland(j,i).eq.2)then f_prob(j,i) = 0 end if end do end do else do j=0,dims2d(0)-1 do i=0,dims2d(1)-1 check = ismissing(grid(j,i)) if (check.eq.False)then if(xland(j,i).eq.1)then f_prob(j,i)= grid(j,i) end if if(xland(j,i).eq.2)then f_prob(j,i)= 0. end if end if end do end do end if if (itime.eq.0)then res5=res1 res5@gsnLeftString = "Accumulated Snow" res5@lbTitleString = "in" ; title string plot(3) = gsn_csm_contour_map(wks,f_prob,res5) ; create plot end if if (itime.eq.1)then res6=res2 res6@gsnLeftString = "Neigh. Prob > 12 in" res6@lbTitleString = "(%)" ; title string plot(7) = gsn_csm_contour_map(wks,f_prob,res6) ; create plot end if if (itime.eq.2)then res7=res3 res7@gsnLeftString = "Neigh. Prob > 18 in" res7@lbTitleString = "(%)" ; title string plot(11) = gsn_csm_contour_map(wks,f_prob,res7) ; create plot end if if (itime.eq.3)then res8=res4 res8@gsnLeftString = "Neigh. Prob > 24 in" res8@lbTitleString = "(%)" ; title string plot(15) = gsn_csm_contour_map(wks,f_prob,res8) ; create plot end if end do gsn_panel(wks,plot,(/4,4/),resP) ; now draw as one plot ; end