[ncl-talk] panel plots with overlays
mberdahl at envsci.rutgers.edu
mberdahl at envsci.rutgers.edu
Mon Jul 27 20:39:43 MDT 2015
Hi Alex,
Thanks for the suggestion. I tried this and now I do not have the same
error message any more, so that is good. However, I am now just getting 3
blank plots of the whole world, even though I'd like plots with overlaid
geopotential height and winds, for just a region over Greenland. These are
followed by a fourth figure that has filled contours of the correct
region, but no map or wind vectors. I'll attach the figures as an
attachment here. They are not panel plots as I want, but each figure is
on its own page. I'll copy the updated code below.
Any thoughts?
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) = 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(contour_zf_hi,winds_hi)
overlay(contour_zf_lo,winds_lo)
overlay(contour_zf_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
> Mira,
>
> Try this:
>
>
> 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(contour_zf_hi,winds_hi)
> overlay(contour_zf_lo,winds_lo)
> overlay(contour_diff,winds_diff)
>
> I don't believe overlay returns anything and you are trying to save it as
> a
> graphic. The code above would replace the lines in your script.
>
> Hope that helps,
>
> -Alex
>
>
>
>
> On Sunday, July 26, 2015, <mberdahl at envsci.rutgers.edu
> <javascript:_e(%7B%7D,'cvml','mberdahl at envsci.rutgers.edu');>> wrote:
>
>> 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
>>
>>
>> _______________________________________________
>> ncl-talk mailing list
>> ncl-talk at ucar.edu
>> List instructions, subscriber options, unsubscribe:
>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>
>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Panel_NAO_z_500.pdf
Type: application/pdf
Size: 201033 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20150727/5843798b/attachment-0001.pdf
More information about the ncl-talk
mailing list