[ncl-talk] compute frequency power density power spectra for selected regions
Sri nandini
bax8609 at uni-hamburg.de
Mon Mar 8 04:22:19 MST 2021
Hello dear ncl-users,
Would anyone know how to remove this red-noise spectrum so i can show
the peaks of high variance for my plot?
Best
sri
On 03.03.21 15:03, Dennis Shea wrote:
> In my opinion, the spectrum reveals a classic red-noise spectrum. Only
> low frequency variance is present.
>
> On Tue, Mar 2, 2021 at 5:53 AM Sri nandini via ncl-talk
> <ncl-talk at mailman.ucar.edu <mailto:ncl-talk at mailman.ucar.edu>> wrote:
>
> Hello dear ncl-users,
>
> I have a question regarding computing frequency power spectra density
> for my selected regions, e.g. Northwestern Atlantic SSH anomalies.
>
> 1) I calculated the frequency spectra of SSH anomalies by (1)
> extracting
> region and subtracting the ensemble mean, (2) temporal removal of
> mean
> and detrending was applied, as well as smoothing (periodogram
> estimates
> using modified Daniel) and percent tapering of the anomalies, (3) the
> function for spectral analysis returns the degrees of freedom
> -spectra
> for each ensemble, so i used a loop over all ensemble members to
> get the
> ensemble mean spectra (4) the variables returned are variance/(unit
> frequency interval) (y-axis) and frequency (cycles/months) for
> x-axis==see attached plot for power Spectra for historical
> Northwestern
> Atlantic SSH anomalies. I don't have a reference for this figure to
> compare with but it doesn't look meaningful, but do you think this is
> correct?
>
> Would be grateful, below is my code
>
> Sri
>
>
> f = addfile ("Hist_monthly_corrected_pi_trend.nc
> <http://trend.nc>", "r")
> hist_anom = f->hist_trend(1632:1871,:,{-30:60},{0:-40})
> hist_anom_NY1=dim_avg_n_Wrap(hist_anom,(/2,3/))
> printVarSummary(hist_anom_NY1)
> hist_anom_NY1=hist_anom_NY1*100
> dimx = dimsizes(hist_anom_NY1)
> ntim = dimx(0) ; 240
> nens = dimx(1) ;100
> hist_ensavg=dim_avg_n(hist_anom_NY1,1)
> printVarSummary(hist_ensavg)
> hist_anom1 = hist_anom_NY1
> do n=0,nens-1
> hist_anom1(:,n) = hist_anom_NY1(:,n) - hist_ensavg
> end do
>
> nLen = ntim/2
> spec = new ( nLen, typeof(hist_anom1))
> spec = 0.0
> z1 = 0.0
> xtsVar= 0.0
> printVarSummary(nLen)
> printVarSummary(spec)
> printVarSummary(hist_anom1)
> ;-------------------------------------;
> ;set function agurments
> ;-------------------------------------;
>
> iopt = 1 ; detrending opt: 0=>remove
> mean
> 1=>remove mean and detrend
> jave = 7 ; Average 7 periodogram
> estimates
> using modified Daniel
> pct = 0.1 ; percent tapered: (0.0 <=
> pct <=
> 1.0) 0.10 common.
> ;-------------------------------------;
> ; calculate mean spectrum:specx_anal returns the degrees of
> freedom as
> a scalar.
> ; calculate confidence interval [here 5 and 95%] and lag1 auto cor.
> ; return 4 curves to be plotted
> ; 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
> printVarSummary(z1)
> 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
> printVarSummary(spec)
> printVarSummary(z1)
> printVarSummary(xvari)
> printVarSummary(sdof)
> ;-----------------------------------------------------------------------------;
> ; create a new variable
> ;-----------------------------------------------------------------------------;
>
> ss = (/ sdof /)
> ss at dof = 2*nens ; # degrees freedom
> ss at bw = sdof at bw ; representing the
> spectral
> band width.
> ss at spcx = spec ;One-dimensional array of
> length N/2.units are variance/(unit frequency interval).
> ss at xvari = xvari ; representing the
> variance
> of the x series on input.
> ss at frq = sdof at frq ; A one-dimensional
> array of
> length N/2 representing frequency (cycles/time).
> ss at xlag1 = (exp(2*z1)-1)/(exp(2*z1)+1) ; match specx_anal
> printVarSummary(ss)
> splt = specx_ci(ss, 0.05, 0.95) ; calc confidence
> interval
> printVarSummary(splt)
> ;-----------------------------------------------------------------------------;
> ; plotting parameters
> ;-----------------------------------------------------------------------------;
> wks = gsn_open_wks("png","PSD_wNA_Monthly_Hist")
>
> res = True ;
> plot
> mods desired
> res at gsnMaximize = True ; blow up plot
> res at gsnDraw = False
> res at gsnFrame = False
>
> res at vpXF = 0.1
> res at vpYF = 0.7
> res at vpWidthF = 0.7
> res at vpHeightF = 0.25
>
> res at trXMinF = 0
> res at trXMaxF = .50
> res at trYMinF = 20
> res at trYMaxF = 0
>
> res at gsnStringFontHeightF = 0.018
> res at gsnRightString = " "
> res at tiXAxisString = "Frequency (cycles/month)" ;the
> fraction (or
> use Period(months))
> res at tiXAxisFontHeightF = 0.018
>
> res at tiYAxisString = "PSD (cm^2/month^-1)" ;"Variance/frq_interval
> res at tiYAxisFontHeightF = 0.018
>
> res at tmXBLabelFontHeightF = 0.016
> res at tmYLLabelFontHeightF = 0.016
>
> ;-------------------------------------;
> ; 4-2. Timeseries plots
> ;-------------------------------------;
>
> colors = (/2, 71, 48/)
> colors_em = (/15, 84, 60/)
>
> res at gsnLeftString = "a) Historical: Northwest Atlantic, SSHa
> Frequency Spectra"
> res at xyLineThicknesses = (/2.,1.,1.,1./) ; Define line
> thicknesses
> res at xyDashPatterns = (/0,0,1,1/) ; Dash patterns
>
> res at xyLineColor = colors(0)
> plot = gsn_csm_xy(wks,ss at frq, ss at spcx,res)
>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at mailman.ucar.edu <mailto:ncl-talk at mailman.ucar.edu>
> List instructions, subscriber options, unsubscribe:
> https://mailman.ucar.edu/mailman/listinfo/ncl-talk
> <https://mailman.ucar.edu/mailman/listinfo/ncl-talk>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://mailman.ucar.edu/pipermail/ncl-talk/attachments/20210308/dc9c52f7/attachment.html>
More information about the ncl-talk
mailing list