气候统计界的瑞士军刀——sacpy

用户11172986

发布于 2024-08-14 21:35:21

1080

发布于 2024-08-14 21:35:21

文章被收录于专栏：气python风雨

气候统计界的瑞士军刀——sacpy

前言

最近收到读者的一封信，表达了对Python在气象数据分析领域应用的兴趣，尤其是气象归因分析方面。我认为读者可能是指气象数据的统计分析和模型验证，其中“归因分析”通常指的是确定气候变化原因的科学方法。考虑到这一点，我深入探索了相关的资源，并在众多开源社区中发现了一个引人注目的项目——sacpy。

sacpy 库简介

sacpy是一个由南京信息工程大学的孟子路同学精心打造的Python库，旨在为气候科学研究提供便利。这个库不仅包含了丰富的统计工具，还特别针对气候数据分析进行了优化，使其成为处理复杂气象数据的理想选择。

项目主页：sacpy GitHub Repository

sacpy库的核心优势在于其封装了一系列广泛应用于气候统计学的功能，包括但不限于：

• EOF（经验正交函数）分析：用于识别空间模式及其随时间的变化。
• SVD（奇异值分解）：帮助揭示不同数据集之间的相关性或共同变化模式。
• T-检验：评估两组数据之间是否存在显著差异。
• 功率谱和交叉谱分析：用于频域中的信号分析，识别周期性和频率成分。

这些工具对于理解气候系统的内在动态、检测长期趋势以及评估极端事件的影响至关重要。sacpy的出现，无疑为气候科学家提供了一个强大且用户友好的工具箱，简化了复杂的统计计算过程，使研究者能够更加专注于科学问题本身，而非繁琐的数据处理和算法实现。

如果您对使用Python进行气象数据分析感兴趣，sacpy绝对值得您深入了解和尝试。无论是学术研究还是专业工作，它都有潜力成为您的得力助手。

可点击原文链接去和鲸社区在线运行

安装

!pip install sacpy -i https://pypi.mirrors.ustc.edu.cn/simple/

!pip install xmca -i https://pypi.mirrors.ustc.edu.cn/simple/

EOF

import sacpy as scp
import numpy as np
import matplotlib.pyplot as plt
# get data
sst = scp.load_sst()["sst"].loc[:, -20:30, 150:275]
ssta = scp.get_anom(sst)
# EOF
eof = scp.EOF(np.array(ssta))
eof.solve()
# get spartial pattern and pc
pc = eof.get_pc(npt=2)
pt = eof.get_pt(npt=2)
# plot
import cartopy.crs as ccrs
import sacpy.Map
lon , lat = np.array(ssta.lon) , np.array(ssta.lat)
fig = plt.figure(figsize=[15,10])
ax = fig.add_subplot(221,projection=ccrs.PlateCarree(central_longitude=180))
m1 = ax.scontourf(lon,lat,pt[0,:,:],cmap='RdBu_r',levels=np.linspace(-0.75,0.75,15),extend="both")
ax.scontour(m1,colors="black")
ax.init_map(ysmall=2.5)
# plt.colorbar(m1)
ax2 = fig.add_subplot(222)
ax2.plot(sst.time,pc[0])
ax3 = fig.add_subplot(223,projection=ccrs.PlateCarree(central_longitude=180))
m2 = ax3.scontourf(lon,lat,pt[1,:,:],cmap='RdBu_r',levels=np.linspace(-0.75,0.75,15),extend="both")
ax3.scontour(m2,colors="black")
ax3.init_map(ysmall=2.5)
ax4 = fig.add_subplot(224)
ax4.plot(sst.time,pc[1])
cb_ax = fig.add_axes([0.1,0.06,0.4,0.02])
fig.colorbar(m1,cax=cb_ax,orientation="horizontal")
plt.show()

SVD

import sacpy as scp
import xarray as xr
import matplotlib.pyplot as plt
import numpy as np
from xmca import array
import sacpy.Map
import cartopy.crs as ccrs

# load data
sst = scp.load_sst()['sst'].loc["1991":"2021", -20:30, 150:275]
ssta = scp.get_anom(sst)
u = scp.load_10mwind()['u']
v = scp.load_10mwind()['v']

