【气象常用】EOF分解

最新推荐文章于 2024-06-03 13:31:43 发布
最新推荐文章于 2024-06-03 13:31:43 发布
阅读量584 收藏 2
点赞数 22
分类专栏： 气象常用 文章标签： python 数据分析
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/weixin_46921265/article/details/138540648
版权

效果图：

主要步骤：

1. 数据准备：我用的1950-2023年夏季的温度数据（更多数据下载方法可看：【数据下载】ERA5地表月平均数据下载_era5月平均数据哪里下-CSDN博客

【数据下载】CMIP6数据下载-CSDN博客

【数据下载】ERA5 各高度层月平均数据下载-CSDN博客）

2. 数据处理：计算了每年的夏季平均温度，及温度距平，并进行了数据淹没

3. 图像绘制：绘制了eof的前四个模态

详细代码：着急的直接拖到最后有完整代码

步骤一：导入库包及图片存储路径并设置中文字体为宋体，西文为新罗马（没有的库包要先下好奥）

###############################################################################
# 导入库及文件
from netCDF4 import Dataset
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import netCDF4 as nc
from cnmaps import get_adm_maps, draw_maps
from matplotlib import rcParams
import cartopy.feature as cfeature
from eofs.standard import Eof
from cartopy.mpl import ticker
config = {
            "font.family": 'serif',
            "font.size": 12,
            "mathtext.fontset": 'stix',
            "font.serif": ['SimSun'],
         }
rcParams.update(config)
rcParams['axes.unicode_minus']=False
###############################################################################

步骤二：读入数据

###############################################################################
# 读入文件及数据
datapath = r'H:/00.csdn/01data/'
figpath = r'H:/00.csdn/02fig/'
shppath = r'H:/00.csdn/04shp/cn_shp/Province_9/Province_9.shp'

data_tem = Dataset(datapath + 't2m_cn_jja_1950-2023.nc', mode='r')
# print(data_tem)
lons = data_tem.variables['longitude'][:]
lats = data_tem.variables['latitude'][:]
lon, lat = np.meshgrid(lons, lats)
time = data_tem.variables['time']
time_new = nc.num2date(time[:], time.units)
# print(time_new)
tem = data_tem.variables['t2m'][:]

data_tem.close()

步骤三：计算夏季温度平均及温度距平，并对数据进行掩膜

##############################################################################
# 计算每年的夏季平均
t_mean = np.empty((74, 601, 901))
for i in range(0, 74, 1):
    t_mean[i, :, :] = np.mean(tem[i*3 : (i+1)*3, : :], 0)
    
# 计算温度距平
dt = t_mean - np.mean(t_mean, 0)

# 掩膜数据
SC = get_adm_maps(province='四川省', only_polygon=True, record='first')
maskout_t = SC.maskout(lon, lat, dt)

步骤四：创建eof分解器（只需在Eof(maskout_t)中将maskout_t换成自己的数据就好啦）

###############################################################################
# 创建EOF分解器
solver = Eof(maskout_t)
eof = solver.eofsAsCorrelation(neofs=4)
pc = solver.pcs(npcs=4, pcscaling=1)
var = solver.varianceFraction(neigs=4)

步骤五：绘制eof的前四个空间模态

###############################################################################
# 绘制图像
fig=plt.figure(figsize=(15, 15))
proj=ccrs.PlateCarree(central_longitude=180)
lonW,lonE,latS,latN=(97, 109, 25, 35)

# 绘制空间模态
def sp_plot(value, location, title, var):
    fig_ax = fig.add_axes(location, projection=proj)
    fig_ax.set_extent([lonW, lonE, latS, latN], crs=ccrs.PlateCarree())
    fig_ax.set_xticks(list(range(lonW, lonE+1, 3)), crs=ccrs.PlateCarree())
    fig_ax.set_yticks(list(range(latS, latN+1, 3)), crs=ccrs.PlateCarree())
    fig_ax.add_feature(cfeature.OCEAN, edgecolor='black')
    fig_ax.add_feature(cfeature.LAKES, alpha=0.5)
    fig_ax.add_feature(cfeature.COASTLINE, lw=1)
    fig_ax.xaxis.set_major_formatter(ticker.LongitudeFormatter())
    fig_ax.yaxis.set_major_formatter(ticker.LatitudeFormatter())

    draw_maps(get_adm_maps(province='四川省', level='省'), color='k', linewidth=1)
    fig_ax.set_title(title, loc='left', fontsize =15)
    fig_ax.set_title('%.2f%%' % (var*100), loc='right', fontsize =15)
    c = fig_ax.contourf(lons, lats, value, levels=np.arange(-1, 1, 0.05), 
                        zorder=0, extend='both', transform=ccrs.PlateCarree(), cmap=plt.cm.RdBu_r)
    return c

