- 积分
- 2366
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2015-3-27
- 最后登录
- 1970-1-1
|
楼主 |
发表于 2018-1-29 09:48:41
|
显示全部楼层
为啥我没设置要贡献还是需要。。。
脚本如下：
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"
load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/run_cor.ncl"
latS_U = -20
latN_U = 20
lonL_U = 40
lonR_U = 120
; latS_P = -20
; latN_P = 20
; lonL_P = 40
; lonR_P = 120
yrStrt = 1974
yrLast = 2016
season = "JJA"
years = ispan(yrStrt,yrLast,1)
nyears = dimsizes(years)
print(nyears)
;================================================
;Variables List:
; var1--------------sst
; var2--------------olr
;================================================
;initial data
sst_in = addfile("/Users/zhangguangli/Desktop/work/data/sst/HadISST_sst.nc", "r")
time_sst = sst_in->time
XXXX = cd_calendar(time_sst,-1 )
t_index_start1 = ind(XXXX.eq.197401)
t_index_end1 = ind(XXXX.eq.201612)
sst = sst_in->sst(t_index_start1:t_index_end1,:,:)
sst_sm = month_to_season(sst, "JJA")
sst_summer = dim_standardize_n(sst_sm, 0, 0)
copy_VarCoords(sst_sm, sst_summer)
printVarSummary(sst)
delete(sst)
sst = dim_standardize_n(sst_summer, 0, 0)
copy_VarCoords(sst_summer, sst)
olr_in = addfile("/Users/zhangguangli/Desktop/work/data/olr.mon.mean_197406_201701.nc", "r")
olr1 = olr_in->olr
olr = short2flt(olr1)
printVarSummary(olr)
olr_0 = dim_avg_n_Wrap(olr(0:2,:,:), 0)
printVarSummary(olr_0)
time_olr = olr_in->time
XXXX1 = cd_calendar(time_olr,-1 )
t_index_start11 = ind(XXXX1.eq.197501)
t_index_end11 = ind(XXXX1.eq.201612)
olr_2 = olr_in->olr(t_index_start11:t_index_end11,:,:)
olr_22= short2flt(olr_2)
printVarSummary(olr_22)
olr_2_sm = month_to_season(olr_22, "JJA")
printVarSummary(olr_2_sm)
olr_sm = new((/43,73,144/),float)
olr_sm(0,:,:) = olr_0(:,:)
olr_sm(1:3,:,:) = olr_2_sm(0:2,:,:)
olr_sm(4,:,:) = dim_avg_n_Wrap(olr_2_sm, 0)
olr_sm(5:42,:,:) = olr_2_sm(4:41,:,:)
copy_VarCoords(olr_0,olr_sm(0,:,:))
olr_sm!0 = "time"
printVarSummary(olr_sm)
delete(olr)
olr = dim_standardize_n(olr_sm, 0, 0)
copy_VarCoords(olr_sm, olr)
printVarSummary(sst)
printVarSummary(olr)
nmca = 2 ; how many MCA we need
var1_region = sst(:,{latS_U:latN_U},{lonL_U:lonR_U});把sst数据reshape成二维
var1_size = dimsizes(var1_region)
n_var1_size = var1_size(1)*var1_size(2)
homlft = new((/nmca,n_var1_size/),float)
hetlft = new((/nmca,n_var1_size/),float)
var1_ano_line = reshape(var1_region,(/var1_size(0),n_var1_size/))
var1_ano_line!0 = "time"
var1_ano_line!1 = "pts"
printVarSummary(var1_ano_line)
printMinMax(var1_ano_line, False)
var2_region = olr(:,{latS_U:latN_U},{lonL_U:lonR_U});把olr数据reshape成二维
var2_size = dimsizes(var2_region)
n_var2_size = var2_size(1)*var2_size(2)
homrgt = new((/nmca,n_var2_size/),float)
hetrgt = new((/nmca,n_var2_size/),float)
var2_ano_line = reshape(var2_region,(/var2_size(0),n_var2_size/))
var2_ano_line!0 = "time"
var2_ano_line!1 = "pts"
printVarSummary(var2_ano_line)
printMinMax(var2_ano_line, False)
ntime = nyears ; # time steps
ncols = n_var1_size ; # columns (stations or grid pts) for S
ncolz = var2_size ; # columns (stations or grid pts) for Z
nsvd = 2 ; # svd patterns to calculate
; [nsvd <= min(ncols, ncolz) ]
; xmsg = -999.9 ; missing value
x = svdcov(var1_ano_line(pts|:,time|:),var2_ano_line(pts|:,time|:),nsvd,homlft,hetlft,homrgt,hetrgt)
print("svdcov: percent variance= " + x)
printVarSummary(x)
ak = onedtond(x@ak,(/nsvd,ntime/))
bk = onedtond(x@bk,(/nsvd,ntime/))
ak!0 = "sv"
ak!1 = "time"
bk!0 = "sv"
bk!1 = "time"
ccr = escorc(ak(0,:), bk(0,:))
print(ccr)
print(ak)
print(bk)
ak_std = dim_standardize_Wrap(ak,1)
bk_std = dim_standardize_Wrap(bk,1)
printVarSummary(ak_std)
printVarSummary(bk_std)
data = new((/2,nyears/),"float")
data(0,:) = ak_std(0,:)
data(1,:) = bk_std(0,:)
ccr_ak = escorc_n(ak(0,:), olr, 0, 0)*-1
copy_VarCoords(olr(0,:,:), ccr_ak)
printVarSummary(ccr_ak)
ccr_bk = escorc_n(bk(0,:), sst, 0, 0)*-1
copy_VarCoords(sst(0,:,:), ccr_bk)
printVarSummary(ccr_bk)
test_ak = rtest(ccr_ak, nyears, 0)
; test_ak = where(ccr_ak.eq.ccr_ak@_FillValue, test_ak@_FillValue, test_ak)
copy_VarCoords(ccr_ak, test_ak)
; test_bk = new((/180,360/),float)
; test_bk = sqrt(nyears-2)*ccr_bk/sqrt(1-ccr_bk*ccr_bk)
; test_bk=(/student_t(test_bk, nyears-2/)
test_bk = rtest(ccr_bk, nyears, 0)
test_bk = where(ccr_bk.eq.ccr_bk@_FillValue, test_bk@_FillValue, test_bk)
copy_VarCoords(ccr_bk, test_bk)
test_bk_reverse = 1 - test_bk
copy_VarCoords(test_bk, test_bk_reverse)
printVarSummary(test_ak)
printMinMax(test_ak, False)
printVarSummary(test_bk)
printMinMax(test_bk, False)
; ==============================================================
; Set the figure parameters
; ==============================================================
wks = gsn_open_wks("png","olr&sst_kongjian")
plot = new(2,"graphic") ;shading
plot_prob = new(2,"graphic") ;shading
gsn_define_colormap(wks,"BlueWhiteOrangeRed")
res = True
;res@cnFillDrawOrder = "PreDraw"
res@gsnDraw = False
res@gsnFrame = False ;don't advance frame
res@gsnAddCyclic = False ;EOF data is not cyclic
res@cnFillOn = True ; turn on color fill
res@cnLinesOn = False
res@cnLineLabelsOn = False ;turn off cn line labels
res@cnInfoLabelOn = False ;turn off contour information label
res@mpFillOn = True ; turn off map fill
res@mpMinLatF = latS_U
res@mpMaxLatF = latN_U
res@mpMinLonF = lonL_U
res@mpMaxLonF = lonR_U
res@mpCenterLonF = (lonL_U+lonR_U)/2
res@lbLabelBarOn = False
; res@lbOrientation = "Vertical"
res@pmLabelBarHeightF = 0.08
res@pmLabelBarWidthF = 0.55
res@lbLabelFontHeightF = 0.018
res@pmLabelBarOrthogonalPosF = 0.14
res@tmXBLabelFontHeightF = 0.022 ;font height of tick labels
res@tmYLLabelFontHeightF = 0.022
res@tmXBTickSpacingF = 10. ;label spacing
res@tmYLTickSpacingF = 10. ;label spacing
res@gsnLeftStringFontHeightF = 0.022
res@gsnRightStringFontHeightF = 0.022
res@gsnLeftString = "(a)SST"
res@gsnRightString = season
res@cnLevelSelectionMode = "ManualLevels"
res@cnMinLevelValF = -1
res@cnMaxLevelValF = 1
res@cnLevelSpacingF = 0.1
plot(0) = gsn_csm_contour_map(wks,ccr_bk,res)
res@tmYLTickSpacingF = 10.
res@tmXBTickSpacingF = 10. ;label spacing
res@pmLabelBarOrthogonalPosF = 0.25
res@mpMinLatF = latS_U
res@mpMaxLatF = latN_U
res@mpMinLonF = lonL_U
res@mpMaxLonF = lonR_U
res@mpCenterLonF = (lonL_U+lonR_U)/2
res@gsnLeftString = "(b)OLR"
res@cnLevelSelectionMode = "ManualLevels"
res@cnMinLevelValF = -1
res@cnMaxLevelValF = 1
res@cnLevelSpacingF = 0.1
plot(1) = gsn_csm_contour_map(wks,ccr_ak,res) ; create a default plot
res2 = True
res2@gsnDraw = False
res2@gsnFrame = False
res2@gsnDraw = False
res2@gsnFrame = False ;don't advance frame
res2@gsnAddCyclic = False ;EOF data is not cyclic
res2@cnFillOn = True ; turn on color fill
res2@cnLinesOn = False
res2@cnLineLabelsOn = False ;turn off cn line labels
res2@cnInfoLabelOn = False ;turn off contour information label
res2@cnFillOn = True
res2@cnLinesOn = False
res2@cnLineLabelsOn = False
res2@cnInfoLabelOn = False
res2@lbLabelBarOn = False
res2@cnMonoFillPattern = False
res2@cnLevelSelectionMode = "ExplicitLevels"
res2@cnLevels = (/0.1/) ;; set to significance level
res2@cnFillPatterns = (/17,-1/)
res2@cnFillColors = (/1,0/)
res2@cnFillDotSizeF = 0.003
res2@gsnLeftString = ""
res2@gsnAddCyclic = True
plot_prob(0) = gsn_csm_contour(wks, test_bk, res2)
; res2@cnLevelSelectionMode = "ExplicitLevels"
; res2@cnLevels = (/0.05/) ;; set to significance level
; res2@cnFillPatterns = (/-1,17/)
; res2@cnFillColors = (/1,0/)
plot_prob(1) = gsn_csm_contour(wks, test_ak, res2)
do i = 0, 1
overlay(plot(i), plot_prob(i))
end do
;************************************************
; panel plot only resources
;************************************************
resP = True
resP@gsnPanelLabelBar = True
; resP@lbOrientation = "Vertical"
resP@gsnPanelMainString = "JJA IO SST&OLR SVD"
gsn_panel(wks,plot,(/2,1/),resP) ; now draw as one plot
;************************************************
wks_ts = gsn_open_wks("png","olr&sst_shijian")
rts = True
rts@gsnDraw = False ; don't draw yet
rts@gsnFrame = False ; don't advance frame yet
rts@vpXF = 0.15
rts@vpWidthF = 0.8
rts@vpHeightF= 0.35
rts@tiYAxisString = "Standardized" ; y-axis label
; rts@tiYAxisString = " " ; y-axis label
rts@tmXBLabelFontHeightF = 0.02 ;font height of tick labels
rts@tmYLLabelFontHeightF = 0.02
rts@gsnLeftStringFontHeightF = 0.02
rts@gsnRightStringFontHeightF = 0.02
rts@gsnYRefLine = 0. ; reference line
rts@xyLineColors = (/"blue","red"/) ; colors chosen
rts@xyLineThicknesses = (/3.0,3.0/) ; line thicknesses
rts@xyDashPatterns = (/0.,0./) ; make all lines solid
rts@trXMinF = yrStrt ; leave a margin for legend
rts@trXMaxF = yrLast
;rts@tmXBMode = "Manual" ; Define own tick mark labels.
rts@tmXBTickSpacingF = 5
rts@trYMinF = -3.0 ; min value on x-axis
rts@trYMaxF = 3.0 ; max value on x-axis
rts@gsnLeftString = "Time Series of SVD1 (R = "+sprintf("%5.2f",ccr)+")"
rts@gsnRightString = sprintf("%5.2f",x(0))+"%"
plot_ts = gsn_csm_xy(wks_ts,years,data*-1,rts)
;---------------------------Add plot legend-----------------------------------
res_lines = True ; polyline mods desired
res_lines@gsLineDashPattern = 0. ; solid line
res_lines@gsLineThicknessF = 3. ; line thicker
res_text = True ; text mods desired
res_text@txFontHeightF = 0.015 ; change text size
res_text@txJust = "CenterLeft" ; text justification
res_lines@gsLineColor = "blue"
yy = (/2.5,2.5/)
xx = (/2007,2010/)
text_pos = 2011
dum11 = gsn_add_polyline(wks_ts,plot_ts,xx,yy,res_lines) ; add polyline
dum12 = gsn_add_text(wks_ts,plot_ts,"SST",text_pos,yy(0),res_text); add text
res_lines@gsLineColor = "red"
yy = yy - 0.6
dum21 = gsn_add_polyline(wks_ts,plot_ts,xx,yy,res_lines) ; add polyline
dum22 = gsn_add_text(wks_ts,plot_ts,"OLR",text_pos,yy(0),res_text); add text
draw(plot_ts)
frame(wks)
frame(wks_ts)
|
|
[复制链接]
窥视卡
雷达卡
发表于 2018-1-29 09:46:15
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
抢沙发
千斤顶
显身卡
发表于 2018-8-12 16:21:11
{:5_213:}{:5_213:}{:5_213:}