;---------------------------------------------------------------------- ; wrf_gsn_10.ncl ;---------------------------------------------------------------------- ; Concepts illustrated: ; - Subsetting a WRF grid three different ways to take an average of points in a lat/lon box ; - Using opacity to emphasize or subdue overlain features ; - Setting the correct WRF map projection using wrf_map_resources ; - Illustrating the use of "where" to subset data ; - Illustrating the use of "getind_latlon2d" to subset data ; - Illustrating the use of "wrf_user_ll_to_xy" to subset data ;---------------------------------------------------------------------- ; 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. ;---------------------------------------------------------------------- ; 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" begin ;---Open WRF output file and read off the lat/lon grid a = addfile("wrfout_sample.nc","r") nt = 0 ; 0 = first time step t2 = wrf_user_getvar(a, "T2", nt) lat2d = wrf_user_getvar(a,"XLAT", nt) lon2d = wrf_user_getvar(a,"XLONG",nt) ;---Lat/lon box of interest where we want to average the data min_lat = 30 max_lat = 40 min_lon = -80 max_lon = -75 ; ; The purpose of the code below is to show two ways of subsetting ; a WRF lat/lon grid: one using the "where" statement, and one ; using "wrf_user_ll_to xy". The important thing to note is that ; you get very different results using these two methods. ; ;---Use "where" statement to get subset of lat/lon range t2@_FillValue = default_fillvalue(typeof(t2)) lon2d@_FillValue = default_fillvalue(typeof(lon2d)) lat2d@_FillValue = default_fillvalue(typeof(lat2d)) lat_sub1 = where(lat2d.lt.min_lat.or.lat2d.gt.max_lat .or. \ lon2d.lt.min_lon.or.lon2d.gt.max_lon, \ lat2d@_FillValue, lat2d) lon_sub1 = where(lat2d.lt.min_lat.or.lat2d.gt.max_lat .or. \ lon2d.lt.min_lon.or.lon2d.gt.max_lon, \ lon2d@_FillValue, lon2d) t2_sub1 = where(lat2d.lt.min_lat.or.lat2d.gt.max_lat .or. \ lon2d.lt.min_lon.or.lon2d.gt.max_lon, \ t2@_FillValue, t2) ;---Use "wrf_user_ll_to_xy" to get subset of lat/lon range lats = (/ min_lat, max_lat /) lons = (/ min_lon, max_lon /) loc = wrf_user_ll_to_xy(a, lons, lats, True) ; Added in NCL V6.6.0 ; Pre NCL V6.6.0, use wrf_user_ll_to_ij. You must subtract one from the ; indexes since this function returns 1-based indexes. ; loc = wrf_user_ll_to_ij(a, lons, lats, True) ; loc = loc-1 ; wrf function returns 1 to n, NCL needs 0 to n-1 ilt1 = loc(1,0) ilt2 = loc(1,1) iln1 = loc(0,0) iln2 = loc(0,1) lat_sub2 = lat2d(ilt1:ilt2 , iln1:iln2) lon_sub2 = lon2d(ilt1:ilt2 , iln1:iln2) t2_sub2 = t2 ; make a copy t2_sub2 = t2_sub2@_FillValue ; set to all missing t2_sub2(ilt1:ilt2,iln1:iln2) = t2(ilt1:ilt2, iln1:iln2) ; fill in subset ; ; Use "getind_latlon2d"; should produce same results as wrf_user_ll_to_xy. ; ; This function can be used on any curvilinear data, whereas ; wrf_user_ll_to_xy is tailored for WRF data. ; nm = getind_latlon2d (lat2d,lon2d, lats, lons) ilt1 = nm(0,0) ilt2 = nm(1,0) iln1 = nm(0,1) iln2 = nm(1,1) lat_sub3 = lat2d(ilt1:ilt2 , iln1:iln2) lon_sub3 = lon2d(ilt1:ilt2 , iln1:iln2) t2_sub3 = t2 ; make a copy t2_sub3 = t2_sub3@_FillValue ; set to all missing t2_sub3(ilt1:ilt2,iln1:iln2) = t2(ilt1:ilt2, iln1:iln2) ; fill in subset ;---Take an average of all the points, and of the two subsetted sets of points t2_avg = avg(t2) t2_sub1_avg = avg(t2_sub1) t2_sub2_avg = avg(t2_sub2) t2_sub3_avg = avg(t2_sub3) ;---Print this information ; printVarSummary(t2) ; printVarSummary(t2_sub1) ; printVarSummary(t2_sub2) ; printVarSummary(t2_sub3) print("--------------------------------------------------") print(" # non-msg points in t2 = " + num(.not.ismissing(t2))) print(" Average of t2 = " + t2_avg) print("--------------------------------------------------") print(" Using 'where' function") print(" # non-msg points in t2_sub1 = " + num(.not.ismissing(t2_sub1))) print(" Average of t2_sub1 = " + t2_sub1_avg) print("--------------------------------------------------") print(" Using 'wrf_user_ll_to_xy' function") print(" # non-msg points in t2_sub2 = " + num(.not.ismissing(t2_sub2))) print(" Average of t2_sub2 = " + t2_sub2_avg) print("--------------------------------------------------") print(" Using 'getind_latlon2d' function") print(" # non-msg points in t2_sub3 = " + num(.not.ismissing(t2_sub3))) print(" Average of t2_sub3 = " + t2_sub3_avg) ;---Start the graphics wtype = "png" if(wtype.eq."png") then wtype@wkWidth = 2000 wtype@wkHeight = 2000 end if wks = gsn_open_wks(wtype,"wrf_gsn") ;---Set some resources for a WRF map res = True res@gsnMaximize = True res@gsnDraw = False res@gsnFrame = False res@gsnLeftString = "min lat/lon = ("+min_lat+","+min_lon+")" res@gsnRightString = "max lat/lon = ("+max_lat+","+max_lon+")" res@mpGridAndLimbOn = True res@mpGridLineColor = "gray50" res@cnFillOn = True res@cnLinesOn = False res@cnLineLabelsOn = False res@cnFillOpacityF = 0.5 ; Fade the colors so we can see lat/lon box better. res@cnLevelSelectionMode = "ManualLevels" mnmxint = nice_mnmxintvl( min(t2), max(t2), 18, False) res@cnMinLevelValF = mnmxint(0) res@cnMaxLevelValF = mnmxint(1) res@cnLevelSpacingF = mnmxint(2) res@lbLabelBarOn = False ; will add in panel ;---Necessary for overlay to work properly in native WRF grid. res = wrf_map_resources(a,res) res@tfDoNDCOverlay = True ; This will always work ; res@tfDoNDCOverlay = "NDCViewport" ; NCL V6.5.0 or later res@gsnAddCyclic = False ;---Create a map to hold all lat/lon points res@tiMainString = "All lat/lon points: avg = " + t2_avg plot_all = gsn_csm_contour_map(wks,t2,res) ;---Create a map to hold subset of points using "where" res@tiMainString = "Using 'where' statement: avg = " + t2_sub1_avg plot1 = gsn_csm_contour_map(wks,t2_sub1,res) ;---Create a map to hold subset of points using "wrf_user_ll_to_i" res@tiMainString = "Using wrf_user_ll_to_xy: avg = " + t2_sub2_avg plot2 = gsn_csm_contour_map(wks,t2_sub2,res) ;---Create a map to hold subset of points using "getind_latlon2d" res@tiMainString = "Using getind_latlon2d: avg = " + t2_sub3_avg plot3 = gsn_csm_contour_map(wks,t2_sub3,res) ;---Add the lat/lon points to each plot mkres = True mkres@gsMarkerColor = "black" mkres@gsMarkerIndex = 16 mkres@gsMarkerSizeF = 0.002 dum0 = gsn_add_polymarker(wks,plot_all,lon2d, lat2d, mkres) dum1 = gsn_add_polymarker(wks,plot1, lon_sub1,lat_sub1,mkres) dum2 = gsn_add_polymarker(wks,plot2, lon_sub2,lat_sub2,mkres) dum3 = gsn_add_polymarker(wks,plot3, lon_sub3,lat_sub3,mkres) ;---Add a box showing the area of interest lnres = True lnres@gsLineColor = "ForestGreen" lnres@gsLineThicknessF = 5.0 dum4 = gsn_add_polyline(wks,plot_all,(/min_lon,max_lon,max_lon,\ min_lon,min_lon/),\ (/min_lat,min_lat,max_lat,\ max_lat,min_lat/),lnres) dum5 = gsn_add_polyline(wks,plot1,(/min_lon,max_lon,\ max_lon,min_lon,min_lon/),\ (/min_lat,min_lat,max_lat,\ max_lat,min_lat/),lnres) dum6 = gsn_add_polyline(wks,plot2,(/min_lon,max_lon,max_lon,\ min_lon,min_lon/),\ (/min_lat,min_lat,max_lat,\ max_lat,min_lat/),lnres) dum7 = gsn_add_polyline(wks,plot3,(/min_lon,max_lon,max_lon,\ min_lon,min_lon/),\ (/min_lat,min_lat,max_lat,\ max_lat,min_lat/),lnres) ;---Add the two points showing min/max lat/lon box corners mkres@gsMarkerSizeF = 0.007 mkres@gsMarkerColor = "NavyBlue" dum8 = gsn_add_polymarker(wks,plot_all,(/min_lon,max_lon/),\ (/min_lat,max_lat/),mkres) dum9 = gsn_add_polymarker(wks,plot1,(/min_lon,max_lon/), \ (/min_lat,max_lat/),mkres) dum10= gsn_add_polymarker(wks,plot2,(/min_lon,max_lon/), \ (/min_lat,max_lat/),mkres) dum11= gsn_add_polymarker(wks,plot3,(/min_lon,max_lon/), \ (/min_lat,max_lat/),mkres) ;---Panel the three plots pres = True pres@gsnMaximize = True pres@gsnPanelLabelBar = True pres@pmLabelBarWidthF = 0.8 pres@lbLabelFontHeightF = 0.015 gsn_panel(wks,(/plot_all,plot1,plot2,plot3/),(/2,2/),pres) end