iter: v1.5

backend.py

@@ -1,26 +1,42 @@
 from gevent import pywsgi, monkey
 monkey.patch_all()
 
+import cosmic
 import tidi
 from utils import *
 import saber
 import radar
 import balloon
-from flask import Flask
+from flask import Flask, request
 from flask_cors import CORS
 from typing import get_args
 import sys
 import matplotlib.font_manager as fm
 
 app = Flask(__name__)
-fm.fontManager.addfont("./SimHei.ttf")
+CORS(app)
 
+fm.fontManager.addfont("./SimHei.ttf")
+@app.before_request
+def auth():
+    # check for method
+    # if it is OPTIONS, do not check for auth
+    if request.method == "OPTIONS":
+        return
+    code = request.headers.get("Authorization")
+    print(code)
+    if code != "0101":
+        return "Unauthorized", 401
+
+@app.route("/ping")
+def ping():
+    return "pong"
 
 app.register_blueprint(balloon.balloon_module, url_prefix="/balloon")
 app.register_blueprint(radar.radar_module, url_prefix="/radar")
 app.register_blueprint(saber.saber_module, url_prefix="/saber")
 app.register_blueprint(tidi.tidi_module, url_prefix="/tidi")
-
+app.register_blueprint(cosmic.cosmic_module, url_prefix="/cosmic")
 # allow cors
 CORS(app)
 
@@ -36,4 +52,6 @@ if __name__ == '__main__':
         server.serve_forever()
 
     elif 'debug' in args:
-        app.run("0.0.0.0",debug=True)
+        app.run("0.0.0.0",port=18200,debug=True)
+    else:
+        raise Exception("Invalied")

balloon/__init__.py

@@ -52,6 +52,10 @@ def list_ballon():
 def list_ballon_year_modes():
     return get_all_modes()
 
+@balloon_module.route("/metadata/stations")
+def list_stations():
+    return []
+
 
 @balloon_module.route("/render/year")
 def render_full_year():

cosmic/__init__.py

@@ -0,0 +1,22 @@
+from io import BytesIO
+from flask import Blueprint, request, send_file
+from matplotlib import pyplot as plt
+
+from cosmic.temp_render import temp_render
+
+
+cosmic_module = Blueprint("Cosmic", __name__)
+
+@cosmic_module.route('/metadata')
+def get_meta():
+    return []
+
+@cosmic_module.route('/temp_render')
+def render():
+    T = request.args.get("T", 16)
+    temp_render(T=int(T))
+    buf = BytesIO()
+    
+    plt.savefig(buf, format="png")
+    buf.seek(0)
+    return send_file(buf, mimetype="image/png")

cosmic/temp_render.py

@@ -0,0 +1,118 @@
+#此代码是对数据处理后的txt数据进行行星波参数提取绘图
+
+import pandas as pd
+import numpy as np
+from scipy.optimize import curve_fit
+import matplotlib.pyplot as plt
+
+# 定义拟合函数
+# 解决绘图中中文不能显示的问题
+import matplotlib
+# 设置中文显示和负号正常显示
+matplotlib.rcParams['font.sans-serif'] = ['SimHei']  # 显示中文
+matplotlib.rcParams['axes.unicode_minus'] = False  # 正常显示负号
+# 读取一年的数据文件
+
+
+def temp_render(path:str = "./cosmic/cosmic.txt", T = 16):
+    def u_func(x, *params):
+        a1, b1, a2, b2, a3, b3, a4, b4, a5, b5, a6, b6, a7, b7, a8, b8, a9, b9 = params
+        return (
+                a1 * np.sin((2 * np.pi / T) * t - 4 * x + b1)
+                + a2 * np.sin((2 * np.pi / T) * t - 3 * x + b2)
+                + a3 * np.sin((2 * np.pi / T) * t - 2 * x + b3)
+                + a4 * np.sin((2 * np.pi / T) * t - x + b4)
+                + a5 * np.sin((2 * np.pi / T) * t + b5)
+                + a6 * np.sin((2 * np.pi / T) * t + x + b6)
+                + a7 * np.sin((2 * np.pi / T) * t + 2 * x + b7)
+                + a8 * np.sin((2 * np.pi / T) * t + 3 * x + b8)
+                + a9 * np.sin((2 * np.pi / T) * t + 4 * x + b9)
+        )
+        
+    df = pd.read_csv(path,  sep='\s+')
+    # 删除有 NaN 的行
+    df = df.dropna(subset=['Temperature'])
+    # 设置初始参数
+    # initial_guess = [1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0]  # v0, a1, b1, a2, b2, a3, b3
+    initial_guess = [0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,0.5, 0.5, 0.5, 0.5,0.5, 0.5, 0.5, 0.5]  # 9个 a 和 9个 b 参数
+
+    # 设置参数界限
+    bounds = (
+    [0, -np.inf, 0, -np.inf, 0, -np.inf, 0, -np.inf, 0, -np.inf, 0, -np.inf, 0, -np.inf, 0, -np.inf, 0, -np.inf],  # 下界
+    [np.inf, np.inf, np.inf, np.inf, np.inf, np.inf, np.inf, np.inf, np.inf, np.inf, np.inf, np.inf, np.inf, np.inf, np.inf,
+    np.inf, np.inf, np.inf])  # 上界
+
+
+
+
+    # 用于存储拟合参数结果的列表
+    all_fit_results = []
+
+    # 设置最小数据量的阈值
+    min_data_points = 36
+
+    # 进行多个时间窗口的拟合
+    #T应该取5、10、16
+    if T not in [5, 10, 16]:
+        raise ValueError("T should be 5, 10, or 16")
+    for start_day in range(0, 365-3*T):  # 最后一个窗口为[351, 366]
+        end_day = start_day + 3 * T  # 每个窗口的结束时间为 start_day + 3*T
+
+        # 选择当前窗口的数据
+        df_8 = df[(df['Time'] >= start_day) & (df['Time'] <= end_day)]
+
+        # 检查当前窗口的数据量
+        if len(df_8) < min_data_points:
+            # 输出数据量不足的警告
+            print(f"数据量不足，无法拟合：{start_day} 到 {end_day}，数据点数量：{len(df_8)}")
+            # 将拟合参数设置为 NaN
+            all_fit_results.append([np.nan] * 18)
+            continue  # 跳过当前时间窗口，继续下一个窗口
+
+        # 提取时间、经度、温度数据
+        t = np.array(df_8['Time'])  # 时间
+        x = np.array(df_8['Longitude_Radians'])  # 经度弧度制
+        temperature = np.array(df_8['Temperature'])  # 温度，因变量
+
+        # 用T进行拟合
+        popt, pcov = curve_fit(u_func, x, temperature, p0=initial_guess, bounds=bounds, maxfev=50000)
+
+        # 将拟合结果添加到列表中
+        all_fit_results.append(popt)
+
+    # 将结果转换为DataFrame
+    columns = ['a1', 'b1', 'a2', 'b2', 'a3', 'b3', 'a4', 'b4', 'a5', 'b5', 'a6', 'b6', 'a7', 'b7', 'a8', 'b8', 'a9', 'b9']
+    fit_df = pd.DataFrame(all_fit_results, columns=columns)  # fit_df即为拟合的参数汇总
+
+    # -------------------------------画图----------------------------
+    #a1-a9,对应波数-4、-3、-2、-1、0、1、2、3、4的行星波振幅
+    a_columns = ['a1', 'a2', 'a3', 'a4', 'a5', 'a6', 'a7', 'a8', 'a9']
+    k_values = list(range(-4, 5))  # 从 -4 到 4
+
+    # 创建一个字典映射 k 值到 a_columns
+    k_to_a = {f'k={k}': a for k, a in zip(k_values, a_columns)}
+
+    # 获取索引并转换为 numpy 数组
+    x_values = fit_df.index.to_numpy()
+
+    # 对每一列生成独立的图
+    for k, col in k_to_a.items():
+        plt.figure(figsize=(8, 6))  # 创建新的图形
+        plt.plot(x_values, fit_df[col].values)
+        plt.title(f'{k} 振幅图')
+        plt.xlabel('Day')
+        plt.ylabel('振幅')
+
+        # 设置横坐标的动态调整
+        adjusted_x_values = x_values + (3 * T + 1) / 2
+        if len(adjusted_x_values) > 50:
+            step = 30
+            tick_positions = adjusted_x_values[::step]  # 选择每30个点
+            tick_labels = [f'{int(val)}' for val in tick_positions]
+        else:
+            tick_positions = adjusted_x_values
+            tick_labels = [f'{int(val)}' for val in tick_positions]
+
+        plt.xticks(ticks=tick_positions, labels=tick_labels)
+
+        plt.show()  # 显示图形

saber/__init__.py

@@ -25,13 +25,13 @@ def extract_payload():
     return send_file(buffer, mimetype="image/png")
 
 
-all_saber_files = glob.glob("./saber/data/**/**.nc", recursive=True)
+_all_saber_files = glob.glob("./saber/data/**/**.nc", recursive=True)
 
 
 @saber_module.route("/metadata")
 def get_files():
     # normalizing the path, and replace \\ with /
-    all_saber_files = [path.replace("\\", "/") for path in all_saber_files]
+    all_saber_files = [path.replace("\\", "/") for path in _all_saber_files]
     return all_saber_files
 
 

staged/cosmic行星波参数全年逐日绘图.py

