[ncl-talk] Issues with period of power spectrum

[Dennis Shea]

It is not clear to me what:
"obtain a power spectrum of daily rainfall data for 24 years  (for a
specific season having 122 days/year)"

[1] 24 yearly segments of 122 days ?
[2] soi  = in->PAVE
     Are the 'PAVE' input as 122*25 days?
     As always, ncl-talk asks that users include the output from (here)
     printVarSummary(soi)

--
Note: You should  NOT perform a single spectrum on a series of 24*122 days.
          You must perform segment averaging.
          (a) compute the raw spectrum of each 122 day segment and then
average
          (b) generally people do not due spectral analysis on series which
contain many zeros like you may have with daily precipitation.




On Mon, Oct 19, 2020 at 12:00 AM Nishtha Agrawal via ncl-talk <
ncl-talk at mailman.ucar.edu> wrote:

> Hi,
>
> I am trying to obtain a power spectrum of daily rainfall data for 24
> years  (for a specific season having 122 days/year). I want to obtain the
> spectra with respect to period instead of frequency. For this I modified
> the available code that uses specx_anal and is given below. I am facing a
> strange issue where the values of period start from 60. I don't get any
> value before the 60 day period. I don't understand if I am making some
> mistake in the following code
>
>    fn  = "data.nc" ; define filename
>    in  = addfile(fn,"r")                       ; open netcdf file
>    soi  = in->PAVE                             ; get data
>   d   = 1    ; detrending opt: 0=>remove mean 1=>remove mean + detrend
>   sm  = 21   ; smooth: should be at least 3 and odd
>   pct = 0.10 ; percent taper: (0.0 <= pct <= 1.0) 0.10 common.
>   sdof = specx_anal(soi,d,sm,pct)
>   splt = specx_ci (sdof, 0.10, 0.90)
>    wks  = gsn_open_wks("png","psd_rainfall")             ; send graphics
> to PNG file
>
>    r = True                                      ; plot mods desired
>    r at tiMainString = "RegCM"                ; title
>    r at tiXAxisString = "Frequency (cycles/day)"  ; xaxis
>    r at tiYAxisString = "Variance/frq_interval"     ; yaxis
>    r at xyLineColors        = (/"foreground","red","blue","blue"/)
>    f = sdof at frq/30.5                              ; for daily data
>    r at trYMinF             = 0.00                 ; manually set lower limit
>    r at trYMaxF             = 3.0                 ;   "          upper
>    r at gsnFrame            = False                ; do not advance frame
>    plot  = gsn_csm_xy(wks,f, splt,r)
>
>    xf   = (/0.40, 0.40+sdof at bw/)                ; Create band width line
>    ys   = 0.75*max(sdof at spcx)                   ; 75% up Y axis
>    yv   = (/ys,ys/)
>    rpl  = True                                  ; resources for polyline
>    rpl at gsLineThicknessF  = 2                    ; Define line thickness
>    gsn_polyline(wks,plot,xf,yv,rpl)             ; Draw BandWidth
>
>    txres= True                                  ; label BW line
>    txres at txFontHeightF = 0.015                  ; font height
>    txres at txJust        = "CenterLeft"           ; Set lable location
>    gsn_text(wks,plot,"BW",0.41+sdof at bw,ys,txres); Label
>    frame (wks)
>
>    r at gsnFrame          = True                   ;reset to default
>    p   = 1/f                                    ; *highly non-linear*
>    p!0 = "f"
>    p&f = f
>    p at long_name = "period"
>    p at units     = "day"
>
>  ;;print("====")
>  ;;print(f+"   "+p+"  "+splt(0,:) )             ; print
>  ;;print("====")
>
>    r at tiXAxisString = "Period (day)"          ; xaxis
>    r at tiYAxisString = "Variance/freq"           ; yaxis
>
>    ip   = ind(p.le.120)                        ; all indices for "short"
> periods
>    plot = gsn_csm_xy(wks,p(ip), splt(:,ip),r)
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at mailman.ucar.edu
> List instructions, subscriber options, unsubscribe:
> https://mailman.ucar.edu/mailman/listinfo/ncl-talk
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://mailman.ucar.edu/pipermail/ncl-talk/attachments/20201020/9e7f157d/attachment.html>