[ncl-talk] panel plots with overlays
mberdahl at envsci.rutgers.edu
mberdahl at envsci.rutgers.edu
Sun Jul 26 21:29:56 MDT 2015
Hi all,
I'm trying to adapt a working code I had that was plotting 3 panel plots
of wind vectors in a certain region. Now, all I want to do is add an
overlay of geopotential height (filled contours) underneath each of these
plots, but the way I'm approaching it is wrong. I've tried simplifying it
to do this without panels, just a single plot, but I still don't get what
I need. The error message I am getting now is:
fatal:syntax error: line 153 in file plotWinds_z_NAO_level.ncl before or
near overlay
plot(0)= overlay
---------------^
My code is copied below. Any suggestions are greatly appreciated.
Thanks,
Mira
;*****************************************
; ; plot average winds overlayed on geopotental filled contours for the
years that are extreme
; the years are found with my matlab script findExtremeYrs.m
; this particlar script does the years that correspond to high correlation
with DJF, SE precip and NAO
;*****************************************
; the original data goes from 1948 January to April 2015.
; I will cut it out to 2014 December, so we have something divisible by 12
so we can do seasonal averages...
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
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 = True ; plot mods desired
res at gsnAddCyclic = False ; has to do with wrapping the longitude at 0/360
;************************************************
; 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
;***********************************************
; ----wind vector plot
;***********************************************
vcres = True
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 mpGeophysicalLineColor = "Navy" ; color of cont. outlines
vcres at mpGeophysicalLineThicknessF = 1.5 ; thickness of outlines
;vcres at gsnLeftString = "DJF High NAO"
winds_hi = gsn_csm_vector_map_ce(wks,uAvgTime_hi,vAvgTime_hi,vcres)
;vcres at gsnLeftString = "DJF Low NAO"
winds_lo = gsn_csm_vector_map_ce(wks,uAvgTime_lo,vAvgTime_lo,vcres)
;vcres at gsnLeftString = "Difference of High - Low"
winds_diff = gsn_csm_vector_map_ce(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
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) = overlay(contour_zf_hi,winds_hi)
plot(1) = overlay(contour_zf_lo,winds_lo)
plot(2) = overlay(contour_diff,winds_diff)
;************************************************
; 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