[ncl-talk] Issues with period of power spectrum

[Nishtha Agrawal]

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)
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://mailman.ucar.edu/pipermail/ncl-talk/attachments/20201019/c182f68e/attachment.html>