MeteoInfoLab脚本示例：EOF分析

MeteoInfo

经验正交函数分析方法(empirical orthogonal function, 缩写为EOF)，也称特征向量分析(eigenvector analysis)，或者主成分分析(principal component analysis,缩写PCA)，是一种分析矩阵数据中的结构特征，提取主要数据特征量的一种方法。Lorenz在1950年代首次将其引入气象和气候研究，现在在地学及其他学科中得到了非常广泛的应用。地学数据分析中通常特征向量对应的是空间样本，所以也称空间特征向量或者空间模态；主成分对应的是时间变化，也称时间系数。因此地学中也将EOF分析称为时空分解。详见此文档： EOFanalysis.pdf (683.36 KB, 下载次数: 139)

2017-10-7 07:25 上传

点击文件名下载附件


论坛里有很多EOF分析的资料和帖子，这里示例用MeteoInfoLab脚本进行EOF分析。示例数据来自此贴：http://bbs.06climate.com/forum.php?mod=viewthread&tid=4753 ，是62年某月的全球温度资料，数据排放按（71*144,62）。由于格点数（71*144）远大于时间序列长度（62），进行时空转换能够极大的提高运算速度。代码里主要进行了矩阵的运算，EOF.T为EOF矩阵的转置，dot函数是矩阵相乘，eig函数获得矩阵的特征值和特征向量（注意运算结果和MatLab略有不同，但本质是一样的）。

1. fn = 'C:/Temp/EOF/test.txt'
   
2. ny = 71
   
3. nx = 144
   
4. m = ny * nx
   
5. n = 62
   
6. ss1 = asciiread(fn, shape=(71,144,n))
   
7. ss1 = ss1[::-1,::1,:]
   
8. X = ss1.reshape(ny * nx, n)
   
9. #Time-space transform
   
10. C = dot(X.T, X)
    
11. E1, EOF1 = linalg.eig(C)
    
12. EOF1 = EOF1.T[:,::-1]
    
13. E = diag(E1[::-1])
    
14. EOFa = dot(X, EOF1)
    
15. EOF = zeros((m,n))
    
16. for i in range(n):
    
17.     EOF[:,i] = EOFa[:,i]/sqrt(E[i,i])
    
18. PC = dot(EOF.T, X)
    
19. PC = PC[::-1,:]
    
20. eof1 = EOF[:,0].reshape(ny, nx)
    
21. 
    
22. #Plot
    
23. axesm()
    
24. lon = linspace(0, 360, nx)
    
25. lat = linspace(-90, 90, ny)
    
26. geoshow('country', edgecolor='k')
    
27. levs = arange(-0.02, 0.021, 0.002)
    
28. layer = contourfm(lon, lat, eof1, levs, smooth=False)
    
29. colorbar(layer)
    
30. title('EOF mode 1')

fn = 'C:/Temp/EOF/test.txt'<br /> ny = 71<br /> nx = 144<br /> m = ny * nx<br /> n = 62<br /> ss1 = asciiread(fn, shape=(71,144,n))<br /> ss1 = ss1[::-1,::1,:]<br /> X = ss1.reshape(ny * nx, n)<br /> #Time-space transform<br /> C = dot(X.T, X)<br /> E1, EOF1 = linalg.eig(C)<br /> EOF1 = EOF1.T[:,::-1]<br /> E = diag(E1[::-1])<br /> EOFa = dot(X, EOF1)<br /> EOF = zeros((m,n))<br /> for i in range(n):<br />     EOF[:,i] = EOFa[:,i]/sqrt(E[i,i])<br /> PC = dot(EOF.T, X)<br /> PC = PC[::-1,:]<br /> eof1 = EOF[:,0].reshape(ny, nx)<br /> <br /> #Plot<br /> axesm()<br /> lon = linspace(0, 360, nx)<br /> lat = linspace(-90, 90, ny)<br /> geoshow('country', edgecolor='k')<br /> levs = arange(-0.02, 0.021, 0.002)<br /> layer = contourfm(lon, lat, eof1, levs, smooth=False)<br /> colorbar(layer)<br /> title('EOF mode 1')

