关于NCL小波分析的求助帖

Yoger

最近在学用NCL做小波分析，用的是老师给我发的单站逐日降水资料（文后有图），运行之后总是提示：warning:ContourPlotInitialize: scalar field is constant; no contour lines will appear; use cnConstFEnableFill to enable fill，而且出现这个warning后程序不再继续运行也不结束，但是资料并不是常数啊，请求各位前辈指导一下。ncl code如下：
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
begin

  fili    = "58238"
  f       = asciiread(fili,(/22158,6/),"float")
  x       = f(:,5)*0.1        ;站点降水
        printVarSummary(x)
  ihp     = 2                             ; band pass
  sigma   = 1.0                           ; Lanczos sigma

  nWgt    = 2001                           ; loose 100 each end                           
  fca     = 1./90.                       ; start freq
  fcb     = 1./10.                        ; last  freq
  wgt     = filwgts_lanczos (nWgt, ihp, fca, fcb, sigma )
  printVarSummary(wgt)
  xBPF    = wgt_runave ( x, wgt, 0 )      ; 10-90 day
  printVarSummary(xBPF)

  time=ispan(1,22158,1)
  N=dimsizes(time)
  ;print(N)
;************************************
; compute wavelet
;************************************
  mother  = 0
  param   = 6.0
  dt      = 0.25    ;timesteps in units of years
  s0      = dt
  dj      = 0.25
  jtot    = 1+floattointeger(((log10(N*dt/s0))/dj)/log10(2.))
  npad    = N
  nadof   = 0
  noise   = 1
  siglvl  = .05
  isigtest= 0

  w = wavelet(xBPF,mother,dt,param,s0,dj,jtot,npad,noise,isigtest,siglvl,nadof)
  printVarSummary(w)
;************************************
; create coodinate arrays for plot
;************************************
  power            = onedtond(w@power,(/jtot,N/))
  power!0          = "period"                        ; Y axis
  power&period     = w@period                      ; convert period to units of years

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

  power@long_name  = "Power Spectrum"
  power@units      = "1/unit-freq"

; compute significance ( >= 1 is significant)
  SIG              = power                            ; transfer meta data
  SIG              = power/conform (power,w@signif,0)
  SIG@long_name    = "Significance"
  SIG@units        = " "
