; https://www.ncl.ucar.edu/Training/Tutorials/WRF_Users_Workshop/WRF_NCL.shtml ; https://www.ncl.ucar.edu/Training/Tutorials/WRF_Users_Workshop/Scripts/wrf_demo_plot_tc_gsn_shapefiles.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 ;---Open WRF output file filename = "wrfout_d02_2018-08-28_00:00:00" dir = "/home/taghizade/IRIMO_NWP/operational/wrfmaps/wrfouts/" f = addfile(dir + filename,"r") ; tc = wrf_user_getvar(f,"tc",0) ; Temperature (C) rh2 = wrf_user_getvar(f,"rh2",0) uv10 = wrf_user_getvar(f,"uvmet10",0) ; u at 10 m, mass point uv10 = uv10*1.94386 ; Turn wind into knots uv10@units = "kts" printVarSummary(uv10) printVarSummary(rh2) ; printVarSummary(tc) wks = gsn_open_wks("png","wrf_demo_plot_tc_uv_rh_gsn_shapefiles") com_res = True com_res@gsnMaximize = True com_res@gsnDraw = False com_res@gsnFrame = False com_res@gsnRightString = "" com_res@gsnLeftString = "" com_res@tfDoNDCOverlay = True ; These two resources are required when overlaying com_res@gsnAddCyclic = False ; plots in native WRF map projection. ; tc_res = com_res rh_res = com_res uv_res = com_res ;---Temperature plot, this will be the base map plot ; tc_res@cnFillOn = True ; tc_res@cnFillPalette = "WhiteBlueGreenYellowRed" ; tc_res@cnLinesOn = False ; tc_res@cnFillOpacityF = 0.75 ; make partially transparent ; tc_res@cnLevelSpacingF = 0.5 ;; tc_res@mpOutlineOn = False ; Turn off NCL's map outlines ; tc_res@tiMainString = "RH2 masked by gsn_contour_shade" ; tc_res@lbOrientation = "Vertical" ; ; tc_res@gsnStringFontHeightF = 0.01 ; tc_res@gsnLeftString = tc@description + " (" + tc@units + " )~C~" + \ ; rh2@description + " (" + rh2@units + " )~C~" + \ ; uv10@description ; ; tc_res = wrf_map_resources(f,tc_res) ; Required for setting native WRF map projection ; ;;---Fix the settings that wrf_map_resources used for map outline color and thicknesses ; tc_res@mpOutlineBoundarySets = "AllBoundaries" ; gives us more outlines ; tc_res@mpDataSetName = "Earth..2" ; better resolution, use "Earth..2" if don't need the resolution ; tc_res@mpGeophysicalLineColor = "white" ; tc_res@mpNationalLineColor = "white" ; tc_res@mpUSStateLineColor = "white" ; tc_res@mpCountyLineColor = "white" ; tc_res@mpGeophysicalLineThicknessF = 5.0 ; tc_res@mpNationalLineThicknessF = 5.0 ; tc_res@mpUSStateLineThicknessF = 5.0 ; tc_res@mpCountyLineThicknessF = 5.0 ; nl = 3 ; level index ; tc_plot = gsn_csm_contour_map(wks,tc(nl,:,:),tc_res) ;---RH line contour plot rh_res@cnLineColor = "magenta" rh_res@cnLineThicknessF = 3.0 ; rh_res@cnLineDashPattern = 3 ; rh_res@cnLineLabelsOn = False ; rh_res@lbLabelBarOn = False rh_res@cnLineLabelFontHeightF = 0.008 ; rh_res@cnInfoLabelFontHeightF = 0.005 rh_res@mpOutlineOn = False ; Turn off NCL's map outlines rh_res@tiMainString = "RH2 masked by gsn_contour_shade" rh_res@gsnStringFontHeightF = 0.01 rh_res@gsnLeftString = rh2@description + " (" + rh2@units + ")~C~" + \ uv10@description + " (" + uv10@units + ")~C~" rh_res = wrf_map_resources(f,rh_res) ; Required for setting native WRF map projection ;---Fix the settings that wrf_map_resources used for map outline color and thicknesses rh_res@mpOutlineBoundarySets = "AllBoundaries" ; gives us more outlines rh_res@mpDataSetName = "Earth..2" ; better resolution, use "Earth..2" if don't need the resolution rh_res@mpGeophysicalLineColor = "white" rh_res@mpNationalLineColor = "white" rh_res@mpUSStateLineColor = "white" rh_res@mpCountyLineColor = "white" rh_res@mpGeophysicalLineThicknessF = 5.0 rh_res@mpNationalLineThicknessF = 5.0 rh_res@mpUSStateLineThicknessF = 5.0 rh_res@mpCountyLineThicknessF = 5.0 ;---Set resources for information label; CMN, CMX, CIU will be automatically filled in by NCL rh_res@cnInfoLabelOn = False rh_res@cnInfoLabelZone = 3 rh_res@cnInfoLabelString = "TopRight" rh_res@cnInfoLabelString = rh2@description + " contours from $CMN$ TO $CMX$ BY $CIU$" rh_res@cnInfoLabelFontHeightF = 0.01 ;---Set all values of rh2 <= 80 to missing. ; rh2@_FillValue = -999 ; rh2 = where(rh2.gt.80,rh2,rh2@_FillValue) rh_plot = gsn_csm_contour_map(wks,rh2,rh_res) opt = True opt@gsnShadeHigh = "magenta" rh_plot = gsn_contour_shade(rh_plot,-999.,80.,opt) ;---Vector plot uv_res@vcMinDistanceF = 0.02 uv_res@vcRefLengthF = 0.02 uv_res@vcRefMagnitudeF = 30. uv_res@vcMinFracLengthF = 0.2 uv_res@vcGlyphStyle = "WindBarb" uv_res@vcLineArrowThicknessF = 3.5 uv_res@vcLineArrowColor = "gray50" uv_res@vcRefAnnoZone = 3 uv_res@vcRefAnnoParallelPosF = 0.0 uv_res@vcRefAnnoJust = "TopLeft" uv_res@vcRefAnnoFontHeightF = 0.01 uv_plot = gsn_csm_vector(wks,uv10(0,:,:),uv10(1,:,:),uv_res) ;--- Add shapefile outlines to the base map plot add_outlines_to_map(wks,rh_plot) add_places_to_map(wks,rh_plot) ; lnres = True ; lnres@gsLineColor = "white" ; lnres@gsLineThicknessF = 3.0 ; default is a little thin ; usa_id = gsn_add_shapefile_polylines(wks,tc_plot,"USA_adm/USA_adm1.shp",lnres) ; mex_id = gsn_add_shapefile_polylines(wks,tc_plot,"MEX_adm/MEX_adm1.shp",lnres) ; cub_id = gsn_add_shapefile_polylines(wks,tc_plot,"CUB_adm_shp/CUB_adm2.shp",lnres) ;---Overlay RH and vector plots on base map plot ; overlay(tc_plot,rh_plot) overlay(rh_plot,uv_plot) draw(rh_plot) ; This draws everything and frame(wks) ; advances the frame end