[ncl-talk] problem with wavelet power calculation

[Yan Yan]

Hi there,



I ran into some problems lately when using wavelet analysis in NCL. Here they are.



I was using wavelet to analyze an EOF PC1, which is in monthly resolution. Two problems stuck out.

One is when I dtrend it, the result time series became constant, -8.520693e+34. The following wavelet gave me a time series with power being all "inf". I couldn't plot the power spectrum. This was obvious the input problem, so I changed the input into the PC1 without detrending, which gave me the same problem, the power of the wavelet time series still showed "inf".



I suspect there is something wrong with this input time series, but I don't know enough to identify what the problem is, and how to deal with it.

PC1 data is attached. Can somebody point me to the right direction on this?



Also, the following is the code I used, which is basically using NCL example code. Since data is monthly, I tried 1. and 1./12 for dt, both gave me the same "inf" power.



wts = eof_ts(0,:)
; wts = dtrend(wts,False)
  mother  = 0
  param   = 6.0
  dt      = 1.
  s0      = dt
  dj      = 0.25
  jtot    = 1 + floattointeger(((log10(N*dt/s0))/dj)/log10(2.))
  npad    = N
  nadof   = new(2,float)
  noise   = 1
  siglvl  = .05
  isigtest= 1



  w = wavelet(wts,mother,dt,param,s0,dj,jtot,npad,noise,isigtest,siglvl,nadof)
;************************************
; create coodinate arrays for plot
;************************************
  power            = onedtond(w at power,(/jtot,N/))
  power!0          = "period"                        ; Y axis
  power&period     = w at period

  power!1          = "time"                          ; X axis
  power&time       = time

  power at long_name  = "Power Spectrum"
  power at units      = "unit^2"

; compute significance ( >= 1 is significant)
  SIG              = power                            ; transfer meta data
  SIG              = power/conform (power,w at signif,0)
  SIG at long_name    = "Significance"
  SIG at units        = " "

;************************************
; initial resource settings
;************************************
  wks = gsn_open_wks("x11","wavelet")             ; send graphics to PNG file
  plot = new(2,graphic)                           ; create graphical array

  res                     = True                  ; plot mods desired
  res at gsnFrame            = False                 ; don't advance frame yet
  res at gsnDraw             = False                 ; don't draw yet

  res at cnFillOn            = True                  ; turn on color
  res at cnFillPalette       = "BlAqGrYeOrReVi200"   ; set color map
  res at cnFillMode          = "RasterFill"          ; turn on raster mode
  res at cnRasterSmoothingOn = True                  ; turn on raster smoothing
  res at cnLinesOn           = False                 ; turn off contour lines
  res at lbOrientation       = "Vertical"            ; vertical label bar

  res at trYReverse          = True                  ; reverse y-axis
  res at tmLabelAutoStride   = True

  res at vpHeightF           = .4                    ; height and width of plot
  res at vpWidthF            = .85


  res at cnLevelSelectionMode = "ManualLevels"       ; set manual contour levels
  res at cnMinLevelValF       = 0.0                  ; set min contour level
  res at cnMaxLevelValF       = 20.0                 ; set max contour level
  res at cnLevelSpacingF      = 2.0                  ; set contour spacing

; power
  res at gsnLeftString       = "Wavelet Power"
  plot(0) = gsn_csm_contour(wks,power({0:35},0:N-1),res)  

; significance
  res at cnMinLevelValF      = 1.0                ; set min contour level
  res at cnMaxLevelValF      = 4.0                ; set max contour level
  res at cnLevelSpacingF     = 3.0                ; set contour spacing
  res at gsnLeftString       = "Wavelet Significance"
  plot(1) = gsn_csm_contour(wks,SIG({0:35},0:N-1),res)


Any help is very much appreciated!



Best,

Yan
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20200205/ab2ea4d3/attachment.html>
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: eofPC1.txt
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20200205/ab2ea4d3/attachment.txt>