uua = scp.get_anom(u)
vua = scp.get_anom(v)
uv = np.concatenate([np.array(uua)[...,np.newaxis],np.array(vua)[...,np.newaxis]],axis=-1)
# calculation
svd = scp.SVD(ssta,uv,complex=False)
svd.solve()
ptl, ptr = svd.get_pt(3)
pcl,pcr = svd.get_pc(3)
upt ,vpt = ptr[...,0] , ptr[...,1]
sst_pt = ptl
## 绘图
import cartopy.crs as ccrs
import sacpy.Map
lon , lat = np.array(ssta.lon) , np.array(ssta.lat)
fig = plt.figure(figsize=[12,8])
ax = fig.add_subplot(221,projection=ccrs.PlateCarree(central_longitude=180))
m1 = ax.scontourf(lon,lat,sst_pt[0],cmap='RdBu_r',levels=np.linspace(-0.1,0.1,13),extend="both")
ax.scontour(m1,colors="black")
ax.squiver(lon,lat,upt[0],vpt[0])
ax.init_map(smally=2.5)
ax2 = fig.add_subplot(222)
ax2.plot(sst.time,pcl[0],label="left")
ax2.plot(sst.time,pcr[0],label="right")
ax2.legend()
ax3 = fig.add_subplot(223,projection=ccrs.PlateCarree(central_longitude=180))
m2 = ax3.scontourf(lon,lat,sst_pt[1],cmap='RdBu_r',levels=np.linspace(-0.1,0.1,13),extend="both")
ax3.squiver(lon,lat,upt[1],vpt[1])
ax3.scontour(m2,colors="black")
ax3.init_map(smally=2.5)
ax4 = fig.add_subplot(224)
ax4.plot(sst.time,pcl[1],label="left")
ax4.plot(sst.time,pcr[1],label="right")
ax4.legend()
cb_ax = fig.add_axes([0.1,0.06,0.4,0.02])
fig.colorbar(m1,cax=cb_ax,orientation="horizontal")
plt.show()

T检验

import sacpy as scp
import numpy as np
import matplotlib.pyplot as plt

sst = scp.load_sst()["sst"]
ssta = scp.get_anom(sst, method=0)
# get Dec Jan Feb SSTA
ssta_djf = scp.XrTools.spec_moth_yrmean(ssta,[12,1,2])
Nino34 = ssta_djf.loc[:, -5:5, 190:240].mean(axis=(1, 2))
# select year of Super El Nino
select = Nino34 >= 1
ssta_sl = ssta_djf[select]
mean, pv = scp.one_mean_test(ssta_sl)
# plot
import sacpy.Map
import cartopy.crs as ccrs
fig = plt.figure(figsize=[12, 8])
ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=180))
lon ,lat = np.array(ssta.lon) , np.array(ssta.lat)
m = ax.scontourf(lon,lat,mean)
n = ax.sig_plot(lon,lat,pv,color="k",marker="..")
ax.init_map(stepx=50, ysmall=2.5)
plt.colorbar(m)
plt.show()

功率谱与交叉谱分析

import sacpy
from sacpy import spectral

Spectral 类 这个类主要用于单个时间序列的功率谱分析。它封装了一系列方法用于计算和分析一个时间序列的频谱特性，如计算功率谱(get_spectra)，拟合红噪声(get_red_noise)，以及获取基于红噪声理论的统计显著性曲线(get_rs_sig)。它适用于想要了解一个信号内部结构的情况，如识别信号中的周期性成分、趋势或是噪声特征。

CrossSpectral 类 继承自 Spectral 类，CrossSpectral 主要用于两个时间序列之间的交叉频谱分析。它提供了额外的方法，如 get_cross_spectra 来获取两个信号的互相关性在频率域的表现，包括协方差、相位和相干度。此外，它也可以计算相干度的显著性(get_coh_sig)。交叉谱分析常用于研究两个信号之间的关系，比如信号间的相位差、同步性和线性依赖程度，特别适用于系统识别和信号耦合分析。

这两个类都提供了一种系统的方式来处理时间序列的频谱分析需求，使得用户能够深入地探索信号的内在特性及其与其他信号的相互作用。

一句话概括，功率谱关注于单一信号的频率成分。交叉谱关注于两个信号之间的频率相关性和相位差异。

import numpy as np
import matplotlib.pyplot as plt
from sacpy.spectral import Spectral

# 创建一个模拟的时间序列
np.random.seed(0)
time = np.linspace(0, 100, 1000)
x = np.sin(2 * np.pi * 1 * time) + np.sin(2 * np.pi * 10 * time) + np.random.normal(size=time.shape)

# 创建 Spectral 分析器实例
spectral_analyzer = Spectral(x=x, M_length=256)

# 计算功率谱
freq, Pxx = spectral_analyzer.get_spectra()

# 绘制功率谱
plt.figure()
plt.semilogy(freq, Pxx)
plt.title('Power Spectral Density')
plt.xlabel('Frequency [Hz]')
plt.ylabel('PSD [V**2/Hz]')
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from sacpy.spectral import CrossSpectral

# 创建两个相关的时间序列
np.random.seed(0)
time = np.linspace(0, 100, 1000)
x = np.sin(2 * np.pi * 1 * time) + np.random.normal(size=time.shape)
y = 0.5 * np.sin(2 * np.pi * 1 * time) + np.random.normal(size=time.shape)

# 创建 CrossSpectral 分析器实例
cross_spectral_analyzer = CrossSpectral(x=x, y=y, M_length=256)
# 获取频率向量
freq, _, _ = cross_spectral_analyzer.get_spectra()

# 计算交叉谱
covariance, phase, coherence = cross_spectral_analyzer.get_cross_spectra()

# 绘制相干度
plt.figure()
plt.plot(freq, coherence)
plt.title('Coherence')
plt.xlabel('Frequency [Hz]')
plt.ylabel('Coherence')
plt.grid(True)
plt.show()

# 绘制相位
plt.figure()
plt.plot(freq, phase)
plt.title('Phase')
plt.xlabel('Frequency [Hz]')
plt.ylabel('Phase [radians]')
plt.grid(True)
plt.show()