;---------------------------------------------------------------------- ; dataonmap_zoom_10.ncl ; ; Concepts illustrated: ; - Plotting WRF data on native grid ; - Plotting WRF data on non-native grid ; - Zooming in on a WRF map ;---------------------------------------------------------------------- ; 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/WRF_contributed.ncl" ;---------------------------------------------------------------------- ; This script was updated Dec 2018 to use wrf_user_ll_to_xy instead of ; wrf_user_ll_to_ij, which will be deprecated in NCL V6.6.0. The ; difference between the two routines is that wrf_user_ll_to_xy returns ; 0-based indexes, and wrf_user_ll_to_ij returns 1-based indexes. ;---------------------------------------------------------------------- begin ;---Read some data off a WRF output file. filename = "wrfout_d03_2019-02-07_12:00:00" dir = "./WRF/test/em_real/" f = addfile(dir + filename,"r") Times = wrf_user_list_times(f) ntimes = dimsizes(Times) ; print(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) pnyl = str_get_cols(Times(0),0,3) pny = str_get_cols(Times(0),2,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) t2 = f->T2 ; TEMP at 2 M (K) [Time | 157] x [south_north | 333] x [west_east | 585] tc2 = t2 - 273.15 copy_VarMeta(t2,tc2) tc2@units = "C" ; ts = 0 ; hgt = wrf_user_getvar(f,"HGT",ts) ; Terrain Height (m) ; ter = wrf_user_getvar(f,"ter",ts) ; Terrain Height (m) ; tc = wrf_user_getvar(f,"tc",ts) ; printVarSummary(tc) ; nl = 0 ; tc2 = tc(nl,:,:) ;---------------------------------------------------------------------- ; This section selects a lat/lon region of interest by specifying ; the lower left and upper right corners of a lat/lon box. ; Use the wrf_user_ll_to_xy function to get the index values ; that *approximately* match the lat/lon area of interest. ;---------------------------------------------------------------------- ;---Set the two lat/lon corners that we want to zoom in on. minlat = 50. maxlat = 51.2 minlon = -125.5 maxlon = -124.5 lats = (/ minlat, maxlat /) lons = (/ minlon, maxlon /) loc = wrf_user_ll_to_xy(f, lons, lats, True) ; loc(0,:) is west-east (x) ; loc(1,:) is south-north (y) ; Use "loc" values to set special resources for zooming in on map. ; wrf_user_ll_to_ij deprecated in NCL V6.6.0. ; ; Index values are returned as 1 to N, so we have to ; subtract 1 for 0 to N-1 indexing that NCL requires. ; loc = wrf_user_ll_to_ij(f, lons, lats, True) ; loc = loc-1 do it = 0, ntimes-1 uv10 = wrf_user_getvar(f,"uvmet10",it) ; wrfvar2 = uvmet10 --> u10,v10 met velocity (m/s) ; printVarSummary(uv10) uv10 = uv10*1.94386 ; Turn wind into knots uv10@units = "kt" u10 = uv10(0,:,:) v10 = uv10(1,:,:) u10@description = "u10 met velocity" v10@description = "v10 met velocity" ;printVarSummary(u10) fa3i = sprinti("%0.3i", it) initLabel = "~C~ Init : "+dowc(0)+" "+Times(0)+" (WRF, 1km, UNBC)" validLabel = " Valid: "+dowc(it)+" "+Times(it)'+initLabel sdwrh = pny + pnm + pnd + pnh + "/w10t2/" system("mkdir -p " + sdwrh) pn = "w10t2"+pnyl+pnm+pnd+pnh+"_WRFf"+fa3i+".png" wks = gsn_open_wks("png", sdwrh+pn) ;---Set resources common to both plots res = True res@gsnMaximize = True ; maximize size res@cnFillOn = True ; color plot desired res@cnFillPalette = "BlGrYeOrReVi200" res@cnLinesOn = False ; turn off contour lines res@cnLineLabelsOn = False ; turn off contour labels ; res@cnLevelSelectionMode = "ExplicitLevels" ; res@cnLevels = ispan(145,345,20) * 0.0001 res@gsnAddCyclic = False ; turn off longitude cyclic point ; res@mpDataBaseVersion = "HighRes" res@pmTickMarkDisplayMode = "Always" ; nicer tickmarks ;---Set common resources for all plots ; res = True res@gsnFrame = False res@gsnDraw = False res@gsnLeftString = "" res@gsnRightString = "" ;---Necessary for contours to be overlaid correctly on WRF projection res@tfDoNDCOverlay = True ; res@tfDoNDCOverlay = "NDCViewport" ; can use this in NCL V6.5.0 or later ;---------------------------------------------------------------------- ; This section shows how to zoom in on WRF data using the native ; projection. ; ; You must set five special resources that will be used by ; wrf_map_resources to calculate the approximate map corners ; need to match the lat/lon box requested. ;---------------------------------------------------------------------- res1 = res ; Make copy of common resource list ; res1@ZoomIn = True ; These five resources are required ; res1@Xstart = loc(0,0) ; when zooming in on WRF data and ; res1@Xend = loc(0,1) ; keeping the same map projection. ; res1@Ystart = loc(1,0) ; res1@Yend = loc(1,1) ; ; Call wrf_map_resources to set various map resources based on the ; five resource values set above, and based on global attributes on ; the file pointed to by "f". ; res1 = wrf_map_resources(f, res1) res1@tfDoNDCOverlay = True ; Tells NCL this is a native projection ; res1@tfDoNDCOverlay = "NDCViewport" ; NCL V6.5.0 or later ;---Overwrite some of the resources set by wrf_map_resources. ; res1@mpUSStateLineColor = "black" ; res1@mpNationalLineColor = "black" ; res1@mpGeophysicalLineColor = "black" ; res1@mpUSStateLineThicknessF = 2.0 ; res1@mpNationalLineThicknessF = 2.0 ; res1@mpGeophysicalLineThicknessF = 2.0 ; ; Set a title and create a native projection WRF plot. ; ; Note that you MUST subset the data here, using the same i,j ; index values you calculated earlier! ; ; res1@tiMainString = "Surface Temperature and 10 Wind Speed" ; plot = gsn_csm_contour_map(wks,tc2(loc(1,0):loc(1,1),loc(0,0):loc(0,1)),res1) plot = gsn_csm_contour(wks,tc2(it,loc(1,0):loc(1,1),loc(0,0):loc(0,1)),res1) ;---Wind vector plot vec_res = res vec_res@vcMinDistanceF = 0.020 vec_res@vcRefLengthF = 0.038 vec_res@vcMinFracLengthF = 0.2 vec_res@vcLineArrowThicknessF = 3.0 vec_res@vcGlyphStyle = "CurlyVector" vec_res@vcRefAnnoOn = False vector = gsn_csm_vector(wks,u10(loc(1,0):loc(1,1),loc(0,0):loc(0,1)),v10(loc(1,0):loc(1,1),loc(0,0):loc(0,1)),vec_res) ;---Map plot map_res = True map_res@gsnFrame = False map_res@gsnDraw = False map_res@tiMainString = "TEMP at 2 M and 10 Wind Speed" map_res@gsnLeftString = tc2@description + " (" + tc2@units + ")" map_res@gsnRightString = "Wind (" + u10@units + ")" map_res@gsnStringFontHeightF = 0.015 map_res@ZoomIn = True ; These five resources are required map_res@Xstart = loc(0,0) ; when zooming in on WRF data and map_res@Xend = loc(0,1) ; keeping the same map projection. map_res@Ystart = loc(1,0) map_res@Yend = loc(1,1) map_res = wrf_map_resources(f,map_res) map_res@mpFillOn = False map_res@mpOutlineOn = True map_res@mpDataBaseVersion = "HighRes" map_res@mpOutlineDrawOrder = "PostDraw" ; map_res@mpUSStateLineColor = "Black" map_res@mpPerimLineColor = "Black" map_res@mpNationalLineColor = "Black" map_res@mpLimbLineColor = "Black" map_res@mpGridLineColor = "Black" map_res@mpGeophysicalLineColor = "Black" map_res@mpUSStateLineThicknessF = 3.0 map_res@mpNationalLineThicknessF = 3.0 map_res@mpGeophysicalLineThicknessF = 3.0 map = gsn_csm_map(wks,map_res) ;-Overlay plots on map and draw. overlay(map,plot) overlay(map,vector) draw(map) ; This will draw all overlaid plots and the map frame(wks) delete(uv10) delete(v10) delete(u10) delete(fa3i) delete(initLabel) delete(validLabel) delete(pn) delete(res) delete(res1) delete(vec_res) delete(map_res) delete(plot) delete(vector) delete(map) delete(wks) end do exit ;---------------------------------------------------------------------- ; This section shows how to zoom in on WRF data using the lat/lon ; arrays read off the file. In this section, we are NOT doing a ; native map projection, but instead are plotting the data on a ; simple lat/lon projection. ;---------------------------------------------------------------------- q2@lat2d = f->XLAT(0,:,:) ; Required for plotting in non-native q2@lon2d = f->XLONG(0,:,:) ; map projection. res2 = res ; Make copy of common resource list res2@mpMinLatF = minlat ; Area we want to zoom in on res2@mpMaxLatF = maxlat res2@mpMinLonF = minlon res2@mpMaxLonF = maxlon res2@pmTitleZone = 4 ; Moves title closer to plot ;---Set a title and create a non-native projection WRF plot. res2@tiMainString = "Zooming in on WRF data using 2D lat/lon arrays" plot = gsn_csm_contour_map(wks,q2(0,:,:),res2) end