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,nlon)
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) ; 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 for 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 ; 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
This archive was generated by hypermail 2.1.8 : Tue Mar 20 2012 - 15:27:15 MDT