NCL脚本利用站点资料做EOF分析

我是小猪

1. ; ==============================================================
   
2. ; eof_0a.ncl
   
3. ;此代码可以做出来！！！更新
   
4. ; Concepts illustrated:
   
5. ;   - Reading a simple ascii file
   
6. ;   - Pretty printing
   
7. ;   - Rearranging data order via named dimension     
   
8. ;   - Weighting the data
   
9. ;   - Calculating EOFs and Principal Components (ie: time series)
   
10. ;   - Reconstructing the original array from EOFs and PCsby unweighting
    
11. ;   - Calculating 'sum-of-square' to verify normalization
    
12. ;   - Calculating cross correlations to verify that each is zero.
    
13. ;     This verifies that they are orthogonal.
    
14. ; =============================================================
    
15. ; John C Davis
    
16. ; Statistics and Data Analysis in Geology
    
17. ; Wiley, 2nd Edition, 1986
    
18. ; Source Data: page 524 , EOF results: page537
    
19. ; =============================================================
    
20. ; Cosine  weighting is performed
    
21. ; =============================================================
    
22. ;;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" ; not needed after 6.0.0
    
23. begin
    
24. ; ================================>  ; PARAMETERS
    
25.   ncol   = 7                           ; # columns (stations or grid points)
    
26.   nrow   = 25                          ; # time steps   
    
27.   neval  = 7                           ; # EOFs to calculate (max=ncol)
    
28.   xmsg   = -999.9                      ; missing value (_FillValue)
    
29.    
    
30.   ntim   = nrow                        ; # time steps   
    
31.   nlat   = ncol                        ; # latitudes  
    
32. ; ================================>  ; READ THE ASCII FILE
    
33.   fname  = "/disk1/dzz/data/test.txt"        ; test data from Davis text book
    
34.                                        ; open "data " as 2-dim array
    
35.   data   = asciiread (fname,(/ntim,nlat/), "float")
    
36.   data!0 = "time"                      ; name the dimensions for subsequent use
    
37.   data!1 = "lat"
    
38.   data@_FillValue = xmsg
    
39. 
    
40. ; ================================>   ; CREATE BOGUS LATITUDES & COSINE WEIGHTING
    
41. ;  lat    = ispan(-60, 60, 20)*1.0   
    
42. ;  lat!0  = "lat"
    
43. ;  lat@units = "degrees_north"
    
44. ;  data&lat  =  lat
    
45. ;  printVarSummary(lat)
    
46.   
    
47. ;  clat   = sqrt(cos(lat*0.01745329))  ; [*]; cosine weighting   
    
48. 
    
49. ; ================================>   ; PRETTY PRINT INPUT DATA
    
50. 
    
51.   opt = True
    
52.   opt@tspace = 5
    
53.   opt@title  = "Davis Input Data"
    
54.   write_matrix (data, (nlat+"f9.3"), opt)
    
55.    
    
56. ; ================================>   ; REORDER TO WHAT EOFUNC EXPECTS
    
57.    
    
58.   x        = data(lat | :,time | :)   ; reorder ... eofunc want 'time' as rightmost dimension
    
59.   printVarSummary(x)
    
60.   print("")
    
61. 
    
62. ;  CLAT     = conform(x, clat, 0)      ; create explicit array (not needed; done for illustration)
    
63. ;  printVarSummary(CLAT)
    
64. ;  print("")
    
65.    
    
66. ; ================================>   ; REORDER TO WHAT EOFUNC EXPECTS
    
67. 
    
68. ;  xw       = x*CLAT                   ; weight the observations; xw -> weighted
    
69. ;  copy_VarCoords(x, xw)
    
70. ;  printVarSummary(xw)
    
71. ;  print("")
    
72.    
    
73. ; ================================>   ; EOFs on weighted observationsS
    
74. 
    
75.   evecv    = eofunc_Wrap    (x,neval,False)
    
