import numpy as np import matplotlib.pyplot as plt import matplotlib.dates as mdates from 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 = "./SABER/data\\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 = "./saber/data" # 初始化某一天、某个纬度、高度范围等参数 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)