[ncl-talk] calculate the ensemble mean power spectrum

[Dennis Shea]

This type of question was recently answered on ncl-talk.

To avoid repeated emails, I extracted the pertinent section.
PLEASE look-at and think-about what is being done!
Add your own print statements to examine any issues.

This has not been tested on your data. There may be errors.
It is your responsibility to fix them.
+++++++++++++++++++++++++++++++++++++++++++++++++
hist_anom1=hist_anom(:,:,{-7},{160})     ;(ntim, ens) ; (240,100)
=====
  nLen = ntim/2
  spec  = new ( nLen, typeof(hist_anom1))
  spec  = 0.0
  z1    = 0.0
  xtsVar= 0.0

  iopt  = 0
  jave  = 1                           ; periodigram for each segment
  pct   = 0.1
                                           ; loop over  every ensemble
member
  do ns=0,nens-1                 ; spectra for each ensemble member
     sdof   = specx_anal(hist_anom1(:.ns),iopt,jave,pct)
     spec   = spec + sdof at spcx      ; sum spectra
     xtsVar = xtsVar + sdof at xvari   ; sum ensemble member variances
                                                      ; sum z-traansformed
lag-1 correlations
     z1     = z1   + 0.5*log((1+sdof at xlag1)/(1-sdof at xlag1)) ; Fischer
z-transform
  end do

  spec  = spec/nens                      ; average spectra ==> ensemble
mean spectra
  z1    = z1/nens                           ; ensemble mean z=transformed
lag1
  xvari = xtsVar/nens                    ; mean variance of ensemble members

                                                    ; create a new variable
  ss           = (/ sdof /)
  ss at dof   = 2*nens                       ; # degrees freedom
  ss at bw    = sdof at bw
  ss at spcx  = spec
  ss at xvari = xvari
  ss at frq     = sdof at frq
  ss at xlag1 = (exp(2*z1)-1)/(exp(2*z1)+1)  ; match specx_anal

  printVarSummary(ss)

;************************************************
; plotting parameters
;************************************************
  pltDir  = "./"
  pltName = "ens_mean"
  pltPath = pltDir+pltName
  pltType = "png"

  wks  = gsn_open_wks(pltType,pltPath)

  r  = True                                                  ; plot mods
desired
  r at gsnMaximize   = True                        ; blow up plot
  r at tiMainString  = "???"                          ; title
  r at tiXAxisString = "Frequency (cycles/day)"    ; xaxis
  r at tiYAxisString = "Variance"                  ; yaxis
  splt = specx_ci(ss, 0.05, 0.95)                ; calc confidence interval
  r at xyLineThicknesses   = (/2.,1.,1.,1./)   ; Define line thicknesses
  r at xyDashPatterns      = (/0,0,1,1/)          ; Dash patterns
  plot = gsn_csm_xy(wks,ss at frq, splt,r)



On Tue, Mar 3, 2020 at 7:08 AM Sri.durgesh Nandini-Weiss via ncl-talk <
ncl-talk at ucar.edu> wrote:

> Dear ncl users,
>
> I am trying to perform power spectral analysis of a time series of
> ensemble data (time*ens==>20years*100ensemble to compute over the year
> and each ensemble member and than calculate the ensemble mean of the
> spectra.
>
> The best way would be to use specx_anal() and write it as a loop but
> this is where i am failing.
>
> I would later on like to know how to add the Markov red noise and the
> upper and lower confidence levels?
>
> would someone advice me on my script below?
>
>
> begin
>
>     f     = addfile ("anom_hist_zo_slr.nc", "r")
>     hist_anom    = f->hist_anom
>     printVarSummary(hist_anom)
>
>     f1     = addfile ("anom_rcp45_zo_slr.nc", "r")
>     rcp45_anom   = f1->rcp45_anom
>     printVarSummary(rcp45_anom)
>
>     dimx = dimsizes(hist_anom)
>     ntim = dimx(0)           ; 240
>     nens = dimx(1)           ; 100
>     nlat = dimx(2)           ; 45
>     nlon = dimx(3)           ; 90
>     lat1=nlat
>     lon1=nlon
>
>     nmos = 12
>     nyrs   = ntim/nmos       ; 20
>     printVarSummary(nyrs)
>
>
> ;==================================================================
> ; print ("Data read in - start calculations")
> ; Region domain lat and lon, ensemble average
> ;==================================================================
>
>    hist_anom1=hist_anom(:,:,{-7},{160}) ;(ntim, ens)
>    printVarSummary(hist_anom1)
>    aveX=dim_avg_n_Wrap(hist_anom1,(/1/))
>    printVarSummary(aveX)     ; [time]
>
>    rcp45_anom1=rcp45_anom(:,:,{-7},{160})                ;(ntim, ens)
>    printVarSummary(rcp45_anom1)
>    aveY=dim_avg_n_Wrap(rcp45_anom1,(/1/))
>    printVarSummary(aveY)      ; [time]
>
>
>    d   = 0                      ; detrending opt: 0=>remove mean
> 1=>remove mean + detrend
>    sm  = 7                    ; smooth periodogram:: should be at least
> 3 and odd
>                                     ; Average 7 periodogram estimates
> using modified Daniell
>    pct = 0.10                 ; percent taper: (0.0 <= pct <= 1.0) 0.10
> common.
>
> ;==================================================================
>    ;i can either use ndtooned and carry out specx_anal on this or make a
> loop.
>
>    aveX = ndtooned(hist_anom1)
>    printVarSummary(aveX)
>
>    aveY = ndtooned(rcp45_anom1)
>    printVarSummary(aveY)
>
>    spec = specx_anal(aveX,d,sm,pct)
>    printVarSummary(spec)
>
>    spec1 = specx_anal(aveY,d,sm,pct)
>    printVarSummary(spec1)
>
>
> Thanx
>
> Sri
>
>
> --
> Dr. Sri Nandini-Weiß
> Research Scientist
>
> Universität Hamburg
> Center for Earth System Research and Sustainability (CEN)
> Cluster of Excellence 'Climate, Climatic Change, and Society' (CLICCS)
>
> Bundesstrasse 53, 20146 Hamburg
> Tel: +49 (0) 40 42838 7472
>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at ucar.edu
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20200303/4124ef42/attachment.html>