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

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 9813|回复: 19

[作图] 关于资料中缺测值的处理存在一个问题

[复制链接]

新浪微博达人勋

发表于 2016-9-4 21:35:36 | 显示全部楼层 |阅读模式

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

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

x
数据用的是hadley的sst资料,下面贴代码

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
  fili    = "/home/lzb/sst.nc"
                                      ; open the netCDF file
  f    = addfile (fili, "r")

  yrStrt = 1870
  yrLast = 2013
; latS   = -60
; latN   =  60
  TIME  = f->time
  YYYY   = cd_calendar(TIME,-1)/100
  iYYYY  = ind(YYYY.ge.yrStrt .and. YYYY.le.yrLast)

  sst    = f->sst(iYYYY,:,:)

;**********************************************

sst = where(sst.eq.-1000,sst,sst@_FillValue)
sst = where(ismissing(sst),-1000,sst)
   
sst@_FillValue = -999
;**********************************************  

;  print(iYYYY)
;  printVarSummary(sst)
  date = iYYYY  
;  print(date)
  dsst = f->sst(iYYYY,:,:)
;  printVarSummary(dsst)
                                      ; compute SEASONAL means
  sst_ts = dsst(time|:,latitude|:,longitude|:)    ;
  sst_ts = runave (sst_ts, 3, 0)      ; 3-mo run ave
;  printVarSummary(sst_ts)

;********************************
; indices for special dates
;*********************************
  ind1 = ind(date.eq.876)       ;
; print(ind1)
  ind2 = ind(date.eq.1271)       ;
  ind3 = ind(date.eq.1536)       ;
  ind4 = ind(date.eq.1727)       ;
;********************************
; compute [2-mon] climatologies
; over specified periods
;********************************
  sstAve1 = clmMonTLL ( sst_ts(ind1:ind2,:,:) )
; printVarSummary(sstAve1)
  sstStd1 = stdMonTLL ( sst_ts(ind1:ind2,:,:) )

  sstAve2 = clmMonTLL ( sst_ts(ind3:ind4,:,:) )
  sstStd2 = stdMonTLL ( sst_ts(ind3:ind4,:,:) )
;********************************
; compute probabilities for means difference
;********************************
  prob = ttest(sstAve2, sstStd2^2, 10 \
              ,sstAve1, sstStd1^2, 10 ,True, False )

  copy_VarCoords (sstAve2, prob)
  prob@long_name = "Probability: difference between means"
                                     ; compute differences  
  difAve = sstAve2-sstAve1
; printVarSummary(difAve)
  copy_VarCoords (sstAve2, difAve)
  difAve@long_name = "98_13-43_75"
  difAve@units     = " C "
;********************************
; create plot
;********************************   
  nmo = 0                                ; for demo, only plot Dec-Feb
  wks   = gsn_open_wks ("pdf", "climHIATUSminus1" ) ; open workstation
  gsn_define_colormap(wks,"BlueRed")
  res = True                             ; plot mods desired

  res@mpCenterLonF         = 180

  res@gsnDraw              = False           ; Do not draw plot
  res@gsnFrame             = False           ; Do not advance frame
  res@gsnSpreadColors      = True         ; spread out color table
  res@cnFillOn             = True         ; turn on color fill
  res@cnLinesOn            = False        ; True is default
  res@cnLineLabelsOn       = False        ; True is default
  res@lbLabelBarOn         = True        ; turn off individual lb's

  res@tiMainString         = "SST Difference: 98_13-43_75"
  res@gsnCenterString      = "5% stippled"
  res@gsnLeftString        = "DJF"

; res@cnLevelSelectionMode = "ManualLevels"  ; set manual contour levels
; res@cnMinLevelValF       = -6             ; set min contour level
; res@cnMaxLevelValF       =  6             ; set max contour level
; res@cnLevelSpacingF      =  0.048            ; set contour spacing

  plot = gsn_csm_contour_map_ce(wks,difAve(nmo,:,:), res)  
  plot = ZeroNegDashLineContour (plot)

  res2 = True                            ; res2 probability plots
  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      = 1.05        ; set max contour level
  res2@cnLevelSpacingF     = 0.05        ; 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.6         ; add extra density

  delete(prob@long_name)
                                         ; add cyclic point
  plot2   = gsn_csm_contour(wks,gsn_add_cyclic_point(prob(nmo,:,:)), res2)
  plot2   = ShadeLtContour(plot2, 0.07, 17) ; shade all areas < 0.07 contour
  overlay (plot, plot2)

  draw (plot)
  frame(wks)

