用Python做SVD分解(包含EOF)，包含源代码与例子（11/06/2021 更新）

蝶恋花

--------------------------------------------------11/06/2021--------------------------------------------------------------

修改了一些typo，加入了最近用xMCA发表的论文， 希望各位大侠可以点个star，谢谢！

-------------------------------------------------10/29/2019---------------------------------------------------------------

今天抽空把student-t 检验加入了xMCA，至此xMCA基本上算是完整了。

后续可能会添加蒙特卡罗检验或bootstrap， 但是这两者检验都需要非常大的计算量（假设全球有一万个格点，每个格点重复一万次，总的计算量就要1亿次了！）。如果添加，应该会把并行用上。

最后如果大家觉得对你平常的研究有帮助，还请大家多多star！谢谢！

代码都托管在github:https://github.com/Yefee/xMCA

API 指南: http://xmca.readthedocs.io/

下面给出一个太平洋海温与美国大陆地表气温的一个SVD例子。该例子来自于UW的大气科学专业的统计课程AS552。例子也在repo中，https://github.com/Yefee/xMCA/tree/master/xMCA/examples

下面是统计检验的栗子。傻瓜式地设置statistical_test为True就可以了。程序会返回lphet, rphet， 分别为左右场的p值。

1. le, re, lphet, rphet = sst_ts.heterogeneousPatterns(n=1, statistical_test=True)
   
2. fig, (ax1, ax2) = plt.subplots(2, 2, figsize=(12, 5))
   
3. le[0].plot(ax=ax1[0])
   
4. re[0].plot(ax=ax1[1])
   
5. 
   
6. # Only plot where p<0.01
   
7. lphet[0].where(lphet[0]<0.01).plot(ax=ax2[0])
   
8. rphet[0].where(rphet[0]<0.01).plot(ax=ax2[1])
   

复制代码

                               
登录/注册后可看大图

                               
登录/注册后可看大图

------------------------------------------------------------------------------------------------------------------------------

各位同学大家好，

最近有用到SVD分解，因此在xarray的基础上写了一个package。现在共享出来供大家使用。SVD也叫MCA，即最大协方差分析。（或者应该反过来说，MCA也叫SVD，个人感觉MCA在数学或者统计上意义更加直观准确）。下一次更新会加入蒙特卡洛检验。

如果大家觉得有帮助希望能够给个star，谢谢大家。

代码都托管在github:https://github.com/Yefee/xMCA

API 指南: http://xmca.readthedocs.io/

下面给出一个太平洋海温与美国大陆地表气温的一个SVD例子。该例子来自于UW的大气科学专业的统计课程AS552。例子也在repo中，https://github.com/Yefee/xMCA/tree/master/xMCA/examples

1. 导入需要用到的包：

1. from xMCA import xMCA
   
2. import xarray as xr
   
3. import matplotlib.pyplot as plt
   
4. %matplotlib inline
   

复制代码

2. 读取美国大陆地表气温以及太平洋海温：

1. usta = xr.open_dataarray('xMCA/examples/data/USTA.nc').transpose(*['time', 'lat', 'lon'])
   
2. usta.name = 'USTA'
   
3. print(usta)
   
4. sstpc = xr.open_dataarray('xMCA/examples/data/SSTPac.nc').transpose(*['time', 'lat', 'lon'])
   
5. sstpc.name = 'SSTPC'
   
6. print(sstpc)
   

复制代码

3. SVD分解，并获取第一与第二模态，其时间扩张系数以及协方差解释的百分比：

1. '''
   
2. decomposition, time should be in the first axis
   
3. lp is for SSTPC
   
4. rp is for USTA
   
5. '''
   
6. 
   
7. sst_ts = xMCA(sstpc, usta)
   
8. sst_ts.solver()
   
9. lp, rp = sst_ts.patterns(n=2) #第一与第二模态， 左场为太平洋海温，右场为气温
   
10. le, re = sst_ts.expansionCoefs(n=2) #时间扩张系数
    