;********************************************************************************
; initial resource settings
;********************************************************************************
  wks = gsn_open_wks("png","wavelet")             ; send graphics to PNG file

  res                     = True                  ; plot mods desired
  res@gsnDraw             = False                 ; Do not draw plot
  res@gsnFrame            = False                 ; Do not advance frome
  res@cnFillOn            = True                  ; turn on color
  res@cnFillPalette       = "BlAqGrYeOrReVi200"   ; set color map
  res@cnFillMode          = "RasterFill"          ; turn on raster mode
  res@cnRasterSmoothingOn = True                  ; turn on raster smoothing
  res@cnLinesOn           = False                 ; turn off contour lines
  res@cnLineLabelsOn      = False
  res@cnInfoLabelOn       = False
  res@trYReverse          = True                  ; reverse y-axis
  res@tmYLMode = "Explicit"
  res@tmYLValues = (/1,2,4,8,16,32,64,128/)
  res@tmYLLabels = (/"1","2","4","8","16","32","64","128"/)
  res@tmLabelAutoStride   = True
  res@vpHeightF           = .4                    ;
  res@vpWidthF            = .7
  res@cnLevelSelectionMode = "ExplicitLevels"       ; set manual contour levels
  res@cnLevels = (/0.5,1.,2.,4./)
  res@gsnRightString       = "Wavelet Power"
  res@gsnLeftString       = "58238 Station Precipitation"
  res@tiYAxisString       = "Period"

  res2 = True                            ; res2 probability plots
  res2@trYReverse          = True
  res2@tmYLMode = "Explicit"
  res2@tmYLValues = (/1,2,4,8,16,32,64,128/)
  res2@tmYLLabels = (/"1","2","4","8","16","32","64","128"/)
  res2@gsnDraw             = False       ; Do not draw plot
  res2@gsnFrame            = False       ; Do not advance frome
  res2@cnLevelSelectionMode = "ManualLevels" ; set manual contour levels
  res2@cnMinLevelValF      = 0.00        ; set min contour level
  res2@cnMaxLevelValF      = 2.00        ; set max contour level
  res2@cnLevelSpacingF     = 1.00        ; set contour spacing
  res2@cnInfoLabelOn       = False
  res2@cnLinesOn           = False       ; do not draw contour lines
  res2@cnLineLabelsOn      = False       ; do not draw contour labels
  res2@cnFillScaleF        = 0.5         ; add extra density
  res2@gsnLeftString = ""
  res2@gsnRightString = ""

  plot = new(2,graphic)
  plot(0) = gsn_csm_contour(wks,power,res)
  iplot = gsn_csm_contour(wks,SIG,res2)
  opt   = True
  opt@gsnShadeFillType = "pattern"
  opt@gsnShadeHigh     = 17
  iplot = gsn_contour_shade(iplot,0, 0.8, opt)
  overlay(plot(0),iplot)                                ; overlay probability plot onto power plot

  scale = w@scale
  Cdelta = w@cdelta
  powernorm = power
  powernorm = power/conform(power,scale,0)
  scaleavg = dj*dt/Cdelta*dim_sum_Wrap(powernorm(time|:,{period|2.:8.}))


  gws = w@gws
  resl = True
  resl@gsnFrame = False
  resl@gsnDraw = False
  resl@trYAxisType = "LogAxis"
  resl@trYReverse          = True                  ; reverse y-axis
  resl@tmYLMode = "Explicit"
  resl@tmYLValues = (/1,2,4,8,16,32,64,128/)
  resl@tmYLLabels = (/"1","2","4","8","16","32","64","128"/)
  plotg = gsn_csm_xy(wks,gws,power&period,resl)

  plotc = gsn_attach_plots(plot(0),plotg,res,resl)

  ress = True
  ress@xyDashPatterns = (/0,1/)
  ress@pmLegendDisplayMode = "Always"
  ress@pmLegendOrthogonalPosF = -1.1
  ress@pmLegendParallelPosF   =  0.2
  ress@pmLegendWidthF           =  0.15
  ress@pmLegendHeightF          =  0.1
  ress@lgPerimOn              = False                  ; turn off box around
  ress@gsnFrame = False
  ress@gsnDraw = False
  ress@vpHeightF           = .4
  ress@vpWidthF            = .7
  ress@xyExplicitLegendLabels = (/"2-8 yr"/)
  ress@tiYAxisString = "variance"
  plot(1) = gsn_csm_xy(wks,power&time,scaleavg,ress)

  pres = True
  pres@gsnMaximize = True
  pres@gsnPaperOrientation = "portrait"
  gsn_panel(wks,plot,(/2,1/),pres)
end

Yoger

oyy. 发表于 2020-3-5 13:18
我在画小波的时候，也遇到这个warning了，但是可以正常出图。所以都是忽略这个warning。

谢谢前辈，已经解决了，过滤后的序列两端是缺测值，剔除就能出图了！

Yoger

顶一顶，有没有前辈能指导一下啊，感谢！{:eb303:}{:eb303:}

oyy.

我在画小波的时候，也遇到这个warning了，但是可以正常出图。所以都是忽略这个warning。

daizhangqi

Yoger 发表于 2020-3-8 10:53
谢谢前辈，已经解决了，过滤后的序列两端是缺测值，剔除就能出图了！

请问如何剔除缺失值呢

EXODUS

请问dt为什么是0.25呢

luoluoluo

oyy. 发表于 2020-3-5 13:18
我在画小波的时候，也遇到这个warning了，但是可以正常出图。所以都是忽略这个warning。

请问一下运行小波分析的程序大概多久时间呢

oyy.

luoluoluo 发表于 2021-4-3 15:55
请问一下运行小波分析的程序大概多久时间呢

得看你电脑的运行内存啊

luoluoluo

oyy. 发表于 2021-4-3 20:51
得看你电脑的运行内存啊

如果是在服务器上呢，因为这个程序运行了几个小时还没结束，我也不知道是正常情况还是程序出了问题