代码收藏家技术教程 2024-11-30

python气象学习笔记（七）：EOF分析运用及代码分析

1.EOF介绍

经验正交函数分析方法（empirical orthogonal function，缩写为EOF），也称特征向量分析（eigenvector analysis），或者主成分分析（principal component annlysis，缩写为PCA），是一种分析矩阵数据中的结构特征，提取主要数据特征量的一种方法。特征向量对应的是空间样本，也称空间特征向量或者空间模态（EOF），在一定程度上反映要素场的空间分布特点；主成分对应的是时间变化，也称时间系数（PC），反映相应空间模态随时间的权重变化。

2.所需要的库

conda install -c conda-forge eofs

库的使用：

from eofs.standard import Eof

3.代码编写

1.导入部分

import numpy as np
import xarray as xr
from eofs.standard import Eof
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import cartopy.io.shapereader as shpreader
import cartopy.mpl.ticker as cticker
import warnings
warnings.filterwarnings("ignore", category=DeprecationWarning)

NumPy (np)：用于数组和数值计算的库。

Xarray (xr)：用于高效处理多维数组数据的库，比如 NetCDF 文件。

eofs.standard.Eof：这个模块用于执行 EOF 分析。

Matplotlib (plt)：用于绘图的基础工具。

make_axes_locatable：用于创建附加轴，比如在地图上绘制小地图。

Cartopy 相关：用来进行地理空间数据可视化，包含坐标参考系统（ccrs）、地理特征（cfeature）、格式化工具（LongitudeFormatter, LatitudeFormatter）、形状文件读取和格式化工具。

Warnings：忽略弃用警告，以保持输出的清晰。

2.数据加载

plt.rcParams['font.sans-serif'] = ['Times New Roman']
plt.rcParams['axes.unicode_minus'] = False

f = xr.open_dataset(r'...\hgt.mon.mean.nc')
hgt_Eu = f.hgt.loc[f.time.dt.month.isin([6]), 850].loc['1948':'2008']
coslat = np.cos(np.deg2rad(hgt_Eu.lat.values))
solver = Eof(hgt_Eu.values, weights=coslat[..., np.newaxis])

设置 Matplotlib 字体和负号：确保图表中能够正确显示负号和选定的字体格式。

数据集加载：从指定路径加载 NetCDF 数据集。

数据选择：获取 1948 年到 2008 年 6 月的特定纬度和高度（850 hPa）的数据。

纬度权重：基于纬度计算权重，提高 EOF 分析的准确性。

EOF 分析：对数据进行 EOF 分析，以提取主要组件。

其中对于下面代码：

hgt_Eu = f.hgt.loc[f.time.dt.month.isin([6]), 850].loc['1948':'2008']

f.hgt: 这里 f 是一个通过 xarray.open_dataset() 打开的数据集，该数据集中的一个变量是 hgt，可能代表大气高度场。
loc[…] 选择器:xarray 提供了类似 pandas 的 loc 索引器，可以过滤数据集中的维度。
f.time.dt.month.isin([6]): f.time 代表时间维度的索引，dt 是日期时间的访问器，用于从时间维度中提取各个时间单位（月）。
通过 .isin([6]) 方法筛选出所有时间维度为6月的数据，这步操作限制了提取出的数据仅为6月份的。
, 850: 这一部分表示选择高度场数据中，压强层为 850 hPa 的数据。这通常用于大气科学领域，850 hPa 是对流层中一个常用的等压面，因为它通常位于天气系统的重要活动层。
.loc[‘1948’:‘2008’]: 这是对时间进行进一步筛选，只提取从 1948 年到 2008 年之间的数据。xarray 支持基于时间序列的索引，这使得可以很容易从数据集中提取某个时间段内的数据。

coslat = np.cos(np.deg2rad(hgt_Eu.lat.values))

