[ncl-talk] Map won't display when using overlays in panel arrangement
mberdahl at envsci.rutgers.edu
mberdahl at envsci.rutgers.edu
Wed Aug 5 19:41:29 MDT 2015
Hi all,
I'm having trouble using overlays within panel plots. I would like to
have three panel plots, each with a map overlaid with a filled contour and
wind vectors. However, when I try to do so, I do not see the map - only
the contours and wind vectors. I can make the plot I'd like (with the
map) when just trying a single figure instead of panel plots. Does anyone
see what might be wrong with my code, copied below?
Thanks in advance,
Mira
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/csm/shea_util.ncl"
;************************************************
begin
;************************************************
; read in netCDF file s
;************************************************
a = addfile("uwnd.mon.mean.alllevels.nc","r") ; u winds
b = addfile("vwnd.mon.mean.alllevels.nc","r") ; v winds
c = addfile("../Geopotential/hgt.mon.mean.nc","r") ; geopotential heights.
;************************************************
; read in zonal [u] and meridional [v] winds (July)
;************************************************
u = a->uwnd(0:803,{500},{45:90},{270:357.5})
v = b->vwnd(0:803,{500},{45:90},{270:357.5}) ; Get u, v, time (1),level
(1000hpa),latitude(-90:90) and longitude(0:360) data.
z = c->hgt(0:803,{500},{45:90},{270:357.5}) ; get geopotenial heights...
printVarSummary(u)
printVarSummary(v)
printVarSummary(z)
; Calculate the seasonal averages.
uDJF = month_to_season(u, "DJF")
vDJF = month_to_season(v, "DJF")
zDJF = month_to_season(z, "DJF")
printVarSummary(uDJF)
printVarSummary(vDJF)
printVarSummary(zDJF)
; from the matlab script i wrote: findExtremeYrs, i pulled out the extreme
years (> or < 1std) that i want to average and plot here.
; for ans = 4 (NAO)
; yearList_hi = 1973 1975 1983 1989 1995
2000 2007 2012
; yearList_lo = 1963 1964 1965 1969 1977
1979 1996 1997 2010 2011
; this data starts at 1948 (this is index 0), so 1953=5, 1963=15 etc.
uDJF_NAO_hi = uDJF((/25,27,35,41,47,52,59,64/),:,:)
uDJF_NAO_lo = uDJF((/15,16,17,21,29,31,48,49,62,63/),:,:)
vDJF_NAO_hi = vDJF((/25,27,35,41,47,52,59,64/),:,:)
vDJF_NAO_lo = vDJF((/15,16,17,21,29,31,48,49,62,63/),:,:)
zDJF_NAO_hi = zDJF((/25,27,35,41,47,52,59,64/),:,:)
zDJF_NAO_lo = zDJF((/15,16,17,21,29,31,48,49,62,63/),:,:)
uAvgTime_hi = dim_avg_n_Wrap(uDJF_NAO_hi,0)
uAvgTime_lo = dim_avg_n_Wrap(uDJF_NAO_lo,0)
printVarSummary(uAvgTime_hi)
printVarSummary(uAvgTime_lo)
vAvgTime_hi = dim_avg_n_Wrap(vDJF_NAO_hi,0)
vAvgTime_lo = dim_avg_n_Wrap(vDJF_NAO_lo,0)
printVarSummary(vAvgTime_hi)
printVarSummary(vAvgTime_lo)
zAvgTime_hi = dim_avg_n_Wrap(zDJF_NAO_hi,0)
zAvgTime_lo = dim_avg_n_Wrap(zDJF_NAO_lo,0)
printVarSummary(zAvgTime_hi)
printVarSummary(zAvgTime_lo)
; dirty way to copy metadata over first.
diff_u = uAvgTime_hi;
diff_v = vAvgTime_hi;
diff_z = zAvgTime_hi;
diff_u = uAvgTime_hi - uAvgTime_lo
diff_v = vAvgTime_hi - vAvgTime_lo
diff_z = zAvgTime_hi - zAvgTime_lo
printVarSummary(diff_u)
printVarSummary(diff_v)
printVarSummary(diff_z)
;************************************************
; create plot
;************************************************
wks = gsn_open_wks("ps","Panel_NAO_z_500") ; open a ps file
gsn_define_colormap(wks,"temp1")
plot = new(3,graphic) ; create a plot array
;---- set common resources for all plots
res = True
res at gsnDraw = False ; dont draw
res at gsnFrame = False ; dont advance frame
res at cnInfoLabelOn = False ; trn off cn info label
res at gsnAddCyclic = False ; has to do with wrapping the longitude at 0/360
res at cnFillPalette = "matlab_jet"
;************************************************
; Choose a subregion
;************************************************
res at mpMaxLatF = 90 ;maximum latitude
res at mpMinLatF = 45 ;minimum latitude
;res at mpMaxLonF = 357.5 ; ;maximum longitude
;res at mpMinLonF = 270 ;minimum longitude
res at mpMaxLonF = 0
res at mpMinLonF = -90
;res at mpFillBoundarySets = AllBoundaries
res at mpOutlineBoundarySets = "National"
res at mpOutlineOn = True
res at mpOutlineDrawOrder = "PostDraw"
mpid = gsn_csm_map(wks,res)
;***********************************************
; ----wind vector plot
;***********************************************
vcres = res
vcres at vcRefAnnoOrthogonalPosF = -1.0 ; move ref vector up
vcres at vcRefMagnitudeF = 10.0 ; define vector ref mag
vcres at vcRefLengthF = 0.045 ; define length of vec ref
vcres at vcGlyphStyle = "CurlyVector" ; turn on curly vectors
vcres at vcMinDistanceF = 0.017
vcres at mpFillOn = False ; turn off gray fill
vcres at mpOutlineBoundarySets = "National" ; turn on country boundaries
;vcres at mpFillBoundarySets = AllBoundaries
vcres at mpGeophysicalLineColor = "Navy" ; color of cont. outlines
vcres at mpGeophysicalLineThicknessF = 1.5 ; thickness of outlines
;vcres at gsnLeftString = "DJF High NAO"
; was previously winds_hi =
gsn_csm_vector_map_ce(wks,uAvgTime_hi,vAvgTime_hi,vcres)
winds_hi = gsn_csm_vector(wks,uAvgTime_hi,vAvgTime_hi,vcres)
;vcres at gsnLeftString = "DJF Low NAO"
winds_lo = gsn_csm_vector(wks,uAvgTime_lo,vAvgTime_lo,vcres)
;vcres at gsnLeftString = "Difference of High - Low"
winds_diff = gsn_csm_vector(wks, diff_u, diff_v,vcres)
;************************************************
;---- geopotential height filled contour plot
;***********************************************
zfres = res
zfres at cnFillOn = True
;zfres at cnLevelSelectionMode = "ExplicitLevels
;zfres at cnLevels = ispan(-20,90,5)
zfres at lbLabelFontHeightF = 0.015
zfres at lbOrientation = "Vertical"
zfres at pmLabelBarOrthogonalPosF = -0.005
zfres at cnFillPalette = "BlWhRe"
contour_zf_hi = gsn_csm_contour(wks,zAvgTime_hi,zfres)
contour_zf_lo = gsn_csm_contour(wks,zAvgTime_lo,zfres)
contour_zf_diff = gsn_csm_contour(wks,diff_z,zfres)
plot(0) = gsn_csm_contour(wks,zAvgTime_hi,zfres)
plot(1) = gsn_csm_contour(wks,zAvgTime_lo,zfres)
plot(2) = gsn_csm_contour(wks,diff_z,zfres)
;overlay(plot(0),winds_hi)
overlay(mpid,plot(0))
overlay(mpid,winds_hi)
;overlay(plot(1),winds_lo)
overlay(mpid,plot(1))
overlay(mpid,winds_lo)
;overlay(plot(2),winds_diff)
overlay(mpid,plot(2))
overlay(mpid,winds_lo)
draw(mpid)
frame(wks)
;************************************************
; create panel
;************************************************
resP = True ; modify the panel plot
resP at txString = "NAO 500mb"
gsn_panel(wks,plot,(/3,1/),resP) ; now draw as one plot;
end
More information about the ncl-talk
mailing list