76.   evecv_ts = eofunc_ts_Wrap (x,evecv,False)
    
77. 
    
78.   print("")
    
79.   printVarSummary(evecv)
    
80.   print("")
    
81.   printVarSummary(evecv_ts)
    
82.   print("")
    
83. 
    
84. ; ================================>   ; PRETTY PRINT OUTPUT EOF RESULTS
    
85. 
    
86.   opt@title  = "Eigenvector components from weighted values: evecv"
    
87.   write_matrix (evecv(lat|:,evn|:), (nlat+"f9.3"), opt) ; reorder to match book
    
88.   print("")
    
89. 
    
90.   opt@title  = "Eigenvector time series from weighted values: evecv_ts"
    
91.   write_matrix (evecv_ts(time|:,evn|:), (nlat+"f9.3"), opt)
    
92.   print(" ")
    
93. 
    
94.   print("evecv_ts@ts_mean="+evecv_ts@ts_mean)
    
95.   print(" ")
    
96. 
    
97. ; ================================>   ; SUM OF THE SQUARES
    
98.                                       ; IF NORMALIZED, THEY SHOULD BE 1
    
99.   sumsqr = dim_sum(evecv^2)
    
100.   print("sum of squares: " + sumsqr)
     
101.   print(" ")
     
102. 
     
103. ; ================================>   ; RECONSTRUCT DATA:
     
104.                                       ; NCL WORKS ON ANOMALIES
     
105.   do n=0,neval-1                        
     
106.      evecv_ts(n,:) = evecv_ts(n,:) + evecv_ts@ts_mean(n)  ; add weighted means   
     
107.   end do                 
     
108. 
     
109.   xRecon = eof2data (evecv,evecv_ts)  ; RECONSTRUCTED WEIGHTED DATA
     
110.   copy_VarCoords(x, xRecon)
     
111. 
     
112.   xRecon = xRecon/CLAT
     
113. 
     
114.   printVarSummary(xRecon)
     
115.   print(" ")
     
116. 
     
117.   opt@title  = "Reconstructed data via NCL: default mode: constant offset"
     
118.   write_matrix (xRecon(time|:,lat|:), (nlat+"f9.3"), opt)
     
119.   print(" ")
     
120. 
     
121. ; ================================>   ; MAX ABSOLUTE DIFFERENCE
     
122.   
     
123.   mxDiff     = max(abs(xRecon-x))
     
124.   print("mxDiff="+mxDiff)
     
125.   print(" ")
     
126. 
     
127. ; ================================>   ; ORTHOGONALITY: DOT PRODUCT
     
128.                                       ; 0=>orthogonal
     
129.   do ne=0,neval-1
     
130.      EVECV = conform(evecv, evecv(ne,:), 1 )
     
131. 
     
132.      print("=====> dotp for col "+ne+"  <=====")
     
133.     do nn = 0, neval - 1
     
134.      dotp = dim_sum(evecv*EVECV)
     
135.      print("dotp patterns: " + dotp)
     
136.      print(" ")
     
137.     end do
     
138.   end do
     
139. 
     
140. ; ================================>   ; CORRELATION
     
141.                                                   
     
142.   do ne = 0, neval - 1
     
143.      corr = escorc(evecv_ts(ne,:),evecv_ts)
     
144.      print("corr patterns: " + corr)
     
145.      print(" ")
     
146.   end do
     
147. end
     
148. 
     

复制代码

