[ncl-talk] panel plots with overlays

mberdahl at envsci.rutgers.edu mberdahl at envsci.rutgers.edu
Wed Jul 29 12:57:46 MDT 2015


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>');>> 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>
>
>




More information about the ncl-talk mailing list