[ncl-talk] panel plots with overlays
Adam Phillips
asphilli at ucar.edu
Mon Aug 3 15:42:12 MDT 2015
Hi Mira,
Try this:
1) Create three distinct maps:
mpid = gsn_csm_map(wks,res)
mpid2 = gsn_csm_map(wks,res)
mpid3 = gsn_csm_map(wks,res)
; .... create the 6 plots as you have been doing:
winds_hi = gsn_csm_vector(wks,uAvgTime_hi,vAvgTime_hi,vcres)
winds_lo = gsn_csm_vector(wks,uAvgTime_lo,vAvgTime_lo,vcres)
winds_diff = gsn_csm_vector(wks, diff_u, diff_v,vcres)
(snip)
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(mpid,plot(0))
overlay(mpid,winds_hi)
overlay(mpid2,plot(1))
overlay(mpid2,winds_lo)
overlay(mpid3,plot(2))
overlay(mpid3,winds_diff)
;************************************************
; create panel
;************************************************
resP = True ; modify the panel plot
resP at txString = "NAO 500mb"
gsn_panel(wks,(/mpid,mpid2,mpid3/),(/3,1/),resP)
Three things to note:
1) If mpid is used as the base plot in overlay (=the first plot indexed),
then that is what you should pass to gsn_panel.
2) You could have called gsn_csm_vector_map instead of gsn_csm_vector to
draw the vectors and the map at the same time if you had wanted.
3) I would recommend creating new resource lists for your zres and vres
resource list (as opposed to creating them from res) to eliminate all the
errors one gets when passing mp* resources to non map plotting functions.
Hope that helps. If not, or if you have any further questions please
respond to the ncl-talk email list.
Adam
On Mon, Aug 3, 2015 at 1:31 PM, <mberdahl at envsci.rutgers.edu> wrote:
> Hi again Alex,
>
> OK - another quick update. When I try to not do panel plots, and just try
> to plot one single figure (I comment the 2 following panels out), then I
> receive the plot that I'm interested in! So it looks like the problem is
> just trying to get this to work as panels, which I'm unclear on how to do.
> The working script for the single plot is copied below.
>
> Best,
> 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
> 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
>
>
>
>
>
>
>
>
>
> Actually, I take that back. I don't get a blank plot (that was a stupid
> error), but I do just get the same thing as before... Overlaid contours
> and vectors, with no base map.
>
> mira
>
> > Hi Alex,
> >
> > Yes! The plot you attached is the map I'd like to have under my wind
> vectors and filled contours. However, when I change the max and min lon
> in my map resources to what you did (0 and -90), I do not receive any
> map - it is a blank page with just a title.
> >
> > I wonder if it has to do with the lines where I extract the longitudes
> in my variables:
> >
> > 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 geopotential
> heights...
> >
> > Do you think it is something to do with this?
> >
> > Mira
> >
> >
> >
> >> Hi Mira,
> >>
> >> I copied these lines from your script and ran them locally:
> >>
> >> 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 mpFillBoundarySets = AllBoundaries
> >> res at mpOutlineBoundarySets = "National"
> >> res at mpOutlineOn = True
> >> res at mpOutlineDrawOrder = "PostDraw"
> >>
> >> mpid = gsn_csm_map(wks,res)
> >>
> >>
> >> I got an error and made a change in your lons. I set them to 0 and -90
> respectively and got this base map. Is this what you were looking for?
> >>
> >> -Alex
> >>
> >>
> >>
> >>
> >>> On Jul 30, 2015, at 9:00 PM, mberdahl at envsci.rutgers.edu wrote:
> >>>
> >>> Hi Alex,
> >>> Thanks again for the help. Setting the workstation colormap worked
> once
> >>> I
> >>> realized I have an older version (6.0.0) of NCL so there is a limited
> set
> >>> of colormaps available.
> >>> I still, however, can't get a base map plotted. I've tried a couple
> things, including what you suggested below, but this doesn't work.
> I've
> >>> tried creating a map and overlaying it with the contour and vectors,
> but
> >>> it doesn't work. I receive the following error and warnings:
> >>>
> >>>
> >>> ^Mfatal:NhlAddOverlay: tranform is already an annotation or overlay:
> 74 ^Mwarning:ContourPlotSetValues: attempt to set overlay member plot
> view ignored
> >>> ^Mwarning:ContourPlotSetValues: attempt to set overlay member plot
> view ignored
> >>> ^Mwarning:ContourPlotSetValues: attempt to set overlay member plot
> view ignored
> >>> ^Mwarning:NhlDraw: cannot draw Plot Member, ID 262, independently
> ^Mwarning:NhlDraw: cannot draw Plot Member, ID 298, independently
> ^Mwarning:NhlDraw: cannot draw Plot Member, ID 334, independently
> >>>
> >>>
> >>>
> >>> My code is copied below.
> >>>
> >>>
> >>> Any thoughts?
> >>> Thanks again,
> >>> 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 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)
> >>>
> >>>
> >>> ;************************************************
> >>> ; 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
> >>>
> >>>
> >>>
> >>>
> >>>
> >>>
> >>>> Hi Mira,
> >>>>
> >>>> As for the color bar, I would suggest either using cnFillColors or
> setting
> >>>> the color table for the workstation, I usually just set the color
> table
> >>>> for workstation.
> >>>>
> >>>> After you open the workstation:
> >>>>
> >>>> wks = gsn_open_wks("ps","Panel_NAO_z_500")
> >>>>
> >>>> Do:
> >>>> gsn_define_colormap(wks, “matlab_jet†)
> >>>>
> >>>> You can choose the color palette that you want and put the name in
> quotes.
> >>>>
> >>>> As for the coastlines, maybe set:
> >>>>
> >>>> mpDataSetName = “Earth..4â€
> >>>>
> >>>> You can also try setting mpFillBoundarySets = AllBoundaries and see
> if that helps.
> >>>>
> >>>> Hope that helps,
> >>>> -Alex
> >>>>
> >>>>
> >>>>
> >>>>
> >>>>> On Jul 29, 2015, at 8:28 PM, mberdahl at envsci.rutgers.edu wrote:
> >>>>>
> >>>>> Hi Alex,
> >>>>>
> >>>>> Great! Now I have three panels, in the correct region, with filled
> contours and overlaid wind vectors. Thanks! I have two outstanding
> problems still: I'm trying to change the color bar for the filled
> contours, as the default in the older version of ncl is not very
> appealing. When I try to do res at cnFillPalette = "matlab_jet", or
> any other choice for that matter, I receive the following warning:
> >>>>>
> >>>>> warning:cnFillPalette is not a valid resource in
> >>>>> Panel_NAO_z_500_vector
> >>>>> at
> >>>>> this time
> >>>>>
> >>>>> I'm not exactly sure why I can't seem to change the contour fill
> palette.
> >>>>>
> >>>>> Secondly, I still don't see a map of coastlines on my figures. Do
> you
> >>>>> know how to add one?
> >>>>>
> >>>>> Thank you very much again,
> >>>>> Mira
> >>>>>
> >>>>>
> >>>>>> Hi Mira,
> >>>>>>
> >>>>>> You arenÂ’t seeing vectors or map because they a buried under the
> filled
> >>>>>> contours. You may want to flip the plots around and do the
> contours as
> >>>>>> a
> >>>>>> map and overlay the vectors on top of the contours.
> >>>>>>
> >>>>>> Upon further analysis of your script, you were doing the overlay on
> contour plots but then using different contour plots for the panel
> plot.
> >>>>>> Below I have altered your script, try it this way and see how it
> goes.
> >>>>>>
> >>>>>> I am not sure if you have made any more recent changes to this
> script
> >>>>>> but
> >>>>>> this is using the last one you sent.
> >>>>>>
> >>>>>> Hope that helps,
> >>>>>> -Alex
> >>>>>>
> >>>>>>
> >>>>>> ;*****************************************
> >>>>>> ; ; 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 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
> >>>>>>
> >>>>>> ;************************************************
> >>>>>> ; Map Options
> >>>>>> ;************************************************
> >>>>>> res at mpOutlineBoundarySets = "National" ; turn on country
> boundaries res at mpGeophysicalLineColor = "Navy" ; color of cont.
> outlines res at mpGeophysicalLineThicknessF = 1.5 ; thickness of
> outlines
> >>>>>>
> >>>>>> ;***********************************************
> >>>>>> ; ----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 gsnLeftString = "DJF High NAO"
> >>>>>> 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
> >>>>>>
> >>>>>> plot(0) = gsn_csm_contour_map_ce(wks,zAvgTime_hi,zfres)
> >>>>>> plot(1) = gsn_csm_contour_map_ce(wks,zAvgTime_lo,zfres)
> >>>>>> plot(2) = gsn_csm_contour_map_ce(wks,diff_z,zfres)
> >>>>>>
> >>>>>> overlay(plot(0),winds_hi)
> >>>>>> overlay(plot(1),winds_lo)
> >>>>>> overlay(plot(2),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
> >>>>>>
> >>>>>>
> >>>>>>> On Jul 29, 2015, at 1:57 PM, mberdahl at envsci.rutgers.edu wrote:
> >>>>>>>
> >>>>>>> Hi Alex,
> >>>>>>> Great - thanks for the suggestions! That does help, but does not
> totally
> >>>>>>> solve my problem. Now I get three panel plots of the correct
> region
> >>>>>>> and
> >>>>>>> all on one page which is great. But only the filled contours
> show. There
> >>>>>>> are no overlaid vectors or map of coastal boundaries. Any
> suggestions
> >>>>>>> would be greatly appreciated.
> >>>>>>> Thanks so much!
> >>>>>>> Mira
> >>>>>>>
> >>>>>>>
> >>>>>>>
> >>>>>>>> Hi Mira,
> >>>>>>>>
> >>>>>>>> I think the issue here is that you set a bunch of resources as
> res but
> >>>>>>>> then set vres to True instead of res. I think switching it to
> vres
> >>>>>>>> =
> >>>>>>>> res
> >>>>>>>> will give you the maps you were expecting. I see you also set
> res =
> >>>>>>>> True
> >>>>>>>> twice, might be helpful to remove the double.
> >>>>>>>>
> >>>>>>>> Hope that helps,
> >>>>>>>> -Alex
> >>>>>>>>
> >>>>>>>>
> >>>>>>>>> On Jul 27, 2015, at 9:39 PM, mberdahl at envsci.rutgers.edu wrote:
> >>>>>>>>>
> >>>>>>>>> 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
> <mailto:mberdahl at envsci.rutgers.edu
> >>>>>>>>>> <mailto:mberdahl at envsci.rutgers.edu
> >>>>>>>>>> <mailto: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
> >>>>>>>>>>>
> >>>>>>>>>>
> >>>>>>>>> <Panel_NAO_z_500.pdf>
> >>>>
> >>>>
> >>>
> >>
> >>
> >
>
>
>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at ucar.edu
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
--
Adam Phillips
Associate Scientist, Climate and Global Dynamics Laboratory, NCAR
www.cgd.ucar.edu/staff/asphilli/ 303-497-1726
<http://www.cgd.ucar.edu/staff/asphilli>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20150803/d738911f/attachment.html
More information about the ncl-talk
mailing list