Why
相关分析我之前也写过,为什么还要有这篇文章?因为我希望将来用模块化的思维使用NCL。计算是一部分,画图是另一部分。这样每次一个routine需要做什么就调用之前写好的template就可以。思路上共3部分:计算;输出;leadlag以及合并。
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
lead_lag_num=12
target_month=13
target_index_name="tp_sh_stat"
target_index_raw = asciiread("~/project/iobm_plus/sst_eof/iobm_lsf/IOBM_monthly_index_"+target_month+".txt", -1, "float")
year = ispan(1979,2014,1)
year_start = 1981
year_end = 2008
year_num = year_end-year_start+1
target_index = target_index_raw
target_index!0 = "year"
target_index&year = year
;index = dim_standardize_Wrap(target_index({year_start:year_end},{target_month}),0)
index = dim_standardize_Wrap(dtrend(target_index({year_start:year_end}),False),0)
copy_VarCoords(target_index({year_start:year_end}),index)
printVarSummary(index)
file_list=systemfunc("ls /uwnd.mon.mean.1958-2014.nc")
a=addfile(file_list,"r")
sst_summer=short2flt(a->uwnd((year_start-1958)*12+target_month-lead_lag_num-1:(year_end-1958)*12+target_month-lead_lag_num-1:12,:,:,:))
printVarSummary(sst_summer)
sst_summer_dtrend=dtrend_n(sst_summer,True,0)
copy_VarCoords(sst_summer,sst_summer_dtrend)
coe=escorc(index,sst_summer_dtrend(lev|:,lat|:,lon|:,time|:))
copy_VarCoords(sst_summer_dtrend(0,:,:,:),coe)
printVarSummary(coe)
printMinMax(coe,False)
coe_ttest=coe
coe_ttest=coe@_FillValue
dim_coe_ttest=dimsizes(coe_ttest)
do i=0,dim_coe_ttest(0)-1
do j=0,dim_coe_ttest(1)-1
do k=0,dim_coe_ttest(2)-1
n = dimsizes(index) ; n=11
df = n-2 ; degree of freedom
t = coe(i,j,k)*sqrt((n-2)/(1-coe(i,j,k)^2))
p = student_t(t, df)
coe_ttest(i,j,k)=p
end do
end do
end do
copy_VarCoords(coe,coe_ttest)
printMinMax(coe_ttest,False)
printVarSummary(coe_ttest)
output_filename=get_script_prefix_name+"_"+target_month+"_"+lead_lag_num+".nc"
system("rm -f "+output_filename)
PDataOut=addfile("./"+get_script_prefix_name+"/"+output_filename,"c")
;***************************************************
; define file options [Version a033]
;***********************************************
setfileoption(PDataOut, "prefill", False)
setfileoption(PDataOut, "suppressclose", True)
setfileoption(PDataOut, "definemode", True)
;***********************************************
; assign file attributes
;***********************************************
;*****************************************************
fAtt = True
fAtt@title = "correlation between GODAS ts and "+target_index_name+"monthly data"
fAtt@source = "GODAS station obs"
fAtt@Conventions = "None"
fAtt@creation_date = systemfunc("date")
fileattdef( PDataOut, fAtt )
;***********************************************
; predefine coordinate information
;***********************************************
dimNames = (/"lev", "lat", "lon" /)
dimSizes = (/dimsizes(coe&lev),dimsizes(coe&lat),dimsizes(coe&lon)/)
dimUnlim = (/False,False,False/)
filedimdef(PDataOut, dimNames,dimSizes,dimUnlim)
filevardef(PDataOut, "lat", typeof(coe&lat), getvardims(coe&lat))
filevarattdef(PDataOut, "lat", coe&lat)
filevardef(PDataOut, "lon", typeof(coe&lon), getvardims(coe&lon) )
filevarattdef(PDataOut, "lon", coe&lon)
filevardef(PDataOut,"lev",typeof(coe&lev), getvardims(coe&lev) )
filevarattdef(PDataOut,"lev", coe&lev)
filevardef(PDataOut,"coe",typeof(coe), getvardims(coe) )
filevarattdef(PDataOut,"coe", coe)
filevardef(PDataOut,"coe_ttest",typeof(coe_ttest), getvardims(coe_ttest) )
filevarattdef(PDataOut,"coe_ttest", coe_ttest)
;***********************************************
; terminate define mode: not necessary but for clarity
;***********************************************
setfileoption(PDataOut, "definemode", False)
;***********************************************
; write data values to predefined locations
; (/ .../) operator transfer values only
;***********************************************
PDataOut->lat = (/ coe&lat /)
PDataOut->lon = (/ coe&lon /)
PDataOut->lev = (/ coe&lev /)
PDataOut->coe = (/ coe /)
PDataOut->coe_ttest = (/ coe_ttest /)
end
接下来用nco的命令ncecat将之前算的文件合并即可。之所以不直接都写在ncl里,一是考虑到高内聚低冗杂的模块化思想,二是考虑到可以降低内存方面的负载。至于为什么不在bash里优雅地使用循环。。好吧,长时间不用已经不想写bash了,觉得bash这种语言可读性不高,简单地句子分分钟明白(写几层循环过几天看整个人都不好了),而且本人Vim玩的很熟,上面地leadlag也就是录个宏搞定地,写下面这个东西也是几个正则表达式就搞定了。
ncecat -u lead_lag "${0%.*}"_1_-12.nc "${0%.*}"_1_-11.nc "${0%.*}"_1_-10.nc "${0%.*}"_1_-9.nc "${0%.*}"_1_-8.nc "${0%.*}"_1_-7.nc "${0%.*}"_1_-6.nc "${0%.*}"_1_-5.nc "${0%.*}"_1_-4.nc "${0%.*}"_1_-3.nc "${0%.*}"_1_-2.nc "${0%.*}"_1_-1.nc "${0%.*}"_1_0.nc "${0%.*}"_1_1.nc "${0%.*}"_1_2.nc "${0%.*}"_1_3.nc "${0%.*}"_1_4.nc "${0%.*}"_1_5.nc "${0%.*}"_1_6.nc "${0%.*}"_1_7.nc "${0%.*}"_1_8.nc "${0%.*}"_1_9.nc "${0%.*}"_1_10.nc "${0%.*}"_1_11.nc "${0%.*}"_1_12.nc "${0%.*}"_1.nc
ncecat -u lead_lag "${0%.*}"_2_-12.nc "${0%.*}"_2_-11.nc "${0%.*}"_2_-10.nc "${0%.*}"_2_-9.nc "${0%.*}"_2_-8.nc "${0%.*}"_2_-7.nc "${0%.*}"_2_-6.nc "${0%.*}"_2_-5.nc "${0%.*}"_2_-4.nc "${0%.*}"_2_-3.nc "${0%.*}"_2_-2.nc "${0%.*}"_2_-1.nc "${0%.*}"_2_0.nc "${0%.*}"_2_1.nc "${0%.*}"_2_2.nc "${0%.*}"_2_3.nc "${0%.*}"_2_4.nc "${0%.*}"_2_5.nc "${0%.*}"_2_6.nc "${0%.*}"_2_7.nc "${0%.*}"_2_8.nc "${0%.*}"_2_9.nc "${0%.*}"_2_10.nc "${0%.*}"_2_11.nc "${0%.*}"_2_12.nc "${0%.*}"_2.nc
ncecat -u lead_lag "${0%.*}"_3_-12.nc "${0%.*}"_3_-11.nc "${0%.*}"_3_-10.nc "${0%.*}"_3_-9.nc "${0%.*}"_3_-8.nc "${0%.*}"_3_-7.nc "${0%.*}"_3_-6.nc "${0%.*}"_3_-5.nc "${0%.*}"_3_-4.nc "${0%.*}"_3_-3.nc "${0%.*}"_3_-2.nc "${0%.*}"_3_-1.nc "${0%.*}"_3_0.nc "${0%.*}"_3_1.nc "${0%.*}"_3_2.nc "${0%.*}"_3_3.nc "${0%.*}"_3_4.nc "${0%.*}"_3_5.nc "${0%.*}"_3_6.nc "${0%.*}"_3_7.nc "${0%.*}"_3_8.nc "${0%.*}"_3_9.nc "${0%.*}"_3_10.nc "${0%.*}"_3_11.nc "${0%.*}"_3_12.nc "${0%.*}"_3.nc
ncecat -u lead_lag "${0%.*}"_4_-12.nc "${0%.*}"_4_-11.nc "${0%.*}"_4_-10.nc "${0%.*}"_4_-9.nc "${0%.*}"_4_-8.nc "${0%.*}"_4_-7.nc "${0%.*}"_4_-6.nc "${0%.*}"_4_-5.nc "${0%.*}"_4_-4.nc "${0%.*}"_4_-3.nc "${0%.*}"_4_-2.nc "${0%.*}"_4_-1.nc "${0%.*}"_4_0.nc "${0%.*}"_4_1.nc "${0%.*}"_4_2.nc "${0%.*}"_4_3.nc "${0%.*}"_4_4.nc "${0%.*}"_4_5.nc "${0%.*}"_4_6.nc "${0%.*}"_4_7.nc "${0%.*}"_4_8.nc "${0%.*}"_4_9.nc "${0%.*}"_4_10.nc "${0%.*}"_4_11.nc "${0%.*}"_4_12.nc "${0%.*}"_4.nc
ncecat -u lead_lag "${0%.*}"_5_-12.nc "${0%.*}"_5_-11.nc "${0%.*}"_5_-10.nc "${0%.*}"_5_-9.nc "${0%.*}"_5_-8.nc "${0%.*}"_5_-7.nc "${0%.*}"_5_-6.nc "${0%.*}"_5_-5.nc "${0%.*}"_5_-4.nc "${0%.*}"_5_-3.nc "${0%.*}"_5_-2.nc "${0%.*}"_5_-1.nc "${0%.*}"_5_0.nc "${0%.*}"_5_1.nc "${0%.*}"_5_2.nc "${0%.*}"_5_3.nc "${0%.*}"_5_4.nc "${0%.*}"_5_5.nc "${0%.*}"_5_6.nc "${0%.*}"_5_7.nc "${0%.*}"_5_8.nc "${0%.*}"_5_9.nc "${0%.*}"_5_10.nc "${0%.*}"_5_11.nc "${0%.*}"_5_12.nc "${0%.*}"_5.nc
ncecat -u lead_lag "${0%.*}"_6_-12.nc "${0%.*}"_6_-11.nc "${0%.*}"_6_-10.nc "${0%.*}"_6_-9.nc "${0%.*}"_6_-8.nc "${0%.*}"_6_-7.nc "${0%.*}"_6_-6.nc "${0%.*}"_6_-5.nc "${0%.*}"_6_-4.nc "${0%.*}"_6_-3.nc "${0%.*}"_6_-2.nc "${0%.*}"_6_-1.nc "${0%.*}"_6_0.nc "${0%.*}"_6_1.nc "${0%.*}"_6_2.nc "${0%.*}"_6_3.nc "${0%.*}"_6_4.nc "${0%.*}"_6_5.nc "${0%.*}"_6_6.nc "${0%.*}"_6_7.nc "${0%.*}"_6_8.nc "${0%.*}"_6_9.nc "${0%.*}"_6_10.nc "${0%.*}"_6_11.nc "${0%.*}"_6_12.nc "${0%.*}"_6.nc
ncecat -u lead_lag "${0%.*}"_7_-12.nc "${0%.*}"_7_-11.nc "${0%.*}"_7_-10.nc "${0%.*}"_7_-9.nc "${0%.*}"_7_-8.nc "${0%.*}"_7_-7.nc "${0%.*}"_7_-6.nc "${0%.*}"_7_-5.nc "${0%.*}"_7_-4.nc "${0%.*}"_7_-3.nc "${0%.*}"_7_-2.nc "${0%.*}"_7_-1.nc "${0%.*}"_7_0.nc "${0%.*}"_7_1.nc "${0%.*}"_7_2.nc "${0%.*}"_7_3.nc "${0%.*}"_7_4.nc "${0%.*}"_7_5.nc "${0%.*}"_7_6.nc "${0%.*}"_7_7.nc "${0%.*}"_7_8.nc "${0%.*}"_7_9.nc "${0%.*}"_7_10.nc "${0%.*}"_7_11.nc "${0%.*}"_7_12.nc "${0%.*}"_7.nc
ncecat -u lead_lag "${0%.*}"_8_-12.nc "${0%.*}"_8_-11.nc "${0%.*}"_8_-10.nc "${0%.*}"_8_-9.nc "${0%.*}"_8_-8.nc "${0%.*}"_8_-7.nc "${0%.*}"_8_-6.nc "${0%.*}"_8_-5.nc "${0%.*}"_8_-4.nc "${0%.*}"_8_-3.nc "${0%.*}"_8_-2.nc "${0%.*}"_8_-1.nc "${0%.*}"_8_0.nc "${0%.*}"_8_1.nc "${0%.*}"_8_2.nc "${0%.*}"_8_3.nc "${0%.*}"_8_4.nc "${0%.*}"_8_5.nc "${0%.*}"_8_6.nc "${0%.*}"_8_7.nc "${0%.*}"_8_8.nc "${0%.*}"_8_9.nc "${0%.*}"_8_10.nc "${0%.*}"_8_11.nc "${0%.*}"_8_12.nc "${0%.*}"_8.nc
ncecat -u lead_lag "${0%.*}"_9_-12.nc "${0%.*}"_9_-11.nc "${0%.*}"_9_-10.nc "${0%.*}"_9_-9.nc "${0%.*}"_9_-8.nc "${0%.*}"_9_-7.nc "${0%.*}"_9_-6.nc "${0%.*}"_9_-5.nc "${0%.*}"_9_-4.nc "${0%.*}"_9_-3.nc "${0%.*}"_9_-2.nc "${0%.*}"_9_-1.nc "${0%.*}"_9_0.nc "${0%.*}"_9_1.nc "${0%.*}"_9_2.nc "${0%.*}"_9_3.nc "${0%.*}"_9_4.nc "${0%.*}"_9_5.nc "${0%.*}"_9_6.nc "${0%.*}"_9_7.nc "${0%.*}"_9_8.nc "${0%.*}"_9_9.nc "${0%.*}"_9_10.nc "${0%.*}"_9_11.nc "${0%.*}"_9_12.nc "${0%.*}"_9.nc
ncecat -u lead_lag "${0%.*}"_10_-12.nc "${0%.*}"_10_-11.nc "${0%.*}"_10_-10.nc "${0%.*}"_10_-9.nc "${0%.*}"_10_-8.nc "${0%.*}"_10_-7.nc "${0%.*}"_10_-6.nc "${0%.*}"_10_-5.nc "${0%.*}"_10_-4.nc "${0%.*}"_10_-3.nc "${0%.*}"_10_-2.nc "${0%.*}"_10_-1.nc "${0%.*}"_10_0.nc "${0%.*}"_10_1.nc "${0%.*}"_10_2.nc "${0%.*}"_10_3.nc "${0%.*}"_10_4.nc "${0%.*}"_10_5.nc "${0%.*}"_10_6.nc "${0%.*}"_10_7.nc "${0%.*}"_10_8.nc "${0%.*}"_10_9.nc "${0%.*}"_10_10.nc "${0%.*}"_10_11.nc "${0%.*}"_10_12.nc "${0%.*}"_10.nc
ncecat -u lead_lag "${0%.*}"_11_-12.nc "${0%.*}"_11_-11.nc "${0%.*}"_11_-10.nc "${0%.*}"_11_-9.nc "${0%.*}"_11_-8.nc "${0%.*}"_11_-7.nc "${0%.*}"_11_-6.nc "${0%.*}"_11_-5.nc "${0%.*}"_11_-4.nc "${0%.*}"_11_-3.nc "${0%.*}"_11_-2.nc "${0%.*}"_11_-1.nc "${0%.*}"_11_0.nc "${0%.*}"_11_1.nc "${0%.*}"_11_2.nc "${0%.*}"_11_3.nc "${0%.*}"_11_4.nc "${0%.*}"_11_5.nc "${0%.*}"_11_6.nc "${0%.*}"_11_7.nc "${0%.*}"_11_8.nc "${0%.*}"_11_9.nc "${0%.*}"_11_10.nc "${0%.*}"_11_11.nc "${0%.*}"_11_12.nc "${0%.*}"_11.nc
ncecat -u lead_lag "${0%.*}"_12_-12.nc "${0%.*}"_12_-11.nc "${0%.*}"_12_-10.nc "${0%.*}"_12_-9.nc "${0%.*}"_12_-8.nc "${0%.*}"_12_-7.nc "${0%.*}"_12_-6.nc "${0%.*}"_12_-5.nc "${0%.*}"_12_-4.nc "${0%.*}"_12_-3.nc "${0%.*}"_12_-2.nc "${0%.*}"_12_-1.nc "${0%.*}"_12_0.nc "${0%.*}"_12_1.nc "${0%.*}"_12_2.nc "${0%.*}"_12_3.nc "${0%.*}"_12_4.nc "${0%.*}"_12_5.nc "${0%.*}"_12_6.nc "${0%.*}"_12_7.nc "${0%.*}"_12_8.nc "${0%.*}"_12_9.nc "${0%.*}"_12_10.nc "${0%.*}"_12_11.nc "${0%.*}"_12_12.nc "${0%.*}"_12.nc
ncecat -u month "${0%.*}"_1.nc "${0%.*}"_2.nc "${0%.*}"_3.nc "${0%.*}"_4.nc "${0%.*}"_5.nc "${0%.*}"_6.nc "${0%.*}"_7.nc "${0%.*}"_8.nc "${0%.*}"_9.nc "${0%.*}"_10.nc "${0%.*}"_11.nc "${0%.*}"_12.nc "${0%.*}".nc
最后生成地文件大概就是这样:
Variable: f
Type: file
filename: coe_uwnd_iobm
path: coe_uwnd_iobm.nc
file global attributes:
creation_date : Wed Feb 22 20:00:53 CST 2017
Conventions : None
source : GODAS station obs
title : correlation between GODAS ts and tp_sh_statmonthly data
history : Fri Feb 24 07:28:23 2017: ncecat -u month coe_uwnd_iobm_1.nc coe_uwnd_iobm_2.nc coe_uwnd_iobm_3.c
Fri Feb 24 07:27:39 2017: ncecat -u lead_lag coe_uwnd_iobm_1_-12.nc coe_uwnd_iobm_1_-11.nc coe_uwnd_iobm_1_-10.nc
nco_openmp_thread_number : 1
dimensions:
month = 12 // unlimited
lat = 145
lon = 288
lev = 37
lead_lag = 25
variables:
float lat ( lat )
long_name : latitude
GridType : Cylindrical Equidistant Projection Grid
units : degrees_north
Dj : 1.25
Di : 1.25
Lo2 : -1.25
La2 : -90
Lo1 : 0
La1 : 90
float lon ( lon )
long_name : longitude
GridType : Cylindrical Equidistant Projection Grid
units : degrees_east
Dj : 1.25
Di : 1.25
Lo2 : -1.25
La2 : -90
Lo1 : 0
La1 : 90
integer lev ( lev )
long_name : isobaric level
units : hPa
float coe ( month, lead_lag, lev, lat, lon )
_FillValue : 9.96921e+36
float coe_ttest ( month, lead_lag, lev, lat, lon )
_FillValue : 9.96921e+36
网友评论