@@ -0,0 +1,115 @@
+#此代码是对数据处理后的txt数据进行行星波参数提取绘图
+
+import pandas as pd
+import numpy as np
+from scipy.optimize import curve_fit
+import matplotlib.pyplot as plt
+
+# 解决绘图中中文不能显示的问题
+import matplotlib
+# 设置中文显示和负号正常显示
+matplotlib.rcParams['font.sans-serif'] = ['SimHei']  # 显示中文
+matplotlib.rcParams['axes.unicode_minus'] = False  # 正常显示负号
+# 读取一年的数据文件
+
+df = pd.read_csv(r'./cosmic.txt',  sep='\s+')
+# 删除有 NaN 的行
+df = df.dropna(subset=['Temperature'])
+# 设置初始参数
+# initial_guess = [1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0]  # v0, a1, b1, a2, b2, a3, b3
+initial_guess = [0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,0.5, 0.5, 0.5, 0.5,0.5, 0.5, 0.5, 0.5]  # 9个 a 和 9个 b 参数
+
+# 设置参数界限
+bounds = (
+[0, -np.inf, 0, -np.inf, 0, -np.inf, 0, -np.inf, 0, -np.inf, 0, -np.inf, 0, -np.inf, 0, -np.inf, 0, -np.inf],  # 下界
+[np.inf, np.inf, np.inf, np.inf, np.inf, np.inf, np.inf, np.inf, np.inf, np.inf, np.inf, np.inf, np.inf, np.inf, np.inf,
+ np.inf, np.inf, np.inf])  # 上界
+
+
+# 定义拟合函数
+def u_func(x, *params):
+    a1, b1, a2, b2, a3, b3, a4, b4, a5, b5, a6, b6, a7, b7, a8, b8, a9, b9 = params
+    return (
+            a1 * np.sin((2 * np.pi / T) * t - 4 * x + b1)
+            + a2 * np.sin((2 * np.pi / T) * t - 3 * x + b2)
+            + a3 * np.sin((2 * np.pi / T) * t - 2 * x + b3)
+            + a4 * np.sin((2 * np.pi / T) * t - x + b4)
+            + a5 * np.sin((2 * np.pi / T) * t + b5)
+            + a6 * np.sin((2 * np.pi / T) * t + x + b6)
+            + a7 * np.sin((2 * np.pi / T) * t + 2 * x + b7)
+            + a8 * np.sin((2 * np.pi / T) * t + 3 * x + b8)
+            + a9 * np.sin((2 * np.pi / T) * t + 4 * x + b9)
+    )
+
+
+# 用于存储拟合参数结果的列表
+all_fit_results = []
+
+# 设置最小数据量的阈值
+min_data_points = 36
+
+# 进行多个时间窗口的拟合
+#T应该取5、10、16
+
+T = 16 # 设置 T，可以动态调整
+for start_day in range(0, 365-3*T):  # 最后一个窗口为[351, 366]
+    end_day = start_day + 3 * T  # 每个窗口的结束时间为 start_day + 3*T
+
+    # 选择当前窗口的数据
+    df_8 = df[(df['Time'] >= start_day) & (df['Time'] <= end_day)]
+
+    # 检查当前窗口的数据量
+    if len(df_8) < min_data_points:
+        # 输出数据量不足的警告
+        print(f"数据量不足，无法拟合：{start_day} 到 {end_day}，数据点数量：{len(df_8)}")
+        # 将拟合参数设置为 NaN
+        all_fit_results.append([np.nan] * 18)
+        continue  # 跳过当前时间窗口，继续下一个窗口
+
+    # 提取时间、经度、温度数据
+    t = np.array(df_8['Time'])  # 时间
+    x = np.array(df_8['Longitude_Radians'])  # 经度弧度制
+    temperature = np.array(df_8['Temperature'])  # 温度，因变量
+
+    # 用T进行拟合
+    popt, pcov = curve_fit(u_func, x, temperature, p0=initial_guess, bounds=bounds, maxfev=50000)
+
+    # 将拟合结果添加到列表中
+    all_fit_results.append(popt)
+
+# 将结果转换为DataFrame
+columns = ['a1', 'b1', 'a2', 'b2', 'a3', 'b3', 'a4', 'b4', 'a5', 'b5', 'a6', 'b6', 'a7', 'b7', 'a8', 'b8', 'a9', 'b9']
+fit_df = pd.DataFrame(all_fit_results, columns=columns)  # fit_df即为拟合的参数汇总
+
+# -------------------------------画图----------------------------
+#a1-a9,对应波数-4、-3、-2、-1、0、1、2、3、4的行星波振幅
+a_columns = ['a1', 'a2', 'a3', 'a4', 'a5', 'a6', 'a7', 'a8', 'a9']
+k_values = list(range(-4, 5))  # 从 -4 到 4
+
+# 创建一个字典映射 k 值到 a_columns
+k_to_a = {f'k={k}': a for k, a in zip(k_values, a_columns)}
+
+# 获取索引并转换为 numpy 数组
+x_values = fit_df.index.to_numpy()
+
+# 对每一列生成独立的图
+for k, col in k_to_a.items():
+    plt.figure(figsize=(8, 6))  # 创建新的图形
+    plt.plot(x_values, fit_df[col].values)
+    plt.title(f'{k} 振幅图')
+    plt.xlabel('Day')
+    plt.ylabel('振幅')
+
+    # 设置横坐标的动态调整
+    adjusted_x_values = x_values + (3 * T + 1) / 2
+    if len(adjusted_x_values) > 50:
+        step = 30
+        tick_positions = adjusted_x_values[::step]  # 选择每30个点
+        tick_labels = [f'{int(val)}' for val in tick_positions]
+    else:
+        tick_positions = adjusted_x_values
+        tick_labels = [f'{int(val)}' for val in tick_positions]
+
+    plt.xticks(ticks=tick_positions, labels=tick_labels)
+
+    plt.show()  # 显示图形

staged/cosmic重力波多天.py