我在NCL官网上找到站点资料的EOF脚本文件，数据也是用的一样的，但是我在xshell4运行之后总是莫名其妙的报错，脚本如下：
原数据如下：
3.760  3.660  0.540  5.275  9.768 13.741  4.782
  8.590  4.990  1.340 10.022  7.500 10.162  2.130
  6.220  6.140  4.520  9.842  2.175  2.732  1.089
  7.570  7.280  7.070 12.662  1.791  2.101  0.822
  9.030  7.080  2.590 11.762  4.539  6.217  1.276
  5.510  3.980  1.300  6.924  5.326  7.304  2.403
  3.270  0.620  0.440  3.357  7.629  8.838  8.389
  8.740  7.000  3.310 11.675  3.529  4.757  1.119
  9.640  9.490  1.030 13.567 13.133 18.519  2.354
  9.730  1.330  1.000  9.871  9.871 11.064  3.704
  8.590  2.980  1.170  9.170  7.851  9.909  2.616
  7.120  5.490  3.680  9.716  2.642  3.430  1.189
  4.690  3.010  2.170  5.983  2.760  3.554  2.013
  5.510  1.340  1.270  5.808  4.566  5.382  3.427
  1.660  1.610  1.570  2.799  1.783  2.087  3.716
  5.900  5.760  1.550  8.388  5.395  7.497  1.973
  9.840  9.270  1.510 13.604  9.017 12.668  1.745
  8.390  4.920  2.540 10.053  3.956  5.237  1.432
  4.940  4.380  1.030  6.678  6.494  9.059  2.807
  7.230  2.300  1.770  7.790  4.393  5.374  2.274
  9.460  7.310  1.040 11.999 11.579 16.182  2.415
  9.550  5.350  4.250 11.742  2.766  3.509  1.054
  4.940  4.520  4.500  8.067  1.793  2.103  1.292
  8.210  3.080  2.420  9.097  3.753  4.657  1.719
  9.410  6.440  5.110 12.495  2.446  3.103  0.914
我运行之后总是报错：
fatal:Variable (evecv) is undefined
fatal:["Execute.c":8565]:Execute: Error occurred at or near line 306
；我的脚本如下：
（不知道我的理解对不对，比如我有全国160站点降水资料的话，用该脚本文件运行之后得出的空间场和时间系数（层次*模态数），然后对时间系数再画xy-plot，对空间场先做NCL插值然后再做contour-plot，是这样的吗？）

1. <div>; ==============================================================
   
2. ; eof_0_640.ncl
   
3. ;
   
4. ; Concepts illustrated:
   
5. ;   - Reading a simple ascii file
   
6. ;   - Pretty printing
   
7. ;   - Rearranging data order via named dimension     
   
8. ;   - Calculating EOFs and Principal Components (ie: time series)
   
9. ;   - Reconstructing the original array from EOFs and PCs
   
10. ;   - Calculating 'sum-of-square' to verify normalization
    
11. ;   - Calculating cross correlations to verify that each is zero.
    
12. ;     This verifies that they are orthogonal.
    
13. ; =============================================================
    
14. ; This script shows how to use the new eofunc_n function added
    
15. ; in NCL V6.4.0 to avoid having to reorder the data before
    
16. ; calculating the EOFs.
    
17. ; =============================================================
    
18. ; John C Davis
    
19. ; Statistics and Data Analysis in Geology
    
20. ; Wiley, 2nd Edition, 1986
    
21. ; Source Data: page 524 , EOF results: page537
    
22. ; =============================================================
    
23. ; No weighting is performed
    
24. ; =============================================================
    
25. ; This file is loaded by default in NCL V6.2.0 and newer
    
26. ; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" ; not needed after 6.0.0
    
27. ; ================================>  ; PARAMETERS
    
28. begin
    
29.   ncol   = 7                           ; # columns (stations or grid points)
    
30.   nrow   = 25                          ; # time steps   
    
31.   neval  = 7                           ; # EOFs to calculate (max=ncol)
    
32.   xmsg   = -999.9                      ; missing value (_FillValue)
    
33.    
    
34.   ntim   = nrow                        ; # time steps   
    
35.   nsta   = ncol                        ; # stations ( or grid points)
    
36. ; ================================>  ; READ THE ASCII FILE
    
37. ;  dir    = "./"
    