f1 = sp_plot( eof[0,:,:], [0.1,0.99,0.47,0.2], '(a) EOF1', var[0])    
f2 = sp_plot( eof[1,:,:], [0.1,0.74,0.47,0.2], '(c) EOF2', var[1])   
f3 = sp_plot( eof[2,:,:], [0.1,0.49,0.47,0.2], '(e) EOF3', var[2])   
f4 = sp_plot( eof[3,:,:], [0.1,0.24,0.47,0.2], '(g) EOF4', var[3])   

cbposition=fig.add_axes([0.19, 0.2, 0.3, 0.015])
fig.colorbar(f4, cax=cbposition, orientation='horizontal', format='%.1f')

步骤六：绘制eof的前四个pc

# 绘制pc
def time_plot(value, location, title):
    fig0 = fig.add_axes(location)
    fig0.set_title(title, loc='left', fontsize = 15)
    fig0.set_ylim(-4, 4)
    fig0.axhline(0, linestyle="--")
    fig0.plot(np.arange(1950, 2024, 1), value, color='blue')
    # return fig0
    
f5 = time_plot(pc[:,0],[0.55, 0.99, 0.47, 0.2], '(b) PC1')
f6 = time_plot(pc[:,1],[0.55, 0.74, 0.47, 0.2], '(d) PC2')
f7 = time_plot(pc[:,2],[0.55, 0.49, 0.47, 0.2], '(f) PC3')
f8 = time_plot(pc[:,3],[0.55, 0.24, 0.47, 0.2], '(h) PC4')

步骤七：保存图像

###########################################################################
# 存图
plt.savefig(figpath+'207 温度 eof分解', dpi=600, bbox_inches = 'tight')
plt.show()

完整代码在这里：

###############################################################################
# 导入库及文件
from netCDF4 import Dataset
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import netCDF4 as nc
from cnmaps import get_adm_maps, draw_maps
from matplotlib import rcParams
import cartopy.feature as cfeature
from eofs.standard import Eof
from cartopy.mpl import ticker
config = {
            "font.family": 'serif',
            "font.size": 12,
            "mathtext.fontset": 'stix',
            "font.serif": ['SimSun'],
         }
rcParams.update(config)
rcParams['axes.unicode_minus']=False
###############################################################################
###############################################################################
# 读入文件及数据
datapath = r'H:/00.csdn/01data/'
figpath = r'H:/00.csdn/02fig/'
shppath = r'H:/00.csdn/04shp/cn_shp/Province_9/Province_9.shp'

data_tem = Dataset(datapath + 't2m_cn_jja_1950-2023.nc', mode='r')
# print(data_tem)
lons = data_tem.variables['longitude'][:]
lats = data_tem.variables['latitude'][:]
lon, lat = np.meshgrid(lons, lats)
time = data_tem.variables['time']
time_new = nc.num2date(time[:], time.units)
# print(time_new)
tem = data_tem.variables['t2m'][:]

data_tem.close()
##############################################################################
# 计算每年的夏季平均
t_mean = np.empty((74, 601, 901))
for i in range(0, 74, 1):
    t_mean[i, :, :] = np.mean(tem[i*3 : (i+1)*3, : :], 0)
    
# 计算温度距平
dt = t_mean - np.mean(t_mean, 0)