hgt_Eu.lat.values:
这是从之前步骤中提取的气象数据的纬度值。xarray 对象的 .values 属性将这些数据提取为 numpy 数组。
np.deg2rad(…):
np.deg2rad() 函数用于将角度从度转换为弧度，因为 numpy 的三角函数如 cos 是基于弧度的。
np.cos(…):
计算弧度值的余弦。余弦函数用于得到纬度的权重。

3.提取信息

eof = solver.eofsAsCorrelation(neofs=3)
pc = solver.pcs(npcs=3, pcscaling=1)
var = solver.varianceFraction()

获取模态场：提取出前三个 EOF 模态。

PC 时间序列：获取前三个主成分的时间序列。

方差解释能力：计算每个 EOF 模态的方差贡献。

4.定义颜色

color1, color2, color3 = [], [], []
for i in range(1948, 2009):
    color1.append('red' if pc[i - 1948, 0] >= 0 else 'blue')
    color2.append('red' if pc[i - 1948, 1] >= 0 else 'blue')
    color3.append('red' if pc[i - 1948, 2] >= 0 else 'blue')

颜色标记：根据 PC 序列的符号为每一年标记不同的颜色（正为红色，负为蓝色），用于后续的折线图。

5.图形设置

fig = plt.figure(figsize=(15, 15), dpi=500)
proj = ccrs.PlateCarree(central_longitude=115)
leftlon, rightlon, lowerlat, upperlat = (70, 140, 15, 55)
lon_formatter = cticker.LongitudeFormatter()
lat_formatter = cticker.LatitudeFormatter()

plt.figure(): Matplotlib 的 figure 函数用于创建一个新的图形窗口。
figsize: 设置图形窗口的大小。本例中为 (15, 15)，单位是英寸。这表示图形的宽高比都是 15 英寸。
dpi: 表示每英寸的点数（dots per inch）。在这里设置为 500，这意味着图形将非常高分辨率（即适合于精细的地理数据可视化），这对于打印或详细查看是非常合适的。
ccrs.PlateCarree(): 这是 Cartopy 中的一个投影类，用于创建等距纬度-经度平面图（Plate Carree 投影）。
central_longitude=115: 设置地图的中心经度为 115°。这意味着地图将在水平方向上居中于 115°经线，通常用来使特定地区（例如中国）在渲染时位于地图的中央部分。
leftlon, rightlon: 分别是地理区域的左边界和右边界经度，范围是从 70°E 到 140°E。
lowerlat, upperlat: 分别是地理区域的下边界和上边界纬度，范围从 15°N 到 55°N。(这些坐标限定了地图绘制的地理范围。对于亚洲尤其是东亚地区的研究，使用这样的范围可以聚焦关注目标地区。)
LongitudeFormatter 和 LatitudeFormatter: Cartopy 提供的格式器，用于将经纬度显示为标准的格式（一般是度，如 120°E, 30°N）。(这些格式器会在后续设置轴标签时使用，以确保地图上显示的经纬度标记格式可读且标准化。)

6.EOF 模态场绘制(（以下代码段的解释可以重复类似理解，因为结构相同，只是涉及的 EOF 不同，末尾会附上完整代码）)