38.   fname  = "/disk1/dzz/data/test.txt"        ; test data from Davis text book
    
39.                                        ; open "data " as 2-dim array
    
40.   data   = asciiread (fname,(/ntim,nsta/), "float")
    
41.   data!0 = "time"                      ; name the dimensions for subsequent use
    
42.   data!1 = "sta"
    
43.   data@_FillValue = xmsg
    
44. 
    
45. ; ================================>   ; PRETTY PRINT INPUT DATA
    
46. 
    
47.   opt = True
    
48.   opt@tspace = 5
    
49.   opt@title  = "Davis Input Data"
    
50.   write_matrix (data, (nsta+"f9.3"), opt)
    
51.    
    
52.   printVarSummary(data)
    
53. 
    
54.   print ("==========> EOFUNC_N <==========")
    
55.   evecv    = eofunc_Wrap  (data,neval,False,0)
    
56.   evecv_ts = eofunc_ts_Wrap (data,evecv,False,0)
    
57. 
    
58.   print("")
    
59.   printVarSummary(evecv)
    
60.   print("")
    
61.   printVarSummary(evecv_ts)
    
62.   print("")
    
63. 
    
64. ; ================================>   ; PRETTY PRINT OUTPUT EOF RESULTS
    
65. 
    
66.   opt@title  = "Eigenvector components: evecv"
    
67.   write_matrix (evecv(sta|,evn|:), (ncol+"f9.3"), opt) ; reorder to match book
    
68.   print("")
    
69. 
    
70.   opt@title  = "Eigenvector time series: evecv_ts"
    
71.   write_matrix (evecv_ts(time|,evn|:), (ncol+"f9.3"), opt)
    
72.   print(" ")
    
73. 
    
74.   print("evecv_ts@ts_mean="+evecv_ts@ts_mean)
    
75.   print(" ")
    
76. 
    
77. ; ================================>   ; SUM OF THE SQUARES
    
78.                                       ; IF NORMALIZED, THEY SHOULD BE 1
    
79.   sumsqr = dim_sum(evecv^2)
    
80.   print("sum of squares:     " + sumsqr)
    
81.   print(" ")
    
82. 
    
83. ; ================================>   ; RECONSTRUCT DATA  
    
84.                                       ; NCL WORKS ON ANOMALIES
    
85.   do n=0,neval-1                        
    
86.      evecv_ts(n,:) = evecv_ts(n,:) + evecv_ts@ts_mean(n)  ; add time series mean
    
87.   end do                 
    
88. 
    
89.   xRecon = eof2data_n (evecv,evecv_ts,0)  ; RECONSTRUCTED DATA
    
90.   copy_VarCoords(data, xRecon)
    
91. 
    
92.   printVarSummary(xRecon)
    
93.   print(" ")
    
94. 
    
95.   opt@title  = "Reconstructed data via NCL: default mode: constant offset"
    
96.   write_matrix (xRecon, (nsta+"f9.3"), opt)
    
97.   print(" ")
    
98. 
    
99. ; ================================>   ; MAX ABSOLUTE DIFFERENCE
    
100.   
     
101.   mxDiff     = max(abs(xRecon-data))
     
102.   print("mxDiff="+mxDiff)
     
103.   print(" ")
     
104. 
     
105. ; ================================>   ; ORTHOGONALITY: DOT PRODUCT
     
106.                                       ; 0=>orthogonal
     
107.   do ne=0,neval-1
     
108.      EVECV = conform(evecv, evecv(ne,:), 1 )
     
109. 
     
110.      print("=====> dotp for col "+ne+"  <=====")
     
111.     do nn = 0, neval - 1
     
112.      dotp = dim_sum(evecv*EVECV)
     
113.      print("dotp patterns: " + dotp)
     
114.      print(" ")
     
115.     end do
     
116.   end do
     
117. 
     
118. ; ================================>   ; CORRELATION
     
119.                                                   
     