# 掩膜数据
SC = get_adm_maps(province='四川省', only_polygon=True, record='first')
maskout_t = SC.maskout(lon, lat, dt)
###############################################################################
# 创建EOF分解器
solver = Eof(maskout_t)
eof = solver.eofsAsCorrelation(neofs=4)
pc = solver.pcs(npcs=4, pcscaling=1)
var = solver.varianceFraction(neigs=4)
###############################################################################
# 绘制图像
fig=plt.figure(figsize=(15, 15))
proj=ccrs.PlateCarree(central_longitude=180)
lonW,lonE,latS,latN=(97, 109, 25, 35)

# 绘制空间模态
def sp_plot(value, location, title, var):
    fig_ax = fig.add_axes(location, projection=proj)
    fig_ax.set_extent([lonW, lonE, latS, latN], crs=ccrs.PlateCarree())
    fig_ax.set_xticks(list(range(lonW, lonE+1, 3)), crs=ccrs.PlateCarree())
    fig_ax.set_yticks(list(range(latS, latN+1, 3)), crs=ccrs.PlateCarree())
    fig_ax.add_feature(cfeature.OCEAN, edgecolor='black')
    fig_ax.add_feature(cfeature.LAKES, alpha=0.5)
    fig_ax.add_feature(cfeature.COASTLINE, lw=1)
    fig_ax.xaxis.set_major_formatter(ticker.LongitudeFormatter())
    fig_ax.yaxis.set_major_formatter(ticker.LatitudeFormatter())

    draw_maps(get_adm_maps(province='四川省', level='省'), color='k', linewidth=1)
    fig_ax.set_title(title, loc='left', fontsize =15)
    fig_ax.set_title('%.2f%%' % (var*100), loc='right', fontsize =15)
    c = fig_ax.contourf(lons, lats, value, levels=np.arange(-1, 1, 0.05), 
                        zorder=0, extend='both', transform=ccrs.PlateCarree(), cmap=plt.cm.RdBu_r)
    return c

f1 = sp_plot( eof[0,:,:], [0.1,0.99,0.47,0.2], '(a) EOF1', var[0])    
f2 = sp_plot( eof[1,:,:], [0.1,0.74,0.47,0.2], '(c) EOF2', var[1])   
f3 = sp_plot( eof[2,:,:], [0.1,0.49,0.47,0.2], '(e) EOF3', var[2])   
f4 = sp_plot( eof[3,:,:], [0.1,0.24,0.47,0.2], '(g) EOF4', var[3])   

cbposition=fig.add_axes([0.19, 0.2, 0.3, 0.015])
fig.colorbar(f4, cax=cbposition, orientation='horizontal', format='%.1f')

# 绘制pc
def time_plot(value, location, title):
    fig0 = fig.add_axes(location)
    fig0.set_title(title, loc='left', fontsize = 15)
    fig0.set_ylim(-4, 4)
    fig0.axhline(0, linestyle="--")
    fig0.plot(np.arange(1950, 2024, 1), value, color='blue')
    # return fig0
    
f5 = time_plot(pc[:,0],[0.55, 0.99, 0.47, 0.2], '(b) PC1')
f6 = time_plot(pc[:,1],[0.55, 0.74, 0.47, 0.2], '(d) PC2')
f7 = time_plot(pc[:,2],[0.55, 0.49, 0.47, 0.2], '(f) PC3')
f8 = time_plot(pc[:,3],[0.55, 0.24, 0.47, 0.2], '(h) PC4')

