[ncl-talk] panel plots with overlays
mberdahl at envsci.rutgers.edu
mberdahl at envsci.rutgers.edu
Thu Jul 30 20:00:18 MDT 2015
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