基於Python進行小波分析

姜颢睿發表於2024-07-29

在氣象學和環境科學的研究中,理解和預測氣象資料的週期性變化至關重要。小波分析作為一種高效的數學工具,近年來在氣象資料的週期性分析中得到了廣泛應用。本文將詳細介紹如何透過Python進行小波分析,以探究氣象資料中的週期性變化。

1 資料來源及下載方式

西北農林科技大學的彭守璋研究員在國家青藏高原科學資料中心公開發布了氣溫、降水、乾燥度等氣象資料,本文所使用的資料為基於其中的 中國1km解析度逐月降水量資料集 轉換得到的某地區歷年降水量資料,可透過FTP進行下載。

Filezilla 桌面版軟體為例,輸入官網中給定的主機號、使用者名稱、密碼、埠等資訊後即可將資料從遠端站點下載至本地。

image

2 程式碼編寫

2.1 匯入相關模組

import openpyxl
import pywt
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib import rcParams

2.2 設定地圖字型及字號

config = {
    "mathtext.fontset": "stix",
    "font.family": "serif",
    "font.serif": ["Times New Roman"],
    "font.size": 24,
    "axes.unicode_minus": False # 處理負號,即-號
}
rcParams.update(config)

2.3 從Excel表讀取資料

# 透過pandas模組讀取Excel表,並提取表格中的年份及降水量資料
file_path = "file_path"	# 此處為Excel檔案路徑
data = pd.read_excel(file_path)
Year_list = list(data["Year"])
Pre_list = list(data["Pre"])

2.4 進行小波變換,繪製小波等值線圖

scales = np.arange(1, 31)	# 設定小波分析的時間尺度
coef, freqs = pywt.cwt(Pre_list, scales, "morl")	# 對降水量資料進行Morlet小波變換

fig, ax = plt.subplots(figsize=(15, 10))

# 繪製小波係數影像
im = plt.imshow(abs(coef), extent=[Year_list[0], Year_list[-1], 30, 1],
                interpolation="bilinear", cmap="gray", aspect="auto",
                vmax=abs(coef).max(), vmin=-abs(coef).max())
plt.colorbar(im, ax=ax)  # 新增色帶圖例

# 繪製小波係數等值線
contour = plt.contour(Year_list, scales, coef.real, colors="black", linewidths=1)

ax.invert_yaxis()   # 反轉Y軸,使時間尺度從大到小排列
ax.set_xticks(np.arange(1901, 2024, 20))
ax.set_yticks(np.arange(5, 31, 5))
ax.set_xlabel("Year")
ax.set_ylabel("Scale")

2.5 計算小波方差,繪製小波方差圖

variance = np.var(coef, axis=1)	# 計算小波方差
variance1 = variance / 10000	# 數值較大,將數值縮小10000倍後顯示在結果圖中


fig, ax = plt.subplots(figsize=(15, 10))

ax.plot(scales, variance1, "k-")	# 繪製方差曲線
ax.set_xlabel("Scale (a)")
# \times為LaTeX語法中的叉乘號
# 在Python中,反斜槓"\"為跳脫字元,此處需同時輸入兩個反斜槓
ax.set_ylabel("Variance $\mathrm{(\\times 10^4)}$")
ax.set_xticks(np.arange(0, 31, 5))
ax.set_yticks(np.arange(0, 9, 2))
fig.savefig("小波方差圖.jpg", dpi=600)

3 等值線及方差圖示例

image

image

相關文章