请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 24604|回复: 21

[作图] ncl 小波分析,打点问题

[复制链接]

新浪微博达人勋

发表于 2017-6-3 10:40:06 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
用ncl做小波分析,想在通过显著性检验的区域打点,可是总感觉打点的位置有问题,以下是我小波以及画图的程序,希望有大神可以指错


load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
begin
f =addfile("NPO_index.nc","r") ; 打开文件
npo=f->npo
N    = dimsizes(npo)
; compute wavelet (Missing values not allowed due to use of FFT)
;************************************  
  mother  = 0
  param   = 6.0
  dt      = 1.              ;for NPO (timesteps in units of years)
  s0      = 0.25
  dj      = 0.25
  jtot    = 1 + floattointeger(((log10(N*dt/s0))/dj)/log10(2.))
; print(jtot)
  npad    = N
  nadof   = new(2,float)
  noise   = 1
  siglvl  = .05
  isigtest= 0
  

  w = wavelet(npo,mother,dt,param,s0,dj,jtot,npad,noise,isigtest,siglvl,nadof)
;************************************
; create coodinate arrays for plot
;************************************
  power            = onedtond(w@power,(/jtot,N/))
  power!0          = "period"                        ; Y axis
  power&period     = w@period
  power!1          = "year"                          ; X axis
  power&year      = year
  power@long_name  = "Power Spectrum"
  power@units      = "unit^2"
  gws              = w@gws
; compute significance ( >= 1 is significant)
  SIG              = power                            ; transfer meta data
  SIG              = power/conform (power,w@signif,0)
  SIG@long_name    = "Significance"
  SIG@units        = " "

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
wks = gsn_open_wks("eps","NPO_wavelet")              ; open ps file
gsn_define_colormap(wks,"BlAqGrYeOrReVi200")    ; choose colormap
res                     = True                  ; plot mods desired
  res@gsnFrame            = False                 ; don't advance frame yet
  res@gsnDraw             = False                 ; don't draw yet
  res@cnFillOn            = True                  ; turn on color
  res@cnFillMode          = "RasterFill"          ; turn on raster mode
  res@cnRasterSmoothingOn = True                  ; turn on raster smoothing
  res@cnLinesOn           = False                 ; turn off contour lines
  res@gsnSpreadColors     = True                  ; use full colormap
  res@lbOrientation       = "Vertical"            ; vertical label bar
  res@trYReverse          = True                  ; reverse y-axis
  res@tmLabelAutoStride   = True
  res@tmYLMode = "Explicit"
  res@tmYLValues = w@period
  res@tmYLLabels = w@period
plot2 = gsn_csm_contour(wks,power({0:66},:),res)
  res1  = True
  res1@cnFillOn=False
  plot2 = ShadeCOI(wks,plot2,w,year,res1)
  res2 = True
  res2@gsnDraw = False
  res2@gsnFrame = False
  res2@cnFillOn = True
  res2@cnLinesOn = False
  res2@cnLineLabelsOn = False
  res2@cnInfoLabelOn = False
  res2@lbLabelBarOn = False
  res2@cnMonoFillPattern = False
  res2@cnLevelSelectionMode = "ExplicitLevels"
  res2@cnLevels=(/0.05/)
  res2@cnFillPatterns = (/17,-1/)
  res2@cnFillColors = (/1,0/)
  res2@gsnLeftString = ""
  plot3 = gsn_csm_contour(wks,SIG({0:66},:),res2)
  overlay(plot2,plot3)
draw(plot2)

end



微信截图_20170603104714.png
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2017-6-4 12:13:04 | 显示全部楼层
已经自行解决了,希望可以帮助和我有过遇到同样问题的人
之前做小波检验时没有注意到ncl官网中明确表示“SIG大于1为显著”,就按以往回归检验时的打点画图了,这样就会出现我之前的问题。
在画图时设置应把
res2@cnLevelSelectionMode = "ExplicitLevels"
  res2@cnLevels=(/0.05/)
改为
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
注:假如,SIG最大值为4,则
res2@cnMaxLevelValF      = 4.00      ; set max contour level
  res2@cnLevelSpacingF     = 3.00        ; set contour spacing
最大值和间隔要按实际情况修改
密码修改失败请联系微信:mofangbao
回复 支持 4 反对 0

使用道具 举报

新浪微博达人勋

发表于 2017-6-12 15:18:36 | 显示全部楼层
帽帽 发表于 2017-6-12 13:20
我没有太明白你的问题,纵坐标我没调过

不好意思我没表达清楚 就是我看一些matlab的例子和一些参考文献 看到横坐标一般都是年份 纵坐标是周期 但是周期一般是 2 4 8 46 32等等 然后我用ncl做的时候也是和你的这个类似 我就不知道应该怎么调整这个纵坐标的值
密码修改失败请联系微信:mofangbao
回复 支持 1 反对 0

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2017-6-4 09:36:00 | 显示全部楼层
求大神啊  ,,,,,,
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2017-6-4 09:46:51 | 显示全部楼层
{:eb303:}{:eb303:}{:eb303:}
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

新浪微博达人勋

发表于 2017-6-4 19:41:48 | 显示全部楼层
多谢楼主分享{:5_213:}
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2017-6-5 10:31:42 | 显示全部楼层
你好我想请问一下  你的那个纵坐标是怎么弄得啊?
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2017-6-5 10:52:13 来自手机 | 显示全部楼层
看不懂,来水一波经验。zzzz
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2017-6-12 13:20:43 | 显示全部楼层
yeah... 发表于 2017-6-5 10:31
你好我想请问一下  你的那个纵坐标是怎么弄得啊?

我没有太明白你的问题,纵坐标我没调过
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2017-6-12 19:43:22 | 显示全部楼层
我是这样改
res@tmYLValues =(/1,2,4,8,16,32,64/)
res@tmYLLabels = (/"1","2","4","8","16","32","64"/)
你可以试试
密码修改失败请联系微信:mofangbao
您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表