- 积分
- 21884
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2013-10-31
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
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的缺测值,但是资料里没有定义,我尝试了不同的定义缺测值的方法(红字)但是都没有得到想要的结果,以下是我的出图,可以看到色标是有问题的。
之后将脚本发给朋友,他用6.3.0的ncl在没有改动脚本的情况下作出了下图(资料都用的hadley海温)
可以看到的是,除了极区以外,其余的pattern两图之间都是相符的,中间打点的是T检验结果,也是没问题的,所以想问问这种情况究竟怎么处理缺测值。。。新人刚刚接触ncl不足一个月,一切还在学习中
|
|