Re: Changing default X-axis (frequency) ticks for spectrum / cospectrum

Thank alot Dennis, that answers my question.Hyacinth

> Date: Thu, 15 Mar 2012 11:14:40 -0600
> From: shea@ucar.edu
> To: hyacinth.1@hotmail.com
> CC: ncl-talk@ucar.edu
> Subject: Re: Changing default X-axis (frequency) ticks for spectrum / cospectrum
>
> period = 1/f (=frequency)
>
> if your sampling is every year and (say) f=0.2 cycles/year, then
> period = 1/[0.2(cycles/year)*(1 year/12 months)]=12/0.2=60 months
>
> On 3/15/12 6:53 AM, Hyacinth Nnamchi wrote:
> > Hi,
> > I'm using ncl -V 6.0.0 to calculate and plot the spectrum and cospectrum
> > of annual time series. My question is that the frequencies appear in
> > fractions (please see attached scripts), but I'd like to have them in
> > years.
> >
> > Thanks for any suggestions,
> >
> > Hyacinth
> >
> >
> > *********************************************
> >
> > ;*************************************************
> > 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
> >
> > ; ==============================================================
> > ; User defined parameters that specify region of globe and
> > ; ==============================================================
> > ; latS = -50.
> > ; latN = 20.
> > ; lonL = -70.
> > ; lonR = 40.
> >
> > yrStrt = 1950
> > yrLast = 2008
> >
> > season = "JJA" ; choose Dec-Jan-Feb seasonal mean
> >
> > ; ==============================================================
> > ; Open the file: Read only the user specified period
> > ; ==============================================================
> >
> > f = addfile ("/home/nnamchi/work/datasets/sst/sst.mnmean.nc", "r")
> >
> > TIME = f->time
> > YYYY = cd_calendar(TIME,-1)/100 ; entire file
> > iYYYY = ind(YYYY.ge.yrStrt .and. YYYY.le.yrLast)
> > slp = short2flt(f->sst(iYYYY,:,:))
> > printVarSummary(slp)
> >
> > delete(TIME)
> > delete(YYYY)
> > delete(iYYYY) ; this changes size
> >
> > f0 = addfile
> > ("/home/nnamchi/work/datasets/rain/precip.mon.total.v201.nc", "r")
> >
> > TIME = f0->time
> > YYYY = cd_calendar(TIME,-1)/100 ; entire file
> > iYYYY = ind(YYYY.ge.yrStrt .and. YYYY.le.yrLast)
> > slpp = short2flt(f0->precip(iYYYY,:,:))
> > printVarSummary(slpp)
> >
> > delete(TIME)
> > delete(YYYY)
> > delete(iYYYY) ; this changes size
> >
> > ; ==============================================================
> > ; dataset longitudes span 0=>357.5
> > ; Because EOFs of the North Atlantic Oscillation are desired
> > ; use the "lonFlip" (contributed. ncl) to reorder
> > ; longitudes to span -180 to 177.5: facilitate coordinate subscripting
> > ; ==============================================================
> > ;=================================================;
> > ; Remove the annual cycle
> > ;=================================================;
> > xClm = clmMonTLL(slp)
> > printVarSummary(xClm) ; (12,nlat,nlon)
> > slp = calcMonAnomTLL (slp, xClm) ; replace with anonamlies
> > slp@long_name = "ANOMALIES: "+slp@long_name
> >
> >
> > ;===
> > xClm1 = clmMonTLL(slpp)
> > printVarSummary(xClm1) ; (12,nlat,n lon)
> > slpp = calcMonAnomTLL (slpp, xClm1) ; replace with anonamlies
> > slpp@long_name = "ANOMALIES: "+slpp@long_name
> >
> >
> > ;=================================================;
> >
> > slp = lonFlip( slp )
> > printVarSummary(slp) ; note the longitude coord
> > ;==
> > slpp = lonFlip(slpp)
> > printVarSummary(slpp) ; note the longitude coord
> > ;;==
> > ; ==============================================================
> > ; compute desired global seasonal mean: month_to_season (contributed.ncl)
> > ; ==============================================================
> > SLP = month_to_season (slp, season)
> > nyrs = dimsizes(SLP&time)
> > printVarSummary(SLP)
> >
> > SLPP = month_to_season (slpp, season)
> > nyrs = dimsizes(SLPP&time)
> > printVarSummary(SLPP)
> >
> > ;=====================================================================
> > ; Calculate the domain averaged time series for spectrum and cospectrum
> > analysis
> > ; First SST indices and the precipitation indices
> > ;==================================================================
> >
> > ; SST Indices
> > nep = wgt_areaave_Wrap(SLP(time |:, {lon|-20:10}, {lat | -15:0}),1.0,
> > 1.0, 0)
> > swp = wgt_areaave_Wrap(SLP(time |:, {lon|-40:-10}, {lat | -40:-25}),1.0,
> > 1.0, 0)
> > atl3 = wgt_areaave_Wrap(SLP(time |:, {lon|-20:0}, {lat | -3:3}),1.0, 1.0, 0)
> > saod = wgt_areaave_Wrap(SLP(time |:, {lon|-20:10}, {lat | -15:0}),1.0,
> > 1.0, 0)-wgt_areaave_Wrap(SLP(time |:, {lon|-40:-10}, {lat |
> > -40:-25}),1.0, 1.0, 0)
> > ;Precipitation Indices
> >
> > gco = wgt_areaave_Wrap(SLPP(time |:, {lon|-15:20}, {lat | 2:12}),1.0,
> > 1.0, 0)
> > sah = wgt_areaave_Wrap(SLPP(time |:, {lon|-20:20}, {lat | 12:20}),1.0,
> > 1.0, 0)
> > waf = wgt_areaave_Wrap(SLPP(time |:, {lon|-20:20}, {lat | 2:20}),1.0,
> > 1.0, 0)
> >
> > ;Normalize each time series by its standard deviation
> > nep = (nep/stddev(nep))
> > swp = (swp/stddev(swp))
> > saod = (saod/stddev(saod))
> > atl3 = (atl3/stddev(atl3))
> >
> > gco = (gco/stddev(gco))
> > sah = (sah/stddev(sah) )
> > waf = (waf/stddev(waf))
> > ;************************************************
> > ; set function arguments for..
> > ; And calculate spectrum and co-specturm
> > ;************************************************
> > ; detrending opt: 0=>remove mean 1=>remove mean and detrend
> > d = 0
> > ; smoothing periodogram: (0 <= sm <= ??.) should be at least 3 and odd
> > sm = 3
> > ; percent tapered: (0.0 <= pct <= 1.0) 0.10 common.
> > pct = 0.10
> > ;************************************************
> > ; calculate spectrum
> > ;************************************************
> >
> > spec = specxy_anal(saod,gco,d,sm,pct)
> > ;************************************************
> > ; plotting parameters that remain constant
> > ;************************************************
> > plot=new(3,graphic) &nbs p; ; creates a graphic array
> >
> > wks = gsn_open_wks("pdf","plot_spec") ; Opens a ps file
> >
> > res = True ; no plot mods desired
> > res@tiXAxisString = "Frequency (cycles/year)" ; xaxis
> > res@gsnFrame = False ; required for panel plots
> > res@gsnDraw = False ; required f or panel plots
> > res@trYMinF = 0.0 ; min value on y-axis
> > res@trYMaxF = 10.0 ; max value on y-axis
> > res@trXMinF = 0.0
> > ;************************************************
> > ; create plot of SAOD spectrum
> > ;************************************************
> > res@tiYAxisString = "SAOD var" ; yaxis
> > plot(0)=gsn_csm_xy(wks,spec@frq,spec@spcx,res); create plot
> > ;***********************************************
> > ; create plot of GCO spectrum
> > ;************************************************
> > res@tiYAxisString = "GCO var" ; ; yaxis
> > plot(1)=gsn_csm_xy(wks,spec@frq,spec@spcy,res); create plot
> >
> > ;***********************************************
> > ; create plot of SAOD/GCO co-spectrum
> > ;************************************************
> > res@tiYAxisString = "GCO var" ; yaxis
> > plot(2)=gsn_csm_xy(wks,spec@frq,spec@cospc,res); create plot
> >
> >
> > ;***********************************************
> > ; Panel the plots
> > ;************************************************
> > res_P = True ; panel mods desired
> > res_P@gsnMaximize = True &nb sp; ; blow up plot
> > gsn_panel(wks,plot,(/3,1/),res_P) ; create panel plot
> >
> > end
> >
> >
> > _______________________________________________
> > ncl-talk mailing list
> > List instructions, subscriber options, unsubscribe:
> > http://mailman.ucar.edu/mailman/listinfo/ncl-talk
                                               

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Fri Mar 16 08:47:26 2012