diff --git a/backend.py b/backend.py index b2436db..13140f6 100644 --- a/backend.py +++ b/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") diff --git a/balloon/__init__.py b/balloon/__init__.py index 8fb8d09..33e5b1f 100644 --- a/balloon/__init__.py +++ b/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(): diff --git a/cosmic/__init__.py b/cosmic/__init__.py new file mode 100644 index 0000000..ddeda83 --- /dev/null +++ b/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") \ No newline at end of file diff --git a/cosmic/temp_render.py b/cosmic/temp_render.py new file mode 100644 index 0000000..0ec7ded --- /dev/null +++ b/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() # 显示图形 \ No newline at end of file diff --git a/saber/__init__.py b/saber/__init__.py index 0a838a1..50777ba 100644 --- a/saber/__init__.py +++ b/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 diff --git a/staged/cosmic行星波参数全年逐日绘图.py b/staged/cosmic行星波参数全年逐日绘图.py new file mode 100644 index 0000000..3cd3e6f --- /dev/null +++ b/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() # 显示图形 \ No newline at end of file diff --git a/staged/cosmic重力波多天.py b/staged/cosmic重力波多天.py new file mode 100644 index 0000000..f829d0a --- /dev/null +++ b/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() \ No newline at end of file diff --git a/tidi/__init__.py b/tidi/__init__.py index 07123b7..d529038 100644 --- a/tidi/__init__.py +++ b/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") \ No newline at end of file + 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") diff --git a/tidi/plot.py b/tidi/plot.py index 32d9c3e..05e905f 100644 --- a/tidi/plot.py +++ b/tidi/plot.py @@ -1,878 +1,903 @@ -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() +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=2015): + 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: + # raise ValueError("Insufficient data for more than 6 longitude groups in the specified latitude band.") + + # 如果代码运行到这里，说明所有分组的数据量都足够或者不足的组数不超过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)) + + +# 初始化一个空的DataFrame来存储所有天的结果 + + + +# ------------------------------------------------------------------------------------------- +# --------meridional------------------------------------------------------------------------- +def process_vmeridional_day(day, year=2015): + 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: + # raise ValueError( + # "Insufficient data for more than 6 longitude groups in the specified latitude band.") + + # 如果代码运行到这里，说明所有分组的数据量都足够或者不足的组数不超过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)) + + +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]}' + +class TidiPlotv2: + def __init__(self, year): + self.year = year + + + + all_days_vzonal_results = pd.DataFrame() + + # 循环处理每一天的数据 + for day in range(1, 365): + u2 = process_vzonal_day(day, year) + all_days_vzonal_results[rf"{day:02d}"] = u2 + + # 将结果按列拼接 + all_days_vzonal_results.columns = [f"{day:02d}" for day in range(1, 365)] + # 初始化一个空的DataFrame来存储所有天的结果 + all_days_vmeridional_results = pd.DataFrame() + + # 循环处理每一天的数据 + for day in range(1, 365): + v2 = process_vmeridional_day(day, year) + all_days_vmeridional_results[rf"{day:02d}"] = v2 + + # 将结果按列拼接 + all_days_vmeridional_results.columns = [f"{day:02d}" for day in range(1, 365)] + + self.all_days_vzonal_results = all_days_vzonal_results + self.all_days_vmeridional_results = all_days_vmeridional_results + + + + # --------------------------------------------------------------------------------------------------- + # --------经纬向风平方和计算动能-------------------------------------------------------------------------------- + + # 使用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 以确保高度从下往上递增 + self.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天 + self.data0_reversed = data0_reversed + + self.HP = HP + self.dates = dates + self.months = [day_to_month(day) for day in dates] + + def plot_v1(self): + h_reversed = self.h_reversed + data0_reversed = self.data0_reversed + dates = self.dates + months = self.months + '''绘图''' + 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() + + def plot_month(self): + HP = self.HP + # --------------绘制月统计图------------------------------------------------------------------- + # 获取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) 形状 + + # 如果大于12， 就把最后一个给删了，保证能转换为(1, 12) + if len(monthly_average) > 12: + monthly_average = monthly_average[:-1] + + 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() diff --git a/tidi/staged/plot_old.py b/tidi/staged/plot_old.py new file mode 100644 index 0000000..32d9c3e --- /dev/null +++ b/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()