积分 2471
贡献
精华
在线时间 小时
注册时间 2016-3-29
最后登录 1970-1-1
30 金钱
做秋雨指数与850hpa温度的相关分析,秋雨指数处理成每年一个数,共计60年,一个是风场的显著性,一个是温度的显著性,显著性区域都成方块状,有人知道是为什么吗?
begin
latS = 0
latN = 90
lonL = 0
lonR = 150
siglvl = 0.05
;;-----------------------------------------------------------------------
f1 = addfile("pre2.nc", "r") ;降水
pre = f1->pre(:,{latS:latN},{lonL:lonR})
index = dim_avg_n_Wrap(pre, (/1,2/));指数,只剩时间维度,后面才能做相关; [time | 60]
;printVarSummary(index)
;;---------------------------------------------------------------------
f2 = addfile("air.mean.nc", "r") ;(都是每年一个数)
airtmp = f2->air(:,{850},{latS:latN},{lonL:lonR}) ; [time | 60] x [lat | 21] x [lon | 37]
printVarSummary(airtmp)
dimx = dimsizes(airtmp)
ntim = dimx(0)
;;-----------------------------------------------------------------------
mxlag = 0 ;同时相关
rairtmp = esccr(index, airtmp({lat|:},{lon|:},{time|:}), mxlag); [21] x [37] x [1]
printVarSummary(rairtmp)
rairtmp1 :=(/rairtmp(:,:,0)/)
copy_VarMeta(airtmp(0,:,:), rairtmp1);把经纬度元数据赋值给rv1
;;---------------------------------------------------------------------------------
n = ntim
df = n-2
t = rairtmp1*sqrt((n-2)/(1-rairtmp1^2)) ;t检验
p = student_t(t, df) ;两侧student-t分布概率值
printVarSummary(p) ; [21] x [37]
psig = 0.05
copy_VarMeta(airtmp(0,:,:), p)
;;----------------------------------------------------------------------------
wks = gsn_open_wks("png", "corr.airtmp")
gsn_define_colormap(wks, "GMT_panoply")
res = True
res@gsnDraw = False
res@gsnFrame = False
res@gsnMaximize = True
res@gsnLeftString = ""
res@gsnRightString = ""
res@pmTitleDisplayMode = "Always"
res@pmTickMarkDisplayMode = "Always"
;;set map;;
mpres = res
mpres@mpDataSetName = "Earth..4"
mpres@mpDataBaseVersion = "MediumRes"
mpres@mpOutlineOn = True
mpres@mpOutlineSpecifiers = (/"China:states","Taiwan"/)
mpres@mpGeophysicalLineThicknessF = 2
mpres@mpNationalLineThicknessF = 2
mpres@mpFillDrawOrder = "PostDraw"
mpres@mpFillOn = False
;;set area;;
mpres@mpMinLatF = 0
mpres@mpMaxLatF = 90
mpres@mpMinLonF = 0
mpres@mpMaxLonF = 150
map = gsn_csm_map(wks,mpres)
;;--------------------------------------------------------------------------相关系数线
cnres = res
cnres@cnFillDrawOrder = "PreDraw"
cnres@cnFillOn = False
cnres@cnLinesOn = True
cnres@cnLineLabelsOn = True
cnres@cnLineThicknesses = 3. ;等值线粗细
cnres@cnSmoothingOn = True
cnres@cnSmoothingDistanceF = 0.001
cnres@cnSmoothingTensionF = -2
cnres@cnInfoLabelOn = False ;右下方等值线信息标签
cnres@gsnContourNegLineDashPattern = 0
cnres@gsnContourPosLineDashPattern = 2
contour = gsn_csm_contour(wks,rairtmp1,cnres)
contour = ColorNegDashZeroPosContour(contour,"black","black","black")
;;-------------------------------------------------------------------------------------
sres = True ; set up a second resource list
sres@gsnDraw = False ; do not draw the plot
sres@gsnFrame = False ; do not advance the frame
sres@cnLineLabelsOn = False ; do not use line labels
sres@cnLinesOn = False ; do not draw contour lines
sres@cnInfoLabelOn = False
sres@cnFillOn = False ; color fill
sres@lbLabelBarOn = False
sres@cnInfoLabelOn = False
sres@gsnLeftString = ""
sres@gsnRightString = ""
sres@cnFillDotSizeF = 0.005
;;---------------------------------------------------------------
opt = True
opt@gsnShadeFillType = "color"
opt@gsnShadeMid = "gray"
Rxt_test = where(p.le.0.05, rairtmp1, rairtmp1@_FillValue)
copy_VarMeta(rairtmp1(:,:), Rxt_test(:,:))
printVarSummary(Rxt_test)
Rxt_max = max(Rxt_test)
Rxt_min = min(Rxt_test)
plot1 = gsn_csm_contour(wks,Rxt_test,sres)
plot1 = gsn_contour_shade(plot1,Rxt_min,Rxt_max,opt)
overlay(map,plot1)
overlay(map,contour)
draw(map)
frame(wks)
end
我来回答