[ncl-talk] panel plots with overlays

mberdahl at envsci.rutgers.edu mberdahl at envsci.rutgers.edu
Wed Jul 29 19:28:37 MDT 2015


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