###########################################################################
# 存图
plt.savefig(figpath+'207 温度 eof分解', dpi=600, bbox_inches = 'tight')
plt.show()
就是99
关注
  • 22
    点赞
  • 踩
  • 2
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
专栏目录
tidemain_eo分解MATLAB_EOF分解_
10-01
对原始数据场进行EOF分解的函数，可以输出EOF分解所需要的各项数据
eof.f90.rar_EOF分解_eof_eof fortran_eof分析气象_降水
09-21
eof分析，关于气象方面的经验正交函数分析，eof.f90 作降水的eof分解
eof.rar_EOF分解_eof_气象
09-20
气象学进行EOF分解代码
python绘制梅尔谱图_推荐Python气象数据处理和绘图系列教程
weixin_39989033的博客
10-21 1051
“作者从NCL转向Python，期间发布了一些列Python教程，从相关工具的安装到气象常规数据处理和分析，由浅入深，通俗易懂。非常适合气象专业的人入门Python！”从NCL转向Python一年多了，首先感谢期间帮助我的大佬们，相对来说Python要比NCL的自学方便很多，因为使用人群要广，然而针对大气专业的一些特殊方法仍需自己摸索，也是因为疫情宅在家里，就把这一年多的学习经验总结整理...
气象数据可视化：NCL数据分析处理与图形绘制
weixin_49857990的博客
12-10 3494
NCL在数据的可视化方面具有世界一流水平，提供了丰富的图形接口，具有高度的可定制性。
气象数据分析之EOF分析以及python的实现
domotokoi的博客
01-11 2万+
import numpy as np import cartopy.crs as ccrs import cartopy.feature as cfeature import matplotlib.pyplot as plt from eofs.multivariate.standard import MultivariateEof from eofs.standard import Eof #计算权重：以纬度cos值得平方作为权重 coslat = np.cos(np.deg2rad(lat)) w
气象中EOF的理解
热门推荐
留下的，留不下的
06-17 2万+
  EOF方法在气象分析中经常用到，接触好长时间，也没有弄清楚他到底在干嘛。为什么这么做是可行的？这里面的计算都有什么意义呢？   看了一些矩阵理解方面的内容后，对EOF算是有了一些粗浅的理解，现在试着去阐述一下其中的院里，好厘清自己对EOF理解。 文章目录几篇参考文章EOF分析——Empirical Orthogonal Function analysisEOF的计算过程EOF的理解内积与投影线...
