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

From: Hyacinth Nnamchi <hyacinth.1_at_nyahnyahspammersnyahnyah>
Date: Thu Mar 15 2012 - 06:53:30 MDT

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

Received on Thu Mar 15 06:53:42 2012

This archive was generated by hypermail 2.1.8 : Tue Mar 20 2012 - 15:27:15 MDT