; ncl WRF_pcp_tv.ncl 'plnm2="WRF_pcp_tv_"' 'filename="wrfout_d01_2018-12-19_0:00:00"' ; cp tvprecp.rgb $NCARG_ROOT/lib/ncarg/colormaps/ ; https://www.ncl.ucar.edu/Applications/Scripts/WRF_pcp_2.ncl ; This script plots Daily precipitation from WRF outputs. ; Written by: Ehsan Taghizadeh ; ; This is a modified version of another script. This script shows ; how to layer plots in order to get different map components drawn ; in the correct order. Here are the various layers: ; ; 1. Terrain plot with filled contours on a map ; 2. Rain total plot with filled and line contours ; 3. Map plot with shapefile outlines and markers/text at city locations ;---------------------------------------------------------------------- ; This script shows how to plot WRF rain totals on a WRF terrain map, ; using transparency for all rain totals below a certain level. ;---------------------------------------------------------------------- ; This script was originally from one contributed by Xiao-Ming Hu (xhu@ou.edu) ; Center for Analysis and Prediction of Storms ; University of Oklahoma ;---------------------------------------------------------------------- ; These files are loaded by default in NCL V6.2.0 and newer ; 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/wrf/WRFUserARW.ncl" ;---------------------------------------------------------------------- ; Create a simple map plot based on the map projection defined on the ; given WRF file. It can be used to lines, markers and text to it. ;---------------------------------------------------------------------- undef("create_wrf_map") function create_wrf_map(wks,wrf_file) local res begin res = True res = wrf_map_resources(wrf_file, res) res@gsnFrame = False res@gsnDraw = False res@pmTickMarkDisplayMode = "never" ; turn off tickmarks and their labels res@mpFillOn = False res@mpOutlineOn = False return(gsn_csm_map(wks,res)) end ;---------------------------------------------------------------------- ; Add shapefile outlines to the given plot. ;---------------------------------------------------------------------- undef("add_outlines_to_map") procedure add_outlines_to_map(wks,plot) local fnames, colors, lnres, n begin ; fnames = (/"/home/taghizade/IRIMO_NWP/shpfiles/shp_ostan/iran_province"/) + ".shp" fnames = "gadm36_IRN_shp/gadm36_IRN_1.shp" colors = (/"Brown"/) ; water, administrative lnres = True ; resources for polylines lnres@gsLineThicknessF = 2.0 ; 2x as thick ; Ehsan: 3.0 ;***Following loop is unnecessary do n=0,dimsizes(fnames)-1 dumstr = unique_string("primitive") lnres@gsLineColor = colors(n) plot@$dumstr$ = gsn_add_shapefile_polylines(wks, plot, fnames(n), lnres) end do end ;---------------------------------------------------------------------- ; Attach city text and markers to the given map plot. ;---------------------------------------------------------------------- undef("add_places_to_map") procedure add_places_to_map(wks,plot) local f, stids, stnames, sticaos, stlats, stlons, stelvs, stdims, xres, mkres begin ;---Read shapefile with "places" information. f = addfile("/home/taghizade/IRIMO_NWP/shpfiles/shp_stations/stations1.shp","r") stids = f->stid stnames = f->stnm sticaos = f->ICAOC stlons = f->stlon stlats = f->stlat stelvs = f->stelv stdims = dimsizes(stnames) ;---Set up text resources to label random areas of interest. txres = True txres@txFontHeightF = 0.010 txres@txJust = "topcenter" ; "topright" ;---Set up marker resources to mark random areas of interest. mkres = True mkres@gsMarkerIndex = 16 mkres@gsMarkerColor = "brown" mkres@gsMarkerSizeF = 7. icao_of_interest = (/ "OITR","OITT","OITL","OICS","OITZ","OIGG", \ "OICC","OIHH","OIIK","OICI","OICK","MAAA", \ "ALKK","MASA","OIII","OIAW","OIFS","OIFM", \ "OIIS","OING","OISY","OIBB","OISS","OIYY", \ "OIMB","OIMM","OIMN","OIKB","OIKK","OIZH" /) ;---Looking for interested icao codes icao_of_interest@_FillValue = "missing" do i=0,stdims-1 if(any(sticaos(i).eq.icao_of_interest)) then ii = ind(icao_of_interest.eq.sticaos(i)) icao_of_interest(ii) = "missing" dumstr = unique_string("primitive") ;---Attach the marker to the map. plot@$dumstr$ = gsn_add_polymarker(wks, plot, stlons(i), stlats(i), mkres) dumstr = unique_string("primitive") ; ;---Attach the text string to the map. plot@$dumstr$ = gsn_add_text(wks, plot, " " + sticaos(i), \ stlons(i), stlats(i), txres) delete(ii) end if end do end ;---------------------------------------------------------------------- ; Create a filled contour plot over a map of terrain. ;---------------------------------------------------------------------- undef("create_terrain_plot") function create_terrain_plot(wks,hgt,wrf_file) local res_ter begin ;---Set resources for terrain plot res_ter = True ; plot mods desired res_ter@gsnFrame = False res_ter@gsnDraw = False res_ter@cnFillOn = True ; color plot desired cmap_ter = read_colormap_file("cmocean_turbid") res_ter@cnFillPalette = cmap_ter(3:,:) res_ter@cnLinesOn = False ; turn off contour lines res_ter@cnLineLabelsOn = False ; turn off contour labels res_ter@cnFillMode = "RasterFill" res_ter@cnFillOpacityF = 1. ; 0.0 (completely transparent) to 1.0 (fully opaque) res_ter@cnMissingValFillColor = "cyan" res_ter@lbOrientation = "Vertical" res_ter@lbBoxEndCapStyle = "TriangleBothEnds" res_ter@pmLabelBarWidthF = 0.06 ; Make labelbar less thick res_ter@lbLabelFontHeightF = 0.014 res_ter@pmLabelBarOrthogonalPosF = 0.01 res_ter@gsnRightString = "" res_ter@gsnLeftString = "" ;--These three lines necessary for plotting in native WRF map projection res_ter = wrf_map_resources(wrf_file, res_ter) ; set map resources to match those on WRF file res_ter@tfDoNDCOverlay = "NDCViewport" ; necessary for correct overlay on map res_ter@gsnAddCyclic = False res_ter@mpOutlineBoundarySets = "AllBoundaries" res_ter@mpDataSetName = "Earth..4" ; Gives us provincial boundaries res_ter@mpNationalLineThicknessF = 2.5 res_ter@mpGeophysicalLineThicknessF = 2.5 res_ter@mpGeophysicalLineColor = "Black" res_ter@mpNationalLineColor = "Black" res_ter@mpFillOn = True res_ter@mpOceanFillColor = "SkyBlue" res_ter@mpInlandWaterFillColor = "SkyBlue" res_ter@mpLandFillColor = "Transparent" ; don't fill land res_ter@mpFillDrawOrder = "PostDraw" ; make sure these are drawn res_ter@pmTickMarkDisplayMode = "Always" ; turn on nicer tickmarks res_ter@tmXBLabelFontHeightF = 0.018 res_ter@tmYLLabelFontHeightF = 0.018 res_ter@tmYLLabelStride = 2 ; label every other tickmark res_ter@tmXBLabelStride = 2 ;---Point the tickmarks inward res_ter@tmYRMajorOutwardLengthF = 0 res_ter@tmYLMajorOutwardLengthF = 0 res_ter@tmXBMajorOutwardLengthF = 0 res_ter@tmXBMinorOutwardLengthF = 0 res_ter@tmXTOn = True res_ter@tmYROn = True res_ter@tmYRLabelsOn = False res_ter@tmXTLabelsOn = False return(gsn_csm_contour_map(wks,hgt,res_ter)) end ;---------------------------------------------------------------------- ; Set resources to be used for a filled contour plot of rain total. ;---------------------------------------------------------------------- undef("set_rain_total_resources") function set_rain_total_resources(wks) local res_tot begin res_tot = True res_tot@gsnFrame = False res_tot@gsnDraw = False cmap = (/(/255,255,255,255/),(/41,249,130,255/),(/72,240,254,255/),\ (/53,152,252,255/),(/19,129,190,255/),(/238,71,237,255/),\ (/228,65,106,255/)/) / 255. cmap(0,:) = (/0,0,0,0/) ; make first color fully transparent res_tot@cnFillOn = True res_tot@cnFillMode = "RasterFill" res_tot@cnFillPalette = cmap ; "precip2_17lev" ; res_tot@cnLineLabelsOn = False ; turn off contour labels res_tot@cnFillOpacityF = 1. ; .85 res_tot@cnInfoLabelOn = True res_tot@cnInfoLabelString = "Min= $ZMN$ Max= $ZMX$" res_tot@cnInfoLabelOrthogonalPosF = -0.05 ; move info label into plot res_tot@cnLevelSelectionMode = "ExplicitLevels" res_tot@cnLevels = (/.1,4,7,10,18,27/) res_tot@lbBoxEndCapStyle = "TriangleBothEnds" res_tot@pmLabelBarHeightF = 0.06 ; Make labelbar less thick res_tot@lbLabelFontHeightF = 0.014 res_tot@pmLabelBarOrthogonalPosF = 0.05 ; Move away from plot res_tot@gsnStringFont = 25 res_tot@gsnStringFontHeightF = 0.01 res_tot@gsnTickMarksOn = False ; Turn off tickmarks res_tot@tiMainString = "" res_tot@gsnRightString = "" return(res_tot) end ;---------------------------------------------------------------------- ; Given an array of plots, make sure they are all the same size and ; then draw them in the order they appear in the array. ;---------------------------------------------------------------------- procedure stack_plots_and_draw(wks,plots) local nplots begin nplots = dimsizes(plots) getvalues plots(0) "vpXF" : vpx "vpYF" : vpy "vpWidthF" : vpw "vpHeightF" : vph end getvalues draw(plots(0)) do np=1,nplots-1 setvalues plots(np) "vpXF" : vpx "vpYF" : vpy "vpWidthF" : vpw "vpHeightF" : vph end setvalues draw(plots(np)) end do frame(wks) end ;---------------------------------------------------------------------- begin ; idomain = 2 ; 3,3 ;---Rain data dir = "/home/taghizade/IRIMO_NWP/operational/wrfmaps/wrfouts/" print(filename) f0 = addfile(dir + filename,"r") ; f0 = addfile("wrfout_d0"+idomain+"_ncrcat_hourly.nc","r") Times = wrf_user_list_times(f0) ; Times = f0->Times ntimes = dimsizes(Times) years = str_get_cols(Times,0,3) months = str_get_cols(Times,5,6) days = str_get_cols(Times,8,9) dow = day_of_week(stringtointeger(years),stringtointeger(months),stringtointeger(days)) dowc = new(ntimes,string) pny = str_get_cols(Times(0),0,3) pnm = str_get_cols(Times(0),5,6) pnd = str_get_cols(Times(0),8,9) pnh = str_get_cols(Times(0),11,12) ;---Extracting day of week days_of_week = (/"Sun","Mon","Tue","Wed","Thu","Fri","Sat"/) dowc = days_of_week(dow) RAINC = f0->RAINC(0,:,:) ; wrfvar1(0,:,:) ; RAINC(0,:,:) RAINNC = f0->RAINNC(0,:,:) ; wrfvar2(0,:,:) ; RAINNC(0,:,:) RainTotal = RAINC(:,:) ;---Terrain data ; fmap = addfile("wrfinput_d0"+idomain+".nc","r") HGT = f0->HGT(0,:,:) ; wrf_user_getvar(a,"HGT",0) ; printMinMax(HGT,0) ; Terrain Height (m) : min=-84.4633 max=3712 ; HGT =(/HGT/1000./) xlat = f0->XLAT(0,:,:) xlong = f0->XLONG(0,:,:) printMinMax(xlong,0) printMinMax(xlat,0) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; Read/open NCL'ls crude (1x1) land-sea mask ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; ; a2 = addfile("$NCARG_ROOT/lib/ncarg/data/cdf/landsea.nc","r") ; lsdata = a2->LSMASK ; lsm = landsea_mask(lsdata, xlat, xlong) ; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; ; HGT2 := HGT ; replicate for plot ; HGT2@_FillValue = -9999 ; HGT2 = where(lsm.eq.0 , HGT2@_FillValue , HGT2) ; HGT2 = where(lsm.eq.2 , HGT2@_FillValue , HGT2) ; HGT2 = where(lsm.eq.4 , HGT2@_FillValue , HGT2) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; do ihour = 4, ntimes(0)-1,4 ; For domain 1 istart = ihour-4 ; For domain 1 iend = ihour wks = gsn_open_wks("png" ,plnm2+pny+pnm+pnd+pnh+"_"+ihour*6); For domain 1 RainTotal = (/(f0->RAINC(iend,:,:) + f0->RAINNC(iend,:,:) - \ (f0->RAINC(istart,:,:) + f0->RAINNC(istart,:,:)))/1. /) res_tot = set_rain_total_resources(wks) res_tot@gsnLeftString = "Terrain Height (" + HGT@units + ")~C~" + \ "Daily Precipitation "+"("+RainTotal@units+ \ ") from "+dowc(istart)+" "+Times(istart)+ \ " to "+dowc(iend)+" "+Times(iend)+"~C~" + \ "Run Time: "+dowc(0)+" "+Times(0)+" Valid Time: +"+ihour+"h"+ \ "~C~(WRF, 9km, IRIMO)~C~" ;---Create three plots terrain_plot = create_terrain_plot(wks,HGT,f0) raintot_plot = gsn_csm_contour(wks,RainTotal,res_tot) wrf_map_plot = create_wrf_map(wks,f0) ;---Add stuff to the WRF map plot add_outlines_to_map(wks,wrf_map_plot) add_places_to_map(wks,wrf_map_plot) ;---Draw the plots in one frame stack_plots_and_draw(wks,(/terrain_plot,raintot_plot,wrf_map_plot/)) end do ; ihour end