python--对站点数据做EOF并做插值绘制填色图
【更新于公众号 | 简谱学记】
09-23 2877
使用python对站点数据进行EOF处理，并使用两种插值方法到规则网格上。
（保姆式教程：从下数据到画图）python如何利用EOF分析SSTA海温异常现象并画图
shayuxing的博客
01-03 1万+
保姆式教程，教你如何利用python分析1979-2004年的太平洋地区的海温异常SSTA，并且绘制出空间模态和时间序列的图片。
EOF.zip_EOF Python_EOF方法_eof_eof分解python_python实现eof
07-14
EOF经验正交分解，是一种主成分分析方法用python语言实现，EOF计算得到了几个模态和对应的时间系数
EOF_信号分解_eof_主成分分析_
09-28
对具有时间维度的数据进行经验正交函数分解（EOF），也称主成分分析（PCA analysis）。对一段时间序列之内的信号进行分解，查看各类主要成分的占比。
InSAR处理及绘图常用GMT命令、bash、csh、matlab语法及样例
river__C的博客
04-18 4420
一两年的科研成果刚被正式接收，在此简单回顾一篇水paper所学和所做的一些东西，也为了以后方便翻阅查找。主要用到了GMTSAR进行SBAS-InSAR处理，中间bash及matlab处理数据，后期GMT成图。以下展示仅为数据处理及部分成图脚本，InSAR处理流程移步隔壁。文章链接以下内容都是基于明确处理目的简单学习结果，做为初学者写的脚本或许都很繁琐，望海涵。
matlab相关性分析及检验,科学网—EOF分析及Mann-Kendall检验相关的Matlab函数 - 刘永和的博文...
weixin_39743414的博客
03-16 516
本人长期以来使用的几个函数。1) EOF(正交经验函数)分解的Matlab函数：function [V,PC]=EOF(X)%V eigen vector (空间函数), V必须是pXp的，p为变量个数%PC 时间函数 Y=V'X%X 原始变量，列对应原始变量类型，行对应样本 原始变量的值都需要标准化处理，即减去均值，除以std%n=3904;[V,D]=eig(X*X'); % V特征向量 ...
如何mysql数据导入到mongdb
codemami的博客
05-30 1500
由于MySQL和MongoDB的数据模型不同（例如，MySQL使用关系模型，而MongoDB使用文档模型），你可能需要转换数据的格式。使用MongoDB驱动程序：你也可以使用MongoDB的官方驱动程序（如Python的pymongo）来编写脚本，将数据直接插入到MongoDB中。注意：如果你的JSON文件包含多个文档，并且它们不是作为数组的一部分（即每个文档都在其自己的行上），则需要使用--jsonArray选项。手动转换：对于小型数据集，你可以手动编辑SQL或CSV文件，将其转换为JSON格式。
专属编程笔记
u014397365的博客
06-01 976
在软件开发中，Utils（或Utilities）目录通常用于存放一些。这些工具类或函数为整个应用程序或多个模块提供便利的功能支持，使得代码更加模块化、易于维护和重用。UtilsUtils目录的设计目的在于避免代码重复，并确保在整个项目中以一致的方式处理共通任务。这样做可以提高开发效率，降低维护成本，并确保代码的整洁性和可读性。
在linux服务器上使用tensorboard，错误记录
wwwwzm的博客
05-30 1198
1. 使用tensorboard命令时，不是从虚拟环境中找tensorboard，而是从（全局路径）中找（/home/ljx/.local/lib/python3.9/site-packages/tensorboard）是一个在 Unix-like 系统（包括 Linux 和 macOS）的命令行界面（如 Bash shell）中使用的命令。2.使用which命令， 查看使用的tensorboard的路径，发现使用的是全局路径，不是虚拟环境路径。是一个特殊的变量，它定义了操作系统搜索可执行文件的目录。
基于NANO 9K 开发板加载PICORV32软核，并建立交叉编译环境
pocean2012的博客
05-30 1194
ISP的实现逻辑是出厂“芯片”（我们这里是软核）自带了BootLoader程序，即出厂引导程序，通过BootLoader可以将程序从串口（上位机）下载到Flash中，实际的时序是通过RST来区分正常启动还是烧录状态，然后上位机的烧录脚本来控制串口时序， 用户程序实际上是通过串口最终下载到了FLASH中，然后程序从flash启动。NANO 9K低成本体验FPGA开发，还是挺香的，官方例子有加载PICORV32软核的例子，可以建立简单的ISP编程环境，通过串口就可以加载程序，很方便。
Python自动化脚本：文件分类
svygh123的专栏
06-01 418
要使用Python根据文件类型进行分类，可以遍历指定目录下的所有文件，检查每个文件的扩展名，然后根据扩展名将文件移动到对应的文件夹中。
Python热重载调试新利器
最新发布
weixin_53707653的博客
06-03 421
Reloading是一个Python工具库，它让我们可以在每次迭代之前从源代码中重新加载（或函数）而不丢失任何当前已执行过程。该工具对于在深度学习模型训练期间编辑源代码尤其好用，它让我们可以添加日志记录、打印统计数据或保存模型，而无需重新启动训练。
eof分解python
01-28
EOF分解（Empirical Orthogonal Function decomposition）是一种主成分分析方法，用于分解多维数据集中的空间和时间变化模态。在气象和海洋学中，EOF分解常用于分析海表面温度、气压等大气和海洋场的变化模态。 以下是使用Python进行EOF分解的示例代码： ```python import numpy as np from scipy.linalg import svd # 假设有一个二维数据集，shape为 (n, m)，其中 n 是时间步数，m 是空间点数 data = np.random.rand(n, m) # 计算数据集的空间平均值 mean = np.mean(data, axis=0) # 将数据集减去空间平均值，得到去平均后的数据集 data_anomaly = data - mean # 计算数据集的协方差矩阵 covariance_matrix = np.cov(data_anomaly.T) # 对协方差矩阵进行奇异值分解 U, s, V = svd(covariance_matrix) # 提取前 k 个模态 k = 3 modes = U[:, :k] # 计算每个时间步的时间系数 time_coefficients = np.dot(data_anomaly, modes) # 合成前 k 个模态 reconstructed_data = np.dot(time_coefficients, modes.T) + mean # 打印结果 print("EOF modes:") print(modes) print("Time coefficients:") print(time_coefficients) print("Reconstructed data:") print(reconstructed_data) ``` 这段代码首先对数据集进行了去平均处理，然后计算了数据集的协方差矩阵，并对其进行了奇异值分解。接着，根据指定的模态数量 k，提取了前 k 个模态，并计算了每个时间步的时间系数。最后，根据时间系数和模态，合成了重构数据集。

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值