[ncl-talk] panel plots with overlays
mberdahl at envsci.rutgers.edu
mberdahl at envsci.rutgers.edu
Mon Aug 3 13:31:33 MDT 2015
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>
>>>>
>>>>
>>>
>>
>>
>
More information about the ncl-talk
mailing list