; 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 ;---------------------------------------------------------------------- ; WRF_pcp_2.ncl ; ; Concepts illustrated: ; - Plotting WRF data ; - Overlaying WRF precipitation on terrain map using gsn_csm functions ; - Setting the correct WRF map projection using wrf_map_resources ; - Creating two contour plots with two sets of filled contours ; - Explicitly setting contour levels ; - Drawing fully transparent filled contours ; - Turning the tickmarks inward on the X and Y axes ; - Moving the contour informational label into the plot ; - Customizing the contour informational label ;---------------------------------------------------------------------- ; 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 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" ;---------------------------------------------------------------------- undef("add_outlines_to_map") procedure add_outlines_to_map(wks,plot) local fnames, colors, lnres, n begin ; ; These two shapefiles contain waterways and administrative ; geographic info. We will add them to the existing map using ; polylines. ; fnames = (/"/home/taghizade/IRIMO_NWP/shpfiles/shp_ostan/iran_province"/) + ".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 ;---------------------------------------------------------------------- ; Add markers for cities to 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 ;---------------------------------------------------------------------- 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 do it=0,ntimes-1 if (dow(it) .eq. 0) then dowc(it) = "Sun" else if (dow(it) .eq. 1) then dowc(it) = "Mon" else if (dow(it) .eq. 2) then dowc(it) = "Tue" else if (dow(it) .eq. 3) then dowc(it) = "Wed" else if (dow(it) .eq. 4) then dowc(it) = "Thu" else if (dow(it) .eq. 5) then dowc(it) = "Fri" else dowc(it) = "Sat" end if end if end if end if end if end if end do 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) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;---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") ; ("OceanLakeLandSnow") res_ter@cnFillPalette = cmap_ter(3:,:) ; cnFillPalette : https://www.ncl.ucar.edu/Document/Graphics/color_table_gallery.shtml ; res_ter@cnFillPalette = "gsltod" ; "MPL_BuGn" ; "topo_15lev" ; "OceanLakeLandSnow" ; "WhiteGreen" ; ; "MPL_Greens" ; "gsltod" ; Select grayscale colormap 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" ; "darkslategray1" ; "cyan" ; make _FillValue 'stand out' ; res_ter@lbLabelBarOn = True res_ter@lbOrientation = "Vertical" res_ter@lbBoxEndCapStyle = "TriangleBothEnds" res_ter@gsnRightString = "" res_ter@gsnStringFont = 25 res_ter@gsnStringFontHeightF = 0.015 res_ter@gsnLeftString = "Terrain Height "+"("+HGT@units+")" res_ter@gsnLeftStringOrthogonalPosF = .11 res_ter@pmLabelBarWidthF = 0.06 ; Make labelbar less thick res_ter@lbLabelFontHeightF = 0.014 res_ter@pmLabelBarOrthogonalPosF = 0.01 res_ter = wrf_map_resources(f0, res_ter) ; set map resources to match those on WRF file res_ter@tfDoNDCOverlay = "NDCViewport" ; necessary for correct overlay on map ; res_ter@tfDoNDCOverlay = True ; old method res_ter@mpOutlineBoundarySets = "AllBoundaries" res_ter@mpDataSetName = "Earth..4" ; Gives us provincial boundaries res_ter@mpNationalLineThicknessF = 2.5 res_ter@mpGeophysicalLineThicknessF = 2.5 ; thickness of map outlines ; Ehsan: Waters boundaries res_ter@mpGeophysicalLineColor = "Black" ; Ehsan: Waters boundaries res_ter@mpNationalLineColor = "Black" ; Ehsan: Countries boundaries ; res_ter@mpProvincialLineThicknessF = 2. ; Ehsan: Doesn't work for me ; res_ter@mpProvincialLineColor = "black" ; Ehsan: Doesn't work for me 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 ;---to fill the ocean in some solid color res_ter@mpFillOn = True res_ter@mpOceanFillColor = "SkyBlue" ; use whatever color you want res_ter@mpInlandWaterFillColor = "SkyBlue" ; use whatever color you want res_ter@mpLandFillColor = "Transparent" ; to make sure land doesn't get filled res_ter@cnFillDrawOrder = "Draw" ; Also can try "Predraw" res_ter@mpFillDrawOrder = "PostDraw" res_ter@mpOutlineDrawOrder = "PostDraw" ;---Set resources for rain total contour plot res_tot = True res_tot@gsnFrame = False res_tot@gsnDraw = False cmap := read_colormap_file("tvprecp.rgb") ; ("cyclic") ; ("wh-bl-gr-ye-re") ; ("NCV_jet") ; ("precip2_17lev") ; ("WhiteBlue") ; ("BlAqGrYeOrReVi200") 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@cnLinesOn = False ; turn off contour lines res_tot@cnLineLabelsOn = False ; turn off contour labels res_tot@cnLineDrawOrder = "PostDraw" res_tot@lbBoxEndCapStyle = "TriangleBothEnds" res_tot@cnFillOpacityF = 1. ; .85 res_tot@cnFillDrawOrder = "PostDraw" res_tot@tfDoNDCOverlay = "NDCViewport" ; necessary for correct overlay on map ; res_tot@tfDoNDCOverlay = True ; old method res_tot@cnLevelSelectionMode = "ExplicitLevels" res_tot@cnLevels = (/.1,4,7,10,18,27/) ; (/.1,1,2,5,10,15,20,30,40,50,60,80,100,120/) ; "mm" ;res_tot@cnLevelSelectionMode = "ManualLevels" ;res_tot@cnMaxLevelValF = 42 ;res_tot@cnMinLevelValF = 2 ;res_tot@cnLevelSpacingF = 4 ;res_tot@ContourParameters = (/2,42,4/) res_tot@pmLabelBarHeightF = 0.06 ; Make labelbar less thick res_tot@lbLabelFontHeightF = 0.014 res_tot@pmLabelBarOrthogonalPosF = -0.008 res_tot@cnInfoLabelOn = True res_tot@cnInfoLabelString = "Min= $ZMN$ Max= $ZMX$" res_tot@cnInfoLabelOrthogonalPosF = -0.104 ; move info label into plot ; res_tot@tiMainFont = "Helvetica-bold" ; res_tot@tiMainFontHeightF = 0.018 res_tot@gsnRightString = "" res_tot@gsnStringFont = 25 res_tot@gsnStringFontHeightF = 0.012 res_tot@gsnLeftStringOrthogonalPosF = 0.02 do ihour = 4, ntimes(0)-1,4 ; For domain 1 istart = ihour-4 ; For domain 1 ; do ihour = 24, ntimes(0)-1,24 ; For domain 2 ; istart = ihour-24 ; For domain 2 iend = ihour ;---Start the graphics wks = gsn_open_wks("png" ,plnm2+pny+pnm+pnd+pnh+"_"+ihour*6); For domain 1 ; wks = gsn_open_wks("png" ,plnm2+pny+pnm+pnd+pnh+"_"+ihour) ; For domain 2 ; wks = gsn_open_wks("png" ,plnm2+Times(0)+"_"+ihour) RainTotal = (/(f0->RAINC(iend,:,:) + f0->RAINNC(iend,:,:) - \ (f0->RAINC(istart,:,:) + f0->RAINNC(istart,:,:)))/1. /) res_tot@tiMainString = "";"Terrain Height "+"("+HGT@units+")" ;RAINC+RAINNC " +chartostring(f0->Times(iend,:)) ; res_tot@tiMainFont = 25 ; res_tot@tiMainFontHeightF = 0.015 res_tot@gsnLeftString = "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)" plot_terrain = gsn_csm_contour_map(wks,HGT,res_ter) add_outlines_to_map(wks,plot_terrain) ; OK add_places_to_map(wks,plot_terrain) ; OK plot_raintot = gsn_csm_contour(wks,RainTotal,res_tot) overlay(plot_terrain, plot_raintot) draw(plot_terrain) frame(wks) end do ; ihour end