@@ -0,0 +1,547 @@
+import os
+import numpy as np
+import pandas as pd
+from scipy.interpolate import interp1d
+from scipy.optimize import curve_fit
+import netCDF4 as nc
+import matplotlib.pyplot as plt
+import seaborn as sns
+# 设置支持中文的字体
+plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']  # 设置字体为微软雅黑
+plt.rcParams['axes.unicode_minus'] = False            # 正常显示负号
+
+# 定义处理单个文件的函数
+def process_single_file(base_folder_path, i):
+    # 构建当前文件夹的路径
+    if i < 10:
+        folder_name = f"atmPrf_repro2021_2008_00{i}"  # 一位数，前面加两个0
+    elif i < 100:
+        folder_name = f"atmPrf_repro2021_2008_0{i}"  # 两位数，前面加一个0
+    else:
+        folder_name = f"atmPrf_repro2021_2008_{i}"   # 三位数，不加0
+    folder_path = os.path.join(base_folder_path, folder_name)
+    # 检查文件夹是否存在
+    if os.path.exists(folder_path):
+        dfs = []
+        # 遍历文件夹中的文件
+        for file_name in os.listdir(folder_path):
+            if file_name.endswith('.0390_nc'):
+                finfo = os.path.join(folder_path, file_name)
+                print(f"正在处理文件: {finfo}")
+                try:
+                    dataset = nc.Dataset(finfo, 'r')
+                    # 提取变量数据
+                    temp = dataset.variables['Temp'][:]
+                    altitude = dataset.variables['MSL_alt'][:]
+                    lat = dataset.variables['Lat'][:]
+                    lon = dataset.variables['Lon'][:]
+                    # 创建DataFrame
+                    df = pd.DataFrame({
+                        'Longitude': lon,
+                        'Latitude': lat,
+                        'Altitude': altitude,
+                        'Temperature': temp
+                    })
+                    dataset.close()
+                    # 剔除高度大于60的行
+                    df = df[df['Altitude'] <= 60]
+                    # 对每个文件的数据进行插值
+                    alt_interp = np.linspace(df['Altitude'].min(), df['Altitude'].max(), 3000)
+                    f_alt = interp1d(df['Altitude'], df['Altitude'], kind='linear', fill_value="extrapolate")
+                    f_lon = interp1d(df['Altitude'], df['Longitude'], kind='linear', fill_value="extrapolate")
+                    f_lat = interp1d(df['Altitude'], df['Latitude'], kind='linear', fill_value="extrapolate")
+                    f_temp = interp1d(df['Altitude'], df['Temperature'], kind='linear', fill_value="extrapolate")
+                    # 计算插值结果
+                    interpolated_alt = f_alt(alt_interp)
+                    interpolated_lon = f_lon(alt_interp)
+                    interpolated_lat = f_lat(alt_interp)
+                    interpolated_temp = f_temp(alt_interp)
+                    # 创建插值后的DataFrame
+                    interpolated_df = pd.DataFrame({
+                        'Altitude': interpolated_alt,
+                        'Longitude': interpolated_lon,
+                        'Latitude': interpolated_lat,
+                        'Temperature': interpolated_temp
+                    })
+                    # 将插值后的DataFrame添加到列表中
+                    dfs.append(interpolated_df)
+                except Exception as e:
+                    print(f"处理文件 {finfo} 时出错: {e}")
+        # 按行拼接所有插值后的DataFrame
+        final_df = pd.concat(dfs, axis=0, ignore_index=True)
+        # 获取 DataFrame 的长度
+        num_rows = len(final_df)
+        # 生成一个每3000个数从0到2999的序列并重复
+        altitude_values = np.tile(np.arange(3000), num_rows // 3000 + 1)[:num_rows]
+        # 将生成的值赋给 DataFrame 的 'Altitude' 列
+        final_df['Altitude'] = altitude_values
+        # 摄氏度换算开尔文
+        final_df['Temperature'] = final_df['Temperature'] + 273.15
+
+        # 筛选出纬度在30到40度之间的数据
+        latitude_filtered_df = final_df[(final_df['Latitude'] >= 30) & (final_df['Latitude'] <= 40)]
+
+        # 划分经度网格，20°的网格
+        lon_min, lon_max = latitude_filtered_df['Longitude'].min(), latitude_filtered_df['Longitude'].max()
+        lon_bins = np.arange(lon_min, lon_max + 20, 20)  # 创建经度网格边界
+        # 将数据分配到网格中
+        latitude_filtered_df['Longitude_Grid'] = np.digitize(latitude_filtered_df['Longitude'], lon_bins) - 1
+
+        # 对相同高度的温度取均值，忽略NaN
+        altitude_temperature_mean = latitude_filtered_df.groupby('Altitude')['Temperature'].mean().reset_index()
+        # 重命名列，使其更具可读性
+        altitude_temperature_mean.columns = ['Altitude', 'Mean_Temperature']
+
+        # 定义高度的范围（这里从0到最短段）
+        altitude_range = range(0, 3000)
+        all_heights_mean_temperature = []  # 用于存储所有高度下的温度均值结果
+        for altitude in altitude_range:
+            # 筛选出当前高度的所有数据
+            altitude_df = latitude_filtered_df[latitude_filtered_df['Altitude'] == altitude]
+            # 对Longitude_Grid同一区间的温度取均值
+            temperature_mean_by_grid = altitude_df.groupby('Longitude_Grid')['Temperature'].mean().reset_index()
+            # 重命名列，使其更具可读性
+            temperature_mean_by_grid.columns = ['Longitude_Grid', 'Mean_Temperature']
+            # 添加高度信息列，方便后续区分不同高度的结果
+            temperature_mean_by_grid['Altitude'] = altitude
+            # 将当前高度的结果添加到列表中
+            all_heights_mean_temperature.append(temperature_mean_by_grid)
+        # 将所有高度的结果合并为一个DataFrame
+        combined_mean_temperature_df = pd.concat(all_heights_mean_temperature, ignore_index=True)
+
+        # 基于Altitude列合并两个DataFrame，只保留能匹配上的行
+        merged_df = pd.merge(combined_mean_temperature_df, altitude_temperature_mean, on='Altitude', how='inner')
+        # 计算差值（减去wn0的扰动）
+        merged_df['Temperature_Difference'] = merged_df['Mean_Temperature_x'] - merged_df['Mean_Temperature_y']
+
+        # 按Altitude分组
+        grouped = merged_df.groupby('Altitude')
+
+        def single_harmonic(x, A, phi):
+            return A * np.sin(2 * np.pi / (18 / k) * x + phi)
+        # 初始化存储每个高度的最佳拟合参数、拟合曲线、残差值以及背景温度的字典
+        fit_results = {}
+        fitted_curves = {}
+        residuals = {}
+        background_temperatures = {}
+        for altitude, group in grouped:
+            y_data = group['Temperature_Difference'].values
+            x_data = np.arange(len(y_data))
+            wn0_data = group['Mean_Temperature_y'].values  # 获取同一高度下的wn0数据
+            # 检查Temperature_Difference列是否全部为NaN
+            if np.all(np.isnan(y_data)):
+                fit_results[altitude] = {'A': [np.nan] * 5, 'phi': [np.nan] * 5}
+                fitted_curves[altitude] = [np.nan * x_data] * 5
+                residuals[altitude] = np.nan * x_data
+                background_temperatures[altitude] = np.nan * x_data
+            else:
+                # 替换NaN值为非NaN值的均值
+                y_data = np.where(np.isnan(y_data), np.nanmean(y_data), y_data)
+                # 初始化存储WN参数和曲线的列表
+                wn_params = []
+                wn_curves = []
+                # 计算wn0（使用Mean_Temperature_y列数据）
+                wn0 = wn0_data
+
+                # 对WN1至WN5进行拟合
+                for k in range(1, 6):
+                    # 更新单谐波函数中的k值
+                    harmonic_func = lambda x, A, phi: single_harmonic(x, A, phi)
+                    # 使用curve_fit进行拟合
+                    popt, pcov = curve_fit(harmonic_func, x_data, y_data, p0=[np.nanmax(y_data) - np.nanmin(y_data), 0])
+                    A_fit, phi_fit = popt
+                    # 存储拟合结果
+                    wn_params.append({'A': A_fit, 'phi': phi_fit})
+                    # 使用拟合参数生成拟合曲线
+                    WN = harmonic_func(x_data, A_fit, phi_fit)
+                    wn_curves.append(WN)
+                    # 计算残差值
+                    y_data = y_data - WN  # 使用残差值作为下一次拟合的y_data
+                # 存储结果
+                fit_results[altitude] = wn_params
+                fitted_curves[altitude] = wn_curves
+                residuals[altitude] = y_data
+                # 计算同一高度下的背景温度（wn0 + wn1 + wn2 + wn3 + wn4 + wn5）
+                wn_sum = np.sum([wn0] + wn_curves, axis=0)
+                background_temperatures[altitude] = wn_sum
+
+        # 将每个字典转换成一个 DataFrame
+        df = pd.DataFrame(residuals)
+        # 使用前向填充（用上一个有效值填充 NaN）
+        df.ffill(axis=1, inplace=True)
+        # 初始化一个新的字典来保存处理结果
+        result = {}
+        # 定义滤波范围
+        lambda_low = 2  # 2 km
+        lambda_high = 15  # 15 km
+        f_low = 2 * np.pi / lambda_high
+        f_high = 2 * np.pi / lambda_low
+        # 循环处理df的每一行（每个高度）
+        for idx, residuals_array in df.iterrows():
+            # 提取有效值
+            valid_values = np.ma.masked_array(residuals_array, np.isnan(residuals_array))
+            compressed_values = valid_values.compressed()  # 去除NaN值后的数组
+            N = len(compressed_values)  # 有效值的数量
+            # 如果有效值为空（即所有值都是NaN），则将结果设置为NaN
+            if N == 0:
+                result[idx] = np.full_like(residuals_array, np.nan)
+            else:
+                # 时间序列和频率
+                dt = 0.02  # 假设的时间间隔
+                n = np.arange(N)
+                f = n / (N * dt)
+                # 傅里叶变换
+                y = np.fft.fft(compressed_values)
+                # 频率滤波
+                yy = y.copy()
+                freq_filter = (f >= f_low) & (f <= f_high)
+                yy[~freq_filter] = 0
+                # 逆傅里叶变换
+                perturbation_after = np.real(np.fft.ifft(yy))
+                # 将处理结果插回到result字典中
+                result[idx] = perturbation_after
+
+        # 处理背景温度和扰动温度数据格式
+        heights = list(background_temperatures.keys())
+        data_length = len(next(iter(background_temperatures.values())))
+        background_matrix = np.zeros((data_length, len(heights)))
+        for idx, height in enumerate(heights):
+            background_matrix[:, idx] = background_temperatures[height]
+
+        heights = list(result.keys())
+        data_length = len(next(iter(result.values())))
+        perturbation_matrix = np.zeros((data_length, len(heights)))
+        for idx, height in enumerate(heights):
+            perturbation_matrix[:, idx] = result[height]
+        perturbation_matrix = perturbation_matrix.T
+
+        # 计算 Brunt-Väisälä 频率和势能
+        heights_for_calc = np.linspace(0, 60, 3000) * 1000
+
+        def brunt_vaisala_frequency(g, BT_z, c_p, heights):
+            # 计算位温随高度的变化率
+            dBT_z_dz = np.gradient(BT_z, heights)
+            # 计算 Brunt-Väisälä 频率，根号内取绝对值
+            frequency_squared = (g / BT_z) * ((g / c_p) + dBT_z_dz)
+            frequency = np.sqrt(np.abs(frequency_squared))
+            return frequency
+
+        # 计算势能
+        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  # 比热容
+        N_z_matrix = []
+        PT_z_matrix = []
+        for i in range(background_matrix.shape[0]):
+            BT_z = np.array(background_matrix[i])
+            PT_z = np.array(perturbation_matrix[i])
+            N_z = brunt_vaisala_frequency(g, BT_z, c_p, heights_for_calc)
+            PW = calculate_gravitational_potential_energy(g, BT_z, N_z, PT_z)
+            N_z_matrix.append(N_z)
+            PT_z_matrix.append(PW)
+        ktemp_Nz = np.vstack(N_z_matrix)
+        ktemp_Ptz = np.vstack(PT_z_matrix)
+        mean_ktemp_Nz = np.mean(ktemp_Nz, axis=0)
+        mean_ktemp_Ptz = np.mean(ktemp_Ptz, axis=0)
+
+        return mean_ktemp_Nz, mean_ktemp_Ptz
+    else:
+        print(f"文件夹 {folder_path} 不存在。")
+        return None, None
+
+# 主循环，处理1到365个文件
+base_folder_path = r"E:\COSMIC\2008"
+all_mean_ktemp_Nz = []
+all_mean_ktemp_Ptz = []
+for file_index in range(101, 104):
+    try:
+        mean_ktemp_Nz, mean_ktemp_Ptz = process_single_file(base_folder_path, file_index)
+        if mean_ktemp_Nz is not None and mean_ktemp_Ptz is not None:
+            all_mean_ktemp_Nz.append(mean_ktemp_Nz)
+            all_mean_ktemp_Ptz.append(mean_ktemp_Ptz)
+    except ValueError as e:
+        print(f"Error processing file index {file_index}: {e}, skipping this file.")
+        continue
+
+# 转换每个数组为二维形状
+final_mean_ktemp_Nz = np.vstack([arr.reshape(1, -1) for arr in all_mean_ktemp_Nz])
+final_mean_ktemp_Ptz = np.vstack([arr.reshape(1, -1) for arr in all_mean_ktemp_Ptz])
+# 使用条件索引替换大于50的值为NaN
+final_mean_ktemp_Ptz[final_mean_ktemp_Ptz > 50] = np.nan
+# heights 为每个高度的值
+heights = np.linspace(0, 60, 3000)
+df_final_mean_ktemp_Ptz = pd.DataFrame(final_mean_ktemp_Ptz)
+df_final_mean_ktemp_Nz = pd.DataFrame(final_mean_ktemp_Nz)
+#-------------------------------------------------绘制年统计图------------------------------------
+#-----------绘制浮力频率年统计图-----------------------
+data=df_final_mean_ktemp_Nz.T
+# 对每个元素进行平方（计算N2）
+data= data ** 2
+data=data*10000#(绘图好看)
+#将大于 10 的值替换为 NaN（个别异常值）
+data[data > 10] = np.nan
+# 绘制热力图的函数
+def plot_heatmap(data, heights, title):
+    plt.figure(figsize=(10, 8))
+    # 绘制热力图，数据中的行代表高度，列代表天数
+    sns.heatmap(data, cmap='coolwarm', xticklabels=1, yticklabels=50, cbar_kws={'label': 'Value'})
+    plt.xlabel('Day')
+    plt.ylabel('Height (km)')
+    plt.title(title)
+    # 设置x轴的刻度，使其每30天显示一个标签
+    num_days = data.shape[1]
+    x_tick_positions = np.arange(0, num_days, 30)
+    x_tick_labels = np.arange(0, num_days, 30)
+    plt.xticks(x_tick_positions, x_tick_labels)
+    # 设置y轴的刻度，使其显示对应的高度
+    plt.yticks(np.linspace(0, data.shape[0] - 1, 6), np.round(np.linspace(heights[0], heights[-1], 6), 2))
+    # 反转 y 轴，使 0 在底部
+    plt.gca().invert_yaxis()
+    plt.show()
+# 调用函数绘制热力图
+plot_heatmap(data, heights, 'Heatmap of final_mean_ktemp_Nz（10^(-4)）')
+#-----------------------------------------------------------------------------
+#-------------绘制重力势能年统计图------------------------------------------------
+data1=df_final_mean_ktemp_Ptz.T
+# 绘制热力图的函数
+def plot_heatmap(data, heights, title):
+    plt.figure(figsize=(10, 8))
+    # 绘制热力图，数据中的行代表高度，列代表天数
+    sns.heatmap(data, cmap='coolwarm', xticklabels=1, yticklabels=50, cbar_kws={'label': 'Value'})
+    plt.xlabel('Day')
+    plt.ylabel('Height (km)')
+    plt.title(title)
+    # 设置x轴的刻度，使其每30天显示一个标签
+    num_days = data.shape[1]
+    x_tick_positions = np.arange(0, num_days, 30)
+    x_tick_labels = np.arange(0, num_days, 30)
+    plt.xticks(x_tick_positions, x_tick_labels)
+    # 设置y轴的刻度，使其显示对应的高度
+    plt.yticks(np.linspace(0, data.shape[0] - 1, 6), np.round(np.linspace(heights[0], heights[-1], 6), 2))
+    # 反转 y 轴，使 0 在底部
+    plt.gca().invert_yaxis()
+    plt.show()
+# 调用函数绘制热力图
+plot_heatmap(data1, heights, 'Heatmap of final_mean_ktemp_Ptz(J/kg)')
+#------------------------绘制月统计图---------------------------------------------------------------------------------
+#----------绘制浮力频率月统计图-------------------------------------------------
+# 获取总列数
+num_columns = data.shape[1]
+# 按30列分组计算均值
+averaged_df = []
+# 逐步处理每30列
+for i in range(0, num_columns, 30):
+    # 获取当前范围内的列，并计算均值
+    subset = data.iloc[:, i:i+30]  # 获取第i到i+29列
+    mean_values = subset.mean(axis=1)  # 对每行计算均值
+    averaged_df.append(mean_values)  # 将均值添加到列表
+# 将结果转化为一个新的 DataFrame
+averaged_df = pd.DataFrame(averaged_df).T
+# 1. 每3000行取一个均值
+# 获取总行数
+num_rows = averaged_df.shape[0]
+# 创建一个新的列表来存储每3000行的均值
+averaged_by_rows_df = []
+# 逐步处理每3000行
+for i in range(0, num_rows, 3000):
+    # 获取当前范围内的行
+    subset = averaged_df.iloc[i:i+3000, :]  # 获取第i到i+99行
+    mean_values = subset.mean(axis=0)  # 对每列计算均值
+    averaged_by_rows_df.append(mean_values)  # 将均值添加到列表
+# 将结果转化为一个新的 DataFrame
+averaged_by_rows_df = pd.DataFrame(averaged_by_rows_df)
+# 绘制折线图
+plt.figure(figsize=(10, 6))  # 设置图形的大小
+plt.plot(averaged_by_rows_df.columns, averaged_by_rows_df.mean(axis=0), marker='o', color='b', label='平均值')
+# 添加标题和标签
+plt.title('每月平均N^2的折线图')
+plt.xlabel('月份')
+plt.ylabel('N^2(10^-4)')
+plt.legend()
+# 显示图形
+plt.grid(True)
+plt.xticks(rotation=45)  # 让x轴标签（月份）倾斜，以便更清晰显示
+plt.tight_layout()
+plt.show()
+#------------重力势能的月统计-----------------------------------
+# 获取总列数
+num_columns = data1.shape[1]
+# 按30列分组计算均值
+averaged_df = []
+# 逐步处理每30列
+for i in range(0, num_columns, 30):
+    # 获取当前范围内的列，并计算均值
+    subset = data1.iloc[:, i:i+30]  # 获取第i到i+29列
+    mean_values = subset.mean(axis=1)  # 对每行计算均值
+    averaged_df.append(mean_values)  # 将均值添加到列表
+# 将结果转化为一个新的 DataFrame
+averaged_df = pd.DataFrame(averaged_df).T
+# 1. 每3000行取一个均值
+# 获取总行数
+num_rows = averaged_df.shape[0]
+# 创建一个新的列表来存储每3000行的均值
+averaged_by_rows_df = []
+# 逐步处理每3000行
+for i in range(0, num_rows, 3000):
+    # 获取当前范围内的行
+    subset = averaged_df.iloc[i:i+3000, :]  # 获取第i到i+99行
+    mean_values = subset.mean(axis=0)  # 对每列计算均值
+    averaged_by_rows_df.append(mean_values)  # 将均值添加到列表
+# 将结果转化为一个新的 DataFrame
+averaged_by_rows_df = pd.DataFrame(averaged_by_rows_df)
+# 绘制折线图
+plt.figure(figsize=(10, 6))  # 设置图形的大小
+plt.plot(averaged_by_rows_df.columns, averaged_by_rows_df.mean(axis=0), marker='o', color='b', label='平均值')
+# 添加标题和标签
+plt.title('每月平均重力势能的折线图')
+plt.xlabel('月份')
+plt.ylabel('重力势能（J/Kg）')
+plt.legend()
+# 显示图形
+plt.grid(True)
+plt.xticks(rotation=45)  # 让x轴标签（月份）倾斜，以便更清晰显示
+plt.tight_layout()
+plt.show()
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+# 获取总列数
+total_columns = data.shape[1]
+# 用于存储每一组30列计算得到的均值列数据（最终会构成新的DataFrame）
+mean_columns = []
+# 分组序号，用于生成列名时区分不同的均值列，从1开始
+group_index = 1
+# 按照每30列一组进行划分（不滑动）
+for start_col in range(0, total_columns, 30):
+    end_col = start_col + 30
+    if end_col > total_columns:
+        end_col = total_columns
+    # 选取当前组的30列（如果不足30列，按实际剩余列数选取）
+    group_data = data.iloc[:, start_col:end_col]
+    # 按行对当前组的列数据求和
+    sum_per_row = group_data.sum(axis=1)
+    # 计算平均（每一组的平均，每行都有一个平均结果）
+    mean_per_row = sum_per_row / (end_col - start_col)
+    # 生成新的列名，格式为'Mean_分组序号'，例如'Mean_1'、'Mean_2'等
+    new_column_name = f'Mean_{group_index}'
+    group_index += 1
+    # 将当前组计算得到的均值列添加到列表中
+    mean_columns.append(mean_per_row)
+# 将所有的均值列合并为一个新的DataFrame（列方向合并）
+new_mean_df = pd.concat(mean_columns, axis=1)
+# 按行对new_mean_df所有列的数据进行求和，得到一个Series，索引与new_mean_df的索引一致，每个元素是每行的总和
+row_sums = new_mean_df.sum(axis=1)
+# 计算所有行总和的均值
+mean_value = row_sums.mean()
+# 设置中文字体为黑体，解决中文显示问题（Windows系统下），如果是其他系统或者有其他字体需求可适当调整
+plt.rcParams['font.sans-serif'] = ['SimHei']
+# 解决负号显示问题
+plt.rcParams['axes.unicode_minus'] = False
+# 提取月份作为x轴标签（假设mean_value的索引就是月份信息）
+months = mean_value.index.tolist()
+# 提取均值数据作为y轴数据
+energy_values = mean_value.tolist()
+# 创建图形和坐标轴对象
+fig, ax = plt.subplots(figsize=(10, 6))
+# 绘制折线图
+ax.plot(months, energy_values, marker='o', linestyle='-', color='b')
+# 设置坐标轴标签和标题
+ax.set_xlabel('月份')
+ax.set_ylabel('平均浮力频率')
+ax.set_title('每月浮力频率变化趋势')
+# 设置x轴刻度，让其旋转一定角度以便更好地显示所有月份标签，避免重叠
+plt.xticks(rotation=45)
+# 显示网格线，增强图表可读性
+ax.grid(True)
+# 显示图形
+plt.show()
+#--------------------------------绘制重力势能月统计图------------------------------
+# 获取总列数
+total_columns = data1.shape[1]
+# 用于存储每一组30列计算得到的均值列数据（最终会构成新的DataFrame）
+mean_columns = []
+# 分组序号，用于生成列名时区分不同的均值列，从1开始
+group_index = 1
+# 按照每30列一组进行划分（不滑动）
+for start_col in range(0, total_columns, 30):
+    end_col = start_col + 30
+    if end_col > total_columns:
+        end_col = total_columns
+    # 选取当前组的30列（如果不足30列，按实际剩余列数选取）
+    group_data = data1.iloc[:, start_col:end_col]
+    # 按行对当前组的列数据求和
+    sum_per_row = group_data.sum(axis=1)
+    # 计算平均（每一组的平均，每行都有一个平均结果）
+    mean_per_row = sum_per_row / (end_col - start_col)
+    # 生成新的列名，格式为'Mean_分组序号'，例如'Mean_1'、'Mean_2'等
+    new_column_name = f'Mean_{group_index}'
+    group_index += 1
+    # 将当前组计算得到的均值列添加到列表中
+    mean_columns.append(mean_per_row)
+# 将所有的均值列合并为一个新的DataFrame（列方向合并）
+new_mean_df = pd.concat(mean_columns, axis=1)
+# 按行对new_mean_df所有列的数据进行求和，得到一个Series，索引与new_mean_df的索引一致，每个元素是每行的总和
+row_sums = new_mean_df.sum(axis=1)
+# 计算所有行总和的均值
+mean_value = row_sums.mean()
+# 设置中文字体为黑体，解决中文显示问题（Windows系统下），如果是其他系统或者有其他字体需求可适当调整
+plt.rcParams['font.sans-serif'] = ['SimHei']
+# 解决负号显示问题
+plt.rcParams['axes.unicode_minus'] = False
+# 提取月份作为x轴标签（假设mean_value的索引就是月份信息）
+months = mean_value.index.tolist()
+# 提取均值数据作为y轴数据
+energy_values = mean_value.tolist()
+# 创建图形和坐标轴对象
+fig, ax = plt.subplots(figsize=(10, 6))
+# 绘制折线图
+ax.plot(months, energy_values, marker='o', linestyle='-', color='b')
+# 设置坐标轴标签和标题
+ax.set_xlabel('月份')
+ax.set_ylabel('平均浮力频率')
+ax.set_title('每月浮力频率变化趋势')
+# 设置x轴刻度，让其旋转一定角度以便更好地显示所有月份标签，避免重叠
+plt.xticks(rotation=45)
+# 显示网格线，增强图表可读性
+ax.grid(True)
+# 显示图形
+plt.show()

tidi/__init__.py

@@ -3,6 +3,7 @@ from io import BytesIO
 from flask import Blueprint, request, send_file
 from matplotlib import pyplot as plt
 
+from tidi.plot import TidiPlotv2
 from tidi.staged.plot import tidi_render
 
 
@@ -32,4 +33,28 @@ def render_wave():
     plt.savefig(buffer, format="png")
     buffer.seek(0)
     
-    return send_file(buffer, mimetype="image/png")    
+    return send_file(buffer, mimetype="image/png")    
+
+@tidi_module.route('/render/month_stats_v1')
+def render_stats_v1():
+    year = request.args.get('year')
+    year = int(year)
+    
+    plotter = TidiPlotv2(year)
+    plotter.plot_v1()
+    buffer = BytesIO()
+    plt.savefig(buffer, format="png")
+    buffer.seek(0)
+    return send_file(buffer, mimetype="image/png")
+
+@tidi_module.route('/render/month_stats_v2')
+def render_stats_v2():
+    year = request.args.get('year')
+    year = int(year)
+    
+    plotter = TidiPlotv2(year)
+    plotter.plot_month()
+    buffer = BytesIO()
+    plt.savefig(buffer, format="png")
+    buffer.seek(0)
+    return send_file(buffer, mimetype="image/png")

tidi/staged/plot_old.py

@@ -0,0 +1,878 @@
+import pandas as pd
+import numpy as np
+from scipy.io import loadmat
+from scipy.optimize import curve_fit
+import matplotlib.pyplot as plt
+import seaborn as sns
+# ---------------------------------------------------------------------------------------
+# -----vzonal----------------------------------------------------------------------------
+
+
+def process_vzonal_day(day, year, ):
+    try:
+        # 读取数据
+        height_data = loadmat(rf"./tidi/data\{year}\{day:03d}_Height.mat")
+        lat_data = loadmat(rf"./tidi/data\{year}\{day:03d}_Lat.mat")
+        lon_data = loadmat(rf"./tidi/data\{year}\{day:03d}_Lon.mat")
+        vmeridional_data = loadmat(
+            rf"./tidi/data\{year}\{day:03d}_VMerdional.mat")
+        vzonal_data = loadmat(rf"./tidi/data\{year}\{day:03d}_Vzonal.mat")
+
+        # 将数据转换为DataFrame
+        height_df = pd.DataFrame(height_data['Height'])
+        lat_df = pd.DataFrame(lat_data['Lat'])
+        lon_df = pd.DataFrame(lon_data['Lon'])
+        vmeridional_df = pd.DataFrame(vmeridional_data['VMerdional'])
+        vzonal_df = pd.DataFrame(vzonal_data['Vzonal'])
+        # 将经纬度拼接为两列并添加到对应的DataFrame中
+        lon_lat_df = pd.concat([lon_df, lat_df], axis=1)
+        lon_lat_df.columns = ['Longitude', 'Latitude']
+        # 筛选出10到30度纬度范围的数据
+        lat_filter = (lat_df.values >= 0) & (lat_df.values <= 20)
+        # 使用纬度范围过滤数据
+        vmeridional_filtered = vmeridional_df.iloc[:, lat_filter.flatten()]
+        vzonal_filtered = vzonal_df.iloc[:, lat_filter.flatten()]
+        lon_lat_filtered = lon_lat_df.iloc[lat_filter.flatten(), :]
+        # 接着对lon_lat_filtered的经度进行分组，0到360度每30度一个区间
+        bins = range(0, 361, 30)
+        group_labels = [f"{i}-{i + 29}" for i in range(0, 360, 30)]
+
+        lon_lat_filtered['Longitude_Group'] = pd.cut(
+            lon_lat_filtered['Longitude'], bins=bins, labels=group_labels)
+
+        # 获取所有唯一的经度分组标签并按照数值顺序排序
+        unique_groups = sorted(lon_lat_filtered['Longitude_Group'].unique(
+        ), key=lambda x: int(x.split('-')[0]))
+        # 按照经度分组获取每个区间对应的vzonal_filtered、vmeridional_filtered数据
+        grouped_data = {}
+        insufficient_data_count = 0  # 用于计数数据不足的组数
+
+        for group in unique_groups:
+            mask = lon_lat_filtered['Longitude_Group'] == group
+            grouped_data[group] = {
+                'vzonal_filtered': vzonal_filtered.loc[:, mask],
+                'vmeridional_filtered': vmeridional_filtered.loc[:, mask],
+                'lon_lat_filtered': lon_lat_filtered.loc[mask]
+            }
+
+            # 计算有效值数量
+            vzonal_count = grouped_data[group]['vzonal_filtered'].notna(
+            ).sum().sum()
+            vmeridional_count = grouped_data[group]['vmeridional_filtered'].notna(
+            ).sum().sum()
+
+            if vzonal_count <= 20 or vmeridional_count <= 20:
+                insufficient_data_count += 1
+
+        # 如果超过6组数据不足，则抛出错误
+        if insufficient_data_count > 6:
+            expected_length = 21
+            return pd.Series(np.nan, index=range(expected_length))
+
+        # 如果代码运行到这里，说明所有分组的数据量都足够或者不足的组数不超过6
+        print("所有分组的数据量都足够")
+        # -----------计算w0------------------------------------------------------------------------------------------
+        # 定义期望的12个区间的分组标签
+        expected_groups = [f"{i}-{i + 29}" for i in range(0, 360, 30)]
+        # 初始化一个空DataFrame来存储所有区间的均值廓线，列名设置为期望的分组标签
+        W0_profiles_df = pd.DataFrame(columns=expected_groups)
+
+        # 遍历grouped_data字典中的每个组
+        for group, data in grouped_data.items():
+            # 提取当前组的vzonal_filtered数据
+            vzonal_filtered = data['vzonal_filtered']
+            # 计算有效数据的均值廓线，跳过NaN值
+            mean_profile = vzonal_filtered.mean(axis=1, skipna=True)
+            # 将当前组的均值廓线作为一列添加到W0_profiles_df DataFrame中
+            W0_profiles_df[group] = mean_profile
+
+        # 检查并填充缺失的区间列，将缺失的列添加并填充为NaN
+        for group in expected_groups:
+            if group not in W0_profiles_df.columns:
+                W0_profiles_df[group] = pd.Series(
+                    [float('NaN')] * len(W0_profiles_df))
+
+        # 打印拼接后的DataFrame以验证
+        print("Concatenated mean profiles for all longitude groups:\n", W0_profiles_df)
+        # 计算每个高度的均值
+        height_mean_profiles = W0_profiles_df.mean(axis=1)
+        # 将每个高度的均值作为新的一行添加到DataFrame中，All_Heights_Mean就是wn0
+        W0_profiles_df['All_Heights_Mean'] = height_mean_profiles
+        wn0_df = W0_profiles_df['All_Heights_Mean']
+        # -------计算残余量--------------------------------------------------------------------------------------
+        # 计算每个经度区间的残余值 (即每个区间的值减去该高度的All_Heights_Mean)
+        residuals_df = W0_profiles_df.drop(columns='All_Heights_Mean').subtract(
+            W0_profiles_df['All_Heights_Mean'], axis=0)
+
+        # --------wn1-------------------------
+        def single_harmonic(x, A, phi, C):
+            return A * np.sin(2 * np.pi / 12 * x + phi) + C
+        # 用于存储每个高度拟合后的参数结果
+        fit_results = []
+        for index, row in residuals_df.iterrows():
+            # 检查该行是否存在NaN值，如果有则跳过拟合，直接将参数设为0
+            if row.isnull().any():
+                fit_results.append((0, 0, 0))
+                continue
+            x = np.arange(12)  # 对应12个位置作为自变量
+            y = row.values
+            try:
+                # 进行曲线拟合
+                popt, _ = curve_fit(single_harmonic, x, y)
+                fit_results.append(popt)
+            except RuntimeError:
+                # 如果拟合过程出现问题（例如无法收敛等），也将参数设为0
+                fit_results.append((0, 0, 0))
+        # 将拟合结果转换为DataFrame
+        fit_results_df = pd.DataFrame(fit_results, columns=['A', 'phi', 'C'])
+        print(fit_results_df)
+        # 用于存储每个高度的拟合值
+        wn1_values = []
+        for index, row in fit_results_df.iterrows():
+            A, phi, C = row
+            x = np.arange(12)  # 同样对应12个位置作为自变量
+            wn1 = single_harmonic(x, A, phi, C)
+            wn1_values.append(wn1)
+        # 将拟合值转换为DataFrame
+        wn1_df = pd.DataFrame(wn1_values, columns=[
+                              f'wn1_{i}' for i in range(12)])
+        print(wn1_df)
+        # 如果wn1_df全为0，则跳过下面的计算，直接令该天的day_log_gwresult全部为NaN
+        if (wn1_df == 0).all().all():
+            return pd.Series(np.nan, index=range(21))
+        # ------------计算temp-wn0-wn1---------------------------------------------------------
+        temp_wn0_wn1 = residuals_df.values - wn1_df.values
+        # 将结果转为 DataFrame
+        temp_wn0_wn1_df = pd.DataFrame(
+            temp_wn0_wn1, columns=residuals_df.columns, index=residuals_df.index)
+
+        # -------wn2--------------------------------------------------------------------------
+        def single_harmonic(x, A, phi, C):
+            return A * np.sin(2 * np.pi / 6 * x + phi) + C
+
+        # 用于存储每个高度拟合后的参数结果
+        fit_results2 = []
+        for index, row in temp_wn0_wn1_df.iterrows():
+            # 检查该行是否存在NaN值，如果有则跳过拟合，直接将参数设为0
+            if row.isnull().any():
+                fit_results2.append((0, 0, 0))
+                continue
+            x = np.arange(12)  # 对应12个位置作为自变量
+            y = row.values
+            try:
+                # 进行曲线拟合
+                popt, _ = curve_fit(single_harmonic, x, y)
+                fit_results2.append(popt)
+            except RuntimeError:
+                # 如果拟合过程出现问题（例如无法收敛等），也将参数设为0
+                fit_results2.append((0, 0, 0))
+        # 将拟合结果转换为DataFrame
+        fit_results2_df = pd.DataFrame(fit_results2, columns=['A', 'phi', 'C'])
+        print(fit_results2_df)
+        # 用于存储每个高度的拟合值
+        wn2_values = []
+        for index, row in fit_results2_df.iterrows():
+            A, phi, C = row
+            x = np.arange(12)  # 同样对应12个位置作为自变量
+            wn2 = single_harmonic(x, A, phi, C)
+            wn2_values.append(wn2)
+        # 将拟合值转换为DataFrame
+        wn2_df = pd.DataFrame(wn2_values, columns=[
+                              f'wn2_{i}' for i in range(12)])
+        print(wn2_df)
+        # ---------计算temp-wn0-wn1-wn2------------------------------------------------------
+        temp_wn0_wn1_wn2 = temp_wn0_wn1_df.values - wn2_df.values
+        # 转换为 DataFrame
+        temp_wn0_wn1_wn2_df = pd.DataFrame(
+            temp_wn0_wn1_wn2, columns=temp_wn0_wn1_df.columns)
+
+        # -------wn3-----------------------------------------------------------------------
+        def single_harmonic(x, A, phi, C):
+            return A * np.sin(2 * np.pi / 4 * x + phi) + C
+
+        # 用于存储每个高度拟合后的参数结果
+        fit_results3 = []
+        for index, row in temp_wn0_wn1_wn2_df.iterrows():
+            # 检查该行是否存在NaN值，如果有则跳过拟合，直接将参数设为0
+            if row.isnull().any():
+                fit_results3.append((0, 0, 0))
+                continue
+            x = np.arange(12)  # 对应12个位置作为自变量
+            y = row.values
+            try:
+                # 进行曲线拟合
+                popt, _ = curve_fit(single_harmonic, x, y)
+                fit_results3.append(popt)
+            except RuntimeError:
+                # 如果拟合过程出现问题（例如无法收敛等），也将参数设为0
+                fit_results3.append((0, 0, 0))
+        # 将拟合结果转换为DataFrame
+        fit_results3_df = pd.DataFrame(fit_results3, columns=['A', 'phi', 'C'])
+        print(fit_results3_df)
+        # 用于存储每个高度的拟合值
+        wn3_values = []
+        for index, row in fit_results3_df.iterrows():
+            A, phi, C = row
+            x = np.arange(12)  # 同样对应12个位置作为自变量
+            wn3 = single_harmonic(x, A, phi, C)
+            wn3_values.append(wn3)
+        # 将拟合值转换为DataFrame
+        wn3_df = pd.DataFrame(wn3_values, columns=[
+                              f'wn3_{i}' for i in range(12)])
+        print(wn3_df)
+        # ---------计算temp-wn0-wn1-wn2-wn3------------------------------------------------------
+        temp_wn0_wn1_wn2_wn3 = temp_wn0_wn1_wn2_df.values - wn3_df.values
+        # 转换为 DataFrame
+        temp_wn0_wn1_wn2_wn3_df = pd.DataFrame(
+            temp_wn0_wn1_wn2_wn3, columns=temp_wn0_wn1_wn2_df.columns)
+        # -------wn4 - ----------------------------------------------------------------------
+
+        def single_harmonic(x, A, phi, C):
+            return A * np.sin(2 * np.pi / 3 * x + phi) + C
+
+        # 用于存储每个高度拟合后的参数结果
+        fit_results4 = []
+        for index, row in temp_wn0_wn1_wn2_wn3_df.iterrows():
+            # 检查该行是否存在NaN值，如果有则跳过拟合，直接将参数设为0
+            if row.isnull().any():
+                fit_results4.append((0, 0, 0))
+                continue
+            x = np.arange(12)  # 对应12个位置作为自变量
+            y = row.values
+            try:
+                # 进行曲线拟合
+                popt, _ = curve_fit(single_harmonic, x, y)
+                fit_results4.append(popt)
+            except RuntimeError:
+                # 如果拟合过程出现问题（例如无法收敛等），也将参数设为0
+                fit_results4.append((0, 0, 0))
+
+        fit_results4_df = pd.DataFrame(fit_results4, columns=['A', 'phi', 'C'])
+        print(fit_results4_df)
+        # 用于存储每个高度的拟合值
+        wn4_values = []
+        for index, row in fit_results4_df.iterrows():
+            A, phi, C = row
+            x = np.arange(12)  # 同样对应12个位置作为自变量
+            wn4 = single_harmonic(x, A, phi, C)
+            wn4_values.append(wn4)
+        # 将拟合值转换为DataFrame
+        wn4_df = pd.DataFrame(wn4_values, columns=[
+                              f'wn4_{i}' for i in range(12)])
+        print(wn4_df)
+        # ---------计算temp-wn0-wn1-wn2-wn3------------------------------------------------------
+        temp_wn0_wn1_wn2_wn3_wn4 = temp_wn0_wn1_wn2_wn3_df.values - wn4_df.values
+        # 转换为 DataFrame
+        temp_wn0_wn1_wn2_wn3_wn4_df = pd.DataFrame(
+            temp_wn0_wn1_wn2_wn3_wn4, columns=temp_wn0_wn1_wn2_wn3_df.columns)
+
+        # -------wn5-----------------------------------------------------------------------
+        def single_harmonic(x, A, phi, C):
+            return A * np.sin(2 * np.pi / 2.4 * x + phi) + C
+
+        # 用于存储每个高度拟合后的参数结果
+        fit_results5 = []
+        for index, row in temp_wn0_wn1_wn2_wn3_wn4_df.iterrows():
+            # 检查该行是否存在NaN值，如果有则跳过拟合，直接将参数设为0
+            if row.isnull().any():
+                fit_results5.append((0, 0, 0))
+                continue
+            x = np.arange(12)  # 对应12个位置作为自变量
+            y = row.values
+            try:
+                # 进行曲线拟合
+                popt, _ = curve_fit(single_harmonic, x, y)
+                fit_results5.append(popt)
+            except RuntimeError:
+                # 如果拟合过程出现问题（例如无法收敛等），也将参数设为0
+                fit_results5.append((0, 0, 0))
+        # 将拟合结果转换为DataFrame
+        fit_results5_df = pd.DataFrame(fit_results5, columns=['A', 'phi', 'C'])
+        print(fit_results5_df)
+        # 用于存储每个高度的拟合值
+        wn5_values = []
+        for index, row in fit_results5_df.iterrows():
+            A, phi, C = row
+            x = np.arange(12)  # 同样对应12个位置作为自变量
+            wn5 = single_harmonic(x, A, phi, C)
+            wn5_values.append(wn5)
+        # 将拟合值转换为DataFrame
+        wn5_df = pd.DataFrame(wn5_values, columns=[
+                              f'wn5_{i}' for i in range(12)])
+        print(wn5_df)
+        # ---------计算temp-wn0-wn1-wn2-wn3------------------------------------------------------
+        temp_wn0_wn1_wn2_wn3_wn4_wn5 = temp_wn0_wn1_wn2_wn3_wn4_df.values - wn5_df.values
+        # 转换为 DataFrame
+        temp_wn0_wn1_wn2_wn3_wn4_wn5_df = pd.DataFrame(temp_wn0_wn1_wn2_wn3_wn4_wn5,
+                                                       columns=temp_wn0_wn1_wn2_wn3_wn4_df.columns)
+
+        # ------计算背景温度=wn0+wn1+wn2+wn3+wn4+wn5---------------------------------------------------
+        background = wn5_df.values + wn4_df.values + \
+            wn3_df.values + wn2_df.values + wn1_df.values
+        # wn0只有一列单独处理相加
+        # 使用 np.isnan 和 np.where 来判断是否为 NaN 或 0，避免这些值参与相加
+        for i in range(21):
+            wn0_value = wn0_df.iloc[i]
+            # 只有当 wn0_value 既不是 NaN 也不是 0 时才加到 background 上
+            if not np.isnan(wn0_value) and wn0_value != 0:
+                background[i, :] += wn0_value
+        # 扰动
+        perturbation = temp_wn0_wn1_wn2_wn3_wn4_wn5_df
+        # ---------傅里叶变换----------------------------------------------------------------------
+        # 初始化一个新的DataFrame来保存处理结果
+        result = pd.DataFrame(
+            np.nan, index=perturbation.index, columns=perturbation.columns)
+        # 定义滤波范围
+        lambda_low = 2  # 2 km
+        lambda_high = 15  # 15 km
+        f_low = 2 * np.pi / lambda_high
+        f_high = 2 * np.pi / lambda_low
+
+        # 循环处理perturbation中的每一列
+        for col in perturbation.columns:
+            x = perturbation[col]
+            # 提取有效值
+            valid_values = x.dropna()
+            N = len(valid_values)  # 有效值的数量
+
+            # 找到第一个有效值的索引
+            first_valid_index = valid_values.index[0] if not valid_values.index.empty else None
+            height_value = height_df.loc[first_valid_index] if first_valid_index is not None else None
+
+            # 如果有效值为空，则跳过该列
+            if N == 0 or height_value is None:
+                continue
+
+            # 时间序列和频率
+            dt = 0.25
+            n = np.arange(N)
+            t = height_value.values + n * dt
+            f = n / (N * dt)
+
+            # 傅里叶变换
+            y = np.fft.fft(valid_values.values)
+
+            # 频率滤波
+            yy = y.copy()
+            freq_filter = (f < f_low) | (f > f_high)
+            yy[freq_filter] = 0  # 过滤掉指定频段
+
+            # 逆傅里叶变换
+            perturbation_after = np.real(np.fft.ifft(yy))
+
+            # 将处理结果插回到result矩阵中
+            result.loc[valid_values.index, col] = perturbation_after
+            u2 = result ** 2
+            u2 = u2.mean(axis=1)
+        return u2
+    except FileNotFoundError:
+        # 如果文件不存在，返回全NaN的Series
+        expected_length = 21
+        return pd.Series(np.nan, index=range(expected_length))
+
+
+# -------------------------------------------------------------------------------------------
+# --------meridional-------------------------------------------------------------------------
+def process_vmeridional_day(day, year):
+    try:
+        # 读取数据
+        height_data = loadmat(rf"./tidi/data\{year}\{day:03d}_Height.mat")
+        lat_data = loadmat(rf"./tidi/data\{year}\{day:03d}_Lat.mat")
+        lon_data = loadmat(rf"./tidi/data\{year}\{day:03d}_Lon.mat")
+        vmeridional_data = loadmat(
+            rf"./tidi/data\{year}\{day:03d}_VMerdional.mat")
+        vzonal_data = loadmat(rf"./tidi/data\{year}\{day:03d}_Vzonal.mat")
+
+        # 将数据转换为DataFrame
+        height_df = pd.DataFrame(height_data['Height'])
+        lat_df = pd.DataFrame(lat_data['Lat'])
+        lon_df = pd.DataFrame(lon_data['Lon'])
+        vmeridional_df = pd.DataFrame(vmeridional_data['VMerdional'])
+        vzonal_df = pd.DataFrame(vzonal_data['Vzonal'])
+        # 将经纬度拼接为两列并添加到对应的DataFrame中
+        lon_lat_df = pd.concat([lon_df, lat_df], axis=1)
+        lon_lat_df.columns = ['Longitude', 'Latitude']
+        # 筛选出10到30度纬度范围的数据
+        lat_filter = (lat_df.values >= 0) & (lat_df.values <= 20)
+        # 使用纬度范围过滤数据
+        vmeridional_filtered = vmeridional_df.iloc[:, lat_filter.flatten()]
+        vzonal_filtered = vzonal_df.iloc[:, lat_filter.flatten()]
+        lon_lat_filtered = lon_lat_df.iloc[lat_filter.flatten(), :]
+        # 接着对lon_lat_filtered的经度进行分组，0到360度每30度一个区间
+        bins = range(0, 361, 30)
+        group_labels = [f"{i}-{i + 29}" for i in range(0, 360, 30)]
+
+        lon_lat_filtered['Longitude_Group'] = pd.cut(
+            lon_lat_filtered['Longitude'], bins=bins, labels=group_labels)
+
+        # 获取所有唯一的经度分组标签并按照数值顺序排序
+        unique_groups = sorted(lon_lat_filtered['Longitude_Group'].unique(
+        ), key=lambda x: int(x.split('-')[0]))
+        # 按照经度分组获取每个区间对应的vzonal_filtered、vmeridional_filtered数据
+        grouped_data = {}
+        insufficient_data_count = 0  # 用于计数数据不足的组数
+
+        for group in unique_groups:
+            mask = lon_lat_filtered['Longitude_Group'] == group
+            grouped_data[group] = {
+                'vzonal_filtered': vzonal_filtered.loc[:, mask],
+                'vmeridional_filtered': vmeridional_filtered.loc[:, mask],
+                'lon_lat_filtered': lon_lat_filtered.loc[mask]
+            }
+
+            # 计算有效值数量
+            vzonal_count = grouped_data[group]['vzonal_filtered'].notna(
+            ).sum().sum()
+            vmeridional_count = grouped_data[group]['vmeridional_filtered'].notna(
+            ).sum().sum()
+
+            if vzonal_count <= 20 or vmeridional_count <= 20:
+                insufficient_data_count += 1
+
+        # 如果超过6组数据不足，则抛出错误
+        if insufficient_data_count > 6:
+            expected_length = 21
+            return pd.Series(np.nan, index=range(expected_length))
+
+        # 如果代码运行到这里，说明所有分组的数据量都足够或者不足的组数不超过6
+        print("所有分组的数据量都足够")
+        # -----------计算w0------------------------------------------------------------------------------------------
+        # 定义期望的12个区间的分组标签
+        expected_groups = [f"{i}-{i + 29}" for i in range(0, 360, 30)]
+        # 初始化一个空DataFrame来存储所有区间的均值廓线，列名设置为期望的分组标签
+        W0_profiles_df = pd.DataFrame(columns=expected_groups)
+
+        # 遍历grouped_data字典中的每个组
+        for group, data in grouped_data.items():
+            # 提取当前组的vzonal_filtered数据
+            vmeridional_filtered = data['vmeridional_filtered']
+            # 计算有效数据的均值廓线，跳过NaN值
+            mean_profile = vmeridional_filtered.mean(axis=1, skipna=True)
+            # 将当前组的均值廓线作为一列添加到W0_profiles_df DataFrame中
+            W0_profiles_df[group] = mean_profile
+
+        # 检查并填充缺失的区间列，将缺失的列添加并填充为NaN
+        for group in expected_groups:
+            if group not in W0_profiles_df.columns:
+                W0_profiles_df[group] = pd.Series(
+                    [float('NaN')] * len(W0_profiles_df))
+
+        # 打印拼接后的DataFrame以验证
+        print("Concatenated mean profiles for all longitude groups:\n", W0_profiles_df)
+        # 计算每个高度的均值
+        height_mean_profiles = W0_profiles_df.mean(axis=1)
+        # 将每个高度的均值作为新的一行添加到DataFrame中，All_Heights_Mean就是wn0
+        W0_profiles_df['All_Heights_Mean'] = height_mean_profiles
+        wn0_df = W0_profiles_df['All_Heights_Mean']
+        # -------计算残余量--------------------------------------------------------------------------------------
+        # 计算每个经度区间的残余值 (即每个区间的值减去该高度的All_Heights_Mean)
+        residuals_df = W0_profiles_df.drop(columns='All_Heights_Mean').subtract(
+            W0_profiles_df['All_Heights_Mean'], axis=0)
+
+        # --------wn1-------------------------
+        def single_harmonic(x, A, phi, C):
+            return A * np.sin(2 * np.pi / 12 * x + phi) + C
+        # 用于存储每个高度拟合后的参数结果
+        fit_results = []
+        for index, row in residuals_df.iterrows():
+            # 检查该行是否存在NaN值，如果有则跳过拟合，直接将参数设为0
+            if row.isnull().any():
+                fit_results.append((0, 0, 0))
+                continue
+            x = np.arange(12)  # 对应12个位置作为自变量
+            y = row.values
+            try:
+                # 进行曲线拟合
+                popt, _ = curve_fit(single_harmonic, x, y)
+                fit_results.append(popt)
+            except RuntimeError:
+                # 如果拟合过程出现问题（例如无法收敛等），也将参数设为0
+                fit_results.append((0, 0, 0))
+        # 将拟合结果转换为DataFrame
+        fit_results_df = pd.DataFrame(fit_results, columns=['A', 'phi', 'C'])
+        print(fit_results_df)
+        # 用于存储每个高度的拟合值
+        wn1_values = []
+        for index, row in fit_results_df.iterrows():
+            A, phi, C = row
+            x = np.arange(12)  # 同样对应12个位置作为自变量
+            wn1 = single_harmonic(x, A, phi, C)
+            wn1_values.append(wn1)
+        # 将拟合值转换为DataFrame
+        wn1_df = pd.DataFrame(wn1_values, columns=[
+                              f'wn1_{i}' for i in range(12)])
+        print(wn1_df)
+        # 如果wn1_df全为0，则跳过下面的计算，直接令该天的day_log_gwresult全部为NaN
+        if (wn1_df == 0).all().all():
+            return pd.Series(np.nan, index=range(21))
+        # ------------计算temp-wn0-wn1---------------------------------------------------------
+        temp_wn0_wn1 = residuals_df.values - wn1_df.values
+        # 将结果转为 DataFrame
+        temp_wn0_wn1_df = pd.DataFrame(
+            temp_wn0_wn1, columns=residuals_df.columns, index=residuals_df.index)
+
+        # -------wn2--------------------------------------------------------------------------
+        def single_harmonic(x, A, phi, C):
+            return A * np.sin(2 * np.pi / 6 * x + phi) + C
+
+        # 用于存储每个高度拟合后的参数结果
+        fit_results2 = []
+        for index, row in temp_wn0_wn1_df.iterrows():
+            # 检查该行是否存在NaN值，如果有则跳过拟合，直接将参数设为0
+            if row.isnull().any():
+                fit_results2.append((0, 0, 0))
+                continue
+            x = np.arange(12)  # 对应12个位置作为自变量
+            y = row.values
+            try:
+                # 进行曲线拟合
+                popt, _ = curve_fit(single_harmonic, x, y)
+                fit_results2.append(popt)
+            except RuntimeError:
+                # 如果拟合过程出现问题（例如无法收敛等），也将参数设为0
+                fit_results2.append((0, 0, 0))
+        # 将拟合结果转换为DataFrame
+        fit_results2_df = pd.DataFrame(fit_results2, columns=['A', 'phi', 'C'])
+        print(fit_results2_df)
+        # 用于存储每个高度的拟合值
+        wn2_values = []
+        for index, row in fit_results2_df.iterrows():
+            A, phi, C = row
+            x = np.arange(12)  # 同样对应12个位置作为自变量
+            wn2 = single_harmonic(x, A, phi, C)
+            wn2_values.append(wn2)
+        # 将拟合值转换为DataFrame
+        wn2_df = pd.DataFrame(wn2_values, columns=[
+                              f'wn2_{i}' for i in range(12)])
+        print(wn2_df)
+        # ---------计算temp-wn0-wn1-wn2------------------------------------------------------
+        temp_wn0_wn1_wn2 = temp_wn0_wn1_df.values - wn2_df.values
+        # 转换为 DataFrame
+        temp_wn0_wn1_wn2_df = pd.DataFrame(
+            temp_wn0_wn1_wn2, columns=temp_wn0_wn1_df.columns)
+
+        # -------wn3-----------------------------------------------------------------------
+        def single_harmonic(x, A, phi, C):
+            return A * np.sin(2 * np.pi / 4 * x + phi) + C
+
+        # 用于存储每个高度拟合后的参数结果
+        fit_results3 = []
+        for index, row in temp_wn0_wn1_wn2_df.iterrows():
+            # 检查该行是否存在NaN值，如果有则跳过拟合，直接将参数设为0
+            if row.isnull().any():
+                fit_results3.append((0, 0, 0))
+                continue
+            x = np.arange(12)  # 对应12个位置作为自变量
+            y = row.values
+            try:
+                # 进行曲线拟合
+                popt, _ = curve_fit(single_harmonic, x, y)
+                fit_results3.append(popt)
+            except RuntimeError:
+                # 如果拟合过程出现问题（例如无法收敛等），也将参数设为0
+                fit_results3.append((0, 0, 0))
+        # 将拟合结果转换为DataFrame
+        fit_results3_df = pd.DataFrame(fit_results3, columns=['A', 'phi', 'C'])
+        print(fit_results3_df)
+        # 用于存储每个高度的拟合值
+        wn3_values = []
+        for index, row in fit_results3_df.iterrows():
+            A, phi, C = row
+            x = np.arange(12)  # 同样对应12个位置作为自变量
+            wn3 = single_harmonic(x, A, phi, C)
+            wn3_values.append(wn3)
+        # 将拟合值转换为DataFrame
+        wn3_df = pd.DataFrame(wn3_values, columns=[
+                              f'wn3_{i}' for i in range(12)])
+        print(wn3_df)
+        # ---------计算temp-wn0-wn1-wn2-wn3------------------------------------------------------
+        temp_wn0_wn1_wn2_wn3 = temp_wn0_wn1_wn2_df.values - wn3_df.values
+        # 转换为 DataFrame
+        temp_wn0_wn1_wn2_wn3_df = pd.DataFrame(
+            temp_wn0_wn1_wn2_wn3, columns=temp_wn0_wn1_wn2_df.columns)
+        # -------wn4 - ----------------------------------------------------------------------
+
+        def single_harmonic(x, A, phi, C):
+            return A * np.sin(2 * np.pi / 3 * x + phi) + C
+
+        # 用于存储每个高度拟合后的参数结果
+        fit_results4 = []
+        for index, row in temp_wn0_wn1_wn2_wn3_df.iterrows():
+            # 检查该行是否存在NaN值，如果有则跳过拟合，直接将参数设为0
+            if row.isnull().any():
+                fit_results4.append((0, 0, 0))
+                continue
+            x = np.arange(12)  # 对应12个位置作为自变量
+            y = row.values
+            try:
+                # 进行曲线拟合
+                popt, _ = curve_fit(single_harmonic, x, y)
+                fit_results4.append(popt)
+            except RuntimeError:
+                # 如果拟合过程出现问题（例如无法收敛等），也将参数设为0
+                fit_results4.append((0, 0, 0))
+
+        fit_results4_df = pd.DataFrame(fit_results4, columns=['A', 'phi', 'C'])
+        print(fit_results4_df)
+        # 用于存储每个高度的拟合值
+        wn4_values = []
+        for index, row in fit_results4_df.iterrows():
+            A, phi, C = row
+            x = np.arange(12)  # 同样对应12个位置作为自变量
+            wn4 = single_harmonic(x, A, phi, C)
+            wn4_values.append(wn4)
+        # 将拟合值转换为DataFrame
+        wn4_df = pd.DataFrame(wn4_values, columns=[
+                              f'wn4_{i}' for i in range(12)])
+        print(wn4_df)
+        # ---------计算temp-wn0-wn1-wn2-wn3------------------------------------------------------
+        temp_wn0_wn1_wn2_wn3_wn4 = temp_wn0_wn1_wn2_wn3_df.values - wn4_df.values
+        # 转换为 DataFrame
+        temp_wn0_wn1_wn2_wn3_wn4_df = pd.DataFrame(
+            temp_wn0_wn1_wn2_wn3_wn4, columns=temp_wn0_wn1_wn2_wn3_df.columns)
+
+        # -------wn5-----------------------------------------------------------------------
+        def single_harmonic(x, A, phi, C):
+            return A * np.sin(2 * np.pi / 2.4 * x + phi) + C
+
+        # 用于存储每个高度拟合后的参数结果
+        fit_results5 = []
+        for index, row in temp_wn0_wn1_wn2_wn3_wn4_df.iterrows():
+            # 检查该行是否存在NaN值，如果有则跳过拟合，直接将参数设为0
+            if row.isnull().any():
+                fit_results5.append((0, 0, 0))
+                continue
+            x = np.arange(12)  # 对应12个位置作为自变量
+            y = row.values
+            try:
+                # 进行曲线拟合
+                popt, _ = curve_fit(single_harmonic, x, y)
+                fit_results5.append(popt)
+            except RuntimeError:
+                # 如果拟合过程出现问题（例如无法收敛等），也将参数设为0
+                fit_results5.append((0, 0, 0))
+        # 将拟合结果转换为DataFrame
+        fit_results5_df = pd.DataFrame(fit_results5, columns=['A', 'phi', 'C'])
+        print(fit_results5_df)
+        # 用于存储每个高度的拟合值
+        wn5_values = []
+        for index, row in fit_results5_df.iterrows():
+            A, phi, C = row
+            x = np.arange(12)  # 同样对应12个位置作为自变量
+            wn5 = single_harmonic(x, A, phi, C)
+            wn5_values.append(wn5)
+        # 将拟合值转换为DataFrame
+        wn5_df = pd.DataFrame(wn5_values, columns=[
+                              f'wn5_{i}' for i in range(12)])
+        print(wn5_df)
+        # ---------计算temp-wn0-wn1-wn2-wn3------------------------------------------------------
+        temp_wn0_wn1_wn2_wn3_wn4_wn5 = temp_wn0_wn1_wn2_wn3_wn4_df.values - wn5_df.values
+        # 转换为 DataFrame
+        temp_wn0_wn1_wn2_wn3_wn4_wn5_df = pd.DataFrame(temp_wn0_wn1_wn2_wn3_wn4_wn5,
+                                                       columns=temp_wn0_wn1_wn2_wn3_wn4_df.columns)
+
+        # ------计算背景温度=wn0+wn1+wn2+wn3+wn4+wn5---------------------------------------------------
+        background = wn5_df.values + wn4_df.values + \
+            wn3_df.values + wn2_df.values + wn1_df.values
+        # wn0只有一列单独处理相加
+        # 使用 np.isnan 和 np.where 来判断是否为 NaN 或 0，避免这些值参与相加
+        for i in range(21):
+            wn0_value = wn0_df.iloc[i]
+            # 只有当 wn0_value 既不是 NaN 也不是 0 时才加到 background 上
+            if not np.isnan(wn0_value) and wn0_value != 0:
+                background[i, :] += wn0_value
+        # 扰动
+        perturbation = temp_wn0_wn1_wn2_wn3_wn4_wn5_df
+        # ---------傅里叶变换----------------------------------------------------------------------
+        # 初始化一个新的DataFrame来保存处理结果
+        result = pd.DataFrame(
+            np.nan, index=perturbation.index, columns=perturbation.columns)
+        # 定义滤波范围
+        lambda_low = 2  # 2 km
+        lambda_high = 15  # 15 km
+        f_low = 2 * np.pi / lambda_high
+        f_high = 2 * np.pi / lambda_low
+
+        # 循环处理perturbation中的每一列
+        for col in perturbation.columns:
+            x = perturbation[col]
+            # 提取有效值
+            valid_values = x.dropna()
+            N = len(valid_values)  # 有效值的数量
+
+            # 找到第一个有效值的索引
+            first_valid_index = valid_values.index[0] if not valid_values.index.empty else None
+            height_value = height_df.loc[first_valid_index] if first_valid_index is not None else None
+
+            # 如果有效值为空，则跳过该列
+            if N == 0 or height_value is None:
+                continue
+
+            # 时间序列和频率
+            dt = 0.25
+            n = np.arange(N)
+            t = height_value.values + n * dt
+            f = n / (N * dt)
+
+            # 傅里叶变换
+            y = np.fft.fft(valid_values.values)
+
+            # 频率滤波
+            yy = y.copy()
+            freq_filter = (f < f_low) | (f > f_high)
+            yy[freq_filter] = 0  # 过滤掉指定频段
+
+            # 逆傅里叶变换
+            perturbation_after = np.real(np.fft.ifft(yy))
+
+            # 将处理结果插回到result矩阵中
+            result.loc[valid_values.index, col] = perturbation_after
+        v2 = result ** 2
+        v2 = v2.mean(axis=1)
+        return v2
+    except FileNotFoundError:
+        # 如果文件不存在，返回全NaN的Series
+        expected_length = 21
+        return pd.Series(np.nan, index=range(expected_length))
+# 初始化一个空的DataFrame来存储所有天的结果
+
+
+# 初始化一个空的DataFrame来存储所有天的结果
+all_days_vzonal_results = pd.DataFrame()
+
+# 循环处理每一天的数据
+for day in range(1, 365):
+    u2 = process_vzonal_day(day, 2019)
+    if u2 is not None:
+        all_days_vzonal_results[rf"{day:02d}"] = u2
+
+# 将结果按列拼接
+# all_days_vzonal_results.columns = [f"{day:02d}" for day in range(1, 365)]
+
+all_days_vmeridional_results = pd.DataFrame()
+
+# 循环处理每一天的数据
+for day in range(1, 365):
+    v2 = process_vmeridional_day(day, 2019)
+    if v2 is not None:
+        all_days_vmeridional_results[rf"{day:02d}"] = v2
+
+# 将结果按列拼接
+# all_days_vmeridional_results.columns = [f"{day:02d}" for day in range(1, 365)]
+
+# ---------------------------------------------------------------------------------------------------
+# --------经纬向风平方和计算动能--------------------------------------------------------------------------------
+
+# 使用numpy.where来检查两个表格中的对应元素是否都不是NaN
+sum_df = np.where(
+    pd.notna(all_days_vmeridional_results) & pd.notna(all_days_vzonal_results),
+    all_days_vmeridional_results + all_days_vzonal_results,
+    np.nan
+)
+HP = 1/2*all_days_vmeridional_results+1/2*all_days_vzonal_results
+heights = [70.0, 72.5, 75.0, 77.5, 80.0, 82.5, 85.0, 87.5, 90.0, 92.5,
+           95.0, 97.5, 100.0, 102.5, 105.0, 107.5, 110.0, 112.5, 115.0, 117.5, 120.0]
+HP.index = heights
+# # 将 DataFrame 保存为 Excel 文件
+# HP.to_excel('HP_data.xlsx')
+# ----------绘年统计图------------------------------------------------------------------------------------------------------------
+data = HP
+# 使用 reset_index() 方法将索引变为第一列
+data = data.reset_index()
+h = data.iloc[:, 0].copy()  # 高度，保留作为纵坐标
+dates = list(range(1, data.shape[1]))  # 日期，作为横坐标
+data0 = data.iloc[:, 1:].copy()  # 绘图数据
+'''数据处理'''
+# 反转 h 以确保高度从下往上递增
+h_reversed = h[::-1].reset_index(drop=True)
+data0_reversed = data0[::-1].reset_index(drop=True)
+# 将数值大于20的数据点替换为nan
+data0_reversed[data0_reversed > 20] = float('nan')
+# 转换成月份，365天
+days_in_month = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
+# 将日期转换为英文月份
+
+
+def day_to_month(day):
+    # 累积每个月的天数，找到对应的月份
+    cumulative_days = 0
+    for i, days in enumerate(days_in_month):
+        cumulative_days += days
+        if day <= cumulative_days:
+            return f'{["January", "February", "March", "April", "May", "June", "July", "August", "September", "October", "November", "December"][i]}'
+
+
+months = [day_to_month(day) for day in dates]
+
+
+'''绘图'''
+plt.rcParams['font.family'] = 'SimSun'  # 宋体
+plt.rcParams['font.size'] = 12           # 中文字号
+plt.rcParams['axes.unicode_minus'] = False  # 正确显示负号
+plt.rcParams['font.sans-serif'] = 'Times New Roman'  # 新罗马
+plt.rcParams['axes.labelsize'] = 14       # 坐标轴标签字号
+plt.rcParams['xtick.labelsize'] = 12      # x轴刻度字号
+plt.rcParams['ytick.labelsize'] = 12      # y轴刻度字号
+plt.rcParams['legend.fontsize'] = 16       # 图例字号
+plt.rcParams['axes.unicode_minus'] = False  # 正确显示负号
+plt.figure(figsize=(10, 6))  # 设置图像大小
+# 绘制热力图，设置 x 和 y 轴的标签
+sns.heatmap(data0_reversed, annot=False, cmap='YlGnBu', linewidths=0.5,
+            yticklabels=h_reversed, xticklabels=months, cbar_kws={'label': ' Gravitational potential energy'})
+
+# 横坐标过长，设置等间隔展示
+interval = 34  # 横坐标显示间隔
+plt.xticks(ticks=range(0, len(dates), interval),
+           labels=months[::interval], rotation=45)  # rotation旋转可不加
+
+# 添加轴标签
+plt.xlabel('Month')  # X轴标签
+plt.ylabel('Height')  # Y轴标签
+
+# 显示图形
+plt.show()
+# --------------绘制月统计图-------------------------------------------------------------------
+# 获取HP的列数
+num_cols = HP.shape[1]
+# 用于存储按要求计算出的均值列数据
+mean_cols = []
+start = 0
+while start < num_cols:
+    end = start + 30
+    if end > num_cols:
+        end = num_cols
+    # 提取每30列（或不满30列的剩余部分）的数据
+    subset = HP.iloc[:, start:end]
+    # 计算该部分数据每一行的均值，得到一个Series，作为新的均值列
+    mean_series = subset.mean(axis=1)
+    mean_cols.append(mean_series)
+    start = end
+# 将所有的均值列合并成一个新的DataFrame
+result_df = pd.concat(mean_cols, axis=1)
+# 对result_df中的每一个元素取自然对数
+result_df_log = result_df.applymap(lambda x: np.log(x))
+# 通过drop方法删除第一行，axis=0表示按行操作，inplace=True表示直接在原DataFrame上修改（若不想修改原DataFrame可设置为False）
+result_df_log.drop(70, axis=0, inplace=True)
+# 计算每个月的平均值
+monthly_average = result_df_log.mean(axis=0)
+# 将结果转换为 (1, 12) 形状
+monthly_average = monthly_average.values.reshape(1, 12)
+monthly_average = monthly_average.ravel()
+
+# 生成x轴的月份标签
+months = ["Jan", "Feb", "Mar", "Apr", "May", "Jun",
+          "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"]
+
+# 绘制折线图
+plt.plot(months, monthly_average, marker='o', linestyle='-', color='b')
+
+# 添加标题和标签
+plt.title("Monthly Average ENERGY(log)")
+plt.xlabel("Month")
+plt.ylabel("Average Energy")
+# 显示图表
+plt.xticks(rotation=45)  # 让月份标签更清晰可读
+plt.grid(True)
+plt.tight_layout()
+# 显示图形
+plt.show()