[ncl-talk] panel plots with overlays
Alex Schaefer
alexander.schaefer at mines.sdsmt.edu
Wed Jul 29 20:46:40 MDT 2015
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>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20150729/e4c6acae/attachment.html
More information about the ncl-talk
mailing list