(保姆式教程:从下数据到画图)python如何利用EOF分析SSTA海温异常现象并画图
最近,在学习如何利用python中的EOF 对太平洋附近的1979-2004年出现的海温异常进行分析。
EOF分析是气象分析中常见的一种分析方法,也被称为经验正交函数。经过EOF分析,可以将几十年的海温数据变成几个空间模态和时间序列,这样就可以通过空间模态大致分析一些变化趋势,话不多说,接下来我们就开始看如何对SSTA进行EOF分解吧!
首先我们需要分析的数据是SSTA,我选取的是1979—2004年的海温数据,下载的网站是
https://www.metoffice.gov.uk/hadobs/hadisst/index.html
选取首页的main data page,进入后有很多可以选择的下载数据。
这里选取第一个文件下载,进入python进行运行。
首先我们读取这个nc文件。
#读取数据
path='C:\\Users\\user\\Desktop\\data\\SST\\HadISST_sst.nc'
SST=xr.open_dataset(path)
查看SST的基本信息:
我们可以看出该数据的时间范围为1840—2021年,但是我们需要分析的是1979—2004年的数据(此处为何选择1979为起点,是因为其实1979年之前的数据准确度都不太够,所以一般分析的时候选取1979作为起点分析。)
其中我们也可以发现我们的经度范围为-180-180,我们此处的分析范围是太平洋地区的厄尔尼诺和拉尼娜现象,所以我们的经纬度范围一定要足够准确才可以,此处我选取 latitude=slice(30,-30),longitude=slice(100,300),那么就会有一个事情需要做,就是我们需要使用cdo对该nc文件进行一个处理,将其中的经度范围从-180-180改为0-360。
关于cdo的内容可以学习这篇文章:cdo常用命令介绍
一些操作可以看这篇文章:如何解决 cdo转换经度-180~180 为0~360
但是操作过程中会发现一些问题,比如我这个数据集即使是用上述方法依然会报错,此时我们去看一下这个nc文件的基本信息
cdo infos HadISST_sst.nc
我们可以发现,在文章中需要将generic转化为lonlat的步骤在这里根本不需要,因为我们本来就有一层是lonlat ,所以我们只要将这一层lonlat取出来作为一个新的nc文件进行转化即可。
cdo selgrid,lonlat HadISST_sst.nc sst2.nc
这就取出来啦,此时进行上述文章中的操作:
cdo sellonlatbox,0,360,-90,90 sst2.nc sst3.nc
这样就成功转化啦!
接下来我们做一些进行EOF分析的准备工作:
首先进行EOF分析必须要安装eof的模块:
conda install -c conda-forge eofs
对数据进行一些处理:
path='C:\\Users\\user\\Desktop\\data\\SST\\sst1.nc'
SST=xr.open_dataset(path).sel(latitude=slice(30,-30),longitude=slice(100,300),time=slice("1979","2004"))
sst1=SST.sst[:]
sst2=np.array(sst1)
lat=SST.latitude[:]
lon=SST.longitude[:]
此处需要将sst转化为array格式才能进行下一步的矩阵运算。
关于为什么要做一定的计算呢?
因为我们需要分析的是海温异常,就需要分析与平均值不同的异常,所以需要将原来的数据和平均值做一个差值,通过差值的大小来判断海温异常的趋势以及分布。
sst=np.array(sst1)
ano=sst1.groupby('time.month')-sst1.groupby('time.month').mean('time', skipna=True)
ano1=np.array(ano)
得到的ano1就是我们要用来做EOF分析的数据集啦!
东西都准备好了,接下来就是我们的主要工作啦!
#计算纬度权重
lat=np.array(lat)
coslat=np.cos(np.deg2rad(lat))
wgts = np.sqrt(coslat)[..., np.newaxis]
#创建EOF分解器
solver=Eof(ano1,weights=wgts)
eof=solver.eofsAsCorrelation(neofs=4)
#此处的neofs的值是我们需要的空间模态数,比如这里我们打算画四个模态
pc=solver.pcs(npcs=4,pcscaling=1)#方差
var=solver.varianceFraction(neigs=4)
分析结束!
接下来就是画图啦!此处我们需要在一张图上画八个子图,左侧为空间模态,右边为时间序列。
fig=plt.figure(figsize=(15,15))#设置画布
proj=ccrs.PlateCarree(central_longitude=180)
leftlon,rightlon,lowerlat,upperlat=(100,290,-30,30)#设置经纬度范围
lon_formatter=ticker.LongitudeFormatter()
lat_formatter=ticker.LatitudeFormatter()
绘制第一模态
fig_ax1=fig.add_axes([0.1,0.95,0.5,0.3],projection=proj)
fig_ax1.set_extent([leftlon,rightlon,lowerlat,upperlat],crs=ccrs.PlateCarree())
fig_ax1.add_feature(cfeature.OCEAN,edgecolor='black')
fig_ax1.add_feature(cfeature.LAKES,alpha=0.5)
fig_ax1.add_feature(cfeature.COASTLINE,lw=1)
fig_ax1.set_xticks(np.arange(leftlon,rightlon,20),crs=ccrs.PlateCarree())
fig_ax1.set_yticks(np.arange(lowerlat,upperlat+5,5),crs=ccrs.PlateCarree())
fig_ax1.xaxis.set_major_formatter(lon_formatter)
fig_ax1.yaxis.set_major_formatter(lat_formatter)
rivers_110m = cfeature.NaturalEarthFeature('physical', 'rivers_lake_centerlines', '110m')
fig_ax1.set_title('(a) EOF1(HadISSTA from 1979-2004)',loc='left',fontsize =15)
fig_ax1.set_title( '%.2f%%' % (var[0]*100),loc='right',fontsize =15)
c1=fig_ax1.contourf(lon,lat, eof[0,:,:], levels=np.arange(-0.9,1.0,0.1 ), zorder=0, extend = 'both',transform=ccrs.PlateCarree(), cmap=plt.cm.RdBu_r)
绘制第一个时间序列:
fig_ax5=fig.add_axes([0.65,0.99,0.47,0.2])
fig_ax5.set_title('(b) PC1',loc='left',fontsize = 15)
fig_ax5.set_ylim(-3.5,3.5)
fig_ax5.axhline(0,linestyle="--")
fig_ax5.plot(np.arange(1979,2005,1/12),pc[:,0],color='blue')
绘制colorbar:
cbposition=fig.add_axes([0.1, 0.2, 0.5, 0.015])
fig.colorbar(c1,cax=cbposition,orientation='horizontal',format='%.1f')
其他模态则重复这些操作四次:(文章最末放出完整代码)
我们看一下成果图:
效果还不错,此时我想要把PC1、PC2、PC3绘制在一张图上,并且以三种不同的线条展示:
ax.plot(np.arange(1979,2005,1/12),pc[:,0],linewidth=1,linestyle='-',color='r',label='PC1')
bar=ax.bar(np.arange(1979,2005,1/12),height=pc[:,1],color='blue',align="center",width=0.1,linewidth=0.1,bottom=None,edgecolor='black',label='PC2')
ax.plot(np.arange(1979,2005,1/12),pc[:,2],linestyle='--',linewidth=1,color='black',label='PC3')
ax.set_ylim(-4,4)
ax.set_title("PC")
ax.set_xlabel("Time")
ax.set_ylabel("y")
plt.legend()
plt.grid()
plt.show()
完美,此时我们就可以对这些图像进行其他分析啦!
完整的代码如下:
import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
from eofs.standard import Eof
import xarray as xr
from cartopy.mpl import ticker#读取数据
path='C:\\Users\\user\\Desktop\\data\\SST\\sst1.nc'SST=xr.open_dataset(path).sel(latitude=slice(30,-30),longitude=slice(100,300),time=slice("1979","2004"))
sst1=SST.sst[:]
sst2=np.array(sst1)
lat=SST.latitude[:]
lon=SST.longitude[:]sst=np.array(sst1)
ano=sst1.groupby('time.month')-sst1.groupby('time.month').mean('time', skipna=True)
ano1=np.array(ano)#计算纬度权重
lat=np.array(lat)
coslat=np.cos(np.deg2rad(lat))
wgts = np.sqrt(coslat)[..., np.newaxis]
#创建EOF分解器
solver=Eof(ano1,weights=wgts)
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)
leftlon,rightlon,lowerlat,upperlat=(100,290,-30,30)
lon_formatter=ticker.LongitudeFormatter()
lat_formatter=ticker.LatitudeFormatter()
# 绘制第一模态
fig_ax1=fig.add_axes([0.1,0.95,0.5,0.3],projection=proj)
fig_ax1.set_extent([leftlon,rightlon,lowerlat,upperlat],crs=ccrs.PlateCarree())
fig_ax1.add_feature(cfeature.OCEAN,edgecolor='black')
fig_ax1.add_feature(cfeature.LAKES,alpha=0.5)
fig_ax1.add_feature(cfeature.COASTLINE,lw=1)
fig_ax1.set_xticks(np.arange(leftlon,rightlon,20),crs=ccrs.PlateCarree())
fig_ax1.set_yticks(np.arange(lowerlat,upperlat+5,5),crs=ccrs.PlateCarree())
fig_ax1.xaxis.set_major_formatter(lon_formatter)
fig_ax1.yaxis.set_major_formatter(lat_formatter)
rivers_110m = cfeature.NaturalEarthFeature('physical', 'rivers_lake_centerlines', '110m')
fig_ax1.set_title('(a) EOF1(HadISSTA from 1979-2004)',loc='left',fontsize =15)
fig_ax1.set_title( '%.2f%%' % (var[0]*100),loc='right',fontsize =15)
c1=fig_ax1.contourf(lon,lat, eof[0,:,:], levels=np.arange(-0.9,1.0,0.1 ), zorder=0, extend = 'both',transform=ccrs.PlateCarree(), cmap=plt.cm.RdBu_r)fig_ax2=fig.add_axes([0.1,0.7,0.5,0.3],projection=proj)
fig_ax2.set_extent([leftlon,rightlon,lowerlat,upperlat],crs=ccrs.PlateCarree())
fig_ax2.add_feature(cfeature.OCEAN,edgecolor='black')
fig_ax2.add_feature(cfeature.LAKES,alpha=0.5)
fig_ax2.add_feature(cfeature.COASTLINE,lw=1)
fig_ax2.set_xticks(np.arange(leftlon,rightlon,20),crs=ccrs.PlateCarree())
fig_ax2.set_yticks(np.arange(lowerlat,upperlat+5,5),crs=ccrs.PlateCarree())
fig_ax2.xaxis.set_major_formatter(lon_formatter)
fig_ax2.yaxis.set_major_formatter(lat_formatter)
rivers_110m = cfeature.NaturalEarthFeature('physical', 'rivers_lake_centerlines', '110m')
fig_ax2.set_title('(c) EOF2(HadISSTA from 1979-2004)',loc='left',fontsize =15)
fig_ax2.set_title( '%.2f%%' % (var[1]*100),loc='right',fontsize =15)
c2=fig_ax2.contourf(lon,lat, eof[1,:,:], levels=np.arange(-0.9,1.0,0.1), zorder=0, extend = 'both',transform=ccrs.PlateCarree(), cmap=plt.cm.RdBu_r)fig_ax3=fig.add_axes([0.1,0.45,0.5,0.3],projection=proj)
fig_ax3.set_extent([leftlon,rightlon,lowerlat,upperlat],crs=ccrs.PlateCarree())
fig_ax3.add_feature(cfeature.OCEAN,edgecolor='black')
fig_ax3.add_feature(cfeature.LAKES,alpha=0.5)
fig_ax3.add_feature(cfeature.COASTLINE,lw=1)
fig_ax3.set_xticks(np.arange(leftlon,rightlon,20),crs=ccrs.PlateCarree())
fig_ax3.set_yticks(np.arange(lowerlat,upperlat+5,5),crs=ccrs.PlateCarree())
fig_ax3.xaxis.set_major_formatter(lon_formatter)
fig_ax3.yaxis.set_major_formatter(lat_formatter)
rivers_110m = cfeature.NaturalEarthFeature('physical', 'rivers_lake_centerlines', '110m')
fig_ax3.set_title('(e) EOF3(HadISSTA from 1979-2004)',loc='left',fontsize =15)
fig_ax3.set_title( '%.2f%%' % (var[2]*100),loc='right',fontsize =15)
c3=fig_ax3.contourf(lon,lat, eof[2,:,:], levels=np.arange(-0.9,1.0,0.1), zorder=0, extend = 'both', transform=ccrs.PlateCarree(), cmap=plt.cm.RdBu_r)fig_ax4=fig.add_axes([0.1,0.2,0.5,0.3],projection=proj)
fig_ax4.set_extent([leftlon,rightlon,lowerlat,upperlat],crs=ccrs.PlateCarree())
fig_ax4.add_feature(cfeature.OCEAN,edgecolor='black')
fig_ax4.add_feature(cfeature.LAKES,alpha=0.5)
fig_ax4.add_feature(cfeature.COASTLINE,lw=1)
fig_ax4.set_xticks(np.arange(leftlon,rightlon,20),crs=ccrs.PlateCarree())
fig_ax4.set_yticks(np.arange(lowerlat,upperlat+5,5),crs=ccrs.PlateCarree())
fig_ax4.xaxis.set_major_formatter(lon_formatter)
fig_ax4.yaxis.set_major_formatter(lat_formatter)
rivers_110m = cfeature.NaturalEarthFeature('physical', 'rivers_lake_centerlines', '110m')
fig_ax4.set_title('(g) EOF4(HadISSTA from 1979-2004)',loc='left',fontsize =15)
fig_ax4.set_title( '%.2f%%' % (var[3]*100),loc='right',fontsize =15)
c4=fig_ax4.contourf(lon,lat, eof[3,:,:], levels=np.arange(-0.9,1.0,0.1), zorder=0, transform=ccrs.PlateCarree(), cmap=plt.cm.RdBu_r)cbposition=fig.add_axes([0.1, 0.2, 0.5, 0.015])
fig.colorbar(c1,cax=cbposition,orientation='horizontal',format='%.1f')fig_ax5=fig.add_axes([0.65,0.99,0.47,0.2])
fig_ax5.set_title('(b) PC1',loc='left',fontsize = 15)
fig_ax5.set_ylim(-3.5,3.5)
fig_ax5.axhline(0,linestyle="--")
fig_ax5.plot(np.arange(1979,2005,1/12),pc[:,0],color='blue')fig_ax6 = fig.add_axes([0.65, 0.74, 0.47, 0.2])
fig_ax6.set_title('(d) PC2',loc='left',fontsize = 15)
fig_ax6.set_ylim(-3.5,3.5)
fig_ax6.axhline(0,linestyle="--")
fig_ax6.plot(np.arange(1979,2005,1/12),pc[:,1],color='blue')fig_ax7 = fig.add_axes([0.65, 0.49, 0.47, 0.2])
fig_ax7.set_title('(f) PC3',loc='left',fontsize = 15)
fig_ax7.set_ylim(-3.5,3.5)
fig_ax7.axhline(0,linestyle="--")
fig_ax7.plot(np.arange(1979,2005,1/12),pc[:,2],color='blue')fig_ax8 = fig.add_axes([0.65, 0.24, 0.47, 0.2])
fig_ax8.set_title('(h) PC4',loc='left',fontsize = 15)
fig_ax8.set_ylim(-3.5,3.5)
fig_ax8.axhline(0,linestyle="--")
fig_ax8.plot(np.arange(1979,2005,1/12),pc[:,3],color='blue')
plt.show()fig=plt.figure(figsize=(10,6))
ax=fig.add_axes([0,0,1,1])ax.plot(np.arange(1979,2005,1/12),pc[:,0],linewidth=1,linestyle='-',color='r',label='PC1')
bar=ax.bar(np.arange(1979,2005,1/12),height=pc[:,1],color='blue',align="center",width=0.1,linewidth=0.1,bottom=None,edgecolor='black',label='PC2')
ax.plot(np.arange(1979,2005,1/12),pc[:,2],linestyle='--',linewidth=1,color='black',label='PC3')
ax.set_ylim(-4,4)
ax.set_title("PC")
ax.set_xlabel("Time")
ax.set_ylabel("y")
plt.legend()
plt.grid()
plt.show()
初次见面,请多关照!希望能解决你的一点小烦恼哦!
一个也也也也在努力学习python的ocean小菜鸟!
水平有限,欢迎指正!!!
欢迎评论、收藏、点赞、转发、关注。
关注我不后悔,记录学习进步的过程~~
(保姆式教程:从下数据到画图)python如何利用EOF分析SSTA海温异常现象并画图相关推荐
- Xilinx FPGA平台DDR3设计保姆式教程(3)MIG IP核使用教程及DDR读写时序
干货来了,用DDR搬砖,只需要会用IP就好,Xilinx官方YYDS! ---------------------------------------------------------------- ...
- Xilinx FPGA平台DDR3设计保姆式教程(1)DDR3基础简介
如果我们只是拿来用ddr搬砖,那么它就简单,知道IP怎么使用就好,但是要想知其所以然,理论知识是必备的,这也是我们初学者所欠缺的东西,慢慢修炼吧! 汇总篇: Xilinx平台DDR3设计保姆式教程(汇 ...
- Node.js的完全卸载与下载安装及各种npm、nvm、nrm配置(保姆式教程---提供全套安装包)---node.js的安装与配置(0)
Node.js的完全卸载与下载安装及各种npm.nvm.nrm配置(保姆式教程-提供全套安装包)-node.js的安装与配置(0) node的卸载 1.打开控制面板 我的电脑右键--->属性-- ...
- Node.js下载安装及各种npm、cnpm、nvm、nrm配置(保姆式教程—提供全套安装包)—nvm的安装与配置(4)
Node.js下载安装及各种npm.cnpm.nvm.nrm配置(保姆式教程-提供全套安装包)-cnpm的安装与配置(3) 五.nvm的下载安装 1.下载 nvm官网下载地址: https://git ...
- Node.js下载安装及各种npm、cnpm、nvm、nrm配置(保姆式教程---提供全套安装包)---npm的安装与配置(2)
Node.js下载安装及各种npm.cnpm.nvm.nrm配置(保姆式教程-提供全套安装包)-node.js的安装与配置(1) 三.配置npm安装的全局模块 需要配置的进行配置(不用C盘的配置,用C ...
- 在Oracle Linux上部署Yunzai Bot v3保姆式教程/甲骨文云/云崽Bot/原神
去我的博客查看本文:在Oracle Linux上部署Yunzai Bot v3保姆式教程 – 肚 (iocky.com) 本文也在Github与gitee可用. 初始配置 直接注册最低配置的就ok了, ...
- 圣诞节快到了,用Python给好友做一个圣诞树小程序吧【保姆式教程】
圣诞节快到了,用Python给好友做一个圣诞树小程序吧[保姆式教程] 马上圣诞节了,一个人的圣诞节可能会有些孤独,我来教你怎么用代码写一棵超级治愈的圣诞树. 话不多说,下面来看具体怎么实现吧! 文章目 ...
- Node.js下载安装及各种npm、cnpm、nvm、nrm配置(保姆式教程—提供全套安装包)—nrm的安装与配置(5)
Node.js下载安装及各种npm.cnpm.nvm.nrm配置(保姆式教程-提供全套安装包)-nvm的安装与配置(4) 一.nrm安装与使用 1.管理员运行cmd,输入如下,全局安装nrm: npm ...
- Node.js下载安装及各种npm、cnpm、nvm、nrm配置(保姆式教程—提供全套安装包)—cnpm的安装与配置(3)
Node.js下载安装及各种npm.cnpm.nvm.nrm配置(保姆式教程-提供全套安装包)-npm的安装与配置(2) 四.安装cnpm 1.管理员身份运行cmd,输入如下命令 npm instal ...
- Node.js下载安装及各种npm、nvm、nrm配置(保姆式教程---提供全套安装包)---node.js的安装与配置(1)
Node.js下载安装及各种npm.nvm.nrm配置(保姆式教程-提供全套安装包)-node.js的安装与配置(1) Node.js的完全卸载与下载安装及各种npm.nvm.nrm配置(保姆式教程- ...
最新文章
- Nginx之反向代理、日志格式、集群、缓存、压缩、URl 重写,读写分离配置
- 怎样把android应用部署到手机上
- linux maven .m2文件夹,Maven .m2文件夹创建(示例代码)
- 112页数学知识整理!机器学习-数学基础回顾.pptx
- 如何处理CRM_ORGMAN 300 error message
- html如何插入swf视频,Html插入SWF方法
- LeetCode 1860. 增长的内存泄露(等差数列)
- LeetCode 1630. 等差子数组
- [转载] Java中的命名参数
- 越狱解决iphone4s外放无声音
- 直接内存访问 (Direct Memory Access, DMA)
- init()函数何时运行?
- Swift的函数嵌套和返回内部函数
- 简要介绍无刷电机的基础知识
- 实战:模拟登录知乎网站(添加cookie)
- matlab矩阵 代表什么,matlab中矩阵AB是什么意思
- MCAL-GTM之时钟管理CMU
- 蚂蚁研究员玉伯:我的技术人生答案
- 爷青结!微软凌晨宣布“再见”!
- Python环境搭建与输入输出