cleanup
This commit is contained in:
parent
a832042b54
commit
f567196561
GPG Key ID:
35AA0AA3DC402D5C
@ -1,175 +0,0 @@
|
||||
# 探空气球只有重力波
|
||||
import matplotlib.pyplot as plt
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
import numpy as np
|
||||
from modules.balloon.extract_wave import calculate_wave, is_terrain_wave
|
||||
|
||||
|
||||
def plot_once(data: pd.DataFrame, path: str, lat: float, g: float):
|
||||
wave = calculate_wave(
|
||||
data[(data["alt"] >= 15) & (data["alt"] <= 25)], lat, g)
|
||||
if len(wave) == 0:
|
||||
return []
|
||||
|
||||
c = is_terrain_wave(data, lat, g)
|
||||
a, b, omega_upper, w_f, λ_z, λ_h, c_x, c_y, c_z, Ek, E_p, MFu1, MFv1, params_u, params_v, params_T = wave[
|
||||
:16]
|
||||
ucmp, vcmp, temp, height, poly_ucmp, poly_vcmp, _, poly_temp, residual_ucmp, residual_vcmp, residual_temp, u_fit, v_fit, T_fit, uh, _ = wave[
|
||||
16:]
|
||||
|
||||
plt.figure(figsize=(16, 10))
|
||||
|
||||
# 二阶多项式拟合曲线
|
||||
|
||||
# u 风扰动量
|
||||
plt.subplot(2, 6, 1)
|
||||
plt.plot(ucmp, height, label="", linestyle="-")
|
||||
plt.plot(poly_ucmp(height), height, label="")
|
||||
plt.ylabel("Height(km)")
|
||||
plt.xlabel("Zonal wind (m/s)")
|
||||
plt.grid(True)
|
||||
|
||||
# v 风扰动量
|
||||
plt.subplot(2, 6, 2)
|
||||
plt.plot(vcmp, height, label="", linestyle="-")
|
||||
plt.plot(poly_vcmp(height), height, label="")
|
||||
plt.ylabel("Height(km)")
|
||||
plt.xlabel("Meridional wind (m/s)")
|
||||
plt.grid(True)
|
||||
|
||||
plt.title("观测的二阶多项式拟合", fontsize=16)
|
||||
|
||||
# 温度扰动量
|
||||
plt.subplot(2, 6, 3)
|
||||
plt.plot(temp, height, label="", linestyle="-")
|
||||
plt.plot(poly_temp(height), height, label="")
|
||||
plt.ylabel("Height(km)")
|
||||
plt.xlabel("Temperature(K)")
|
||||
plt.grid(True)
|
||||
|
||||
# 绘制正弦拟合图
|
||||
|
||||
# u 风扰动量
|
||||
plt.subplot(2, 6, 4)
|
||||
plt.plot(residual_ucmp, height, label="", marker="o", linestyle="None")
|
||||
plt.plot(u_fit, height, label="")
|
||||
plt.ylabel("Height(km)")
|
||||
plt.xlabel("Zonal wind (m/s)")
|
||||
plt.grid(True)
|
||||
|
||||
# v 风扰动量
|
||||
plt.subplot(2, 6, 5)
|
||||
plt.plot(residual_vcmp, height, label="", marker="o", linestyle="None")
|
||||
plt.plot(v_fit, height, label="")
|
||||
plt.ylabel("Height(km)")
|
||||
plt.xlabel("Meridional wind (m/s)")
|
||||
plt.grid(True)
|
||||
|
||||
plt.title("扰动分量的正弦波拟合", fontsize=16)
|
||||
|
||||
# 温度扰动量
|
||||
plt.subplot(2, 6, 6)
|
||||
plt.plot(residual_temp, height, label="", marker="o", linestyle="None")
|
||||
plt.plot(T_fit, height, label="")
|
||||
plt.ylabel("Height(km)")
|
||||
plt.xlabel("Temperature(K)")
|
||||
plt.grid(True)
|
||||
|
||||
# 标记3个特定高度
|
||||
specified_heights = [15, 15.05, 15.1]
|
||||
markers = ["*", "D", "^"] # 星号，菱形，三角形
|
||||
|
||||
# a. U-V矢量曲线
|
||||
plt.subplot(2, 2, 3)
|
||||
plt.plot(u_fit, v_fit) # 绘制U-V曲线
|
||||
for h, marker in zip(specified_heights, markers):
|
||||
index = np.abs(height - h).argmin() # 找到离指定高度最近的索引
|
||||
plt.scatter(u_fit[index], v_fit[index],
|
||||
marker=marker, s=100, label=f"{h} km")
|
||||
|
||||
# 设置坐标轴范围和比例
|
||||
plt.xlim(-8, 8)
|
||||
plt.ylim(-4, 4)
|
||||
# plt.gca().set_aspect("equal") # 确保横纵坐标刻度长短一致
|
||||
plt.axvline(0, color="gray", linestyle="--")
|
||||
plt.axhline(0, color="gray", linestyle="--")
|
||||
plt.xlabel("Zonal Wind (m/s)")
|
||||
plt.ylabel("Meridional Wind (m/s)")
|
||||
plt.legend() # 显示图例
|
||||
plt.grid(True) # 显示网格线
|
||||
plt.text(0.05, 1, "(a)", transform=plt.gca().transAxes,
|
||||
fontsize=14, verticalalignment="top")
|
||||
plt.title("径向风-纬向风矢量图")
|
||||
|
||||
# b. 水平风-温度的矢端曲线
|
||||
plt.subplot(2, 2, 4)
|
||||
plt.plot(uh, T_fit) # 绘制水平风-温度曲线
|
||||
for h, marker in zip(specified_heights, markers):
|
||||
index = np.abs(height - h).argmin() # 找到离指定高度最近的索引
|
||||
plt.scatter(uh[index], T_fit[index],
|
||||
marker=marker, s=100, label=f"{h} km")
|
||||
|
||||
# 设置坐标轴范围和比例
|
||||
plt.xlim(-4, 4)
|
||||
plt.ylim(-2, 2)
|
||||
# plt.gca().set_aspect("equal") # 确保横纵坐标刻度长短一致
|
||||
plt.axvline(0, color="gray", linestyle="--")
|
||||
plt.axhline(0, color="gray", linestyle="--")
|
||||
plt.xlabel("Horizontal wind (m/s)")
|
||||
plt.ylabel("Temp (K)")
|
||||
plt.legend()
|
||||
plt.grid(True)
|
||||
plt.text(0.05, 1, "(b)", transform=plt.gca().transAxes,
|
||||
fontsize=14, verticalalignment="top")
|
||||
plt.title("温度-水平风矢量图")
|
||||
|
||||
# 设置中文显示和负号正常显示
|
||||
plt.rcParams["font.sans-serif"] = ["SimHei"] # 显示中文
|
||||
plt.rcParams["axes.unicode_minus"] = False # 正常显示负号
|
||||
|
||||
plt.tight_layout()
|
||||
# plt.savefig(path, transparent=True)
|
||||
plt.savefig(path)
|
||||
plt.close()
|
||||
|
||||
"""
|
||||
地形重力波 c 是否是地形重力波
|
||||
垂直传播方向 a 1上
|
||||
水平传播方向 b degree
|
||||
本征（固有）频率 round(omega_upper, 6) rad/s
|
||||
本征周期 round(2 * np.pi / omega_upper / 3600, 2) h
|
||||
本征频率与固有频率的比值（即长短轴之比） round(w_f, 2)
|
||||
垂直波长 round(λ_z, 2) km
|
||||
水平波长 round(λ_h, 2) km
|
||||
纬向本征相速度 round(c_x, 2) m/s
|
||||
经向本征相速度 round(c_y, 2) m/s
|
||||
垂直本征相速度 round(c_z, 3) m/s
|
||||
波动能 round(Ek, 4) J/kg
|
||||
势能 round(E_p, 4) J/kg
|
||||
纬向动量通量 round(MFu1, 4) mPa
|
||||
经向动量通量 round(MFv1, 4) mPa
|
||||
纬向风扰动振幅 round(params_u, 2) m/s
|
||||
经向风扰动振幅 round(params_v, 2) m/s
|
||||
温度扰动振幅 round(params_T, 2) K
|
||||
"""
|
||||
return [
|
||||
c,
|
||||
a,
|
||||
b,
|
||||
round(omega_upper, 6),
|
||||
round(2 * np.pi / omega_upper / 3600, 2),
|
||||
round(w_f, 2),
|
||||
round(λ_z, 2),
|
||||
round(λ_h, 2),
|
||||
round(c_x, 2),
|
||||
round(c_y, 2),
|
||||
round(c_z, 3),
|
||||
round(Ek, 4),
|
||||
round(E_p, 4),
|
||||
round(MFu1, 4),
|
||||
round(MFv1, 4),
|
||||
round(params_u, 2),
|
||||
round(params_v, 2),
|
||||
round(params_T, 2),
|
||||
]
|
||||
@ -1,388 +0,0 @@
|
||||
# 探空气球只有重力波
|
||||
|
||||
import matplotlib.pyplot as plt
|
||||
import matplotlib.gridspec as gridspec
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
import seaborn as sns
|
||||
from windrose import WindroseAxes
|
||||
|
||||
|
||||
def plot_year(data: pd.DataFrame, path: str, lat: float, g: float):
|
||||
if data.size == 0:
|
||||
return False
|
||||
|
||||
data.loc[:, "date"] = data["file_name"].str[:15]
|
||||
filtered_df = data[[
|
||||
"date"] + [col for col in data.columns if col != "file_name" and col != "date"]]
|
||||
filtered_df.reset_index(drop=True, inplace=True)
|
||||
|
||||
filtered_df = filtered_df.drop_duplicates(
|
||||
subset="date", keep="last") # 使用 drop_duplicates 函数，保留每组重复中的最后一个记录
|
||||
filtered_df.reset_index(drop=True, inplace=True)
|
||||
|
||||
filtered_df = filtered_df[filtered_df["w_f"] < 10] # 筛选 w_f 值小于 10 的数据
|
||||
|
||||
# todo:1-删除不合理的数据
|
||||
# 1.先剔除明显异常的值
|
||||
for column in filtered_df.columns[1:]: # 不考虑第一列日期列
|
||||
filtered_df = filtered_df[(
|
||||
filtered_df[column] >= -9999) & (filtered_df[column] <= 9999)] # 460
|
||||
|
||||
# 2.再用四分位数法，适合所有数据集
|
||||
def remove_outliers_iqr(df):
|
||||
for column in df.columns[9:]: # 从第10列开始
|
||||
Q1 = df[column].quantile(0.25)
|
||||
Q3 = df[column].quantile(0.75)
|
||||
IQR = Q3 - Q1
|
||||
lower_bound = Q1 - 5 * IQR
|
||||
upper_bound = Q3 + 5 * IQR
|
||||
|
||||
# 过滤掉异常值
|
||||
df = df[(df[column] >= lower_bound) & (df[column] <= upper_bound)]
|
||||
return df
|
||||
|
||||
filtered_df = remove_outliers_iqr(filtered_df) # 408
|
||||
|
||||
# 画图
|
||||
|
||||
fig = plt.figure(figsize=(36, 20))
|
||||
gs = gridspec.GridSpec(4, 12)
|
||||
|
||||
ax1 = fig.add_subplot(gs[0, 0:3])
|
||||
ax2 = fig.add_subplot(gs[0, 3:6])
|
||||
ax3 = fig.add_subplot(gs[0, 6:9])
|
||||
ax4 = fig.add_subplot(gs[0, 9:12])
|
||||
|
||||
ax5_1 = fig.add_subplot(gs[1, 0:2])
|
||||
ax5_2 = fig.add_subplot(gs[1, 2:4])
|
||||
ax5_3 = fig.add_subplot(gs[1, 4:6])
|
||||
ax6_1 = fig.add_subplot(gs[1, 6:8])
|
||||
ax6_2 = fig.add_subplot(gs[1, 8:10])
|
||||
ax6_3 = fig.add_subplot(gs[1, 10:12])
|
||||
|
||||
ax7_1 = fig.add_subplot(gs[2, 0:2])
|
||||
ax7_2 = fig.add_subplot(gs[2, 2:4])
|
||||
ax8 = fig.add_subplot(gs[2, 4:8])
|
||||
ax9 = fig.add_subplot(gs[2, 8:12])
|
||||
|
||||
ax10 = []
|
||||
for i in range(0, 10, 3):
|
||||
ax10.append(fig.add_subplot(gs[3, i: i + 3], projection="windrose"))
|
||||
|
||||
sns.set_theme(style="whitegrid", font="SimHei") # 设置绘图样式为白色背景和网格线
|
||||
|
||||
# 1、w/f比值
|
||||
# 计算bins的边界
|
||||
min_val = 1
|
||||
max_val = 10
|
||||
bin_width = 0.5
|
||||
# 创建bins的边界
|
||||
bins = np.arange(min_val, max_val + bin_width, bin_width)
|
||||
|
||||
sns.histplot(filtered_df["w_f"], bins=bins, kde=False, edgecolor="black",
|
||||
stat="percent", ax=ax1) # 加上stat='percent'表示算的是频率
|
||||
# 设置x轴范围
|
||||
ax1.set_xlim(min_val, max_val)
|
||||
# 添加标题和标签
|
||||
ax1.set_title("w/f值统计结果", fontsize=24)
|
||||
ax1.set_xlabel("w/f")
|
||||
ax1.set_ylabel("Occurrence(%)")
|
||||
|
||||
# 2、周期
|
||||
# 计算bins的边界
|
||||
min_val = 1
|
||||
max_val = 10
|
||||
bin_width = 0.5
|
||||
# 创建bins的边界
|
||||
bins = np.arange(min_val, max_val + bin_width, bin_width)
|
||||
|
||||
min_val = np.floor(filtered_df["zhou_qi"].min()) # 向下取整
|
||||
max_val = np.ceil(filtered_df["zhou_qi"].max()) # 向上取整
|
||||
bin_width = 1
|
||||
bins = np.arange(min_val, max_val + bin_width, bin_width)
|
||||
|
||||
sns.histplot(filtered_df["zhou_qi"], bins=bins, kde=False,
|
||||
edgecolor="black", stat="percent", ax=ax2) # 加上stat='percent'表示算的是频率
|
||||
# 设置x轴范围
|
||||
ax2.set_xlim(min_val, max_val)
|
||||
# 添加标题和标签
|
||||
ax2.set_title("周期统计结果", fontsize=24)
|
||||
ax2.set_xlabel("h")
|
||||
ax2.set_ylabel("Occurrence(%)")
|
||||
|
||||
# 3、垂直波长
|
||||
# 创建bins的边界
|
||||
min_val = np.floor(filtered_df["ver_wave_len"].min()) # 向下取整
|
||||
max_val = np.ceil(filtered_df["ver_wave_len"].max()) # 向上取整
|
||||
bin_width = 0.5
|
||||
bins = np.arange(min_val, max_val + bin_width, bin_width)
|
||||
|
||||
sns.histplot(filtered_df["ver_wave_len"], bins=bins, kde=False,
|
||||
edgecolor="black", stat="percent", ax=ax3) # 加上stat='percent'表示算的是频率
|
||||
# 设置x轴范围
|
||||
ax3.set_xlim(min_val, max_val)
|
||||
# 添加标题和标签
|
||||
ax3.set_title("垂直波长分布", fontsize=24)
|
||||
ax3.set_xlabel("Vertical wavelength(km)")
|
||||
ax3.set_ylabel("Occurrence(%)")
|
||||
|
||||
# 4、水平波长
|
||||
# 创建bins的边界
|
||||
min_val = np.floor(filtered_df["hori_wave_len"].min()) # 向下取整
|
||||
max_val = np.ceil(filtered_df["hori_wave_len"].max()) # 向上取整
|
||||
bin_width = 100
|
||||
bins = np.arange(min_val, max_val + bin_width, bin_width)
|
||||
|
||||
sns.histplot(filtered_df["hori_wave_len"], bins=bins, kde=False,
|
||||
edgecolor="black", stat="percent", ax=ax4) # 加上stat='percent'表示算的是频率
|
||||
# 设置x轴范围
|
||||
ax4.set_xlim(min_val, max_val)
|
||||
# 添加标题和标签
|
||||
ax4.set_title("水平波长分布", fontsize=24)
|
||||
ax4.set_xlabel("Horizontal wavelength(km)")
|
||||
ax4.set_ylabel("Occurrence(%)")
|
||||
|
||||
# 5、本征相速度
|
||||
|
||||
# 纬向本征相速度
|
||||
# 计算bins的边界
|
||||
min_val = np.floor(filtered_df["c_x"].min()) # 向下取整
|
||||
max_val = np.ceil(filtered_df["c_x"].max()) # 向上取整
|
||||
bin_width = 10
|
||||
bins = np.arange(min_val, max_val + bin_width, bin_width)
|
||||
# 创建纬向直方图
|
||||
sns.histplot(filtered_df["c_x"], bins=bins, kde=False,
|
||||
edgecolor="black", stat="percent", ax=ax5_1)
|
||||
ax5_1.set_xlim(min_val, max_val)
|
||||
ax5_1.set_title("纬向本征相速度", fontsize=24)
|
||||
ax5_1.set_xlabel("Zonal phase speed (m/s)")
|
||||
ax5_1.set_ylabel("Occurrence (%)")
|
||||
|
||||
# 经向本征相速度
|
||||
# 计算bins的边界
|
||||
min_val = np.floor(filtered_df["c_y"].min()) # 向下取整
|
||||
max_val = np.ceil(filtered_df["c_y"].max()) # 向上取整
|
||||
bin_width = 10
|
||||
bins = np.arange(min_val, max_val + bin_width, bin_width)
|
||||
# 创建经向直方图
|
||||
sns.histplot(filtered_df["c_y"], bins=bins, kde=False,
|
||||
edgecolor="black", stat="percent", ax=ax5_2)
|
||||
ax5_2.set_xlim(min_val, max_val)
|
||||
ax5_2.set_title("经向本征相速度", fontsize=24)
|
||||
ax5_2.set_xlabel("Meridional phase speed (m/s)")
|
||||
ax5_2.set_ylabel("Occurrence (%)")
|
||||
|
||||
# 垂直本征相速度
|
||||
# 计算bins的边界
|
||||
min_val = filtered_df["c_z"].min() # -1.148
|
||||
max_val = filtered_df["c_z"].max() # 0.624
|
||||
bin_width = 0.1
|
||||
bins = np.arange(min_val, max_val + bin_width, bin_width)
|
||||
# 创建垂直直方图
|
||||
sns.histplot(filtered_df["c_z"], bins=bins, kde=False,
|
||||
edgecolor="black", stat="percent", ax=ax5_3)
|
||||
ax5_3.set_xlim(min_val, max_val)
|
||||
ax5_3.set_title("垂直本征相速度", fontsize=24)
|
||||
ax5_3.set_xlabel("Vertical phase speed (m/s)")
|
||||
ax5_3.set_ylabel("Occurrence (%)")
|
||||
|
||||
# 6、扰动振幅
|
||||
|
||||
# 纬向扰动振幅
|
||||
# 计算bins的边界
|
||||
min_val = np.floor(filtered_df["u1"].min()) # 向下取整
|
||||
max_val = np.ceil(filtered_df["u1"].max()) # 向上取整
|
||||
bin_width = 0.5
|
||||
bins = np.arange(min_val, max_val + bin_width, bin_width)
|
||||
# 创建纬向直方图
|
||||
sns.histplot(filtered_df["u1"], bins=bins, kde=False,
|
||||
edgecolor="black", stat="percent", ax=ax6_1)
|
||||
ax6_1.set_xlim(min_val, max_val)
|
||||
ax6_1.set_title(" ", fontsize=24)
|
||||
ax6_1.set_xlabel("Zonal wind amplitude (m/s)")
|
||||
ax6_1.set_ylabel("Occurrence (%)")
|
||||
|
||||
# 经向扰动振幅
|
||||
# 计算bins的边界
|
||||
min_val = np.floor(filtered_df["v1"].min()) # 向下取整
|
||||
max_val = np.ceil(filtered_df["v1"].max()) # 向上取整
|
||||
bin_width = 0.5
|
||||
bins = np.arange(min_val, max_val + bin_width, bin_width)
|
||||
# 创建经向直方图
|
||||
sns.histplot(filtered_df["v1"], bins=bins, kde=False,
|
||||
edgecolor="black", stat="percent", ax=ax6_2)
|
||||
ax6_2.set_xlim(min_val, max_val)
|
||||
ax6_2.set_title("扰动振幅统计结果", fontsize=24)
|
||||
ax6_2.set_xlabel("Meridional wind amplitude (m/s)")
|
||||
ax6_2.set_ylabel("Occurrence (%)")
|
||||
|
||||
# 垂直扰动振幅
|
||||
# 计算bins的边界
|
||||
min_val = np.floor(filtered_df["T1"].min()) # 向下取整
|
||||
max_val = np.ceil(filtered_df["T1"].max()) # 向上取整
|
||||
bin_width = 0.5
|
||||
bins = np.arange(min_val, max_val + bin_width, bin_width)
|
||||
# 创建垂直直方图
|
||||
sns.histplot(filtered_df["T1"], bins=bins, kde=False,
|
||||
edgecolor="black", stat="percent", ax=ax6_3)
|
||||
ax6_3.set_xlim(min_val, max_val)
|
||||
ax6_3.set_title(" ", fontsize=24)
|
||||
ax6_3.set_xlabel("Temperature amplitude (K)")
|
||||
ax6_3.set_ylabel("Occurrence (%)")
|
||||
|
||||
# 7、动量通量
|
||||
# 挑选出向上传的重力波
|
||||
filtered_df1 = filtered_df[filtered_df["a"] == 1]
|
||||
|
||||
# 绘制第一个子图
|
||||
# 计算bins的边界
|
||||
min_val = np.floor(filtered_df1["MFu"].min()) # 向下取整
|
||||
max_val = np.ceil(filtered_df1["MFu"].max()) # 向上取整
|
||||
bin_width = 0.1
|
||||
bins = np.arange(min_val, max_val + bin_width, bin_width)
|
||||
|
||||
sns.histplot(filtered_df1["MFu"], bins=bins, kde=False,
|
||||
edgecolor="black", stat="percent", ax=ax7_1)
|
||||
ax7_1.set_xlim(min_val, max_val) # 设置x轴范围
|
||||
ax7_1.set_title("纬向动量通量统计结果", fontsize=24) # 添加标题
|
||||
ax7_1.set_xlabel("Zonal momentum flux(mPa)") # x轴标签
|
||||
ax7_1.set_ylabel("Occurrence(%)") # y轴标签
|
||||
|
||||
# 绘制第二个子图
|
||||
# 计算bins的边界
|
||||
min_val = np.floor(filtered_df1["MFv"].min()) # 向下取整
|
||||
max_val = np.ceil(filtered_df1["MFv"].max()) # 向上取整
|
||||
bin_width = 0.1
|
||||
bins = np.arange(min_val, max_val + bin_width, bin_width)
|
||||
|
||||
sns.histplot(filtered_df1["MFv"], bins=bins, kde=False,
|
||||
edgecolor="black", stat="percent", ax=ax7_2)
|
||||
ax7_2.set_xlim(min_val, max_val) # 设置x轴范围
|
||||
ax7_2.set_title("经向动量通量统计结果", fontsize=24) # 添加标题
|
||||
ax7_2.set_xlabel("Meridional momentum flux(mPa)") # x轴标签
|
||||
ax7_2.set_ylabel("Occurrence(%)") # y轴标签
|
||||
|
||||
# 10、水平传播方向
|
||||
|
||||
filtered_df["date1"] = filtered_df["date"].str.split(
|
||||
"T").str[0] # 再加一列，只保留日期部分
|
||||
filtered_df["date1"] = pd.to_datetime(
|
||||
filtered_df["date1"], format="%Y%m%d") # 确保 'date1' 列是日期格式
|
||||
|
||||
# 添加季节列
|
||||
def get_season(date):
|
||||
month = date.month
|
||||
if month in [12, 1, 2]:
|
||||
return "Winter"
|
||||
elif month in [3, 4, 5]:
|
||||
return "Spring"
|
||||
elif month in [6, 7, 8]:
|
||||
return "Summer"
|
||||
else:
|
||||
return "Fall"
|
||||
|
||||
filtered_df["season"] = filtered_df["date1"].apply(get_season) # 添加季节列
|
||||
seasons = ["Winter", "Spring", "Summer", "Fall"]
|
||||
|
||||
for ax, season in zip(ax10, seasons):
|
||||
season_data = filtered_df[filtered_df["season"] == season]
|
||||
windrose = WindroseAxes.from_ax(ax)
|
||||
ax.set_title(season, fontsize=18)
|
||||
windrose.bar(season_data["b"], np.ones_like(
|
||||
season_data["b"]), normed=False) # normed=True表示占每个季节的占比
|
||||
|
||||
# # 添加总标题
|
||||
# fig.suptitle("水平传播方向在各个季节的分布变化", fontsize=16, fontweight="bold")
|
||||
|
||||
# 8、垂直传播方向
|
||||
# 设置 日期'date1' 列为索引
|
||||
filtered_df.set_index("date1", inplace=True)
|
||||
|
||||
# 按月份分组并计算百分比
|
||||
monthly_stats_df = (
|
||||
filtered_df.groupby([filtered_df.index.month, filtered_df.index.year])
|
||||
.apply(lambda x: pd.Series({"Upload (%)": (x["a"] == 1).mean() * 100, "Downward (%)": (x["a"] == -1).mean() * 100}))
|
||||
.reset_index(level=1, drop=True)
|
||||
)
|
||||
|
||||
# 按月份汇总这些年的数据
|
||||
monthly_avg_stats_df = monthly_stats_df.groupby(level=0).mean()
|
||||
# 确保索引是 numpy 数组
|
||||
dates = monthly_avg_stats_df.index.to_numpy()
|
||||
|
||||
# 绘制折线图
|
||||
ax8.plot(dates, monthly_avg_stats_df["Upload (%)"].to_numpy(
|
||||
), marker="o", label="Up (%)")
|
||||
ax8.plot(dates, monthly_avg_stats_df["Downward (%)"].to_numpy(
|
||||
), marker="o", label="Down (%)")
|
||||
|
||||
# 设置月份标签
|
||||
ax8.set_xticks(
|
||||
ticks=np.arange(1, 13), labels=["Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"], rotation=0
|
||||
) # 不倾斜，rotation=0表示倾斜45°
|
||||
# 添加图例、标题和坐标轴标签
|
||||
ax8.legend()
|
||||
ax8.set_title("每月上传/下传重力波占比", fontsize=24)
|
||||
ax8.set_xlabel("Month")
|
||||
ax8.set_ylabel("Percentage (%)")
|
||||
|
||||
# 9、动能和势能
|
||||
filtered_df.reset_index(inplace=True) # 取消索引并恢复为默认整数索引
|
||||
# 提取年月
|
||||
filtered_df["year_month"] = filtered_df["date1"].dt.to_period("M")
|
||||
# 计算每个月的动能和势能的平均值
|
||||
monthly_avg = filtered_df.groupby("year_month")[["Ek", "E_p"]].mean()
|
||||
|
||||
# 创建完整的月份范围（因为有的月份没数据）
|
||||
full_range = pd.period_range(
|
||||
start=monthly_avg.index.min(), end=monthly_avg.index.max(), freq="M")
|
||||
|
||||
# 创建一个新的 DataFrame 使用完整的年份月份范围
|
||||
full_monthly_avg = pd.DataFrame(index=full_range)
|
||||
# 将原始数据合并到新的 DataFrame 中
|
||||
full_monthly_avg = full_monthly_avg.join(monthly_avg)
|
||||
|
||||
# 确保 'Ek' 和 'E_p' 列为数值型
|
||||
full_monthly_avg["Ek"] = pd.to_numeric(
|
||||
full_monthly_avg["Ek"], errors="coerce")
|
||||
full_monthly_avg["E_p"] = pd.to_numeric(
|
||||
full_monthly_avg["E_p"], errors="coerce")
|
||||
|
||||
# 只显示每年6月、12月，简化图形
|
||||
|
||||
# 绘制 Ek、E_p
|
||||
ax9.plot(full_monthly_avg.index.values.astype(
|
||||
str), full_monthly_avg["Ek"].values, marker="o", linestyle="-", color="r", label="动能月平均值")
|
||||
ax9.plot(full_monthly_avg.index.values.astype(
|
||||
str), full_monthly_avg["E_p"].values, marker="o", linestyle="-", color="b", label="势能月平均值")
|
||||
|
||||
# 添加标题和标签
|
||||
ax9.set_title("动能和势能分布情况", fontsize=24)
|
||||
ax9.set_xlabel("Month", fontsize=14)
|
||||
ax9.set_ylabel("Wave energy (J/kg)", fontsize=14)
|
||||
|
||||
# 设定横轴标签
|
||||
months = full_monthly_avg.index.values.astype(str)
|
||||
# 获取所有年份的6月和12月的索引
|
||||
june_indices = [i for i, date in enumerate(months) if date.endswith("-06")]
|
||||
december_indices = [i for i, date in enumerate(
|
||||
months) if date.endswith("-12")]
|
||||
selected_indices = june_indices + december_indices
|
||||
|
||||
# 只显示选定的标签
|
||||
ax9.set_xticks(ticks=selected_indices, labels=[
|
||||
months[i] for i in selected_indices], rotation=45)
|
||||
|
||||
# 添加网格和图例
|
||||
# plt.grid()
|
||||
ax9.legend() # 显示图例
|
||||
|
||||
plt.rcParams["font.sans-serif"] = ["SimHei"] # 显示中文
|
||||
plt.rcParams["axes.unicode_minus"] = False # 正常显示负号
|
||||
|
||||
plt.tight_layout()
|
||||
plt.savefig(path)
|
||||
plt.close()
|
||||
|
||||
return True
|
||||
File diff suppressed because it is too large
Load Diff
@ -1,828 +0,0 @@
|
||||
# 推断出这是重力波
|
||||
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
import matplotlib.dates as mdates
|
||||
|
||||
from CONSTANT import DATA_BASEPATH
|
||||
from modules.saber.utils import *
|
||||
# from matplotlib.colors import LinearSegmentedColormap
|
||||
# 设置字体为支持中文的字体
|
||||
plt.rcParams['font.family'] = 'SimHei' # 设置为黑体（需要你的环境中有该字体）
|
||||
plt.rcParams['axes.unicode_minus'] = False # 解决负号'-'显示为方块的问题
|
||||
|
||||
|
||||
# ----------------------------------------------------------------------------------------------------------------------------
|
||||
# 5---main程序环节 -
|
||||
# ----------------------------------------------------------------------------------------------------------------------------
|
||||
# 按日处理资料
|
||||
|
||||
|
||||
def day_process_main(file_path, day_read, latitude_min, latitude_max, altitude_min, altitude_max, lamda_low, lamda_high, lvboin):
|
||||
|
||||
dataset, tplatitude, tplongitude, tpaltitude, ktemp, time, date, date_time = data_nc_load(
|
||||
file_path)
|
||||
# 2018年的第94天 # 4月4日的日期，以年份和年内的第几天表示
|
||||
df = day_data_read(date, day_read, tplatitude)
|
||||
cycles = data_cycle_identify(df, latitude_min, latitude_max)
|
||||
ktemp_cycles, altitude_cycles = data_cycle_generate(
|
||||
cycles, ktemp, tpaltitude, altitude_min, altitude_max)
|
||||
ktemp_wn0 = ktemp_cycles - \
|
||||
np.mean(ktemp_cycles, axis=0) # 观测值-平均温度
|
||||
ktemp_fit_wn1, ktemp_wn1 = fit_wave(ktemp_wn0, 1)
|
||||
ktemp_fit_wn2, ktemp_wn2 = fit_wave(ktemp_wn1, 2)
|
||||
ktemp_fit_wn3, ktemp_wn3 = fit_wave(ktemp_wn2, 3)
|
||||
ktemp_fit_wn4, ktemp_wn4 = fit_wave(ktemp_wn3, 4)
|
||||
ktemp_fit_wn5, ktemp_wn5 = fit_wave(ktemp_wn4, 5)
|
||||
ktemp_fft, ktemp_fft_lvbo, ktemp_ifft = fft_ifft_wave(
|
||||
ktemp_wn5, lamda_low, lamda_high, altitude_min, altitude_max, lvboin)
|
||||
ktemp_Nz, ktemp_Ptz = power_indices(
|
||||
ktemp_cycles, ktemp_wn5, ktemp_ifft, altitude_min, altitude_max)
|
||||
|
||||
return ktemp_cycles, altitude_cycles, ktemp_wn0, ktemp_fit_wn1, ktemp_wn1, ktemp_fit_wn2, ktemp_wn2, ktemp_fit_wn3, ktemp_wn3, ktemp_fit_wn4, ktemp_wn4, ktemp_fit_wn5, ktemp_wn5, ktemp_fft, ktemp_fft_lvbo, ktemp_ifft, ktemp_Nz, ktemp_Ptz
|
||||
# ----------------------------------------------------------------------------------------------------------------------------
|
||||
# 按文件中单日处理资料
|
||||
|
||||
|
||||
def day_process_maing(dataset, tplatitude, tplongitude, tpaltitude, ktemp, time, date, date_time, day_read,
|
||||
latitude_min, latitude_max, altitude_min, altitude_max, lamda_low, lamda_high, lvboin):
|
||||
df = day_data_read(date, day_read, tplatitude)
|
||||
cycles = data_cycle_identify(df, latitude_min, latitude_max)
|
||||
|
||||
if not cycles: # 如果周期列表为空，返回18个None值
|
||||
return (None,) * 18
|
||||
|
||||
ktemp_cycles, altitude_cycles = data_cycle_generate(
|
||||
cycles, ktemp, tpaltitude, altitude_min, altitude_max)
|
||||
|
||||
if ktemp_cycles is None or altitude_cycles is None: # 再次检查周期数据是否为空
|
||||
return (None,) * 18
|
||||
|
||||
ktemp_wn0 = ktemp_cycles - \
|
||||
np.mean(ktemp_cycles, axis=0) # 按照纬向计算平均温度 wn0_temp
|
||||
ktemp_fit_wn1, ktemp_wn1 = fit_wave(ktemp_wn0, 1)
|
||||
ktemp_fit_wn2, ktemp_wn2 = fit_wave(ktemp_wn1, 2)
|
||||
ktemp_fit_wn3, ktemp_wn3 = fit_wave(ktemp_wn2, 3)
|
||||
ktemp_fit_wn4, ktemp_wn4 = fit_wave(ktemp_wn3, 4)
|
||||
ktemp_fit_wn5, ktemp_wn5 = fit_wave(ktemp_wn4, 5)
|
||||
ktemp_fft, ktemp_fft_lvbo, ktemp_ifft = fft_ifft_wave(ktemp_wn5, lamda_low, lamda_high, altitude_min, altitude_max,
|
||||
lvboin)
|
||||
ktemp_Nz, ktemp_Ptz = power_indices(
|
||||
ktemp_cycles, ktemp_wn5, ktemp_ifft, altitude_min, altitude_max)
|
||||
|
||||
return ktemp_cycles, altitude_cycles, ktemp_wn0, ktemp_fit_wn1, ktemp_wn1, ktemp_fit_wn2, ktemp_wn2, ktemp_fit_wn3, ktemp_wn3, ktemp_fit_wn4, ktemp_wn4, ktemp_fit_wn5, ktemp_wn5, ktemp_fft, ktemp_fft_lvbo, ktemp_ifft, ktemp_Nz, ktemp_Ptz
|
||||
|
||||
# ----------------------------------------------------------------------------------------------------------------------------
|
||||
# 按月处理资料
|
||||
|
||||
|
||||
def mon_process_main(file_path, latitude_min, latitude_max, altitude_min, altitude_max, lamda_low, lamda_high, lvboin):
|
||||
# 打开文件并读取数据
|
||||
dataset, tplatitude, tplongitude, tpaltitude, ktemp, time, date, date_time = data_nc_load(
|
||||
file_path)
|
||||
|
||||
ktemp_cycles_mon = []
|
||||
altitude_cycles_mon = []
|
||||
ktemp_wn0_mon = []
|
||||
ktemp_fit_wn1_mon = []
|
||||
ktemp_wn1_mon = []
|
||||
ktemp_fit_wn2_mon = []
|
||||
ktemp_wn2_mon = []
|
||||
ktemp_fit_wn3_mon = []
|
||||
ktemp_wn3_mon = []
|
||||
ktemp_fit_wn4_mon = []
|
||||
ktemp_wn4_mon = []
|
||||
ktemp_fit_wn5_mon = []
|
||||
ktemp_wn5_mon = []
|
||||
ktemp_fft_mon = []
|
||||
ktemp_fft_lvbo_mon = []
|
||||
ktemp_ifft_mon = []
|
||||
ktemp_Nz_mon = []
|
||||
ktemp_Ptz_mon = []
|
||||
|
||||
# 遍历每一天，处理数据
|
||||
for day_read in date_time:
|
||||
print(f"读取日期 {day_read}")
|
||||
# 处理单日数据
|
||||
results = day_process_maing(
|
||||
dataset, tplatitude, tplongitude, tpaltitude, ktemp, time, date, date_time,
|
||||
day_read, latitude_min, latitude_max, altitude_min, altitude_max, lamda_low, lamda_high, lvboin
|
||||
)
|
||||
|
||||
# 检查结果是否包含有效数据
|
||||
if results is not None and len(results) == 18:
|
||||
ktemp_cycles0, altitude_cycles0, ktemp_wn00, ktemp_fit_wn10, ktemp_wn10, ktemp_fit_wn20, ktemp_wn20, ktemp_fit_wn30, ktemp_wn30, ktemp_fit_wn40, ktemp_wn40, ktemp_fit_wn50, ktemp_wn50, ktemp_fft0, ktemp_fft_lvbo0, ktemp_ifft0, ktemp_Nz0, ktemp_Ptz0 = results
|
||||
|
||||
# 将有效数据添加到月度列表中
|
||||
ktemp_cycles_mon.append(ktemp_cycles0)
|
||||
altitude_cycles_mon.append(altitude_cycles0)
|
||||
ktemp_wn0_mon.append(ktemp_wn00)
|
||||
ktemp_fit_wn1_mon.append(ktemp_fit_wn10)
|
||||
ktemp_wn1_mon.append(ktemp_wn10)
|
||||
ktemp_fit_wn2_mon.append(ktemp_fit_wn20)
|
||||
ktemp_wn2_mon.append(ktemp_wn20)
|
||||
ktemp_fit_wn3_mon.append(ktemp_fit_wn30)
|
||||
ktemp_wn3_mon.append(ktemp_wn30)
|
||||
ktemp_fit_wn4_mon.append(ktemp_fit_wn40)
|
||||
ktemp_wn4_mon.append(ktemp_wn40)
|
||||
ktemp_fit_wn5_mon.append(ktemp_fit_wn50)
|
||||
ktemp_wn5_mon.append(ktemp_wn50)
|
||||
ktemp_fft_mon.append(ktemp_fft0)
|
||||
ktemp_fft_lvbo_mon.append(ktemp_fft_lvbo0)
|
||||
ktemp_ifft_mon.append(ktemp_ifft0)
|
||||
ktemp_Nz_mon.append(ktemp_Nz0)
|
||||
ktemp_Ptz_mon.append(ktemp_Ptz0)
|
||||
|
||||
# 返回整个月的数据
|
||||
return (date_time, ktemp_cycles_mon, altitude_cycles_mon, ktemp_wn0_mon, ktemp_fit_wn1_mon, ktemp_wn1_mon,
|
||||
ktemp_fit_wn2_mon, ktemp_wn2_mon, ktemp_fit_wn3_mon, ktemp_wn3_mon, ktemp_fit_wn4_mon, ktemp_wn4_mon,
|
||||
ktemp_fit_wn5_mon, ktemp_wn5_mon, ktemp_fft_mon, ktemp_fft_lvbo_mon, ktemp_ifft_mon, ktemp_Nz_mon,
|
||||
ktemp_Ptz_mon)
|
||||
# 滤波后NZ、PTZ重力波势能指标计算
|
||||
# ----------------------------------------------------------------------------------------------------------------------------
|
||||
# 按年处理资料
|
||||
|
||||
|
||||
def process_yearly_data(year, path, latitude_min, latitude_max, altitude_min, altitude_max, lamda_low, lamda_high, lvboin):
|
||||
# 创建空列表来存储每月的数据
|
||||
date_time_list = []
|
||||
ktemp_cycles_mon_list = []
|
||||
altitude_cycles_mon_list = []
|
||||
ktemp_wn0_mon_list = []
|
||||
ktemp_fit_wn1_mon_list = []
|
||||
ktemp_wn1_mon_list = []
|
||||
ktemp_fit_wn2_mon_list = []
|
||||
ktemp_wn2_mon_list = []
|
||||
ktemp_fit_wn3_mon_list = []
|
||||
ktemp_wn3_mon_list = []
|
||||
ktemp_fit_wn4_mon_list = []
|
||||
ktemp_wn4_mon_list = []
|
||||
ktemp_fit_wn5_mon_list = []
|
||||
ktemp_wn5_mon_list = []
|
||||
ktemp_fft_mon_list = []
|
||||
ktemp_fft_lvbo_mon_list = []
|
||||
ktemp_ifft_mon_list = []
|
||||
ktemp_Nz_mon_list = []
|
||||
ktemp_Ptz_mon_list = []
|
||||
|
||||
# 循环处理每个月的数据
|
||||
for month in range(1, 13):
|
||||
# 获取当前月的文件路径
|
||||
file_path = filename_read(year, month, path)
|
||||
|
||||
try:
|
||||
# 调用 mon_process_main 函数处理文件并获取结果
|
||||
(date_time, ktemp_cycles_mon, altitude_cycles_mon, ktemp_wn0_mon, ktemp_fit_wn1_mon,
|
||||
ktemp_wn1_mon, ktemp_fit_wn2_mon, ktemp_wn2_mon, ktemp_fit_wn3_mon, ktemp_wn3_mon,
|
||||
ktemp_fit_wn4_mon, ktemp_wn4_mon, ktemp_fit_wn5_mon, ktemp_wn5_mon, ktemp_fft_mon,
|
||||
ktemp_fft_lvbo_mon, ktemp_ifft_mon, ktemp_Nz_mon, ktemp_Ptz_mon) = mon_process_main(
|
||||
file_path, latitude_min, latitude_max, altitude_min, altitude_max, lamda_low, lamda_high, lvboin
|
||||
)
|
||||
|
||||
# 将每月的结果添加到相应的列表中
|
||||
date_time_list.extend(date_time)
|
||||
ktemp_cycles_mon_list.extend(ktemp_cycles_mon)
|
||||
altitude_cycles_mon_list.extend(altitude_cycles_mon)
|
||||
ktemp_wn0_mon_list.extend(ktemp_wn0_mon)
|
||||
ktemp_fit_wn1_mon_list.extend(ktemp_fit_wn1_mon)
|
||||
ktemp_wn1_mon_list.extend(ktemp_wn1_mon)
|
||||
ktemp_fit_wn2_mon_list.extend(ktemp_fit_wn2_mon)
|
||||
ktemp_wn2_mon_list.extend(ktemp_wn2_mon)
|
||||
ktemp_fit_wn3_mon_list.extend(ktemp_fit_wn3_mon)
|
||||
ktemp_wn3_mon_list.extend(ktemp_wn3_mon)
|
||||
ktemp_fit_wn4_mon_list.extend(ktemp_fit_wn4_mon)
|
||||
ktemp_wn4_mon_list.extend(ktemp_wn4_mon)
|
||||
ktemp_fit_wn5_mon_list.extend(ktemp_fit_wn5_mon)
|
||||
ktemp_wn5_mon_list.extend(ktemp_wn5_mon)
|
||||
ktemp_fft_mon_list.extend(ktemp_fft_mon)
|
||||
ktemp_fft_lvbo_mon_list.extend(ktemp_fft_lvbo_mon)
|
||||
ktemp_ifft_mon_list.extend(ktemp_ifft_mon)
|
||||
ktemp_Nz_mon_list.extend(ktemp_Nz_mon)
|
||||
ktemp_Ptz_mon_list.extend(ktemp_Ptz_mon)
|
||||
|
||||
except Exception as e:
|
||||
print(f"处理文件 {file_path} 时出错（月份：{month}）：{e}")
|
||||
continue # 出现错误时跳过当前月份，继续处理下个月
|
||||
|
||||
# 返回所有月份的处理结果
|
||||
return {
|
||||
"date_time_list": date_time_list,
|
||||
"ktemp_cycles_mon_list": ktemp_cycles_mon_list,
|
||||
"altitude_cycles_mon_list": altitude_cycles_mon_list,
|
||||
"ktemp_wn0_mon_list": ktemp_wn0_mon_list,
|
||||
"ktemp_fit_wn1_mon_list": ktemp_fit_wn1_mon_list,
|
||||
"ktemp_wn1_mon_list": ktemp_wn1_mon_list,
|
||||
"ktemp_fit_wn2_mon_list": ktemp_fit_wn2_mon_list,
|
||||
"ktemp_wn2_mon_list": ktemp_wn2_mon_list,
|
||||
"ktemp_fit_wn3_mon_list": ktemp_fit_wn3_mon_list,
|
||||
"ktemp_wn3_mon_list": ktemp_wn3_mon_list,
|
||||
"ktemp_fit_wn4_mon_list": ktemp_fit_wn4_mon_list,
|
||||
"ktemp_wn4_mon_list": ktemp_wn4_mon_list,
|
||||
"ktemp_fit_wn5_mon_list": ktemp_fit_wn5_mon_list,
|
||||
"ktemp_wn5_mon_list": ktemp_wn5_mon_list,
|
||||
"ktemp_fft_mon_list": ktemp_fft_mon_list,
|
||||
"ktemp_fft_lvbo_mon_list": ktemp_fft_lvbo_mon_list,
|
||||
"ktemp_ifft_mon_list": ktemp_ifft_mon_list,
|
||||
"ktemp_Nz_mon_list": ktemp_Nz_mon_list,
|
||||
"ktemp_Ptz_mon_list": ktemp_Ptz_mon_list
|
||||
}
|
||||
|
||||
# -----------------------------------------------------------------------------------------------------------------------------
|
||||
# 6-主要统计分析图绘制 -
|
||||
# -----------------------------------------------------------------------------------------------------------------------------
|
||||
# 6-1 示例的逐层滤波效果图---不同波数 曲线图
|
||||
|
||||
|
||||
def day_fit_wave_plot(height_no, ktemp_wn0, ktemp_fit_wn1, ktemp_wn1, ktemp_fit_wn2, ktemp_wn2, ktemp_fit_wn3, ktemp_wn3, ktemp_fit_wn4, ktemp_wn4, ktemp_fit_wn5, ktemp_wn5):
|
||||
|
||||
N = len(ktemp_wn0[:, height_no])
|
||||
# 循环周期索引
|
||||
x = np.arange(N)
|
||||
|
||||
y1_1 = ktemp_wn0[:, height_no]
|
||||
y1_2 = ktemp_fit_wn1[:, height_no]
|
||||
|
||||
y2_1 = ktemp_wn1[:, height_no]
|
||||
y2_2 = ktemp_fit_wn2[:, height_no]
|
||||
|
||||
y3_1 = ktemp_wn2[:, height_no]
|
||||
y3_2 = ktemp_fit_wn3[:, height_no]
|
||||
|
||||
y4_1 = ktemp_wn3[:, height_no]
|
||||
y4_2 = ktemp_fit_wn4[:, height_no]
|
||||
|
||||
y5_1 = ktemp_wn4[:, height_no]
|
||||
y5_2 = ktemp_fit_wn5[:, height_no]
|
||||
|
||||
y6 = ktemp_wn5[:, height_no]
|
||||
|
||||
plt.figure(figsize=(16, 10)) # 调整图形大小
|
||||
# 原始信号的时间序列
|
||||
plt.subplot(2, 3, 1)
|
||||
plt.plot(x, y1_1, label='原始信号')
|
||||
plt.plot(x, y1_2, label='拟合信号', linestyle='--')
|
||||
plt.title('（a）波数k=1')
|
||||
plt.xlabel('Cycles', labelpad=10) # 增加标签间距
|
||||
plt.ylabel('温度 (K)', labelpad=10) # 增加标签间距
|
||||
plt.legend()
|
||||
plt.tight_layout() # 调整子图参数以适应图形区域
|
||||
|
||||
plt.subplot(2, 3, 2)
|
||||
plt.plot(x, y2_1, label='原始信号')
|
||||
plt.plot(x, y2_2, label='拟合信号', linestyle='--')
|
||||
plt.title('（b）波数k=2')
|
||||
plt.xlabel('Cycles', labelpad=10) # 增加标签间距
|
||||
plt.ylabel('温度 (K)', labelpad=10) # 增加标签间距
|
||||
plt.legend()
|
||||
plt.tight_layout() # 调整子图参数以适应图形区域
|
||||
|
||||
plt.subplot(2, 3, 3)
|
||||
plt.plot(x, y3_1, label='原始信号')
|
||||
plt.plot(x, y3_2, label='拟合信号', linestyle='--')
|
||||
plt.title('（c）波数k=3')
|
||||
plt.xlabel('Cycles', labelpad=10) # 增加标签间距
|
||||
plt.ylabel('温度 (K)', labelpad=10) # 增加标签间距
|
||||
plt.legend()
|
||||
plt.tight_layout() # 调整子图参数以适应图形区域
|
||||
|
||||
plt.subplot(2, 3, 4)
|
||||
plt.plot(x, y4_1, label='原始信号')
|
||||
plt.plot(x, y4_2, label='拟合信号', linestyle='--')
|
||||
plt.title('（d）波数k=4')
|
||||
plt.xlabel('Cycles', labelpad=10) # 增加标签间距
|
||||
plt.ylabel('温度 (K)', labelpad=10) # 增加标签间距
|
||||
plt.legend()
|
||||
plt.tight_layout() # 调整子图参数以适应图形区域
|
||||
|
||||
plt.subplot(2, 3, 5)
|
||||
plt.plot(x, y5_1, label='原始信号')
|
||||
plt.plot(x, y5_2, label='拟合信号', linestyle='--')
|
||||
plt.title('（e）波数k=5')
|
||||
plt.xlabel('Cycles', labelpad=10) # 增加标签间距
|
||||
plt.ylabel('温度 (K)', labelpad=10) # 增加标签间距
|
||||
plt.legend()
|
||||
plt.tight_layout() # 调整子图参数以适应图形区域
|
||||
|
||||
plt.subplot(2, 3, 6)
|
||||
plt.plot(x, y6, label='滤波信号')
|
||||
plt.title('（f）滤波1-5后信号')
|
||||
plt.xlabel('Cycles', labelpad=10) # 增加标签间距
|
||||
plt.ylabel('温度 (K)', labelpad=10) # 增加标签间距
|
||||
plt.tight_layout() # 调整子图参数以适应图形区域
|
||||
|
||||
# 调整子图之间的边距
|
||||
plt.subplots_adjust(top=0.8, bottom=0.2, left=0.1,
|
||||
right=0.8, hspace=0.3, wspace=0.2)
|
||||
|
||||
plt.show()
|
||||
|
||||
|
||||
def day_fit_wave_plotg(height_no, ktemp_wn0, ktemp_fit_wn1, ktemp_wn1, ktemp_fit_wn2,
|
||||
ktemp_wn2, ktemp_fit_wn3, ktemp_wn3, ktemp_fit_wn4,
|
||||
ktemp_wn4, ktemp_fit_wn5, ktemp_wn5):
|
||||
|
||||
N = len(ktemp_wn0[:, height_no])
|
||||
x = np.arange(N)
|
||||
|
||||
y1_1, y1_2 = ktemp_wn0[:, height_no], ktemp_fit_wn1[:, height_no]
|
||||
y2_1, y2_2 = ktemp_wn1[:, height_no], ktemp_fit_wn2[:, height_no]
|
||||
y3_1, y3_2 = ktemp_wn2[:, height_no], ktemp_fit_wn3[:, height_no]
|
||||
y4_1, y4_2 = ktemp_wn3[:, height_no], ktemp_fit_wn4[:, height_no]
|
||||
y5_1, y5_2 = ktemp_wn4[:, height_no], ktemp_fit_wn5[:, height_no]
|
||||
y6 = ktemp_wn5[:, height_no]
|
||||
|
||||
plt.figure(figsize=(16, 10))
|
||||
|
||||
y_limits = (min(min(y1_1), min(y2_1), min(y3_1), min(y4_1), min(y5_1), min(y6)),
|
||||
max(max(y1_1), max(y2_1), max(y3_1), max(y4_1), max(y5_1), max(y6)))
|
||||
|
||||
for i, (y1, y2) in enumerate([(y1_1, y1_2), (y2_1, y2_2), (y3_1, y3_2),
|
||||
(y4_1, y4_2), (y5_1, y5_2), (y6, None)]):
|
||||
plt.subplot(2, 3, i + 1)
|
||||
plt.plot(x, y1, label='原始信号')
|
||||
if y2 is not None:
|
||||
plt.plot(x, y2, label='拟合信号', linestyle='--')
|
||||
plt.title(f'（{"abcdef"[i]}）波数k={i + 1 if i < 5 else "滤波0-5"}')
|
||||
plt.xlabel('Cycles', labelpad=10)
|
||||
plt.ylabel('温度 (K)', labelpad=10)
|
||||
plt.legend()
|
||||
plt.xticks(x) # 设置横坐标为整数
|
||||
plt.ylim(y_limits) # 设置统一纵坐标范围
|
||||
plt.tight_layout()
|
||||
|
||||
plt.subplots_adjust(top=0.8, bottom=0.2, left=0.1,
|
||||
right=0.8, hspace=0.3, wspace=0.2)
|
||||
plt.show()
|
||||
|
||||
|
||||
# -----------------------------------------------------------------------------------------------------------------------------
|
||||
# 6-2 示例的高度滤波处理--不同循环周期 曲线图
|
||||
def day_fft_ifft_plot(cycle_no, ktemp_wn5, ktemp_fft, ktemp_ifft, altitude_min, altitude_max, lamda_low, lamda_high):
|
||||
|
||||
N = len(ktemp_wn5[cycle_no, :])
|
||||
# 采样时间间隔，其倒数等于采用频率，以1km为标准尺度等同于1s，假设波的速度为1km/s
|
||||
dt = (altitude_max-altitude_min)/(N-1)
|
||||
# 时间序列索引
|
||||
n = np.arange(N)
|
||||
f = n / (N * dt)
|
||||
t = np.round(np.linspace(altitude_min, altitude_max, N), 2)
|
||||
|
||||
# 原始扰动温度
|
||||
x = ktemp_wn5[cycle_no, :]
|
||||
# 傅里叶变换频谱分析
|
||||
y = ktemp_fft[cycle_no, :]
|
||||
# 滤波后的傅里叶变换频谱分析
|
||||
yy = ktemp_fft_lvbo[cycle_no, :]
|
||||
# 傅里叶逆变换后的扰动温度
|
||||
yyy = ktemp_ifft[cycle_no, :]
|
||||
|
||||
plt.figure(figsize=(15, 10)) # 调整图形大小
|
||||
# 原始信号的时间序列
|
||||
plt.subplot(2, 2, 1)
|
||||
plt.plot(t, x)
|
||||
plt.title('（a）原始信号')
|
||||
plt.xlabel('高度 (km)', labelpad=10) # 增加标签间距
|
||||
plt.ylabel('温度 (K)', labelpad=10) # 增加标签间距
|
||||
# 原始振幅谱
|
||||
plt.subplot(2, 2, 2)
|
||||
plt.plot(f, np.abs(y) * 2 / N)
|
||||
plt.title('（b））原始振幅谱')
|
||||
plt.xlabel('频率/Hz', labelpad=10) # 增加标签间距
|
||||
plt.ylabel('振幅', labelpad=10) # 增加标签间距
|
||||
|
||||
# 通过IFFT回到时间域
|
||||
plt.subplot(2, 2, 3)
|
||||
plt.plot(t, yyy)
|
||||
plt.title('（c））傅里叶逆变换')
|
||||
plt.xlabel('高度 (km)', labelpad=10) # 增加标签间距
|
||||
plt.ylabel('温度 (K)', labelpad=10) # 增加标签间距
|
||||
|
||||
# 滤波后的振幅谱
|
||||
plt.subplot(2, 2, 4)
|
||||
plt.plot(f, np.abs(yy) * 2 / N)
|
||||
plt.title(f'（d）滤除波长 < {lamda_low} km, > {lamda_high} km的波')
|
||||
plt.xlabel('频率/Hz', labelpad=10) # 增加标签间距
|
||||
plt.ylabel('振幅', labelpad=10) # 增加标签间距
|
||||
|
||||
# 调整子图之间的边距
|
||||
plt.subplots_adjust(top=0.8, bottom=0.2, left=0.2,
|
||||
right=0.8, hspace=0.3, wspace=0.2)
|
||||
|
||||
plt.show()
|
||||
|
||||
# 绘制原始信号与滤波后信号
|
||||
plt.figure(figsize=(12, 6)) # 调整图形大小
|
||||
plt.plot(t, x, label='原始信号')
|
||||
plt.plot(t, yyy, label='滤波后信号', linestyle='--')
|
||||
plt.title('信号比较')
|
||||
plt.xlabel('高度 (km)', labelpad=10) # 增加标签间距
|
||||
plt.ylabel('温度 (K)', labelpad=10) # 增加标签间距
|
||||
plt.legend()
|
||||
plt.tight_layout() # 调整子图参数以适应图形区域
|
||||
plt.show()
|
||||
|
||||
|
||||
def day_fft_ifft_plotg(cycle_no, ktemp_wn5, ktemp_fft, ktemp_ifft, altitude_min, altitude_max, lamda_low, lamda_high):
|
||||
|
||||
N = len(ktemp_wn5[cycle_no, :])
|
||||
dt = (altitude_max - altitude_min) / (N - 1) # 采样时间间隔
|
||||
n = np.arange(N) # 时间序列索引
|
||||
f = n / (N * dt)
|
||||
t = np.round(np.linspace(altitude_min, altitude_max, N), 2)
|
||||
|
||||
x = ktemp_wn5[cycle_no, :] # 原始扰动温度
|
||||
y = ktemp_fft[cycle_no, :] # 傅里叶变换频谱分析
|
||||
yy = ktemp_fft_lvbo[cycle_no, :] # 滤波后的傅里叶变换频谱分析
|
||||
yyy = ktemp_ifft[cycle_no, :] # 傅里叶逆变换后的扰动温度
|
||||
|
||||
plt.figure(figsize=(15, 10)) # 调整图形大小
|
||||
|
||||
# 计算纵坐标范围
|
||||
temp_limits = (min(min(x), min(yyy)), max(max(x), max(yyy)))
|
||||
amp_limits = (0, max(np.max(np.abs(y) * 2 / N),
|
||||
np.max(np.abs(yy) * 2 / N)))
|
||||
|
||||
# 原始信号的时间序列
|
||||
plt.subplot(2, 2, 1)
|
||||
plt.plot(t, x)
|
||||
plt.title('（a）原始信号')
|
||||
plt.xlabel('高度 (km)', labelpad=10)
|
||||
plt.ylabel('温度 (K)', labelpad=10)
|
||||
plt.ylim(temp_limits) # 设置统一纵坐标范围
|
||||
|
||||
# 原始振幅谱
|
||||
plt.subplot(2, 2, 2)
|
||||
plt.plot(f, np.abs(y) * 2 / N)
|
||||
plt.title('（b）原始振幅谱')
|
||||
plt.xlabel('频率/Hz', labelpad=10)
|
||||
plt.ylabel('振幅', labelpad=10)
|
||||
plt.ylim(amp_limits) # 设置统一纵坐标范围
|
||||
|
||||
# 通过IFFT回到时间域
|
||||
plt.subplot(2, 2, 3)
|
||||
plt.plot(t, yyy)
|
||||
plt.title('（c）傅里叶逆变换')
|
||||
plt.xlabel('高度 (km)', labelpad=10)
|
||||
plt.ylabel('温度 (K)', labelpad=10)
|
||||
plt.ylim(temp_limits) # 设置统一纵坐标范围
|
||||
|
||||
# 滤波后的振幅谱
|
||||
plt.subplot(2, 2, 4)
|
||||
plt.plot(f, np.abs(yy) * 2 / N)
|
||||
plt.title(f'（d）滤除波长 < {lamda_low} km, > {lamda_high} km的波')
|
||||
plt.xlabel('频率/Hz', labelpad=10)
|
||||
plt.ylabel('振幅', labelpad=10)
|
||||
plt.ylim(amp_limits) # 设置统一纵坐标范围
|
||||
|
||||
# 调整子图之间的边距
|
||||
plt.subplots_adjust(top=0.8, bottom=0.2, left=0.2,
|
||||
right=0.8, hspace=0.3, wspace=0.2)
|
||||
plt.show()
|
||||
|
||||
# 绘制原始信号与滤波后信号
|
||||
plt.figure(figsize=(6, 8)) # 调整图形大小
|
||||
plt.plot(x, t, label='原始信号')
|
||||
plt.plot(yyy, t, label='滤波后信号', linestyle='--')
|
||||
plt.title('信号比较')
|
||||
plt.ylabel('高度 (km)', labelpad=10) # 增加标签间距
|
||||
plt.xlabel('温度 (K)', labelpad=10) # 增加标签间距
|
||||
plt.legend()
|
||||
plt.xlim(temp_limits) # 设置统一纵坐标范围
|
||||
plt.tight_layout()
|
||||
plt.show()
|
||||
|
||||
# -----------------------------------------------------------------------------------------------------------------------------
|
||||
# 6-3 示例的按高度的重力波势能变化曲线图
|
||||
|
||||
|
||||
def day_cycle_power_wave_plot(cycle_no, altitude_min, altitude_max, ktemp_Nz, ktemp_Ptz):
|
||||
|
||||
N = len(ktemp_Nz[cycle_no, :])
|
||||
y = np.round(np.linspace(altitude_min, altitude_max, N), 2)
|
||||
x1 = ktemp_Nz[cycle_no, :]
|
||||
x2 = ktemp_Ptz[cycle_no, :]
|
||||
|
||||
plt.figure(figsize=(12, 10)) # 调整图形大小
|
||||
# 原始信号的时间序列
|
||||
plt.subplot(1, 2, 1)
|
||||
plt.plot(x1[::-1], y, label='原始信号')
|
||||
plt.title('（a）Nz')
|
||||
plt.xlabel('Nz', labelpad=10) # 增加标签间距
|
||||
plt.ylabel('高度 (km)', labelpad=10) # 增加标签间距
|
||||
|
||||
# 原始信号的时间序列
|
||||
plt.subplot(1, 2, 2)
|
||||
plt.plot(x2[::-1], y, label='原始信号')
|
||||
plt.title('（b）Ptz')
|
||||
plt.xlabel('Ptz', labelpad=10) # 增加标签间距
|
||||
plt.ylabel('高度 (km)', labelpad=10) # 增加标签间距
|
||||
|
||||
# 调整子图之间的边距
|
||||
plt.subplots_adjust(top=0.8, bottom=0.2, left=0.1,
|
||||
right=0.8, hspace=0.3, wspace=0.2)
|
||||
plt.tight_layout() # 调整子图参数以适应图形区域
|
||||
|
||||
plt.show()
|
||||
# -----------------------------------------------------------------------------------------------------------------------------
|
||||
# 6-4 按日统计的按周期计算的不同高度的重力波势能 平面图
|
||||
|
||||
|
||||
def day_power_wave_plot(altitude_min, altitude_max, ktemp_Nz, ktemp_Ptz):
|
||||
# 假设 ktemp_Nz 和 ktemp_Ptz 以及 altitude_min, altitude_max 已经定义好
|
||||
x = np.arange(ktemp_Nz.shape[0])
|
||||
y = np.round(np.linspace(altitude_min, altitude_max, ktemp_Nz.shape[1]), 2)
|
||||
|
||||
# 创建一个图形，并指定两个子图
|
||||
fig, axs = plt.subplots(1, 2, figsize=(15, 10))
|
||||
|
||||
# 第一幅图 (a) NZ
|
||||
cax1 = axs[0].imshow(ktemp_Nz.T[::-1], aspect='auto', cmap='viridis',
|
||||
origin='lower', extent=[x[0], x[-1], y[0], y[-1]])
|
||||
fig.colorbar(cax1, ax=axs[0]) # 为第一幅图添加颜色条
|
||||
axs[0].set_title('(a) NZ')
|
||||
axs[0].set_ylabel('Height (km)')
|
||||
axs[0].set_xlabel('Cycles')
|
||||
axs[0].set_yticks(np.linspace(30, 90, 7))
|
||||
axs[0].set_yticklabels(np.round(np.linspace(30, 90, 7), 1))
|
||||
axs[0].set_xticks(np.arange(15))
|
||||
axs[0].set_xticklabels(np.arange(15))
|
||||
|
||||
# 第二幅图 (b) PTZ
|
||||
cax2 = axs[1].imshow(ktemp_Ptz.T[::-1], aspect='auto', cmap='viridis',
|
||||
origin='lower', extent=[x[0], x[-1], y[0], y[-1]])
|
||||
fig.colorbar(cax2, ax=axs[1]) # 为第二幅图添加颜色条
|
||||
axs[1].set_title('(b) PTZ')
|
||||
axs[1].set_ylabel('Height (km)')
|
||||
axs[1].set_xlabel('Cycles')
|
||||
axs[1].set_yticks(np.linspace(30, 90, 7))
|
||||
axs[1].set_yticklabels(np.round(np.linspace(30, 90, 7), 1))
|
||||
axs[1].set_xticks(np.arange(15))
|
||||
axs[1].set_xticklabels(np.arange(15))
|
||||
|
||||
# 调整子图之间的边距
|
||||
plt.subplots_adjust(top=0.9, bottom=0.1, left=0.05,
|
||||
right=0.95, hspace=0.3, wspace=0.3)
|
||||
plt.tight_layout() # 调整布局以避免重叠
|
||||
plt.show()
|
||||
# -----------------------------------------------------------------------------------------------------------------------------
|
||||
# 6-5 按月统计的每日重力波势能随天变化的图
|
||||
|
||||
|
||||
def month_power_wave_plot(altitude_min, altitude_max, ktemp_Nz_mon, ktemp_Ptz_mon):
|
||||
if ktemp_Nz_mon and ktemp_Nz_mon[0] is not None:
|
||||
nz_shape = np.array(ktemp_Nz_mon[0]).shape
|
||||
else:
|
||||
nz_shape = (15, 157)
|
||||
if ktemp_Ptz_mon and ktemp_Ptz_mon[0] is not None:
|
||||
ptz_shape = np.array(ktemp_Ptz_mon[0]).shape
|
||||
else:
|
||||
ptz_shape = (15, 157)
|
||||
y = np.round(np.linspace(altitude_min, altitude_max, nz_shape[1]), 2)
|
||||
x = np.arange(len(date_time.data))
|
||||
# 处理 ktemp_Nz_mon
|
||||
ktemp_Nz_plot = np.array([np.mean(day_data if day_data is not None else np.zeros(
|
||||
nz_shape), axis=0) for day_data in ktemp_Nz_mon])
|
||||
ktemp_Ptz_plot = np.array(
|
||||
[np.mean(day_data if day_data is not None else np.zeros(nz_shape), axis=0) for day_data in ktemp_Ptz_mon])
|
||||
# 处理 ktemp_Ptz_mon（以100为界剔除异常值）
|
||||
# ktemp_Ptz_plot = np.array([np.mean(day_data if day_data is not None and np.all(day_data <= 100) else np.zeros(ptz_shape), axis=0) for day_data in ktemp_Ptz_mon])
|
||||
# 创建一个图形，并指定两个子图
|
||||
fig, axs = plt.subplots(1, 2, figsize=(15, 10))
|
||||
|
||||
# 第一幅图 (a) NZ
|
||||
cax1 = axs[0].imshow(ktemp_Nz_plot.T[::-1], aspect='auto', cmap='rainbow', origin='lower',
|
||||
extent=[x[0], x[-1], y[0], y[-1]])
|
||||
fig.colorbar(cax1, ax=axs[0]) # 为第一幅图添加颜色条
|
||||
axs[0].set_title('(a) NZ')
|
||||
axs[0].set_xlabel('Time')
|
||||
axs[0].set_ylabel('Height')
|
||||
axs[0].set_yticks(np.linspace(30, 100, 8))
|
||||
axs[0].set_yticklabels(np.round(np.linspace(30, 100, 8), 1))
|
||||
axs[0].set_xticks(x)
|
||||
axs[0].set_xticklabels(x)
|
||||
|
||||
# 第二幅图 (b) PTZ
|
||||
cax2 = axs[1].imshow(np.log(ktemp_Ptz_plot.T[::-1]), aspect='auto',
|
||||
cmap='rainbow', origin='lower', extent=[x[0], x[-1], y[0], y[-1]])
|
||||
fig.colorbar(cax2, ax=axs[1]) # 为第二幅图添加颜色条
|
||||
axs[1].set_title('(b) PTZ')
|
||||
axs[1].set_xlabel('Time')
|
||||
axs[1].set_ylabel('Height')
|
||||
axs[1].set_yticks(np.linspace(30, 100, 8))
|
||||
axs[1].set_yticklabels(np.round(np.linspace(30, 100, 8), 1))
|
||||
axs[1].set_xticks(x)
|
||||
axs[1].set_xticklabels(x)
|
||||
|
||||
# 调整子图之间的边距
|
||||
plt.subplots_adjust(top=0.9, bottom=0.1, left=0.05,
|
||||
right=0.95, hspace=0.3, wspace=0.3)
|
||||
plt.tight_layout() # 调整布局以避免重叠
|
||||
plt.show()
|
||||
|
||||
|
||||
# -----------------------------------------------------------------------------------------------------------------------------
|
||||
# 6-6 按年统计的每月重力波势能随月变化的图
|
||||
|
||||
|
||||
def year_power_wave_plot(year, path, latitude_min, latitude_max, altitude_min, altitude_max, lamda_low, lamda_high,
|
||||
lvboin):
|
||||
# 假设我们已经从process_yearly_data函数中获取了一年的Nz和Ptz数据
|
||||
results = process_yearly_data(year, path, latitude_min, latitude_max, altitude_min, altitude_max, lamda_low,
|
||||
lamda_high, lvboin)
|
||||
ktemp_Nz_mon_list = results["ktemp_Nz_mon_list"]
|
||||
ktemp_Ptz_mon_list = results["ktemp_Ptz_mon_list"]
|
||||
ktemp_Ptz_mon_list.pop(0)
|
||||
ktemp_Nz_mon_list.pop(0)
|
||||
|
||||
# 准备日期数据，这里假设date_time_list是一年中的所有日期
|
||||
date_time_list = results["date_time_list"]
|
||||
date_time_list.pop(0)
|
||||
date_nums = mdates.date2num(date_time_list) # 将日期转换为matplotlib可以理解的数字格式
|
||||
|
||||
# 获取date_time_list长度作为横坐标新的依据
|
||||
x_ticks_length = len(date_time_list)
|
||||
x_ticks = np.arange(0, x_ticks_length, 30)
|
||||
x_labels = [date_time_list[i] if i < len(
|
||||
date_time_list) else "" for i in x_ticks]
|
||||
|
||||
# 准备高度数据
|
||||
# 假设高度数据有157个点
|
||||
y = np.round(np.linspace(altitude_min, altitude_max, 157), 2)
|
||||
|
||||
# 创建一个图形，并指定两个子图
|
||||
fig, axs = plt.subplots(1, 2, figsize=(15, 10))
|
||||
|
||||
# 处理 ktemp_Nz_mon
|
||||
ktemp_Nz_plot = np.array(
|
||||
[np.mean(day_data if day_data is not None else np.zeros((15, 157)), axis=0) for day_data in ktemp_Nz_mon_list])
|
||||
# 处理 ktemp_Ptz_mon
|
||||
ktemp_Ptz_plot = np.array(
|
||||
[np.mean(day_data if day_data is not None else np.zeros((15, 157)), axis=0) for day_data in ktemp_Ptz_mon_list])
|
||||
# ktemp_Ptz_plot = np.array(
|
||||
# [np.mean(day_data if day_data is not None and np.all(day_data <= 100) else np.zeros((15, 157)), axis=0) for
|
||||
# day_data in ktemp_Ptz_mon_list])
|
||||
|
||||
# 第一幅图 (a) NZ
|
||||
cax1 = axs[0].imshow(ktemp_Nz_plot.T[::-1], aspect='auto', cmap='rainbow', origin='lower',
|
||||
extent=[0, x_ticks_length - 1, y[0], y[-1]])
|
||||
fig.colorbar(cax1, ax=axs[0]) # 为第一幅图添加颜色条
|
||||
axs[0].set_title('(a) NZ')
|
||||
axs[0].set_xlabel('Time')
|
||||
axs[0].set_ylabel('Height')
|
||||
axs[0].set_yticks(np.linspace(30, 100, 8))
|
||||
axs[0].set_yticklabels(np.round(np.linspace(30, 100, 8), 1))
|
||||
axs[0].set_xticks(x_ticks) # 设置新的横坐标刻度
|
||||
axs[0].set_xticklabels(x_labels, rotation=45)
|
||||
|
||||
# 第二幅图 (b) PTZ
|
||||
cax2 = axs[1].imshow(np.log(ktemp_Ptz_plot.T[::-1]), aspect='auto', cmap='rainbow', origin='lower',
|
||||
extent=[0, x_ticks_length - 1, y[0], y[-1]])
|
||||
fig.colorbar(cax2, ax=axs[1]) # 为第二幅图添加颜色条
|
||||
axs[1].set_title('(b) PTZ')
|
||||
axs[1].set_xlabel('Time')
|
||||
axs[1].set_ylabel('Height')
|
||||
axs[1].set_yticks(np.linspace(30, 100, 8))
|
||||
axs[1].set_yticklabels(np.round(np.linspace(30, 100, 8), 1))
|
||||
axs[1].set_xticks(x_ticks) # 设置新的横坐标刻度
|
||||
axs[1].set_xticklabels(x_labels, rotation=45)
|
||||
|
||||
# 调整子图之间的边距
|
||||
plt.subplots_adjust(top=0.9, bottom=0.1, left=0.05,
|
||||
right=0.95, hspace=0.3, wspace=0.3)
|
||||
plt.tight_layout() # 调整布局以避免重叠
|
||||
plt.show()
|
||||
|
||||
|
||||
# -----------------------------------------------------------------------------------------------------------------------------
|
||||
# 7 主程序，运行数据,并输出主要统计分析图 -
|
||||
# -----------------------------------------------------------------------------------------------------------------------------
|
||||
# 7-1 挑选某一天进行运行主要程序
|
||||
file_path = f"{DATA_BASEPATH.saber}\\2016\\SABER_Temp_O3_April2016_v2.0.nc"
|
||||
# day_read=2018113,
|
||||
day_read = 2016094
|
||||
# 初始化某一天、某个纬度、高度范围等参数
|
||||
latitude_min = 30.0
|
||||
latitude_max = 40.0
|
||||
altitude_min = 30.0
|
||||
altitude_max = 100.0
|
||||
lamda_low = 2
|
||||
lamda_high = 15
|
||||
lvboin = True
|
||||
|
||||
(ktemp_cycles, altitude_cycles,
|
||||
ktemp_wn0,
|
||||
ktemp_fit_wn1, ktemp_wn1,
|
||||
ktemp_fit_wn2, ktemp_wn2,
|
||||
ktemp_fit_wn3, ktemp_wn3,
|
||||
ktemp_fit_wn4, ktemp_wn4,
|
||||
ktemp_fit_wn5, ktemp_wn5,
|
||||
ktemp_fft, ktemp_fft_lvbo, ktemp_ifft,
|
||||
ktemp_Nz, ktemp_Ptz) = day_process_main(
|
||||
file_path,
|
||||
# day_read=2018113,
|
||||
day_read,
|
||||
# 初始化某一天、某个纬度、高度范围等参数
|
||||
latitude_min,
|
||||
latitude_max,
|
||||
altitude_min,
|
||||
altitude_max,
|
||||
lamda_low, lamda_high, lvboin)
|
||||
# -----------------------------------------------------------------------------------------------------------------------------
|
||||
# 7-2 挑选某一个月进行运行主要程序
|
||||
(date_time, ktemp_cycles_mon, altitude_cycles_mon, # 某月份 逐日时间 每日不同循环周期温度 不同循环周期高度数据
|
||||
# 距平扰动温度
|
||||
ktemp_wn0_mon,
|
||||
# 波数为1的拟合温度 滤波后的扰动温度
|
||||
ktemp_fit_wn1_mon, ktemp_wn1_mon,
|
||||
# 波数为2的拟合温度 滤波后的扰动温度
|
||||
ktemp_fit_wn2_mon, ktemp_wn2_mon,
|
||||
# 波数为3的拟合温度 滤波后的扰动温度
|
||||
ktemp_fit_wn3_mon, ktemp_wn3_mon,
|
||||
# 波数为4的拟合温度 滤波后的扰动温度
|
||||
ktemp_fit_wn4_mon, ktemp_wn4_mon,
|
||||
# 波数为5的拟合温度 滤波后的扰动温度
|
||||
ktemp_fit_wn5_mon, ktemp_wn5_mon,
|
||||
# 滤波0-5后扰动温度傅里叶频谱分析 滤除波长为2km-15km内后频谱分析 滤波后的扰动温度
|
||||
ktemp_fft_mon, ktemp_fft_lvbo_mon, ktemp_ifft_mon,
|
||||
# 滤波后NZ、PTZ重力波势能指标计算
|
||||
ktemp_Nz_mon, ktemp_Ptz_mon) = mon_process_main(
|
||||
file_path,
|
||||
# 初始化某一天、某个纬度、高度范围等参数
|
||||
latitude_min,
|
||||
latitude_max,
|
||||
altitude_min,
|
||||
altitude_max,
|
||||
lamda_low, lamda_high, lvboin)
|
||||
|
||||
# -----------------------------------------------------------------------------------------------------------------------------
|
||||
# 7-3挑选某一年进行运行主要程序
|
||||
|
||||
|
||||
def filename_read(year, month_num, path):
|
||||
# 月份映射
|
||||
month_map = {
|
||||
1: "January", 2: "February", 3: "March", 4: "April",
|
||||
5: "May", 6: "June", 7: "July", 8: "August",
|
||||
9: "September", 10: "October", 11: "November", 12: "December"}
|
||||
|
||||
# 获取月份名称
|
||||
month = month_map.get(month_num, "Invalid month number")
|
||||
|
||||
# 构造文件路径：E:\SABER\year\SABER_Temp_O3_MonthYear_v2.0.nc
|
||||
file_path = f"{path}\\{year}\\SABER_Temp_O3_{month}{year}_v2.0.nc"
|
||||
|
||||
# 打印路径
|
||||
print(file_path)
|
||||
|
||||
# 返回文件路径
|
||||
return file_path
|
||||
|
||||
|
||||
year = 2018
|
||||
path = f"{DATA_BASEPATH.saber}"
|
||||
# 初始化某一天、某个纬度、高度范围等参数
|
||||
latitude_min = 30.0
|
||||
latitude_max = 40.0
|
||||
altitude_min = 20.0
|
||||
altitude_max = 105.0
|
||||
lamda_low = 2
|
||||
lamda_high = 15
|
||||
lvboin = True
|
||||
# 调用函数处理所有月份
|
||||
results = process_yearly_data(year, path, latitude_min, latitude_max,
|
||||
altitude_min, altitude_max, lamda_low, lamda_high, lvboin)
|
||||
|
||||
# 解包返回的字典
|
||||
date_time_list = results["date_time_list"]
|
||||
ktemp_cycles_mon_list = results["ktemp_cycles_mon_list"]
|
||||
altitude_cycles_mon_list = results["altitude_cycles_mon_list"]
|
||||
ktemp_wn0_mon_list = results["ktemp_wn0_mon_list"]
|
||||
ktemp_fit_wn1_mon_list = results["ktemp_fit_wn1_mon_list"]
|
||||
ktemp_wn1_mon_list = results["ktemp_wn1_mon_list"]
|
||||
ktemp_fit_wn2_mon_list = results["ktemp_fit_wn2_mon_list"]
|
||||
ktemp_wn2_mon_list = results["ktemp_wn2_mon_list"]
|
||||
ktemp_fit_wn3_mon_list = results["ktemp_fit_wn3_mon_list"]
|
||||
ktemp_wn3_mon_list = results["ktemp_wn3_mon_list"]
|
||||
ktemp_fit_wn4_mon_list = results["ktemp_fit_wn4_mon_list"]
|
||||
ktemp_wn4_mon_list = results["ktemp_wn4_mon_list"]
|
||||
ktemp_fit_wn5_mon_list = results["ktemp_fit_wn5_mon_list"]
|
||||
ktemp_wn5_mon_list = results["ktemp_wn5_mon_list"]
|
||||
ktemp_fft_mon_list = results["ktemp_fft_mon_list"]
|
||||
ktemp_fft_lvbo_mon_list = results["ktemp_fft_lvbo_mon_list"]
|
||||
ktemp_ifft_mon_list = results["ktemp_ifft_mon_list"]
|
||||
ktemp_Nz_mon_list = results["ktemp_Nz_mon_list"]
|
||||
ktemp_Ptz_mon_list = results["ktemp_Ptz_mon_list"]
|
||||
# -----------------------------------------------------------------------------------------------------------------------------
|
||||
# 7-3 绘制不同结果图
|
||||
height_no = 1
|
||||
day_fit_wave_plotg(height_no, ktemp_wn0, ktemp_fit_wn1, ktemp_wn1, ktemp_fit_wn2, ktemp_wn2,
|
||||
ktemp_fit_wn3, ktemp_wn3, ktemp_fit_wn4, ktemp_wn4, ktemp_fit_wn5, ktemp_wn5)
|
||||
cycle_no = 1
|
||||
day_fft_ifft_plotg(cycle_no, ktemp_wn5, ktemp_fft, ktemp_ifft,
|
||||
altitude_min, altitude_max, lamda_low, lamda_high)
|
||||
day_cycle_power_wave_plot(cycle_no, altitude_min,
|
||||
altitude_max, ktemp_Nz, ktemp_Ptz)
|
||||
day_power_wave_plot(altitude_min, altitude_max, ktemp_Nz, ktemp_Ptz)
|
||||
month_power_wave_plot(altitude_min, altitude_max, ktemp_Nz_mon, ktemp_Ptz_mon)
|
||||
year_power_wave_plot(year, path, latitude_min, latitude_max,
|
||||
altitude_min, altitude_max, lamda_low, lamda_high, lvboin)
|
||||
@ -1,842 +0,0 @@
|
||||
# saber 重力波
|
||||
|
||||
from io import BytesIO
|
||||
import netCDF4 as nc
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
import pandas as pd
|
||||
from scipy.optimize import curve_fit
|
||||
# from matplotlib.colors import LinearSegmentedColormap
|
||||
# 设置字体为支持中文的字体
|
||||
plt.rcParams['font.family'] = 'SimHei' # 设置为黑体（需要你的环境中有该字体）
|
||||
plt.rcParams['axes.unicode_minus'] = False # 解决负号'-'显示为方块的问题
|
||||
|
||||
# ----------------------------------------------------------------------------------------------------------------------------
|
||||
# 1---打开文件并读取不同变量数据 -
|
||||
# ----------------------------------------------------------------------------------------------------------------------------
|
||||
|
||||
|
||||
def data_nc_load(file_path):
|
||||
|
||||
dataset = nc.Dataset(file_path, 'r')
|
||||
|
||||
# 纬度数据，二维数组形状为(42820,379) 42820为事件，379则为不同高度
|
||||
tplatitude = dataset.variables['tplatitude'][:, :]
|
||||
tplongitude = dataset.variables['tplongitude'][:, :] # 经度数据
|
||||
tpaltitude = dataset.variables['tpaltitude'][:, :] # 高度，二维数组形状为(42820,379)
|
||||
time = dataset.variables['time'][:, :] # 二维数组形状为(42820,379)
|
||||
date = dataset.variables['date'][:]
|
||||
date_time = np.unique(date) # 输出数据时间信息
|
||||
ktemp = dataset.variables['ktemp'][:, :] # 温度数据，二维数组形状为(42820,379)
|
||||
|
||||
return dataset, tplatitude, tplongitude, tpaltitude, ktemp, time, date, date_time
|
||||
# ----------------------------------------------------------------------------------------------------------------------------
|
||||
|
||||
# ----------------------------------------------------------------------------------------------------------------------------
|
||||
# 2---筛选某一天、某个纬度和高度范围15个不同cycle的温度数据 -
|
||||
# ----------------------------------------------------------------------------------------------------------------------------
|
||||
# 2-1 读取某一天的所有事件及其对应的纬度数据
|
||||
|
||||
|
||||
def day_data_read(date, day_read, tplatitude):
|
||||
|
||||
events = np.where(date == day_read)[0] # 读取筛选天的事件编号 4294-5714位置，从0开始编号
|
||||
time_events = date[date == day_read] # 读取筛选天的事件编号 4294-5714位置，从0开始编号
|
||||
latitudes = tplatitude[events, 189] # 输出每个事件中间位置 即第189个经纬度
|
||||
df = pd.DataFrame([ # 创建一个包含事件编号、纬度的列表，共1421个事件
|
||||
{'time': time, 'event': event, 'latitude': lat}
|
||||
for time, event, lat in zip(time_events, events, latitudes)])
|
||||
|
||||
# print(df.head()) # 打印前几行数据以检查
|
||||
|
||||
return df
|
||||
# ----------------------------------------------------------------------------------------------------------------------------
|
||||
# 2-2 将事件按照卫星轨迹周期进行输出处理，并输出落在纬度范围内的每个周期的事件的行号
|
||||
|
||||
|
||||
def data_cycle_identify(df, latitude_min, latitude_max):
|
||||
|
||||
cycles = [] # 存储每个周期的事件编号列表
|
||||
# 存储当前周期的事件编号
|
||||
current_cycle_events = []
|
||||
prev_latitude = None
|
||||
|
||||
# 遍历DataFrame中的每一行以识别周期和筛选事件
|
||||
for index, row in df.iterrows():
|
||||
current_event = int(row['event'])
|
||||
current_latitude = row['latitude']
|
||||
|
||||
if prev_latitude is not None and prev_latitude < 0 and current_latitude >= 0: # 检查是否是新周期的开始（纬度从负变正，且首次变正）
|
||||
# 重置当前周期的事件编号列表
|
||||
current_cycle_events = []
|
||||
|
||||
if latitude_min <= current_latitude <= latitude_max: # 如果事件的纬度在指定范围内，添加到当前周期的事件编号列表
|
||||
current_cycle_events.append(current_event)
|
||||
|
||||
if prev_latitude is not None and prev_latitude >= 0 and current_latitude < 0: # 检查是否是周期的结束（纬度从正变负）
|
||||
|
||||
# 添加当前周期的事件编号列表到周期列表
|
||||
if current_cycle_events: # 确保当前周期有事件编号
|
||||
cycles.append(current_cycle_events)
|
||||
current_cycle_events = [] # 重置当前周期的事件编号列表
|
||||
prev_latitude = current_latitude
|
||||
|
||||
if current_cycle_events: # 处理最后一个周期，如果存在的话
|
||||
cycles.append(current_cycle_events)
|
||||
|
||||
print(f"一天周期为 {len(cycles)}")
|
||||
for cycle_index, cycle in enumerate(cycles, start=0):
|
||||
# 屏幕显示每个循环周期的事件
|
||||
print(f"周期 {cycle_index} 包含事件个数: {len(cycle)} 具体事件为: {cycle} ")
|
||||
|
||||
return cycles
|
||||
# ----------------------------------------------------------------------------------------------------------------------------
|
||||
# 2-3---按照循环周期合并同周期数据，并输出处理后的温度数据、对应的高度数据
|
||||
|
||||
|
||||
def data_cycle_generate(cycles, ktemp, tpaltitude, altitude_min, altitude_max):
|
||||
|
||||
# 初始化列表存储每个周期的温度
|
||||
ktemp_cycles = []
|
||||
# 初始化每个循环周期的高度数据
|
||||
altitude_cycles = []
|
||||
print(f"周期数为 {len(cycles)}")
|
||||
if cycles is not None: # 检查周期是否有事件编号
|
||||
for event in cycles:
|
||||
# 获取每个周期各个事件的ktemp数据
|
||||
ktemp_cycles_events = np.array(ktemp[event, :])
|
||||
ktemp_cycles_events[np.logical_or(
|
||||
ktemp_cycles_events == -999, ktemp_cycles_events > 999)] = np.NaN # 缺失值处理，避免影响结果
|
||||
ktemp_cycles_mean = np.nanmean(
|
||||
ktemp_cycles_events, axis=0) # 对所有周期的 ktemp 数据取均值
|
||||
|
||||
# altitude_cycles_events =np.array(tpaltitude[event,:]) # 获取每个周期各个事件的tpaltitude数据
|
||||
# altitude_cycles_mean = np.nanmean(altitude_cycles_events, axis=0) # 对所有周期的 altitude 数据取均值
|
||||
# 使用第一个的高度来表征所有的
|
||||
altitude_cycles_mean = tpaltitude[event[0], :]
|
||||
altitude_indices = np.where((altitude_cycles_mean >= altitude_min) & (
|
||||
altitude_cycles_mean <= altitude_max))[0]
|
||||
|
||||
ktemp_cycles.append(np.array(ktemp_cycles_mean[altitude_indices]))
|
||||
altitude_cycles.append(
|
||||
np.array(altitude_cycles_mean[altitude_indices]))
|
||||
|
||||
# 找到最短列表的长度
|
||||
min_length = min(len(arr) for arr in ktemp_cycles)
|
||||
# 创建新的列表，将每个子列表截断为最短长度
|
||||
ktemp_cycles = np.vstack([arr[:min_length] for arr in ktemp_cycles])
|
||||
altitude_cycles = np.vstack([arr[:min_length] for arr in altitude_cycles])
|
||||
|
||||
return ktemp_cycles, altitude_cycles
|
||||
# ----------------------------------------------------------------------------------------------------------------------------
|
||||
|
||||
# ----------------------------------------------------------------------------------------------------------------------------
|
||||
# 4---高度相同下不同循环周期数据的波拟合和滤波处理 -
|
||||
# ----------------------------------------------------------------------------------------------------------------------------
|
||||
|
||||
|
||||
# 对输入数据按照行（即纬向）进行波数为k的滤波,数据为15*157
|
||||
def fit_wave(ktemp_wn0, k):
|
||||
|
||||
wave_fit = []
|
||||
|
||||
def single_harmonic(x, A, phi, C, k):
|
||||
return A * np.sin(2 * np.pi / (15/k) * x + phi) + C
|
||||
|
||||
# 数据转置并对每行进行操作，按照同高度数据进行处理
|
||||
for rtemp in ktemp_wn0.T:
|
||||
# 为当前高度层创建索引数组
|
||||
indices = np.arange(rtemp.size)
|
||||
def fit_temp(x, A, phi, C): return single_harmonic(
|
||||
x, A, phi, C, k) # 定义只拟合 A, phi, C 的 lambda 函数，k 固定
|
||||
params, params_covariance = curve_fit(
|
||||
fit_temp, indices, rtemp) # 使用 curve_fit 进行拟合
|
||||
# 提取拟合参数 A, phi, C
|
||||
A, phi, C = params
|
||||
# 使用拟合参数计算 wn3
|
||||
rtemp1 = single_harmonic(indices, A, phi, C, k)
|
||||
# 存储拟合参数和拟合曲线
|
||||
wave_fit.append(np.array(rtemp1))
|
||||
wave_fit = np.vstack(wave_fit)
|
||||
|
||||
wave_fit = wave_fit.T
|
||||
wave_wn = ktemp_wn0-wave_fit
|
||||
|
||||
return wave_fit, wave_wn
|
||||
# ----------------------------------------------------------------------------------------------------------------------------
|
||||
|
||||
|
||||
# ----------------------------------------------------------------------------------------------------------------------------
|
||||
# 4---同周期下不同高度数据的波拟合和滤波处理 -
|
||||
# ----------------------------------------------------------------------------------------------------------------------------
|
||||
# 对输入数据按照列（即纬向）进行滤波，滤除波长2和15km以内的波, 数据为15*157
|
||||
def fft_ifft_wave(ktemp_wn5, lamda_low, lamda_high, altitude_min, altitude_max, lvboin):
|
||||
|
||||
ktemp_fft = []
|
||||
ktemp_ifft = []
|
||||
ktemp_fft_lvbo = []
|
||||
|
||||
for rtemp in ktemp_wn5:
|
||||
# 采样点数或长度
|
||||
N = len(rtemp)
|
||||
# 采样时间间隔，其倒数等于采用频率，以1km为标准尺度等同于1s，假设波的速度为1km/s
|
||||
dt = (altitude_max-altitude_min)/(N-1)
|
||||
# 时间序列索引
|
||||
n = np.arange(N)
|
||||
# # t = altitude_min + n * dt # 时间向量
|
||||
# t = np.round(np.linspace(altitude_min, altitude_max, N),2)
|
||||
# 频率索引向量
|
||||
f = n / (N * dt)
|
||||
# 对输入信号进行傅里叶变换
|
||||
y = np.fft.fft(rtemp)
|
||||
|
||||
# f_low = 2*np.pi / lamda_high # 定义波长滤波范围（以频率计算） # 高频截止频率
|
||||
# f_high = 2*np.pi / lamda_low
|
||||
|
||||
# 定义波长滤波范围（以频率计算） # 高频截止频率
|
||||
f_low = 1 / lamda_high
|
||||
# 低频截止频率
|
||||
f_high = 1 / lamda_low
|
||||
|
||||
# 创建滤波后的频域信号
|
||||
yy = y.copy()
|
||||
|
||||
# 使用逻辑索引过滤特定频段（未确定）
|
||||
if lvboin:
|
||||
freq_filter = (f > f_low) & (f < f_high) # 创建逻辑掩码
|
||||
else:
|
||||
freq_filter = (f < f_low) | (f > f_high) # 创建逻辑掩码
|
||||
|
||||
yy[freq_filter] = 0 # 过滤掉指定频段
|
||||
yy_ifft = np.real(np.fft.ifft(yy))
|
||||
|
||||
ktemp_fft.append(y)
|
||||
# 存储拟合参数和拟合曲线
|
||||
ktemp_ifft.append(np.array(yy_ifft))
|
||||
ktemp_fft_lvbo.append(yy)
|
||||
|
||||
ktemp_fft = np.vstack(ktemp_fft)
|
||||
ktemp_ifft = np.vstack(ktemp_ifft)
|
||||
ktemp_fft_lvbo = np.vstack(ktemp_fft_lvbo)
|
||||
|
||||
return ktemp_fft, ktemp_fft_lvbo, ktemp_ifft
|
||||
# ----------------------------------------------------------------------------------------------------------------------------
|
||||
# ----------------------------------------------------------------------------------------------------------------------------
|
||||
# 5---同周期下不同高度数据的BT_z背景位等指标计算 -
|
||||
# ----------------------------------------------------------------------------------------------------------------------------
|
||||
# ----------------------------------------------------------------------------------
|
||||
|
||||
|
||||
def power_indices(ktemp_cycles, ktemp_wn5, ktemp_ifft, altitude_min, altitude_max):
|
||||
|
||||
# 定义Brunt-Väisälä频率计算函数
|
||||
def brunt_vaisala_frequency(g, BT_z, c_p, height):
|
||||
# 计算位温随高度的变化率
|
||||
dBT_z_dz = np.gradient(BT_z, height)
|
||||
# 计算 Brunt-Väisälä 频率
|
||||
return np.sqrt((g / BT_z) * ((g / c_p) + dBT_z_dz))
|
||||
|
||||
# 定义势能计算函数
|
||||
def calculate_gravitational_potential_energy(g, BT_z, N_z, PT_z):
|
||||
return 0.5 * ((g / N_z) ** 2) * ((PT_z / BT_z) ** 2)
|
||||
|
||||
# 重力加速度
|
||||
g = 9.81
|
||||
# 空气比热容
|
||||
c_p = 1004.5
|
||||
|
||||
height = np.linspace(altitude_min, altitude_max,
|
||||
ktemp_cycles.shape[1]) * 1000 # 高度
|
||||
background = ktemp_cycles - ktemp_wn5
|
||||
|
||||
# 初始化结果矩阵
|
||||
N_z_matrix = []
|
||||
# 初始化结果矩阵
|
||||
PT_z_matrix = []
|
||||
|
||||
# 循环处理background和filtered_perturbation所有行
|
||||
for i in range(background.shape[0]):
|
||||
BT_z = np.array(background[i])
|
||||
# 滤波后的扰动
|
||||
PT_z = np.array(ktemp_ifft[i])
|
||||
|
||||
# 调用Brunt-Väisälä频率函数
|
||||
N_z = brunt_vaisala_frequency(g, BT_z, c_p, height)
|
||||
PT_z = calculate_gravitational_potential_energy(
|
||||
g, BT_z, N_z, PT_z) # 调用势能函数
|
||||
N_z_matrix.append(N_z)
|
||||
PT_z_matrix.append(PT_z)
|
||||
|
||||
ktemp_Nz = np.vstack(N_z_matrix)
|
||||
ktemp_Ptz = np.vstack(PT_z_matrix)
|
||||
return ktemp_Nz, ktemp_Ptz
|
||||
# ----------------------------------------------------------------------------------------------------------------------------
|
||||
# 5---main程序环节 -
|
||||
# ----------------------------------------------------------------------------------------------------------------------------
|
||||
# 按日处理资料
|
||||
|
||||
|
||||
def day_process_main(file_path, day_read, latitude_min, latitude_max, altitude_min, altitude_max, lamda_low, lamda_high, lvboin):
|
||||
|
||||
dataset, tplatitude, tplongitude, tpaltitude, ktemp, time, date, date_time = data_nc_load(
|
||||
file_path)
|
||||
# 2018年的第94天 # 4月4日的日期，以年份和年内的第几天表示
|
||||
df = day_data_read(date, day_read, tplatitude)
|
||||
cycles = data_cycle_identify(df, latitude_min, latitude_max)
|
||||
ktemp_cycles, altitude_cycles = data_cycle_generate(
|
||||
cycles, ktemp, tpaltitude, altitude_min, altitude_max)
|
||||
# 按照纬向计算平均温度 wn0_temp
|
||||
ktemp_wn0 = ktemp_cycles - np.mean(ktemp_cycles, axis=0)
|
||||
ktemp_fit_wn1, ktemp_wn1 = fit_wave(ktemp_wn0, 1)
|
||||
ktemp_fit_wn2, ktemp_wn2 = fit_wave(ktemp_wn1, 2)
|
||||
ktemp_fit_wn3, ktemp_wn3 = fit_wave(ktemp_wn2, 3)
|
||||
ktemp_fit_wn4, ktemp_wn4 = fit_wave(ktemp_wn3, 4)
|
||||
ktemp_fit_wn5, ktemp_wn5 = fit_wave(ktemp_wn4, 5)
|
||||
ktemp_fft, ktemp_fft_lvbo, ktemp_ifft = fft_ifft_wave(
|
||||
ktemp_wn5, lamda_low, lamda_high, altitude_min, altitude_max, lvboin)
|
||||
ktemp_Nz, ktemp_Ptz = power_indices(
|
||||
ktemp_cycles, ktemp_wn5, ktemp_ifft, altitude_min, altitude_max)
|
||||
|
||||
return ktemp_cycles, altitude_cycles, ktemp_wn0, ktemp_fit_wn1, ktemp_wn1, ktemp_fit_wn2, ktemp_wn2, ktemp_fit_wn3, ktemp_wn3, ktemp_fit_wn4, ktemp_wn4, ktemp_fit_wn5, ktemp_wn5, ktemp_fft, ktemp_fft_lvbo, ktemp_ifft, ktemp_Nz, ktemp_Ptz
|
||||
# ----------------------------------------------------------------------------------------------------------------------------
|
||||
# 按文件中单日处理资料
|
||||
|
||||
|
||||
def day_process_maing(dataset, tplatitude, tplongitude, tpaltitude, ktemp, time, date, date_time, day_read, latitude_min, latitude_max, altitude_min, altitude_max, lamda_low, lamda_high, lvboin):
|
||||
|
||||
# dataset, tplatitude, tplongitude, tpaltitude, ktemp, time, date, date_time = data_nc_load (file_path)
|
||||
# 2018年的第94天 # 4月4日的日期，以年份和年内的第几天表示
|
||||
df = day_data_read(date, day_read, tplatitude)
|
||||
cycles = data_cycle_identify(df, latitude_min, latitude_max)
|
||||
ktemp_cycles, altitude_cycles = data_cycle_generate(
|
||||
cycles, ktemp, tpaltitude, altitude_min, altitude_max)
|
||||
# 按照纬向计算平均温度 wn0_temp
|
||||
ktemp_wn0 = ktemp_cycles - np.mean(ktemp_cycles, axis=0)
|
||||
ktemp_fit_wn1, ktemp_wn1 = fit_wave(ktemp_wn0, 1)
|
||||
ktemp_fit_wn2, ktemp_wn2 = fit_wave(ktemp_wn1, 2)
|
||||
ktemp_fit_wn3, ktemp_wn3 = fit_wave(ktemp_wn2, 3)
|
||||
ktemp_fit_wn4, ktemp_wn4 = fit_wave(ktemp_wn3, 4)
|
||||
ktemp_fit_wn5, ktemp_wn5 = fit_wave(ktemp_wn4, 5)
|
||||
ktemp_fft, ktemp_fft_lvbo, ktemp_ifft = fft_ifft_wave(
|
||||
ktemp_wn5, lamda_low, lamda_high, altitude_min, altitude_max, lvboin)
|
||||
ktemp_Nz, ktemp_Ptz = power_indices(
|
||||
ktemp_cycles, ktemp_wn5, ktemp_ifft, altitude_min, altitude_max)
|
||||
|
||||
return ktemp_cycles, altitude_cycles, ktemp_wn0, ktemp_fit_wn1, ktemp_wn1, ktemp_fit_wn2, ktemp_wn2, ktemp_fit_wn3, ktemp_wn3, ktemp_fit_wn4, ktemp_wn4, ktemp_fit_wn5, ktemp_wn5, ktemp_fft, ktemp_fft_lvbo, ktemp_ifft, ktemp_Nz, ktemp_Ptz
|
||||
|
||||
# ----------------------------------------------------------------------------------------------------------------------------
|
||||
# 按月处理资料
|
||||
|
||||
|
||||
def mon_process_main(file_path, latitude_min, latitude_max, altitude_min, altitude_max, lamda_low, lamda_high, lvboin):
|
||||
|
||||
dataset, tplatitude, tplongitude, tpaltitude, ktemp, time, date, date_time = data_nc_load(
|
||||
file_path)
|
||||
|
||||
ktemp_cycles_mon = []
|
||||
altitude_cycles_mon = []
|
||||
ktemp_wn0_mon = []
|
||||
ktemp_fit_wn1_mon = []
|
||||
ktemp_wn1_mon = []
|
||||
ktemp_fit_wn2_mon = []
|
||||
ktemp_wn2_mon = []
|
||||
ktemp_fit_wn3_mon = []
|
||||
ktemp_wn3_mon = []
|
||||
ktemp_fit_wn4_mon = []
|
||||
ktemp_wn4_mon = []
|
||||
ktemp_fit_wn5_mon = []
|
||||
ktemp_wn5_mon = []
|
||||
ktemp_fft_mon = []
|
||||
ktemp_fft_lvbo_mon = []
|
||||
ktemp_ifft_mon = []
|
||||
ktemp_Nz_mon = []
|
||||
ktemp_Ptz_mon = []
|
||||
|
||||
for day_read in date_time:
|
||||
print(f"读取日期 {day_read}")
|
||||
(ktemp_cycles0, altitude_cycles0,
|
||||
ktemp_wn00,
|
||||
ktemp_fit_wn10, ktemp_wn10,
|
||||
ktemp_fit_wn20, ktemp_wn20,
|
||||
ktemp_fit_wn30, ktemp_wn30,
|
||||
ktemp_fit_wn40, ktemp_wn40,
|
||||
ktemp_fit_wn50, ktemp_wn50,
|
||||
ktemp_fft0, ktemp_fft_lvbo0, ktemp_ifft0,
|
||||
ktemp_Nz0, ktemp_Ptz0) = day_process_maing(
|
||||
dataset, tplatitude, tplongitude, tpaltitude, ktemp, time, date, date_time,
|
||||
day_read, latitude_min, latitude_max, altitude_min, altitude_max, lamda_low, lamda_high, lvboin)
|
||||
|
||||
ktemp_cycles_mon.append(ktemp_cycles0)
|
||||
altitude_cycles_mon.append(altitude_cycles0)
|
||||
ktemp_wn0_mon.append(ktemp_wn00)
|
||||
ktemp_fit_wn1_mon.append(ktemp_fit_wn50)
|
||||
ktemp_wn1_mon.append(ktemp_wn50)
|
||||
ktemp_fit_wn2_mon.append(ktemp_fit_wn50)
|
||||
ktemp_wn2_mon.append(ktemp_wn50)
|
||||
ktemp_fit_wn3_mon.append(ktemp_fit_wn50)
|
||||
ktemp_wn3_mon.append(ktemp_wn50)
|
||||
ktemp_fit_wn4_mon.append(ktemp_fit_wn50)
|
||||
ktemp_wn4_mon.append(ktemp_wn50)
|
||||
ktemp_fit_wn5_mon.append(ktemp_fit_wn50)
|
||||
ktemp_wn5_mon.append(ktemp_wn50)
|
||||
ktemp_fft_mon.append(ktemp_fft0)
|
||||
ktemp_fft_lvbo_mon.append(ktemp_fft_lvbo0)
|
||||
ktemp_ifft_mon.append(ktemp_ifft0)
|
||||
ktemp_Nz_mon.append(ktemp_Nz0)
|
||||
ktemp_Ptz_mon.append(ktemp_Ptz0)
|
||||
|
||||
return [date_time, ktemp_cycles_mon, altitude_cycles_mon, # 某月份 逐日时间 每日不同循环周期温度 不同循环周期高度数据
|
||||
# 距平扰动温度
|
||||
ktemp_wn0_mon,
|
||||
# 波数为1的拟合温度 滤波后的扰动温度
|
||||
ktemp_fit_wn1_mon, ktemp_wn1_mon,
|
||||
# 波数为2的拟合温度 滤波后的扰动温度
|
||||
ktemp_fit_wn2_mon, ktemp_wn2_mon,
|
||||
# 波数为3的拟合温度 滤波后的扰动温度
|
||||
ktemp_fit_wn3_mon, ktemp_wn3_mon,
|
||||
# 波数为4的拟合温度 滤波后的扰动温度
|
||||
ktemp_fit_wn4_mon, ktemp_wn4_mon,
|
||||
# 波数为5的拟合温度 滤波后的扰动温度
|
||||
ktemp_fit_wn5_mon, ktemp_wn5_mon,
|
||||
# 滤波0-5后扰动温度傅里叶频谱分析 滤除波长为2km-15km内后频谱分析 滤波后的扰动温度
|
||||
ktemp_fft_mon, ktemp_fft_lvbo_mon, ktemp_ifft_mon,
|
||||
ktemp_Nz_mon, ktemp_Ptz_mon] # 滤波后NZ、PTZ重力波势能指标计算
|
||||
|
||||
# -----------------------------------------------------------------------------------------------------------------------------
|
||||
# 6-主要统计分析图绘制 -
|
||||
# -----------------------------------------------------------------------------------------------------------------------------
|
||||
# 6-1 示例的逐层滤波效果图---不同波数 曲线图
|
||||
|
||||
|
||||
day_fit_wave_modes = ["k=1", "k=2", "k=3", "k=4", "k=5", "滤波"]
|
||||
|
||||
|
||||
def day_fit_wave_plot(height_no, ktemp_wn0, ktemp_fit_wn1, ktemp_wn1, ktemp_fit_wn2, ktemp_wn2, ktemp_fit_wn3, ktemp_wn3, ktemp_fit_wn4, ktemp_wn4, ktemp_fit_wn5, ktemp_wn5, mode):
|
||||
|
||||
N = len(ktemp_wn0[:, height_no])
|
||||
# 循环周期索引
|
||||
x = np.arange(N)
|
||||
|
||||
y1_1 = ktemp_wn0[:, height_no]
|
||||
y1_2 = ktemp_fit_wn1[:, height_no]
|
||||
|
||||
y2_1 = ktemp_wn1[:, height_no]
|
||||
y2_2 = ktemp_fit_wn2[:, height_no]
|
||||
|
||||
y3_1 = ktemp_wn2[:, height_no]
|
||||
y3_2 = ktemp_fit_wn3[:, height_no]
|
||||
|
||||
y4_1 = ktemp_wn3[:, height_no]
|
||||
y4_2 = ktemp_fit_wn4[:, height_no]
|
||||
|
||||
y5_1 = ktemp_wn4[:, height_no]
|
||||
y5_2 = ktemp_fit_wn5[:, height_no]
|
||||
|
||||
y6 = ktemp_wn5[:, height_no]
|
||||
|
||||
if mode not in ["k=1", "k=2", "k=3", "k=4", "k=5", "滤波"]:
|
||||
raise ValueError(
|
||||
"mode 参数应为 'k=1', 'k=2', 'k=3', 'k=4', 'k=5', 'lvbo' 中的一个")
|
||||
|
||||
plt.figure(figsize=(16, 10)) # 调整图形大小
|
||||
# 原始信号的时间序列
|
||||
if mode == "k=1":
|
||||
plt.plot(x, y1_1, label='原始信号')
|
||||
plt.plot(x, y1_2, label='拟合信号', linestyle='--')
|
||||
plt.title('（a）波数k=1')
|
||||
plt.xlabel('Cycles', labelpad=10) # 增加标签间距
|
||||
plt.ylabel('温度 (°C)', labelpad=10) # 增加标签间距
|
||||
plt.legend()
|
||||
plt.tight_layout() # 调整子图参数以适应图形区域
|
||||
|
||||
if mode == "k=2":
|
||||
plt.plot(x, y2_1, label='原始信号')
|
||||
plt.plot(x, y2_2, label='拟合信号', linestyle='--')
|
||||
plt.title('（b）波数k=2')
|
||||
plt.xlabel('Cycles', labelpad=10) # 增加标签间距
|
||||
plt.ylabel('温度 (°C)', labelpad=10) # 增加标签间距
|
||||
plt.legend()
|
||||
plt.tight_layout() # 调整子图参数以适应图形区域
|
||||
|
||||
if mode == "k=3":
|
||||
plt.plot(x, y3_1, label='原始信号')
|
||||
plt.plot(x, y3_2, label='拟合信号', linestyle='--')
|
||||
plt.title('（c）波数k=3')
|
||||
plt.xlabel('Cycles', labelpad=10) # 增加标签间距
|
||||
plt.ylabel('温度 (°C)', labelpad=10) # 增加标签间距
|
||||
plt.legend()
|
||||
plt.tight_layout() # 调整子图参数以适应图形区域
|
||||
|
||||
if mode == "k=4":
|
||||
plt.plot(x, y4_1, label='原始信号')
|
||||
plt.plot(x, y4_2, label='拟合信号', linestyle='--')
|
||||
plt.title('（d）波数k=4')
|
||||
plt.xlabel('Cycles', labelpad=10) # 增加标签间距
|
||||
plt.ylabel('温度 (°C)', labelpad=10) # 增加标签间距
|
||||
plt.legend()
|
||||
plt.tight_layout() # 调整子图参数以适应图形区域
|
||||
|
||||
if mode == "k=5":
|
||||
plt.plot(x, y5_1, label='原始信号')
|
||||
plt.plot(x, y5_2, label='拟合信号', linestyle='--')
|
||||
plt.title('（e）波数k=5')
|
||||
plt.xlabel('Cycles', labelpad=10) # 增加标签间距
|
||||
plt.ylabel('温度 (°C)', labelpad=10) # 增加标签间距
|
||||
plt.legend()
|
||||
plt.tight_layout() # 调整子图参数以适应图形区域
|
||||
|
||||
if mode == "滤波":
|
||||
plt.plot(x, y6, label='滤波信号')
|
||||
plt.title('（f）滤波1-5后信号')
|
||||
plt.xlabel('Cycles', labelpad=10) # 增加标签间距
|
||||
plt.ylabel('温度 (°C)', labelpad=10) # 增加标签间距
|
||||
plt.tight_layout() # 调整子图参数以适应图形区域
|
||||
|
||||
|
||||
def day_fit_wave_plotg(height_no, ktemp_wn0, ktemp_fit_wn1, ktemp_wn1, ktemp_fit_wn2,
|
||||
ktemp_wn2, ktemp_fit_wn3, ktemp_wn3, ktemp_fit_wn4,
|
||||
ktemp_wn4, ktemp_fit_wn5, ktemp_wn5):
|
||||
|
||||
N = len(ktemp_wn0[:, height_no])
|
||||
x = np.arange(N)
|
||||
|
||||
y1_1, y1_2 = ktemp_wn0[:, height_no], ktemp_fit_wn1[:, height_no]
|
||||
y2_1, y2_2 = ktemp_wn1[:, height_no], ktemp_fit_wn2[:, height_no]
|
||||
y3_1, y3_2 = ktemp_wn2[:, height_no], ktemp_fit_wn3[:, height_no]
|
||||
y4_1, y4_2 = ktemp_wn3[:, height_no], ktemp_fit_wn4[:, height_no]
|
||||
y5_1, y5_2 = ktemp_wn4[:, height_no], ktemp_fit_wn5[:, height_no]
|
||||
y6 = ktemp_wn5[:, height_no]
|
||||
|
||||
plt.figure(figsize=(16, 10))
|
||||
|
||||
y_limits = (min(min(y1_1), min(y2_1), min(y3_1), min(y4_1), min(y5_1), min(y6)),
|
||||
max(max(y1_1), max(y2_1), max(y3_1), max(y4_1), max(y5_1), max(y6)))
|
||||
|
||||
for i, (y1, y2) in enumerate([(y1_1, y1_2), (y2_1, y2_2), (y3_1, y3_2),
|
||||
(y4_1, y4_2), (y5_1, y5_2), (y6, None)]):
|
||||
plt.subplot(2, 3, i + 1)
|
||||
plt.plot(x, y1, label='原始信号')
|
||||
if y2 is not None:
|
||||
plt.plot(x, y2, label='拟合信号', linestyle='--')
|
||||
plt.title(f'（{"abcdef"[i]}）波数k={i + 1 if i < 5 else "滤波"}')
|
||||
plt.xlabel('Cycles', labelpad=10)
|
||||
plt.ylabel('温度 (°C)', labelpad=10)
|
||||
plt.legend()
|
||||
plt.xticks(x) # 设置横坐标为整数
|
||||
plt.ylim(y_limits) # 设置统一纵坐标范围
|
||||
plt.tight_layout()
|
||||
|
||||
plt.subplots_adjust(top=0.8, bottom=0.2, left=0.1,
|
||||
right=0.8, hspace=0.3, wspace=0.2)
|
||||
plt.show()
|
||||
|
||||
|
||||
# -----------------------------------------------------------------------------------------------------------------------------
|
||||
# 6-2 示例的高度滤波处理--不同循环周期 曲线图
|
||||
day_fft_ifft_modes = ["原始信号", "fft", "ifft", "滤波", "比较"]
|
||||
|
||||
|
||||
def day_fft_ifft_plot(cycle_no, ktemp_wn5, ktemp_fft, ktemp_ifft, altitude_min, altitude_max, lamda_low, lamda_high, ktemp_fft_lvbo, mode):
|
||||
|
||||
N = len(ktemp_wn5[cycle_no, :])
|
||||
# 采样时间间隔，其倒数等于采用频率，以1km为标准尺度等同于1s，假设波的速度为1km/s
|
||||
dt = (altitude_max-altitude_min)/(N-1)
|
||||
# 时间序列索引
|
||||
n = np.arange(N)
|
||||
f = n / (N * dt)
|
||||
t = np.round(np.linspace(altitude_min, altitude_max, N), 2)
|
||||
|
||||
# 原始扰动温度
|
||||
x = ktemp_wn5[cycle_no, :]
|
||||
# 傅里叶变换频谱分析
|
||||
y = ktemp_fft[cycle_no, :]
|
||||
# 滤波后的傅里叶变换频谱分析
|
||||
yy = ktemp_fft_lvbo[cycle_no, :]
|
||||
# 傅里叶逆变换后的扰动温度
|
||||
yyy = ktemp_ifft[cycle_no, :]
|
||||
|
||||
if mode not in ["原始信号", "fft", "ifft", "滤波", "比较"]:
|
||||
raise ValueError(
|
||||
"mode 参数应为 '原始信号', 'fft', 'ifft', 'lv bo' 中的一个")
|
||||
|
||||
plt.figure(figsize=(15, 10)) # 调整图形大小
|
||||
# 原始信号的时间序列
|
||||
if mode == "原始信号":
|
||||
plt.plot(t, x)
|
||||
plt.title('（a）原始信号')
|
||||
plt.xlabel('高度 (km)', labelpad=10) # 增加标签间距
|
||||
plt.ylabel('温度 (°C)', labelpad=10) # 增加标签间距
|
||||
# 原始振幅谱
|
||||
if mode == "fft":
|
||||
plt.plot(f, np.abs(y) * 2 / N)
|
||||
plt.title('（b））原始振幅谱')
|
||||
plt.xlabel('频率/Hz', labelpad=10) # 增加标签间距
|
||||
plt.ylabel('振幅', labelpad=10) # 增加标签间距
|
||||
|
||||
# 通过IFFT回到时间域
|
||||
if mode == "ifft":
|
||||
plt.plot(t, yyy)
|
||||
plt.title('（c））傅里叶逆变换')
|
||||
plt.xlabel('高度 (km)', labelpad=10) # 增加标签间距
|
||||
plt.ylabel('温度 (°C)', labelpad=10) # 增加标签间距
|
||||
|
||||
if mode == "滤波":
|
||||
plt.plot(f, np.abs(yy) * 2 / N)
|
||||
plt.title(f'（d）滤除波长 < {lamda_low} km, > {lamda_high} km的波')
|
||||
plt.xlabel('频率/Hz', labelpad=10) # 增加标签间距
|
||||
plt.ylabel('振幅', labelpad=10) # 增加标签间距
|
||||
|
||||
# 调整子图之间的边距
|
||||
|
||||
# 绘制原始信号与滤波后信号
|
||||
if mode == '比较':
|
||||
plt.plot(t, x, label='原始信号')
|
||||
plt.plot(t, yyy, label='滤波后信号', linestyle='--')
|
||||
plt.title('信号比较')
|
||||
plt.xlabel('高度 (km)', labelpad=10) # 增加标签间距
|
||||
plt.ylabel('温度 (°C)', labelpad=10) # 增加标签间距
|
||||
plt.legend()
|
||||
plt.tight_layout() # 调整子图参数以适应图形区域
|
||||
|
||||
|
||||
day_fft_ifft_modes = ["原始信号", "原始振幅谱", "傅里叶逆变换", "滤波后振幅谱", "信号比较"]
|
||||
|
||||
|
||||
def day_fft_ifft_plotg(cycle_no, ktemp_wn5, ktemp_fft, ktemp_ifft, altitude_min, altitude_max, lamda_low, lamda_high, mode):
|
||||
|
||||
N = len(ktemp_wn5[cycle_no, :])
|
||||
dt = (altitude_max - altitude_min) / (N - 1) # 采样时间间隔
|
||||
n = np.arange(N) # 时间序列索引
|
||||
f = n / (N * dt)
|
||||
t = np.round(np.linspace(altitude_min, altitude_max, N), 2)
|
||||
|
||||
x = ktemp_wn5[cycle_no, :] # 原始扰动温度
|
||||
y = ktemp_fft[cycle_no, :] # 傅里叶变换频谱分析
|
||||
yy = ktemp_fft_lvbo[cycle_no, :] # 滤波后的傅里叶变换频谱分析
|
||||
yyy = ktemp_ifft[cycle_no, :] # 傅里叶逆变换后的扰动温度
|
||||
|
||||
plt.figure(figsize=(15, 10)) # 调整图形大小
|
||||
|
||||
# 计算纵坐标范围
|
||||
temp_limits = (min(min(x), min(yyy)), max(max(x), max(yyy)))
|
||||
amp_limits = (0, max(np.max(np.abs(y) * 2 / N),
|
||||
np.max(np.abs(yy) * 2 / N)))
|
||||
|
||||
if mode not in ["原始信号", "原始振幅谱", "傅里叶逆变换", "滤波后振幅谱", "信号比较"]:
|
||||
raise ValueError(
|
||||
"mode 参数应为 '原始信号', '原始振幅谱', '傅里叶逆变换', '滤波后振幅谱', '信号比较' 中的一个")
|
||||
|
||||
# 原始信号的时间序列
|
||||
if mode == "原始信号":
|
||||
plt.plot(t, x)
|
||||
plt.title('（a）原始信号')
|
||||
plt.xlabel('高度 (km)', labelpad=10)
|
||||
plt.ylabel('温度 (°C)', labelpad=10)
|
||||
plt.ylim(temp_limits) # 设置统一纵坐标范围
|
||||
|
||||
# 原始振幅谱
|
||||
if mode == "原始振幅谱":
|
||||
plt.plot(f, np.abs(y) * 2 / N)
|
||||
plt.title('（b）原始振幅谱')
|
||||
plt.xlabel('频率/Hz', labelpad=10)
|
||||
plt.ylabel('振幅', labelpad=10)
|
||||
plt.ylim(amp_limits) # 设置统一纵坐标范围
|
||||
|
||||
# 通过IFFT回到时间域
|
||||
if mode == "傅里叶逆变换":
|
||||
plt.plot(t, yyy)
|
||||
plt.title('（c）傅里叶逆变换')
|
||||
plt.xlabel('高度 (km)', labelpad=10)
|
||||
plt.ylabel('温度 (°C)', labelpad=10)
|
||||
plt.ylim(temp_limits) # 设置统一纵坐标范围
|
||||
|
||||
# 滤波后的振幅谱
|
||||
if mode == "滤波后振幅谱":
|
||||
plt.plot(f, np.abs(yy) * 2 / N)
|
||||
plt.title(f'（d）滤除波长 < {lamda_low} km, > {lamda_high} km的波')
|
||||
plt.xlabel('频率/Hz', labelpad=10)
|
||||
plt.ylabel('振幅', labelpad=10)
|
||||
plt.ylim(amp_limits) # 设置统一纵坐标范围
|
||||
|
||||
if mode == "信号比较":
|
||||
# 绘制原始信号与滤波后信号
|
||||
plt.figure(figsize=(6, 8)) # 调整图形大小
|
||||
plt.plot(x, t, label='原始信号')
|
||||
plt.plot(yyy, t, label='滤波后信号', linestyle='--')
|
||||
plt.title('信号比较')
|
||||
plt.ylabel('高度 (km)', labelpad=10) # 增加标签间距
|
||||
plt.xlabel('温度 (°C)', labelpad=10) # 增加标签间距
|
||||
plt.legend()
|
||||
plt.xlim(temp_limits) # 设置统一纵坐标范围
|
||||
plt.tight_layout()
|
||||
# -----------------------------------------------------------------------------------------------------------------------------
|
||||
# 6-3 示例的按高度的重力波势能变化曲线图
|
||||
|
||||
|
||||
def day_cycle_power_wave_plot(cycle_no, altitude_min, altitude_max, ktemp_Nz, ktemp_Ptz):
|
||||
|
||||
N = len(ktemp_Nz[cycle_no, :])
|
||||
y = np.round(np.linspace(altitude_min, altitude_max, N), 2)
|
||||
x1 = ktemp_Nz[cycle_no, :]
|
||||
x2 = ktemp_Ptz[cycle_no, :]
|
||||
|
||||
plt.figure(figsize=(12, 10)) # 调整图形大小
|
||||
# 原始信号的时间序列
|
||||
plt.subplot(1, 2, 1)
|
||||
plt.plot(x1, y, label='原始信号')
|
||||
plt.title('（a）Nz')
|
||||
plt.xlabel('Nz', labelpad=10) # 增加标签间距
|
||||
plt.ylabel('高度 (km)', labelpad=10) # 增加标签间距
|
||||
|
||||
# 原始信号的时间序列
|
||||
plt.subplot(1, 2, 2)
|
||||
plt.plot(x2, y, label='原始信号')
|
||||
plt.title('（b）Ptz')
|
||||
plt.xlabel('Ptz', labelpad=10) # 增加标签间距
|
||||
plt.ylabel('高度 (km)', labelpad=10) # 增加标签间距
|
||||
|
||||
# 调整子图之间的边距
|
||||
plt.subplots_adjust(top=0.8, bottom=0.2, left=0.1,
|
||||
right=0.8, hspace=0.3, wspace=0.2)
|
||||
plt.tight_layout() # 调整子图参数以适应图形区域
|
||||
|
||||
plt.show()
|
||||
# -----------------------------------------------------------------------------------------------------------------------------
|
||||
# 6-4 按日统计的按周期计算的不同高度的重力波势能 平面图
|
||||
|
||||
|
||||
def day_power_wave_plot(altitude_min, altitude_max, ktemp_Nz, ktemp_Ptz):
|
||||
# 假设 ktemp_Nz 和 ktemp_Ptz 以及 altitude_min, altitude_max 已经定义好
|
||||
x = np.arange(ktemp_Nz.shape[0])
|
||||
y = np.round(np.linspace(altitude_min, altitude_max, ktemp_Nz.shape[1]), 2)
|
||||
|
||||
# 创建一个图形，并指定两个子图
|
||||
fig, axs = plt.subplots(1, 2, figsize=(15, 10))
|
||||
|
||||
# 第一幅图 (a) NZ
|
||||
cax1 = axs[0].imshow(ktemp_Nz.T, aspect='auto', cmap='viridis',
|
||||
origin='lower', extent=[x[0], x[-1], y[0], y[-1]])
|
||||
fig.colorbar(cax1, ax=axs[0]) # 为第一幅图添加颜色条
|
||||
axs[0].set_title('(a) NZ')
|
||||
axs[0].set_ylabel('Height (km)')
|
||||
axs[0].set_xlabel('Cycles')
|
||||
axs[0].set_yticks(np.linspace(30, 90, 7))
|
||||
axs[0].set_yticklabels(np.round(np.linspace(30, 90, 7), 1))
|
||||
axs[0].set_xticks(np.arange(15))
|
||||
axs[0].set_xticklabels(np.arange(15))
|
||||
|
||||
# 第二幅图 (b) PTZ
|
||||
cax2 = axs[1].imshow(ktemp_Ptz.T, aspect='auto', cmap='viridis',
|
||||
origin='lower', extent=[x[0], x[-1], y[0], y[-1]])
|
||||
fig.colorbar(cax2, ax=axs[1]) # 为第二幅图添加颜色条
|
||||
axs[1].set_title('(b) PTZ')
|
||||
axs[1].set_ylabel('Height (km)')
|
||||
axs[1].set_xlabel('Cycles')
|
||||
axs[1].set_yticks(np.linspace(30, 90, 7))
|
||||
axs[1].set_yticklabels(np.round(np.linspace(30, 90, 7), 1))
|
||||
axs[1].set_xticks(np.arange(15))
|
||||
axs[1].set_xticklabels(np.arange(15))
|
||||
|
||||
# 调整子图之间的边距
|
||||
plt.subplots_adjust(top=0.9, bottom=0.1, left=0.05,
|
||||
right=0.95, hspace=0.3, wspace=0.3)
|
||||
plt.tight_layout() # 调整布局以避免重叠
|
||||
plt.show()
|
||||
# -----------------------------------------------------------------------------------------------------------------------------
|
||||
# 6-5 按月统计的每日重力波势能随天变化的图
|
||||
|
||||
|
||||
# -----------------------------------------------------------------------------------------------------------------------------
|
||||
# 6-6 按年统计的每月重力波势能随月变化的图
|
||||
|
||||
|
||||
# -----------------------------------------------------------------------------------------------------------------------------
|
||||
# 7 主程序，运行数据,并输出主要统计分析图 -
|
||||
# -----------------------------------------------------------------------------------------------------------------------------
|
||||
# 7-1 挑选某一天进行运行主要程序
|
||||
def render_saber(file_path, mode, mode_args=None):
|
||||
# day_read=2018113,
|
||||
day_read = 2018094
|
||||
# 初始化某一天、某个纬度、高度范围等参数
|
||||
latitude_min = 30.0
|
||||
latitude_max = 50.0
|
||||
altitude_min = 30.0
|
||||
altitude_max = 90.0
|
||||
lamda_low = 2
|
||||
lamda_high = 15
|
||||
lvboin = True
|
||||
|
||||
(ktemp_cycles, altitude_cycles,
|
||||
ktemp_wn0,
|
||||
ktemp_fit_wn1, ktemp_wn1,
|
||||
ktemp_fit_wn2, ktemp_wn2,
|
||||
ktemp_fit_wn3, ktemp_wn3,
|
||||
ktemp_fit_wn4, ktemp_wn4,
|
||||
ktemp_fit_wn5, ktemp_wn5,
|
||||
ktemp_fft, ktemp_fft_lvbo, ktemp_ifft,
|
||||
ktemp_Nz, ktemp_Ptz) = day_process_main(
|
||||
file_path,
|
||||
# day_read=2018113,
|
||||
day_read,
|
||||
# 初始化某一天、某个纬度、高度范围等参数
|
||||
latitude_min,
|
||||
latitude_max,
|
||||
altitude_min,
|
||||
altitude_max,
|
||||
lamda_low, lamda_high, lvboin)
|
||||
# -----------------------------------------------------------------------------------------------------------------------------
|
||||
# 7-2 挑选某一个月进行运行主要程序
|
||||
(date_time, ktemp_cycles_mon, altitude_cycles_mon, # 某月份 逐日时间 每日不同循环周期温度 不同循环周期高度数据
|
||||
# 距平扰动温度
|
||||
ktemp_wn0_mon,
|
||||
# 波数为1的拟合温度 滤波后的扰动温度
|
||||
ktemp_fit_wn1_mon, ktemp_wn1_mon,
|
||||
# 波数为2的拟合温度 滤波后的扰动温度
|
||||
ktemp_fit_wn2_mon, ktemp_wn2_mon,
|
||||
# 波数为3的拟合温度 滤波后的扰动温度
|
||||
ktemp_fit_wn3_mon, ktemp_wn3_mon,
|
||||
# 波数为4的拟合温度 滤波后的扰动温度
|
||||
ktemp_fit_wn4_mon, ktemp_wn4_mon,
|
||||
# 波数为5的拟合温度 滤波后的扰动温度
|
||||
ktemp_fit_wn5_mon, ktemp_wn5_mon,
|
||||
# 滤波0-5后扰动温度傅里叶频谱分析 滤除波长为2km-15km内后频谱分析 滤波后的扰动温度
|
||||
ktemp_fft_mon, ktemp_fft_lvbo_mon, ktemp_ifft_mon,
|
||||
# 滤波后NZ、PTZ重力波势能指标计算
|
||||
ktemp_Nz_mon, ktemp_Ptz_mon) = mon_process_main(
|
||||
file_path,
|
||||
# 初始化某一天、某个纬度、高度范围等参数
|
||||
latitude_min,
|
||||
latitude_max,
|
||||
altitude_min,
|
||||
altitude_max,
|
||||
lamda_low, lamda_high, lvboin)
|
||||
|
||||
# -----------------------------------------------------------------------------------------------------------------------------
|
||||
# 7-3 绘制不同截过图
|
||||
|
||||
if mode == "day_fit_wave_plot":
|
||||
|
||||
height_no = 1
|
||||
# day_fit_wave_plotg(height_no, ktemp_wn0, ktemp_fit_wn1, ktemp_wn1, ktemp_fit_wn2, ktemp_wn2,
|
||||
# ktemp_fit_wn3, ktemp_wn3, ktemp_fit_wn4, ktemp_wn4, ktemp_fit_wn5, ktemp_wn5)
|
||||
day_fit_wave_plot(height_no, ktemp_wn0, ktemp_fit_wn1, ktemp_wn1, ktemp_fit_wn2, ktemp_wn2,
|
||||
ktemp_fit_wn3, ktemp_wn3, ktemp_fit_wn4, ktemp_wn4, ktemp_fit_wn5, ktemp_wn5, mode_args)
|
||||
if mode == "day_fft_ifft_plot":
|
||||
cycle_no = 1
|
||||
day_fft_ifft_plotg(cycle_no, ktemp_wn5, ktemp_fft, ktemp_ifft,
|
||||
altitude_min, altitude_max, lamda_low, lamda_high, mode_args)
|
||||
if mode == "day_cycle_power_wave_plot":
|
||||
|
||||
day_cycle_power_wave_plot(cycle_no, altitude_min,
|
||||
altitude_max, ktemp_Nz, ktemp_Ptz)
|
||||
if mode == "day_power_wave_plot":
|
||||
day_power_wave_plot(altitude_min, altitude_max, ktemp_Nz, ktemp_Ptz)
|
||||
|
||||
buffer = BytesIO()
|
||||
plt.savefig(buffer, format='png')
|
||||
buffer.seek(0)
|
||||
return buffer
|
||||
|
||||
|
||||
ALL_RENDER_MODES = ["day_fit_wave_plot", "day_fft_ifft_plot",
|
||||
"day_cycle_power_wave_plot", "day_power_wave_plot"]
|
||||
Loading…
x
Reference in New Issue
Block a user