Hello,
I wrote a script that subtracts NOHRSC snow depth data from my WRF modeled
snow depth data and plots the difference in a panel plot. The NOHRSC snow
depth data are only available within US boundaries.
When I plot up my data, values outside of US boundaries appear on the very
high end of my range. I tried to use 'mask' as shown in mask_3.ncl to
exclude those values but it does not appear to be working. The high
values do not show up when I use 'max' to determine the maximum value:
For example, max_xdiff is the maximum value before I implement mask.
Variable: max_xdiff
Type: float
Total Size: 4 bytes
1 values
Number of Dimensions: 1
Dimensions and sizes: [1]
Coordinates:
(0) 1.240401
max_range_only is the maximum value after I implement mask.
Variable: max_range_only
Type: float
Total Size: 4 bytes
1 values
Number of Dimensions: 1
Dimensions and sizes: [1]
Coordinates:
(0) 1.240401
As you can see, the max value does not change, however, the high values
located outside of the NOHRSC domain are still plotting up on the high
range of my scale.
Full script is below.
Thanks,
Liz
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/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl"
;============================================================
; The main code
;============================================================
begin
;---Average the daily NOHRSC data by month
varn = "mslt"
dirn = "/glade/p/work/burakows/nohrsc/mNOHRSC_mslt_200811_200904/"
filn = systemfunc("cd "+dirn+" ; ls "+varn+"*.nc")
nfiln = dimsizes(filn)
;Loop through all of the filenames to retrieve year and month from daily
files
do nfn=0,nfiln-1,32
yyyymm = toint( str_get_field(filn(nfn),2,"_.") ) ; parse
NOHRSC name for time
yyyy = yyyymm/10000 ; current
year as integer
mm = (yyyymm-(yyyy*10000))/100 ; current month
ncocmd = "ncea -O -v z
"+dirn+varn+"_"+yyyy+sprinti("%0.2i",mm)+"*.nc
"+dirn+"avg_"+varn+"_"+yyyy+sprinti("%0.2i",mm)+".nc"
print("ncocmd = "+ncocmd) ; print the nco command
system(ncocmd) ; execute the nco command
end do
print(systemfunc("ls "+dirn+"avg_*.nc"))
;---Remove nfn, nfiln, and filn when loop is complete
delete ([/nfn,nfiln,filn/])
;---Specify interpolation to be used
method = "conserve"
wgtFileDir = "/glade/p/work/burakows/plot/NCL/ESMF_regridding/"
wgtFileName = "NOHRSC_1kmM2_to_WRF."+method+"_wgt.nc"
wgtFilePath = wgtFileDir+wgtFileName
print("wgtFilePath="+wgtFilePath)
;---NOHRSC files [input]
varn = "mslt" ; NOHRSC
name
dirn = "/glade/p/work/burakows/nohrsc/mNOHRSC_mslt_200811_200904/"
filn = systemfunc("cd "+dirn+" ; ls avg_"+varn+"*.nc")
nfiln = dimsizes(filn)
print(filn(0:nfiln-1)) ; 1st 5 file
;---WRF directory
varw = "SNOWH" ; WRF name
dirw =
"/glade/p/work/burakows/plot/albedo_CLMV35current_200811-200904/"
varm = "SNOWH" ; WRF monthly
name
prtFlg = True
wrfo = addfile(dirw+"wrfout_d03_2008-11-01_00:00:00.nc","r") ;
read any wrfout file
x = wrfo->T2 ; read T2 and preserve metadata
x@lat2d = wrfo->XLAT(0,:,:) ; Assign lat
to x
x@lon2d = wrfo->XLONG(0,:,:) ; Assign lon
to x
minlat = min(x@lat2d)
maxlat = max(x@lat2d)
minlon = min(x@lon2d)
maxlon = max(x@lon2d)
bndadd = 0.10 ; spacing around the plot edges
;----------------------------------------------------------------------
; Plotting options section
;----------------------------------------------------------------------
pltType = "ps" ; plot type
pltDir = "./" ; plot directory
pltName = "Diff_SNOWH_WRFCLM_NOHRSC_2009" ; plot name (ps file)
pltPath = pltDir+pltName ; plot path
wks = gsn_open_wks(pltType,pltPath) ; create workstation
for ps file
gsn_define_colormap(wks,"GreenMagenta16") ; define color
table (temperature)
res = True
res@gsnDraw = False
res@gsnFrame = False
res@cnFillOn = True ; color plot desired
res@cnLinesOn = False ; turn off contour
lines
res@cnLineLabelsOn = False ; turn off contour
labels
res@cnInfoLabelOn = False ; turn off info
label (top labels of indvid. plots)
res@cnFillMode = "RasterFill" ; turn raster on
res@cnLevelSelectionMode = "ManualLevels" ; Set contour levels
manually
res@cnMinLevelValF = -.5 ; minimum contour
(degrees C, or mm)
res@cnMaxLevelValF = .5 ; maximum contour
(degrees C, or mm)
res@cnLevelSpacingF = .1 ; contour interval
res@cnNoDataLabelOn = False ; Attempt to mask out the
no data point over ocean and Canada
res@lbLabelBarOn = False ; Will turn on in
panel later
res@lbOrientation = "Horizontal" ; Horizontal label bar
res@mpFillOn = False
res@mpOutlineOn = True
res@mpOutlineBoundarySets = "AllBoundaries"
res@mpProjection = "CylindricalEquidistant"
res@mpLimitMode = "LatLon" ; required
res@mpMinLatF = minlat-bndadd
res@mpMaxLatF = maxlat+bndadd
res@mpMinLonF = minlon-bndadd
res@mpMaxLonF = maxlon+bndadd
res@mpCenterLonF = (minlon + maxlon)*0.5
res@mpCenterLatF = (minlat + maxlat)*0.5
res@gsnLeftString = "" ; Turn off left
subtitle
res@gsnRightString = "" ; Turn off right
subtitle
res@gsnMajorLatSpacing = 2
;res@gsnMajorLonSpacing = 2
res@gsnMinorLonSpacing = 2
res@cnLinesOn = False
res@cnLineLabelsOn = False
res@cnFillOn = True
res@cnFillMode = "RasterFill"
res@gsnAddCyclic = False ; regional grid (changes
central meridian)/xwo
;---Create new ncl variable to hold panel plot contents
plot_diff = new(6,graphic) ; from p.80 of
NCL manual, right before looping (panel1e.ncl)
WrfNOHRSCdiff = new((/6,192,132/),float,1e20) ; create
new variable WrfPr(time,y,x)
labels = new(6,string) ; create new
variable labels for figure labels in panel plot
wrfdim = dimsizes(WrfNOHRSCdiff)
;----------------------------------------------------------------------
; Loop through six NOHRSC and WRF monthly files to calculate difference
;----------------------------------------------------------------------
do nfn=0,nfiln-1
yyyymm = toint( str_get_field(filn(nfn),3,"_.") ) ; parse
NOHRSC name for time
yyyy = yyyymm/100 ; current
year as integer
mm = yyyymm-(yyyy*100) ; current month
fn = addfile(dirn+filn(nfn), "r")
XN = fn->z ; original
NOHRSC variable (all 'z')
xn = ESMF_regrid_with_weights(XN,wgtFilePath,False); create new
variable on WRF grid
xnm = xn/1000
; matching
WRF file
filw = "avg_"+varm+"_"+yyyy+"-"+sprinti("%0.2i",mm)+".nc" ;
average temperature, snowh, snow
if (isfilepresent(dirw+filw)) then
fw = addfile(dirw+filw, "r") ; Open WRF
monthly file (if present)
xw = fw->$varm$ ; WRF
monthly variable
;----------------------------------------------------------------------
; Subtract NOHRSC from WRF
;----------------------------------------------------------------------
fwo = addfile(dirw+"wrfout_d03_2008-11-01_00:00:00.nc","r"); Open any
wrfout daily file to obtain
; lat and lon
;xwo = fwo->SNOWH ; read T2
and preserve metadata
xdiff = xw(0,:,:)-xnm ; Avg
monthly WRF minus Prism for month
max_xdiff = max(xdiff)
print(max_xdiff)
range_only = xdiff
range_only = mask(xdiff,(xdiff.le.5),True) ; Mask out values
less than <-8 (basically NaN, outside domain)
range_only@long_name = "difference" + yyyymm ; Assign
long name to difference
range_only@units = "mm/month" ; Assign
units to difference
range_only@_FillValue = -9999 ; Assign fill value to -9999
range_only@lat2d = fwo->XLAT(0,:,:) ; Assign lat
to range_only
range_only@lon2d = fwo->XLONG(0,:,:) ; Assign lon
to range_only
printVarSummary(range_only)
labels(nfn) = yyyymm ; assign
yyyymm of current plot to labels for panel plot
WrfNOHRSCdiff(nfn,:,:)=range_only ; assign
range_only to WrfNOHRSCdiff variable
plot_diff(nfn) =
gsn_csm_contour_map(wks,WrfNOHRSCdiff(nfn,:,:),res) ; create
plot of month nfn in plot_diff
if (prtFlg) then ; first
time only
printVarSummary(XN)
printVarSummary(xn)
printVarSummary(xw)
prtFlg = False
end if ; nfn
end if ; isfilepresent
end do
printVarSummary(WrfNOHRSCdiff)
;----------------------------------------------------------------------
; Panel Plotting options section
;----------------------------------------------------------------------
pres = True
pres@gsnMaximize = True
pres@gsnPanelLabelBar = True
pres@gsnPanelYWhiteSpace= 5
pres@gsnPanelXWhiteSpace= 5
pres@txString = "Monthly total SNOWH WRFCLM minus NOHRSC"
pres@gsnPanelFigureStrings = labels
gsn_panel(wks,plot_diff,(/2,3/),pres)
end
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Fri Feb 28 15:25:15 2014
This archive was generated by hypermail 2.1.8 : Mon Mar 03 2014 - 14:26:18 MST