120.   do ne = 0, neval - 1
     
121.      corr = escorc(evecv_ts(ne,:),evecv_ts)
     
122.      print("corr patterns: " + corr)
     
123.      print(" ")
     
124.   end do
     
125. end</div>

复制代码

原脚本是eof_0.ncl那个例子
一开始脚本文件用的是eofunc_n_Wrap和eofunc_ts_n_Wrap函数，但是我运行之后报错：
fatal:["Execute.c":8565]:Execute: Error occurred at or near line 177
Undefined identifier: (eofunc_n_Wrap) is undefined, can't continue
所以去查这个函数官网上说其适用于NCL6.4及其之后版本，所以我改用eofunc_Wrap和eofunc_Wrap,运行之后继续报错：
fatal:["Execute.c":8565]:Execute: Error occurred at or near line 556
Variable (evecv) is undefined
继续修改函数，将 evecv    = eofunc_Wrap  (data,neval,False，0)改为： evecv    = eofunc_Wrap  (data,neval,False)  （我没看懂官网上的“0”指的是什么意思），然后报错：
Undefined identifier: (eof2data_n) is undefined, can't continue
fatal:["Execute.c":8565]:Execute: Error occurred at or near line 711
，其实到此已经出现eof2data_n函数之前的结果：
Variable: evecv
Dimensions and sizes: [evn | 7] x [time | 25]
；
Variable: evecv_ts
Number of Dimensions: 2
Dimensions and sizes: [evn | 7] x [sta | 7]

但是奇怪的是，空间场和时间系数按理说应该分别是7*sta(7),和7*time（25），但是输出为什么会反着来呢，而且最终的结果实在老报错输出不来。

小弟才疏学浅，望各位不吝赐教，实在是纠结了好几天啊（如大家找不粗问题的话可以借鉴我的错误的流程）

孤独卡齐莫多

看不懂楼主的问题阿....

我是小猪

孤独卡齐莫多 发表于 2016-7-15 17:38
看不懂楼主的问题阿....

sorry啊，可能说的太多了，主要的问题就是第一个脚本文件导入数据之后出现一系列的问题，怎么也改不好，不知道为什么？

我是小猪

孤独卡齐莫多 发表于 2016-7-15 17:38
看不懂楼主的问题阿....

sorry啊，可能说的太多了，主要的问题就是第一个脚本文件导入数据之后出现一系列的问题，怎么也改不好，不知道为什么？

我是小猪

孤独卡齐莫多 发表于 2016-7-15 17:38
看不懂楼主的问题阿....

sorry啊，可能说的太多了，主要的问题就是第一个脚本文件导入数据之后出现一系列的问题，怎么也改不好，不知道为什么？

123出去玩

只有时间和纬度变量？没有经度可以做出来？

dovely

楼主你用的资料是 降水站点资料吗？那是要插值成格点才能继续EOF吗

我是小猪

dovely 发表于 2017-4-27 10:11
楼主你用的资料是 降水站点资料吗？那是要插值成格点才能继续EOF吗

例子是不需要插值的，不过过来快一年了我再回来看当初发的这个程序，太粗糙了，哈哈，当时还是太年轻，我这几天再重做看看

dovely

我是小猪 发表于 2017-4-27 13:00
例子是不需要插值的，不过过来快一年了我再回来看当初发的这个程序，太粗糙了，哈哈，当时还是太年轻，我 ...

哇！那楼主你看看之后跟我分享一下！我这两天正在学习怎么做站点资料的EOF呢！

我是小猪

dovely 发表于 2017-5-2 16:49
哇！那楼主你看看之后跟我分享一下！我这两天正在学习怎么做站点资料的EOF呢！

恩 我看完发现那个只是单纯的用站点进行EOF分解，分析出来前几个模态的站点之后，你再根据站点经纬度插值就行； 或者你先用站点数据插值再用eofunc函数进行EOF分解也是一样的