end


问题是这样,用ncdump v sst sst.nc |less 检查出除了资料中默认缺测值-1e+30外,在极区存在很多 -1000的缺测值,但是资料里没有定义,我尝试了不同的定义缺测值的方法(红字)但是都没有得到想要的结果,以下是我的出图,可以看到色标是有问题的。
QQ图片20160904212734.png

之后将脚本发给朋友,他用6.3.0的ncl在没有改动脚本的情况下作出了下图(资料都用的hadley海温)

QQ图片20160904212724.png

可以看到的是,除了极区以外,其余的pattern两图之间都是相符的,中间打点的是T检验结果,也是没问题的,所以想问问这种情况究竟怎么处理缺测值。。。新人刚刚接触ncl不足一个月,一切还在学习中

密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2016-9-8 17:31:39 | 显示全部楼层
关于这个问题我之前专门给Met Office发过邮件问过,这是工作人员的回答(不知道怎么上图,直接复制过来了):

Dear Lee,

Thank you for your enquiry.

The SST values of -1000 indicate areas of ocean where the sea ice concentration is 100%. We use -1000 rather than the missing data indicator because there is no defined SST (being fully ice covered, there is no sea surface) but there is data.
I hope that you will find this useful.

Kind Regards

Christina

密码修改失败请联系微信:mofangbao
回复 支持 1 反对 0

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2016-9-4 21:36:46 | 显示全部楼层
欢迎交流,谢谢大家
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2016-9-5 08:32:58 | 显示全部楼层
hadley sst 里-1.8和-1000好像都表示海冰,具体区别我也没搞明白
你只要在读入数据之后将所有小于0的都变成_fillvalue就行了。虽然说海水-4度才结冰,但是一般数据中用0度截一下就可以了
把缺省值改成-999没必要
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2016-9-5 11:03:33 | 显示全部楼层
井中月 发表于 2016-9-5 08:32
hadley sst 里-1.8和-1000好像都表示海冰,具体区别我也没搞明白
你只要在读入数据之后将所有小于0的都变 ...

对的,我也发现有少数的-1.8 ......
我用一个where函数   
sst = where(sst.lt.0,sst,sst@_FillValue)
这样写图还是没有变化,是语句位置有问题吗?
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2016-9-5 11:09:55 | 显示全部楼层
井中月 发表于 2016-9-5 08:32
hadley sst 里-1.8和-1000好像都表示海冰,具体区别我也没搞明白
你只要在读入数据之后将所有小于0的都变 ...

找到错误了、、、多谢!
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2016-9-9 15:43:16 | 显示全部楼层
Adam-Lee 发表于 2016-9-8 17:31
关于这个问题我之前专门给Met Office发过邮件问过,这是工作人员的回答(不知道怎么上图,直接复制过来了) ...

对,师兄师姐说-1000是海冰的值,那-1.8又是代表什么呢
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2016-9-10 11:14:37 | 显示全部楼层
本帖最后由 Adam-Lee 于 2016-9-10 11:18 编辑
南气斑织逍遥 发表于 2016-9-9 15:43
对,师兄师姐说-1000是海冰的值,那-1.8又是代表什么呢

-1.8°C是默认的表面海水凝结的温度,小于-1.8°C的海温值没有意义
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2016-9-10 13:40:17 | 显示全部楼层
Adam-Lee 发表于 2016-9-10 11:14
-1.8°C是默认的表面海水凝结的温度,小于-1.8°C的海温值没有意义

好的,谢谢
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2017-1-7 21:42:06 | 显示全部楼层
南气斑织逍遥 发表于 2016-9-5 11:09
找到错误了、、、多谢!

lz你好,这个问题你是怎么解决的呢?
密码修改失败请联系微信:mofangbao
您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

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

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

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