[ncl-talk] panel plots with overlays
Alex Schaefer
alexander.schaefer at mines.sdsmt.edu
Mon Jul 27 21:12:17 MDT 2015
Hi Mira,
I think the issue here is that you set a bunch of resources as res but then set vres to True instead of res. I think switching it to vres = res will give you the maps you were expecting. I see you also set res = True twice, might be helpful to remove the double.
Hope that helps,
-Alex
> On Jul 27, 2015, at 9:39 PM, mberdahl at envsci.rutgers.edu wrote:
>
> Hi Alex,
> Thanks for the suggestion. I tried this and now I do not have the same
> error message any more, so that is good. However, I am now just getting 3
> blank plots of the whole world, even though I'd like plots with overlaid
> geopotential height and winds, for just a region over Greenland. These are
> followed by a fourth figure that has filled contours of the correct
> region, but no map or wind vectors. I'll attach the figures as an
> attachment here. They are not panel plots as I want, but each figure is
> on its own page. I'll copy the updated code below.
> Any thoughts?
> Thanks!
> Mira
>
>
>
>
> ;*****************************************
> ; ; plot average winds overlayed on geopotental filled contours for the
> years that are extreme
> ; the years are found with my matlab script findExtremeYrs.m
> ; this particlar script does the years that correspond to high correlation
> with DJF, SE precip and NAO
> ;*****************************************
> ; the original data goes from 1948 January to April 2015.
> ; I will cut it out to 2014 December, so we have something divisible by 12
> so we can do seasonal averages...
>
>
> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
> ;************************************************
> begin
> ;************************************************
> ; read in netCDF file s
> ;************************************************
> a = addfile("uwnd.mon.mean.alllevels.nc","r") ; u winds
> b = addfile("vwnd.mon.mean.alllevels.nc","r") ; v winds
>
> c = addfile("../Geopotential/hgt.mon.mean.nc","r") ; geopotential heights.
>
> ;************************************************
> ; read in zonal [u] and meridional [v] winds (July)
> ;************************************************
>
> u = a->uwnd(0:803,{500},{45:90},{270:357.5})
> v = b->vwnd(0:803,{500},{45:90},{270:357.5}) ; Get u, v, time (1),level
> (1000hpa),latitude(-90:90) and longitude(0:360) data.
> z = c->hgt(0:803,{500},{45:90},{270:357.5}) ; get geopotenial heights...
>
> printVarSummary(u)
> printVarSummary(v)
> printVarSummary(z)
>
> ; Calculate the seasonal averages.
> uDJF = month_to_season(u, "DJF")
> vDJF = month_to_season(v, "DJF")
> zDJF = month_to_season(z, "DJF")
>
> printVarSummary(uDJF)
> printVarSummary(vDJF)
> printVarSummary(zDJF)
>
> ; from the matlab script i wrote: findExtremeYrs, i pulled out the extreme
> years (> or < 1std) that i want to average and plot here.
>
> ; for ans = 4 (NAO)
> ; yearList_hi = 1973 1975 1983 1989 1995
> 2000 2007 2012
> ; yearList_lo = 1963 1964 1965 1969 1977
> 1979 1996 1997 2010 2011
>
>
> ; this data starts at 1948 (this is index 0), so 1953=5, 1963=15 etc.
>
> uDJF_NAO_hi = uDJF((/25,27,35,41,47,52,59,64/),:,:)
> uDJF_NAO_lo = uDJF((/15,16,17,21,29,31,48,49,62,63/),:,:)
>
> vDJF_NAO_hi = vDJF((/25,27,35,41,47,52,59,64/),:,:)
> vDJF_NAO_lo = vDJF((/15,16,17,21,29,31,48,49,62,63/),:,:)
>
> zDJF_NAO_hi = zDJF((/25,27,35,41,47,52,59,64/),:,:)
> zDJF_NAO_lo = zDJF((/15,16,17,21,29,31,48,49,62,63/),:,:)
>
>
> uAvgTime_hi = dim_avg_n_Wrap(uDJF_NAO_hi,0)
> uAvgTime_lo = dim_avg_n_Wrap(uDJF_NAO_lo,0)
>
> printVarSummary(uAvgTime_hi)
> printVarSummary(uAvgTime_lo)
>
> vAvgTime_hi = dim_avg_n_Wrap(vDJF_NAO_hi,0)
> vAvgTime_lo = dim_avg_n_Wrap(vDJF_NAO_lo,0)
>
> printVarSummary(vAvgTime_hi)
> printVarSummary(vAvgTime_lo)
>
> zAvgTime_hi = dim_avg_n_Wrap(zDJF_NAO_hi,0)
> zAvgTime_lo = dim_avg_n_Wrap(zDJF_NAO_lo,0)
>
> printVarSummary(zAvgTime_hi)
> printVarSummary(zAvgTime_lo)
>
> ; dirty way to copy metadata over first.
> diff_u = uAvgTime_hi;
> diff_v = vAvgTime_hi;
> diff_z = zAvgTime_hi;
>
> diff_u = uAvgTime_hi - uAvgTime_lo
> diff_v = vAvgTime_hi - vAvgTime_lo
> diff_z = zAvgTime_hi - zAvgTime_lo
>
> printVarSummary(diff_u)
> printVarSummary(diff_v)
> printVarSummary(diff_z)
>
> ;************************************************
> ; create plot
> ;************************************************
> wks = gsn_open_wks("ps","Panel_NAO_z_500") ; open a ps file
> plot = new(3,graphic) ; create a plot array
>
> ;---- set common resources for all plots
> res = True
> res at gsnDraw = False ; dont draw
> res at gsnFrame = False ; dont advance frame
> res at cnInfoLabelOn = False ; trn off cn info label
> res = True ; plot mods desired
> res at gsnAddCyclic = False ; has to do with wrapping the longitude at 0/360
> ;************************************************
> ; Choose a subregion
> ;************************************************
> res at mpMaxLatF = 90 ;maximum latitude
> res at mpMinLatF = 45 ;minimum latitude
> res at mpMaxLonF = 357.5 ; ;maximum longitude
> res at mpMinLonF = 270 ;minimum longitude
>
> ;***********************************************
> ; ----wind vector plot
> ;***********************************************
> vcres = True
> vcres at vcRefAnnoOrthogonalPosF = -1.0 ; move ref vector up
> vcres at vcRefMagnitudeF = 10.0 ; define vector ref mag
> vcres at vcRefLengthF = 0.045 ; define length of vec ref
> vcres at vcGlyphStyle = "CurlyVector" ; turn on curly vectors
> vcres at vcMinDistanceF = 0.017
> vcres at mpFillOn = False ; turn off gray fill
> vcres at mpOutlineBoundarySets = "National" ; turn on country boundaries
> vcres at mpGeophysicalLineColor = "Navy" ; color of cont. outlines
> vcres at mpGeophysicalLineThicknessF = 1.5 ; thickness of outlines
>
>
> ;vcres at gsnLeftString = "DJF High NAO"
> winds_hi = gsn_csm_vector_map_ce(wks,uAvgTime_hi,vAvgTime_hi,vcres)
> ;vcres at gsnLeftString = "DJF Low NAO"
> winds_lo = gsn_csm_vector_map_ce(wks,uAvgTime_lo,vAvgTime_lo,vcres)
> ;vcres at gsnLeftString = "Difference of High - Low"
> winds_diff = gsn_csm_vector_map_ce(wks, diff_u, diff_v,vcres)
> ;************************************************
> ;---- geopotential height filled contour plot
> ;***********************************************
> zfres = res
> zfres at cnFillOn = True
> ;zfres at cnLevelSelectionMode = "ExplicitLevels
> ;zfres at cnLevels = ispan(-20,90,5)
> zfres at lbLabelFontHeightF = 0.015
> zfres at lbOrientation = "Vertical"
> zfres at pmLabelBarOrthogonalPosF = -0.005
>
>
> contour_zf_hi = gsn_csm_contour(wks,zAvgTime_hi,zfres)
> contour_zf_lo = gsn_csm_contour(wks,zAvgTime_lo,zfres)
> contour_zf_diff = gsn_csm_contour(wks,diff_z,zfres)
>
> plot(0) = gsn_csm_contour(wks,zAvgTime_hi,zfres)
> plot(1) = gsn_csm_contour(wks,zAvgTime_lo,zfres)
> plot(2) = gsn_csm_contour(wks,diff_z,zfres)
>
> overlay(contour_zf_hi,winds_hi)
> overlay(contour_zf_lo,winds_lo)
> overlay(contour_zf_diff,winds_diff)
>
>
> ;************************************************
> ; create panel
> ;************************************************
> resP = True ; modify the panel plot
> resP at txString = "NAO 500mb"
> gsn_panel(wks,plot,(/3,1/),resP) ; now draw as one plot;
>
>
> end
>
>
>
>
>
>
>
>
>
>
>
>
>> Mira,
>>
>> Try this:
>>
>>
>> plot(0) = gsn_csm_contour(wks,zAvgTime_hi,zfres)
>> plot(1) = gsn_csm_contour(wks,zAvgTime_lo,zfres)
>> plot(2) = gsn_csm_contour(wks,diff_z,zfres)
>>
>> overlay(contour_zf_hi,winds_hi)
>> overlay(contour_zf_lo,winds_lo)
>> overlay(contour_diff,winds_diff)
>>
>> I don't believe overlay returns anything and you are trying to save it as
>> a
>> graphic. The code above would replace the lines in your script.
>>
>> Hope that helps,
>>
>> -Alex
>>
>>
>>
>>
>> On Sunday, July 26, 2015, <mberdahl at envsci.rutgers.edu
>> <javascript:_e(%7B%7D,'cvml','mberdahl at envsci.rutgers.edu <mailto:mberdahl at envsci.rutgers.edu>');>> wrote:
>>
>>> Hi all,
>>>
>>> I'm trying to adapt a working code I had that was plotting 3 panel plots
>>> of wind vectors in a certain region. Now, all I want to do is add an
>>> overlay of geopotential height (filled contours) underneath each of
>>> these
>>> plots, but the way I'm approaching it is wrong. I've tried simplifying
>>> it
>>> to do this without panels, just a single plot, but I still don't get
>>> what
>>> I need. The error message I am getting now is:
>>>
>>> fatal:syntax error: line 153 in file plotWinds_z_NAO_level.ncl before
>>> or
>>> near overlay
>>> plot(0)= overlay
>>> ---------------^
>>>
>>>
>>> My code is copied below. Any suggestions are greatly appreciated.
>>> Thanks,
>>> Mira
>>>
>>> ;*****************************************
>>> ; ; plot average winds overlayed on geopotental filled contours for the
>>> years that are extreme
>>> ; the years are found with my matlab script findExtremeYrs.m
>>> ; this particlar script does the years that correspond to high
>>> correlation
>>> with DJF, SE precip and NAO
>>> ;*****************************************
>>> ; the original data goes from 1948 January to April 2015.
>>> ; I will cut it out to 2014 December, so we have something divisible by
>>> 12
>>> so we can do seasonal averages...
>>>
>>>
>>> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
>>> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
>>> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
>>> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
>>> ;************************************************
>>> begin
>>> ;************************************************
>>> ; read in netCDF file s
>>> ;************************************************
>>> a = addfile("uwnd.mon.mean.alllevels.nc","r") ; u winds
>>> b = addfile("vwnd.mon.mean.alllevels.nc","r") ; v winds
>>>
>>> c = addfile("../Geopotential/hgt.mon.mean.nc","r") ; geopotential
>>> heights.
>>>
>>> ;************************************************
>>> ; read in zonal [u] and meridional [v] winds (July)
>>> ;************************************************
>>>
>>> u = a->uwnd(0:803,{500},{45:90},{270:357.5})
>>> v = b->vwnd(0:803,{500},{45:90},{270:357.5}) ; Get u, v, time (1),level
>>> (1000hpa),latitude(-90:90) and longitude(0:360) data.
>>> z = c->hgt(0:803,{500},{45:90},{270:357.5}) ; get geopotenial
>>> heights...
>>>
>>> printVarSummary(u)
>>> printVarSummary(v)
>>> printVarSummary(z)
>>>
>>> ; Calculate the seasonal averages.
>>> uDJF = month_to_season(u, "DJF")
>>> vDJF = month_to_season(v, "DJF")
>>> zDJF = month_to_season(z, "DJF")
>>>
>>> printVarSummary(uDJF)
>>> printVarSummary(vDJF)
>>> printVarSummary(zDJF)
>>>
>>> ; from the matlab script i wrote: findExtremeYrs, i pulled out the
>>> extreme
>>> years (> or < 1std) that i want to average and plot here.
>>>
>>> ; for ans = 4 (NAO)
>>> ; yearList_hi = 1973 1975 1983 1989 1995
>>> 2000 2007 2012
>>> ; yearList_lo = 1963 1964 1965 1969 1977
>>> 1979 1996 1997 2010 2011
>>>
>>>
>>> ; this data starts at 1948 (this is index 0), so 1953=5, 1963=15 etc.
>>>
>>> uDJF_NAO_hi = uDJF((/25,27,35,41,47,52,59,64/),:,:)
>>> uDJF_NAO_lo = uDJF((/15,16,17,21,29,31,48,49,62,63/),:,:)
>>>
>>> vDJF_NAO_hi = vDJF((/25,27,35,41,47,52,59,64/),:,:)
>>> vDJF_NAO_lo = vDJF((/15,16,17,21,29,31,48,49,62,63/),:,:)
>>>
>>> zDJF_NAO_hi = zDJF((/25,27,35,41,47,52,59,64/),:,:)
>>> zDJF_NAO_lo = zDJF((/15,16,17,21,29,31,48,49,62,63/),:,:)
>>>
>>>
>>> uAvgTime_hi = dim_avg_n_Wrap(uDJF_NAO_hi,0)
>>> uAvgTime_lo = dim_avg_n_Wrap(uDJF_NAO_lo,0)
>>>
>>> printVarSummary(uAvgTime_hi)
>>> printVarSummary(uAvgTime_lo)
>>>
>>> vAvgTime_hi = dim_avg_n_Wrap(vDJF_NAO_hi,0)
>>> vAvgTime_lo = dim_avg_n_Wrap(vDJF_NAO_lo,0)
>>>
>>> printVarSummary(vAvgTime_hi)
>>> printVarSummary(vAvgTime_lo)
>>>
>>> zAvgTime_hi = dim_avg_n_Wrap(zDJF_NAO_hi,0)
>>> zAvgTime_lo = dim_avg_n_Wrap(zDJF_NAO_lo,0)
>>>
>>> printVarSummary(zAvgTime_hi)
>>> printVarSummary(zAvgTime_lo)
>>>
>>> ; dirty way to copy metadata over first.
>>> diff_u = uAvgTime_hi;
>>> diff_v = vAvgTime_hi;
>>> diff_z = zAvgTime_hi;
>>>
>>> diff_u = uAvgTime_hi - uAvgTime_lo
>>> diff_v = vAvgTime_hi - vAvgTime_lo
>>> diff_z = zAvgTime_hi - zAvgTime_lo
>>>
>>> printVarSummary(diff_u)
>>> printVarSummary(diff_v)
>>> printVarSummary(diff_z)
>>>
>>> ;************************************************
>>> ; create plot
>>> ;************************************************
>>> wks = gsn_open_wks("ps","Panel_NAO_z_500") ; open a ps file
>>> plot = new(3,graphic) ; create a plot array
>>>
>>> ;---- set common resources for all plots
>>> res = True
>>> res at gsnDraw = False ; dont draw
>>> res at gsnFrame = False ; dont advance frame
>>> res at cnInfoLabelOn = False ; trn off cn info label
>>> res = True ; plot mods desired
>>> res at gsnAddCyclic = False ; has to do with
>>> wrapping
>>> the longitude at 0/360
>>> ;************************************************
>>> ; Choose a subregion
>>> ;************************************************
>>> res at mpMaxLatF = 90 ;maximum latitude
>>> res at mpMinLatF = 45 ;minimum latitude
>>> res at mpMaxLonF = 357.5 ; ;maximum longitude
>>> res at mpMinLonF = 270 ;minimum longitude
>>>
>>> ;***********************************************
>>> ; ----wind vector plot
>>> ;***********************************************
>>> vcres = True
>>> vcres at vcRefAnnoOrthogonalPosF = -1.0 ; move ref vector up
>>> vcres at vcRefMagnitudeF = 10.0 ; define vector ref mag
>>> vcres at vcRefLengthF = 0.045 ; define length of vec
>>> ref
>>> vcres at vcGlyphStyle = "CurlyVector" ; turn on curly vectors
>>> vcres at vcMinDistanceF = 0.017
>>> vcres at mpFillOn = False ; turn off gray fill
>>> vcres at mpOutlineBoundarySets = "National" ; turn on country
>>> boundaries
>>> vcres at mpGeophysicalLineColor = "Navy" ; color of cont.
>>> outlines
>>> vcres at mpGeophysicalLineThicknessF = 1.5 ; thickness of outlines
>>>
>>>
>>> ;vcres at gsnLeftString = "DJF High NAO"
>>> winds_hi = gsn_csm_vector_map_ce(wks,uAvgTime_hi,vAvgTime_hi,vcres)
>>> ;vcres at gsnLeftString = "DJF Low NAO"
>>> winds_lo = gsn_csm_vector_map_ce(wks,uAvgTime_lo,vAvgTime_lo,vcres)
>>> ;vcres at gsnLeftString = "Difference of High - Low"
>>> winds_diff = gsn_csm_vector_map_ce(wks, diff_u, diff_v,vcres)
>>> ;************************************************
>>> ;---- geopotential height filled contour plot
>>> ;***********************************************
>>> zfres = res
>>> zfres at cnFillOn = True
>>> ;zfres at cnLevelSelectionMode = "ExplicitLevels
>>> ;zfres at cnLevels = ispan(-20,90,5)
>>> zfres at lbLabelFontHeightF = 0.015
>>> zfres at lbOrientation = "Vertical"
>>> zfres at pmLabelBarOrthogonalPosF = -0.005
>>>
>>>
>>> contour_zf_hi = gsn_csm_contour(wks,zAvgTime_hi,zfres)
>>> contour_zf_lo = gsn_csm_contour(wks,zAvgTime_lo,zfres)
>>> contour_zf_diff = gsn_csm_contour(wks,diff_z,zfres)
>>>
>>> plot(0) = overlay(contour_zf_hi,winds_hi)
>>> plot(1) = overlay(contour_zf_lo,winds_lo)
>>> plot(2) = overlay(contour_diff,winds_diff)
>>>
>>>
>>> ;************************************************
>>> ; create panel
>>> ;************************************************
>>> resP = True ; modify the panel plot
>>> resP at txString = "NAO 500mb"
>>> gsn_panel(wks,plot,(/3,1/),resP) ; now draw as one plot;
>>>
>>>
>>> end
>>>
>>>
>>> _______________________________________________
>>> ncl-talk mailing list
>>> ncl-talk at ucar.edu
>>> List instructions, subscriber options, unsubscribe:
>>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>>
>>
> <Panel_NAO_z_500.pdf>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20150727/5b347ca7/attachment.html
More information about the ncl-talk
mailing list