fig_ax1 = fig.add_axes([0.1, 0.8, 0.5, 0.3],projection = proj)
fig_ax1.set_extent([leftlon, rightlon, lowerlat, upperlat], crs=ccrs.PlateCarree())
fig_ax1.add_feature(cfeature.COASTLINE.with_scale('50m'))
fig_ax1.add_feature(cfeature.LAKES, alpha=0.5)
fig_ax1.set_xticks(np.arange(leftlon,rightlon+10,10), crs=ccrs.PlateCarree())
fig_ax1.set_yticks(np.arange(lowerlat,upperlat+10,10), crs=ccrs.PlateCarree())
fig_ax1.xaxis.set_major_formatter(lon_formatter)
fig_ax1.yaxis.set_major_formatter(lat_formatter)
china = shpreader.Reader(r'D:\shp\bou2_4l.dbf').geometries()
fig_ax1.add_geometries(china, ccrs.PlateCarree(),facecolor='none', edgecolor='black',zorder = 1)
fig_ax1.set_title('(a) EOF1',loc='left',fontsize =15)
fig_ax1.set_title( '%.2f%%' % (var[0]*100),loc='right',fontsize =15)
c1=fig_ax1.contourf(hgt_Eu.lon,hgt_Eu.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.add_axes(): 在之前创建的 fig 图形对象上添加一个新的子轴（子图）。括号内的列表 [0.1, 0.8, 0.5, 0.3] 定义了子图在图形中的相对位置和大小，具体表示为 [左, 下, 宽度, 高度]。这些值都是相对于整个图形尺寸的比例。

projection = proj: 指定使用之前定义的地理投影 proj（PlateCarree 投影）来绘制地图。

set_extent(): 设置子图的地理显示范围。leftlon, rightlon, lowerlat, upperlat 是之前定义的地理范围参数，这里用来指定地图展示的经纬度边界。

crs=ccrs.PlateCarree(): 明确指定当前范围的坐标参考系统是 PlateCarree，确保范围参数被正确解释。

add_feature(): 添加地理特征到地图上。

cfeature.COASTLINE.with_scale(‘50m’): 添加全球海岸线特征，‘50m’ 表示使用中等分辨率的数据。

cfeature.LAKES: 添加湖泊特征，alpha=0.5 将这些湖泊的透明度设为 50%，使得湖泊在背景上显示得半透明。

set_xticks() 和 set_yticks(): 分别设置 X 轴（经度）和 Y 轴（纬度）的刻度位置。

np.arange(): 使用 NumPy 的 arange 函数生成从 leftlon 到 rightlon 以及从 lowerlat 到 upperlat 的刻度，步长为 10 度。

crs=ccrs.PlateCarree(): 指定这些刻度的位置基于 PlateCarree 投影的坐标参考系统。

为 X 轴和 Y 轴设置格式化器，分别使用之前定义的 lon_formatter 和 lat_formatter，这将使得坐标轴上的经纬度标签以合理的格式显示（如度数和方向）。

shpreader.Reader(): 读取指定路径的 shapefile 文件。

.geometries(): 获取形状文件中的几何对象（例如边界）。

add_geometries(): 将读取到的中国边界几何形状添加到地图上，facecolor=‘none’ 意味着无填充色，edgecolor=‘black’ 设置边界颜色为黑色，zorder=1 控制绘制顺序。

设置两个标题：一个在左侧，标识图例编号“(a) EOF1”；一个在右侧，显示 EOF 第一个模态的方差贡献，以百分比形式展示。

fontsize=15 设置标题的字体大小。

contourf(): 绘制等高线填充图，将数据 (EOF 模态) 以不同颜色的区域显示。

levels=np.arange(-0.9, 1.0, 0.1): 定义填充的等高线级别，从 -0.9 到 1.0，每隔0.1 等级。

zorder=0: 设置绘制的顺序层级，与其他地理元素叠放。

extend=‘both’: 指定当数据超出定义的颜色级范围时应如何呈现（此处为两端都超出时都显示）。

transform=ccrs.PlateCarree(): 指定数据的坐标参考系统为PlateCarree。

cmap=plt.cm.RdBu_r: 使用颜色映射 RdBu_r（反转的红蓝色映射）来表示数值高低，通常用来强调正负变化。

7.添加小型地理子图(（以下代码段的解释可以重复类似理解，因为结构相同，只是涉及的 EOF 不同，末尾会附上完整代码）)

fig_ax11 = fig.add_axes([0.525, 0.08, 0.072, 0.15],projection = proj)
fig_ax11.set_extent([105, 125, 0, 25], crs=ccrs.PlateCarree())
fig_ax11.add_feature(cfeature.COASTLINE.with_scale('50m'))
china = shpreader.Reader(r'...\bou2_4l.dbf').geometries()
fig_ax11.add_geometries(china, ccrs.PlateCarree(),facecolor='none', edgecolor='black',zorder = 1)

fig.add_axes(): 在主图上添加一个新的子轴（这里是作为插图或缩略图）。方括号中的四个值 [0.525, 0.08, 0.072, 0.15] 定义了插图在主图中的位置和大小，具体为 [左, 下, 宽度, 高度]，这些值是相对于整个图形尺寸的比例。
其中：0.525 和 0.08：分别表示插图的左边缘和下边缘离图形左边缘和底部的相对距离。0.072 和 0.15：分别表示插图的宽度和高度，其单位也是相对比例。

projection = proj: 使用之前定义的 proj 投影（PlateCarree 投影）来绘制该插图。

set_extent(): 设置插图显示的地理范围，这里展示了一个更加精细聚焦的区域：经度范围从 105°E 到 125°E。纬度范围从 0° 到 25°N。

crs=ccrs.PlateCarree(): 指定该范围是根据 PlateCarree 坐标参考系统来理解的。

add_feature(): 在插图中添加地理特征，这里是海岸线。

cfeature.COASTLINE.with_scale(‘50m’): 指定使用中等分辨率的海岸线数据，以确保在缩略图中能较清晰地显示海岸线。

shpreader.Reader(): 读取指定位置的 shapefile 文件，这里用于获取中国边界的几何信息。

.geometries(): 返回形状文件中的几何图形对象。

add_geometries(): 将读取到的中国边界添加到地图中，facecolor=‘none’ 表示无填充颜色，而 edgecolor=‘black’ 将边界线绘制为黑色。

zorder=1: 指定绘制顺序，确保边界线在地图的其他元素之上绘制，以清晰显示。

8. EOF 分析的第一个主成分（PC1）的时间序列（后面的主成分分析结构一样，只需要修改数值即可，末尾附上完整代码）

fig_ax4 = fig.add_axes([0.65, 0.808, 0.47, 0.285])
fig_ax4.set_title('(b) PC1',loc='left',fontsize = 15)
fig_ax4.set_ylim(min(pc[:,0]),max(pc[:,0]))
fig_ax4.axhline(0,linestyle="--")
fig_ax4.bar(np.arange(1948,2009,1),pc[:,0],color=color1)

fig.add_axes(): 在主图上添加新的子轴（子图），这次是用于显示一个条形图。列表参数 [0.65, 0.808, 0.47, 0.285] 描述了子图在主图形中的位置和大小，分别表示为 [左, 下, 宽度, 高度]：0.65 和 0.808：定义了子图左下角在主图中的位置（以相对比例表示）。0.47 和 0.285：定义了子图的宽度和高度（同样为相对比例）。

set_title(): 为这个子图设置标题。

‘(b) PC1’ 表示该子图是关于第一个主成分（PC1）的图像，同时用字母标识为 (b)。

loc=‘left’ 确保标题在子图左侧。

fontsize=15 设置标题的字体大小。

set_ylim(): 设置 y 轴的显示范围。

min(pc[:,0]) 和 max(pc[:,0]) 分别计算 PC1 时间序列中最小和最大值，以确保条形图覆盖所有数据范围。

axhline(): 在 y=0 的位置添加一条水平线，通常用于强调零值基线。

linestyle=“–” 将该参考线设置为虚线样式，以便于区别数据条（颜色）和基线。

bar(): 绘制条形图，x 轴对应年份，y 轴对应 PC1 的值。

np.arange(1948, 2009, 1) 生成年份的范围，从 1948 年到 2008 年。

pc[:,0] 提供每一年的 PC1 值，以展示它随时间的变化。

color=color1：使用之前定义的颜色数组 color1，其中颜色基于 PC1 的符号（红色表示正值，蓝色表示负值），以视觉地突出正负值之间的变化。

9.显示图形

plt.show()

完整代码一：

import numpy as np
import xarray as xr
from eofs.standard import Eof
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LongitudeFormatter,LatitudeFormatter
import cartopy.io.shapereader as shpreader
import cartopy.mpl.ticker as cticker
import warnings
warnings.filterwarnings("ignore",category=DeprecationWarning)

# 更改字体
plt.rcParams['font.sans-serif'] = ['Times New Roman']
plt.rcParams['axes.unicode_minus'] = False

f = xr.open_dataset(r'...\hgt.mon.mean.nc')  # 读取数据集
hgt_Eu = f.hgt.loc[f.time.dt.month.isin([6]), 850].loc['1948':'2008']  # 中国地区一月高度场，压强为850pa
coslat = np.cos(np.deg2rad(hgt_Eu.lat.values))  # 纬度权重
solver = Eof(hgt_Eu.values, weights=coslat[..., np.newaxis])  # EOF分解器

##获取前三个模态，获取对应的PC序列和解释方差
eof = solver.eofsAsCorrelation(neofs=3)  # 模态场
pc = solver.pcs(npcs=3,pcscaling=1)  # 时间序列
var = solver.varianceFraction()  # 解释方差

#####以上是EOF构建及处理步骤
color1, color2, color3 = [], [], []

for i in range(1948, 2009):
    color1.append('red' if pc[i - 1948, 0] >= 0 else 'blue')
    color2.append('red' if pc[i - 1948, 1] >= 0 else 'blue')
    color3.append('red' if pc[i - 1948, 2] >= 0 else 'blue')


fig = plt.figure(figsize=(15,15),dpi=500)
proj = ccrs.PlateCarree(central_longitude=115)
leftlon, rightlon, lowerlat, upperlat = (70,140,15,55)
lon_formatter = cticker.LongitudeFormatter()
lat_formatter = cticker.LatitudeFormatter()

fig_ax1 = fig.add_axes([0.1, 0.8, 0.5, 0.3],projection = proj)
fig_ax1.set_extent([leftlon, rightlon, lowerlat, upperlat], crs=ccrs.PlateCarree())
fig_ax1.add_feature(cfeature.COASTLINE.with_scale('50m'))
fig_ax1.add_feature(cfeature.LAKES, alpha=0.5)
fig_ax1.set_xticks(np.arange(leftlon,rightlon+10,10), crs=ccrs.PlateCarree())
fig_ax1.set_yticks(np.arange(lowerlat,upperlat+10,10), crs=ccrs.PlateCarree())
fig_ax1.xaxis.set_major_formatter(lon_formatter)
fig_ax1.yaxis.set_major_formatter(lat_formatter)
china = shpreader.Reader(r'...\bou2_4l.dbf').geometries()
fig_ax1.add_geometries(china, ccrs.PlateCarree(),facecolor='none', edgecolor='black',zorder = 1)
fig_ax1.set_title('(a) EOF1',loc='left',fontsize =15)
fig_ax1.set_title( '%.2f%%' % (var[0]*100),loc='right',fontsize =15)
c1=fig_ax1.contourf(hgt_Eu.lon,hgt_Eu.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.45, 0.5, 0.3],projection = proj)
fig_ax2.set_extent([leftlon, rightlon, lowerlat, upperlat], crs=ccrs.PlateCarree())
fig_ax2.add_feature(cfeature.COASTLINE.with_scale('50m'))
fig_ax2.add_feature(cfeature.LAKES, alpha=0.5)
fig_ax2.set_xticks(np.arange(leftlon,rightlon+10,10), crs=ccrs.PlateCarree())
fig_ax2.set_yticks(np.arange(lowerlat,upperlat+10,10), crs=ccrs.PlateCarree())
fig_ax2.xaxis.set_major_formatter(lon_formatter)
fig_ax2.yaxis.set_major_formatter(lat_formatter)
china = shpreader.Reader(r'...\bou2_4l.dbf').geometries()
fig_ax2.add_geometries(china, ccrs.PlateCarree(),facecolor='none', edgecolor='black',zorder = 1)
fig_ax2.set_title('(c) EOF2',loc='left',fontsize =15)
fig_ax2.set_title( '%.2f%%' % (var[1]*100),loc='right',fontsize =15)
c2=fig_ax2.contourf(hgt_Eu.lon,hgt_Eu.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.1, 0.5, 0.3],projection = proj)
fig_ax3.set_extent([leftlon, rightlon, lowerlat, upperlat], crs=ccrs.PlateCarree())
fig_ax3.add_feature(cfeature.COASTLINE.with_scale('50m'))
fig_ax3.add_feature(cfeature.LAKES, alpha=0.5)
fig_ax3.set_xticks(np.arange(leftlon,rightlon+10,10), crs=ccrs.PlateCarree())
fig_ax3.set_yticks(np.arange(lowerlat,upperlat+10,10), crs=ccrs.PlateCarree())
fig_ax3.xaxis.set_major_formatter(lon_formatter)
fig_ax3.yaxis.set_major_formatter(lat_formatter)
china = shpreader.Reader(r'...\bou2_4l.dbf').geometries()
fig_ax3.add_geometries(china, ccrs.PlateCarree(),facecolor='none', edgecolor='black',zorder = 1)
fig_ax3.set_title('(d) EOF3',loc='left',fontsize =15)
fig_ax3.set_title( '%.2f%%' % (var[2]*100),loc='right',fontsize =15)
c2=fig_ax3.contourf(hgt_Eu.lon,hgt_Eu.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_ax11 = fig.add_axes([0.525, 0.08, 0.072, 0.15],projection = proj)
fig_ax11.set_extent([105, 125, 0, 25], crs=ccrs.PlateCarree())
fig_ax11.add_feature(cfeature.COASTLINE.with_scale('50m'))
china = shpreader.Reader(r'...\bou2_4l.dbf').geometries()
fig_ax11.add_geometries(china, ccrs.PlateCarree(),facecolor='none', edgecolor='black',zorder = 1)


fig_ax22 = fig.add_axes([0.525, 0.43, 0.072, 0.15],projection = proj)
fig_ax22.set_extent([105, 125, 0, 25], crs=ccrs.PlateCarree())
fig_ax22.add_feature(cfeature.COASTLINE.with_scale('50m'))
china = shpreader.Reader(r'...\bou2_4l.dbf').geometries()
fig_ax22.add_geometries(china, ccrs.PlateCarree(),facecolor='none', edgecolor='black',zorder = 1)

fig_ax33 = fig.add_axes([0.525, 0.78, 0.072, 0.15],projection = proj)
fig_ax33.set_extent([105, 125, 0, 25], crs=ccrs.PlateCarree())
fig_ax33.add_feature(cfeature.COASTLINE.with_scale('50m'))
china = shpreader.Reader(r'...\bou2_4l.dbf').geometries()
fig_ax33.add_geometries(china, ccrs.PlateCarree(),facecolor='none', edgecolor='black',zorder = 1)

cbposition=fig.add_axes([0.102, 0.04, 0.5, 0.015])
fig.colorbar(c1,cax=cbposition,orientation='horizontal',format='%.1f',)

fig_ax4 = fig.add_axes([0.65, 0.808, 0.47, 0.285])
fig_ax4.set_title('(b) PC1',loc='left',fontsize = 15)
fig_ax4.set_ylim(min(pc[:,0]),max(pc[:,0]))
fig_ax4.axhline(0,linestyle="--")
fig_ax4.bar(np.arange(1948,2009,1),pc[:,0],color=color1)

fig_ax5 = fig.add_axes([0.65, 0.458, 0.47, 0.285])
fig_ax5.set_title('(d) PC2',loc='left',fontsize = 15)
fig_ax5.set_ylim(min(pc[:,1]),max(pc[:,1]))
fig_ax5.axhline(0,linestyle="--")
fig_ax5.bar(np.arange(1948,2009,1),pc[:,1],color=color2)

fig_ax6 = fig.add_axes([0.65, 0.108, 0.47, 0.285])
fig_ax6.set_title('(f) PC3',loc='left',fontsize = 15)
fig_ax6.set_ylim(min(pc[:,2]),max(pc[:,2]))
fig_ax6.axhline(0,linestyle="--")
fig_ax6.bar(np.arange(1948,2009,1),pc[:,2],color=color3)
plt.show()