请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

新浪微博登陆

只需一步, 快速开始

QQ登录

只需一步,快速开始

搜索
查看: 702|回复: 21

[源代码] 用Python做SVD分解(包含EOF),包含源代码与例子

[复制链接] |关注本帖

新浪微博达人勋

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
本帖最后由 蝶恋花 于 2019-6-8 07:57 编辑

各位同学大家好,

最近有用到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. sst_ts = xMCA(sstpc, usta)
  7. sst_ts.solver()
  8. lp, rp = sst_ts.patterns(n=2) #第一与第二模态, 左场为太平洋海温,右场为气温
  9. le, re = sst_ts.expansionCoefs(n=2) #时间扩张系数
  10. frac = sst_ts.covFracs(n=2) # 协方差解释百分比
  11. 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. rp[0].plot(ax=ax2[0])
  5. 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. le[0].plot(ax=ax2[0])
  9. 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. sst_ts = xMCA(sstpc, sstpc.rename('SSTPC_copy')) # EOF!
  7. sst_ts.solver()
  8. lp, _ = sst_ts.patterns(n=2)
  9. le, _ = sst_ts.expansionCoefs(n=2)
  10. frac = sst_ts.covFracs(n=2)
  11. 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)
复制代码

                               
登录/注册后可看大图


                               
登录/注册后可看大图

评分

参与人数 2金钱 +40 收起 理由
rceclx + 20 赞一个!
chongzika + 20 赞一个!

查看全部评分

已有1人关注本帖

zhtv587
密码修改失败请联系qq:937062711

新浪微博达人勋

发表于 2019-6-7 08:35:47 | 显示全部楼层 |取消关注该作者的回复
赞赞赞!!!!!
密码修改失败请联系qq:937062711
回复

使用道具 举报

新浪微博达人勋

发表于 2019-6-7 09:14:05 | 显示全部楼层 |取消关注该作者的回复
密码修改失败请联系qq:937062711
回复

使用道具 举报

新浪微博达人勋

发表于 2019-6-7 10:03:22 | 显示全部楼层 |取消关注该作者的回复
给力啊,学习了
密码修改失败请联系qq:937062711
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2019-6-7 16:26:31 | 显示全部楼层 |取消关注该作者的回复
高手,辟谷期高手,我有一个问题:SVD分解时候,计算的时候陆地温度,但是格点中有包含在海洋区域内的,如何计算?
密码修改失败请联系qq:937062711
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2019-6-7 18:35:05 | 显示全部楼层 |取消关注该作者的回复
学到了 目前我也在用python处理数据 有空可以一起交流一下!(有空可以向您请教一下)
密码修改失败请联系qq:937062711
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2019-6-7 21:26:28 | 显示全部楼层 |取消关注该作者的回复
我是windows下,安装之后为什么出现这个错误啊
求大神指点
微信图片_20190607212509.png
密码修改失败请联系qq:937062711
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2019-6-8 06:30:21 | 显示全部楼层 |取消关注该作者的回复
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.
密码修改失败请联系qq:937062711
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2019-6-8 06:34:01 | 显示全部楼层 |取消关注该作者的回复
乱世一杯茶 发表于 2019-6-7 16:26
高手,辟谷期高手,我有一个问题:SVD分解时候,计算的时候陆地温度,但是格点中有包含在海洋区域内的,如 ...

It doesn't matter. Left and right field can be anything, given the same time dimension.
密码修改失败请联系qq:937062711
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2019-6-8 16:57:08 | 显示全部楼层 |取消关注该作者的回复
young89 发表于 2019-6-7 21:26
我是windows下,安装之后为什么出现这个错误啊
求大神指点

请问你是怎么装到windows下的,我试了好多次都安装不成功
密码修改失败请联系qq:937062711
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表