11. frac = sst_ts.covFracs(n=2) # 协方差解释百分比
    
12. print(frac)
    

复制代码

4. 绘制第一模态及其时间扩张系数：

1. fig, (ax1, ax2) = plt.subplots(2, 2, figsize=(12, 5))
   
2. lp[0].plot(ax=ax1[0])
   
3. le[0].plot(ax=ax1[1])
   
4. 
   
5. rp[0].plot(ax=ax2[0])
   
6. re[0].plot(ax=ax2[1])
   

复制代码

                               
登录/注册后可看大图

5. 第一模态的同场相关与异场相关：

1. # 同场相关的图可见左场第一模态的时间扩展系数为ENSO模态，右场第一模态为北冷南暖的偶极模态
   
2. lh, rh = sst_ts.homogeneousPatterns(n=1)
   
3. # 异场相关可见美国气温在interannual与ENSO有很高的相关性，反之亦然。
   
4. le, re = sst_ts.heterogeneousPatterns(n=1)
   
5. fig, (ax1, ax2) = plt.subplots(2, 2, figsize=(12, 5))
   
6. lh[0].plot(ax=ax1[0])
   
7. rh[0].plot(ax=ax1[1])
   
8. 
   
9. le[0].plot(ax=ax2[0])
   
10. re[0].plot(ax=ax2[1])
    

复制代码

                               
登录/注册后可看大图

下面给出一个太平洋海温EOF的例子。其实很简单，把USTA变成SSTPC就是EOF了。

1. '''
   
2. decomposition, time should be in the first axis
   
3. lp is for SSTPC
   
4. rp is for USTA
   
5. '''
   
6. 
   
7. sst_ts = xMCA(sstpc, sstpc.rename('SSTPC_copy')) # EOF！
   
8. sst_ts.solver()
   
9. lp, _ = sst_ts.patterns(n=2)
   
10. le, _ = sst_ts.expansionCoefs(n=2)
    
11. frac = sst_ts.covFracs(n=2)
    
12. print(frac)

复制代码

第一模态与PC

1. fig, ax1 = plt.subplots(1, 2, figsize=(12, 5))
   
2. lp[0].plot(ax=ax1[0])
   
3. le[0].plot(ax=ax1[1])

复制代码

                               
登录/注册后可看大图

                               
登录/注册后可看大图

将第一个PC与原场做回归

1. lh, _ = sst_ts.homogeneousPatterns(n=1)
   
2. fig, ax1= plt.subplots()
   
3. lh[0].plot(ax=ax1)

复制代码

                               
登录/注册后可看大图

                               
登录/注册后可看大图

蝶恋花

实在不行你就把xMCA那个py文件拷贝到你要用的目录直接import就好了，一样用的

好久不见

赞赞赞！！！！！

木人

强

chongzika

给力啊，学习了

乱世一杯茶

高手，辟谷期高手，我有一个问题：SVD分解时候，计算的时候陆地温度，但是格点中有包含在海洋区域内的，如何计算？

zhtv587

学到了 目前我也在用python处理数据 有空可以一起交流一下！（有空可以向您请教一下）

young89

我是windows下，安装之后为什么出现这个错误啊
求大神指点

蝶恋花

young89 发表于 2019-6-7 21:26
我是windows下，安装之后为什么出现这个错误啊
求大神指点

Not sure if it is fixed.
Sorry for the inconvenience. I tried to fix it, but I don't have a windows machine and cannot test it.

蝶恋花

乱世一杯茶 发表于 2019-6-7 16:26
高手，辟谷期高手，我有一个问题：SVD分解时候，计算的时候陆地温度，但是格点中有包含在海洋区域内的，如 ...

It doesn't matter. Left and right field can be anything, given the same time dimension.

勿忘草

young89 发表于 2019-6-7 21:26
我是windows下，安装之后为什么出现这个错误啊
求大神指点

请问你是怎么装到windows下的，我试了好多次都安装不成功