zephyr/zephyr-backend
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
zephyr-backend/modules/saber/archive/saber_render.py

829 lines
36 KiB
Python
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

# 推断出这是重力波
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)
Reference in New Issue View Git Blame Copy Permalink