From f567196561349f9b338a1ed2ece56b3fbbc597a8 Mon Sep 17 00:00:00 2001 From: Dustella Date: Thu, 20 Feb 2025 13:10:13 +0800 Subject: [PATCH] cleanup --- modules/balloon/archive/plot_once.py | 175 ---- modules/balloon/archive/plot_year.py | 388 ------- modules/saber/archive/legacy.py | 1077 -------------------- modules/saber/archive/saber_render.py | 828 --------------- modules/saber/archive/zjy_wave_fit_plot.py | 842 --------------- 5 files changed, 3310 deletions(-) delete mode 100644 modules/balloon/archive/plot_once.py delete mode 100644 modules/balloon/archive/plot_year.py delete mode 100644 modules/saber/archive/legacy.py delete mode 100644 modules/saber/archive/saber_render.py delete mode 100644 modules/saber/archive/zjy_wave_fit_plot.py diff --git a/modules/balloon/archive/plot_once.py b/modules/balloon/archive/plot_once.py deleted file mode 100644 index fc6a6f5..0000000 --- a/modules/balloon/archive/plot_once.py +++ /dev/null @@ -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), - ] diff --git a/modules/balloon/archive/plot_year.py b/modules/balloon/archive/plot_year.py deleted file mode 100644 index d4b7c7b..0000000 --- a/modules/balloon/archive/plot_year.py +++ /dev/null @@ -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 diff --git a/modules/saber/archive/legacy.py b/modules/saber/archive/legacy.py deleted file mode 100644 index 238ca91..0000000 --- a/modules/saber/archive/legacy.py +++ /dev/null @@ -1,1077 +0,0 @@ -import numpy as np -import matplotlib.pyplot as plt -import netCDF4 as nc -import pandas as pd -from scipy.optimize import curve_fit -import matplotlib.dates as mdates -# 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): - if not cycles: # 如果周期列表为空，跳过当前迭代 - return None, None - - ktemp_cycles = [] # 初始化列表存储每个周期的温度 - altitude_cycles = [] # 初始化每个循环周期的高度数据 - for event in cycles: - ktemp_cycles_events = np.array(ktemp[event, :]) # 获取每个周期各个事件的ktemp数据 - 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_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 = 157 # 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) - 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 = "E:\\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 = "E:\\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) diff --git a/modules/saber/archive/saber_render.py b/modules/saber/archive/saber_render.py deleted file mode 100644 index 7694417..0000000 --- a/modules/saber/archive/saber_render.py +++ /dev/null @@ -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) diff --git a/modules/saber/archive/zjy_wave_fit_plot.py b/modules/saber/archive/zjy_wave_fit_plot.py deleted file mode 100644 index 5942ee7..0000000 --- a/modules/saber/archive/zjy_wave_fit_plot.py +++ /dev/null @@ -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"]