用python 做 REOF

qq492947833

最近鄙人在做自己的论文的时候发现一个问题，python目前似乎没有方法做reof，只能做eof。但是众所周知，在我们分析降水、气温分区时，reof是重要工具。如果不能做reof，那就失去了在做分区研究时候的优势。
  在气象家园python模块输入reof，发现讲reof的并没有多少，只有两个帖子，一个是提问题的帖子，另外一个代码量过多，过于复杂，难以理解。
  
  但是我偶然之间发现了华裔教授郑中华在github上的帖子（网址：https://github.com/zzheng93/pyEOF），发现其编写了reof代码，可以很好的做reof。但是其API文档（网址：https://pyeof.readthedocs.io/en/latest/notebooks/basic_usage.html#EOF-Analysis）有很多多余的代码，不够简洁，而且还有许多代码没有写在API中，会导致报错。于是我发邮件联系了郑中华教授，跟他深入交流了reof程序问题，终于在这几天有了自己的理解，并且做出来了自己的reof图。
  首先来说必要函数库：
  numpy、pandas这两个库都可以通过conda install命令安装，并且大部分朋友已经安装了这两个库。新引入的库有rpy2，因为程序用到了r语言，该库用conda install -c conda-forge rpy2安装。还有pyEOF库，用于做reof，用pip install pyEOF命令安装。当然要自己安装r语言，在https://mirrors.tuna.tsinghua.edu.cn/CRAN/上可以自行下载。
  #安装完之后，需要在python中先运行命令，


import os

os.environ['R_HOME'] = r'D:\R_HOME\R-4.1.0'


#由于我的R语言安装在D盘的R_HOME文件夹下的R-4.1.0文件夹下，所以我这么输入，如果大家的路径不同，需要更改为自己的安装路径。
#不执行该命令的话pyEOF库导入不进来！


#接下来是引入函数库，

import numpy as np

import pandas as pd

import rpy2

from pyEOF import df_eof,get_time_space


#只需要这4个库就可以。
#然后就是代码部分，

latlong = v.shape[1]

lonlong = v.shape[2]


#v是我们的变量，其必须要是xr.Dataarray格式。做这一步的原因是我们后续要把变量还原成二维，方便绘图。


vs = v.to_dataframe([v_name]).reset_index()

#由于郑教授的reof程序只能处理pd.DataFrame格式数据，所以要先把变量v转化为dataframe格式，其中v_name为变量名称，填写字符串类型，如'air'，最好与nc文件中的变量名称一致。

vs.columns = [timename,latname,lonname,v_name]


#这一行规定了vs的每一列名称，其中timename是时间维名称，如'time';latname是纬度维名称，如'lat';lonname是经度维名称，如'lon';v_name是变量名称，与上面的v_name相同。timename，latname和lonname最好也与nc文件中各个维度名称一致。

vss = get_time_space(vs, time_dim = timename, lumped_space_dims = [latname,lonname])


#这一行为了把三维的vs压缩成二维，因为郑教授的reof程序只能做二维reof。其把经纬度压缩了，其中timename，latname和lonname与上面的一致。

pca = df_eof(vss,pca_type="varimax",n_components=eofnum)


eof = np.array(pca.eofs(s=2, n=eofnum))


pc = np.array(pca.pcs(s=2, n=eofnum))  


var = np.array(pca.evf(n=v.shape[0]))


#这里开始做reof，第一行做reof，第二行取空间模态，第三行取时间系数，第四行取方差贡献。其中pca_type='varimax'不用修改，eofnum代表依照前几个模态做旋转，填写整数，如10代表依照前10个模态做旋转。s=2也不需要修改。eof就是空间场，pc是时间系数，var是方差贡献。


eof = eof[eofmodel-1,:]
pc = pc[:,eofmodel-1]

#通过这种方法可以单独把一个模态的空间场和时间系数取出来，其中eofmodel就是你要取的模态，例如eofmodel为1代表取第1模态，以此类推。但是eofmodel不能大于eofnum。


eof = eof.reshape(latlong,lonlong)


#最后通过reshape把eof重新变成二维变量，方便绘图。


print(var[eofmodel-1]/np.nansum(var))

#还可以通过这种方式把方差贡献输出来，其中eofmodel与上述含义相同。


最后得到的eof是一个二维场，具体到某一模态的空间场，pc是一个一维场，具体到某一模态的时间系数（如下图为eofnum取10，eofmodel取5，即以前十模态为依据进行旋转，取第五模态绘图的空间场和时间系数图）。

当然我自己也在写自动绘图函数（网址：http://bbs.06climate.com/forum.php?mod=viewthread&tid=97940&extra=&page=1），也有视频教程（网址：https://www.bilibili.com/video/BV1d341167kM/，大概从21分08秒开始讲的reof从安装库到编写代码到测试），大家感兴趣可以去看看。

                               
登录/注册后可看大图

补充：评论区有网友@zhuyidong 说明了安装pyEOF时不要有中文路径，感谢补充。不过本人现在在用python时基本上不用中文路径，很多中文路径，包括使用xarray读取nc时，如果有中文路径有时也会报错。所以建议大家使用python时尽量减少中文路径的使用，从而避免不必要的报错。

zhuyidong

qq492947833 发表于 2021-8-7 13:45
感谢补充！的确python中对中文路径很不友好，很多nc文件如果有中文路径可能也会报错。

是的，为了做reof，配了r环境，运行时报了这个错，就是中文路径问题

qq492947833

视频还在审核，估计一会就能看到了

一大碗年糕

林夕

Alliris

谢谢楼主分享

看云淡风清

您分享的内容都很实用

zhuyidong

我昨天也试了一下，正好看看结果是不是一样的，感谢楼主的分享
还有大家再安装pyEOF不要有中文路径，不然会报错

kmdqlb

qq492947833 发表于 2021-8-2 21:12
视频还在审核，估计一会就能看到了

老师，你好，谢谢你的python自绘图视频，因为很多小伙伴都是python初学者，希望在百度网盘之类的地方，共享一下程序中所用到所用到数据，这样才可以实际操作老师的程序，也才知道数据应该是的时空维度怎么排列，以便日后工作中使用该函数，再次感谢老师的辛勤付出，谢谢！

VV09

真不戳！~~~！！

岫鹿森涌

牛啊牛啊牛啊，马着学习了