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 Thu Mar 15 11:14:53 2012