给meteo模块增加了eof和varimax函数，eof函数用于做EOF分析，有是否进行时空转换的选项。varimax是进行矩阵varimax旋转的函数，可用于进行REOF（旋转EOF）分析。下面是一个例子，计算结果未进行严格验证，如果有错欢迎指出。

1. fn = 'C:/Temp/EOF/test.txt'
   
2. ny = 71
   
3. nx = 144
   
4. m = ny * nx
   
5. n = 62
   
6. ss1 = asciiread(fn, shape=(71,144,n))
   
7. ss1 = ss1[::-1,::1,:]
   
8. X = ss1.reshape(ny * nx, n)
   
9. 
   
10. #EOF calculation - Time-space transform
    
11. EOF, E, PC = meteo.eof(X, transform=True)
    
12. eof1 = EOF[:,0].reshape(ny, nx)
    
13. 
    
14. #REOF calculation using varimax rotation
    
15. REOF = meteo.varimax(EOF[:,:5])[0]
    
16. reof1 = REOF[:,0].reshape(ny, nx)
    
17. 
    
18. #Plot
    
19. lon = linspace(0, 360, nx)
    
20. lat = linspace(-90, 90, ny)
    
21. #EOF mode 1
    
22. subplot(2,1,1)
    
23. axesm(newaxes=False)
    
24. geoshow('country', edgecolor='k')
    
25. levs = arange(-0.02, 0.021, 0.002)
    
26. layer = contourfm(lon, lat, eof1, levs, smooth=False)
    
27. colorbar(layer)
    
28. title('EOF mode 1')
    
29. #REOF mode 1
    
30. subplot(2,1,2)
    
31. axesm(newaxes=False)
    
32. geoshow('country', edgecolor='k')
    
33. levs = arange(-0.02, 0.021, 0.002)
    
34. layer = contourfm(lon, lat, reof1, levs, smooth=False)
    
35. colorbar(layer)
    
36. title('REOF mode 1')

fn = 'C:/Temp/EOF/test.txt'<br /> ny = 71<br /> nx = 144<br /> m = ny * nx<br /> n = 62<br /> ss1 = asciiread(fn, shape=(71,144,n))<br /> ss1 = ss1[::-1,::1,:]<br /> X = ss1.reshape(ny * nx, n)<br /> <br /> #EOF calculation - Time-space transform<br /> EOF, E, PC = meteo.eof(X, transform=True)<br /> eof1 = EOF[:,0].reshape(ny, nx)<br /> <br /> #REOF calculation using varimax rotation<br /> REOF = meteo.varimax(EOF[:,:5])[0]<br /> reof1 = REOF[:,0].reshape(ny, nx)<br /> <br /> #Plot<br /> lon = linspace(0, 360, nx)<br /> lat = linspace(-90, 90, ny)<br /> #EOF mode 1<br /> subplot(2,1,1)<br /> axesm(newaxes=False)<br /> geoshow('country', edgecolor='k')<br /> levs = arange(-0.02, 0.021, 0.002)<br /> layer = contourfm(lon, lat, eof1, levs, smooth=False)<br /> colorbar(layer)<br /> title('EOF mode 1')<br /> #REOF mode 1<br /> subplot(2,1,2)<br /> axesm(newaxes=False)<br /> geoshow('country', edgecolor='k')<br /> levs = arange(-0.02, 0.021, 0.002)<br /> layer = contourfm(lon, lat, reof1, levs, smooth=False)<br /> colorbar(layer)<br /> title('REOF mode 1')

遥感卫星图像

谢谢楼主分享，能做一个REOF的例子就更好了，期待。谢谢老师

又是那隻貓

王老师，厉害了

river

赶紧来学习，支持王老师！

MeteoInfo

rcdxph 发表于 2017-10-8 12:22
老师，能否做一个 REOF的 例子啊，好像REOF比EOF 更广范啊！

最近比较忙，有空了再弄吧。

板板板板转

万年难题得以解决

板板板板转

万年难题得以解决

lgz

谢谢楼主分享

杜小卜

期待老师做一个REOF 的分析案例。。

气象wzl

谢谢楼主分享