999 lines
36 KiB
Python
999 lines
36 KiB
Python
import os
|
||
import pandas as pd
|
||
import numpy as np
|
||
import matplotlib.pyplot as plt
|
||
from scipy.optimize import curve_fit
|
||
import seaborn as sns
|
||
|
||
# 解决绘图中中文不能显示的问题
|
||
import matplotlib
|
||
# 设置中文显示和负号正常显示
|
||
matplotlib.rcParams['font.sans-serif'] = ['SimHei'] # 显示中文
|
||
matplotlib.rcParams['axes.unicode_minus'] = False # 正常显示负号
|
||
|
||
# 定义所有模型
|
||
def tidal_model(t, v0, a1, b1, a2, b2, a3, b3, a4, b4):
|
||
"""潮汐波模型"""
|
||
T1, T2, T3, T4 = 6, 8, 12, 24 # 周期以小时为单位
|
||
return (v0 + a1 * np.sin((2 * np.pi / T1) * t + b1)
|
||
+ a2 * np.sin((2 * np.pi / T2) * t + b2)
|
||
+ a3 * np.sin((2 * np.pi / T3) * t + b3)
|
||
+ a4 * np.sin((2 * np.pi / T4) * t + b4))
|
||
|
||
def planetary_model_2day(t, v0, a, b):
|
||
"""2日行星波模型"""
|
||
T = 48 # 周期为2天（48小时）
|
||
return v0 + a * np.sin((2 * np.pi / T) * t + b)
|
||
|
||
def planetary_model_5day(t, v0, a, b):
|
||
"""5日行星波模型"""
|
||
T = 120 # 周期为5天（120小时）
|
||
return v0 + a * np.sin((2 * np.pi / T) * t + b)
|
||
|
||
def planetary_model_10day(t, v0, a, b):
|
||
"""10日行星波模型"""
|
||
T = 240 # 周期为10天（240小时）
|
||
return v0 + a * np.sin((2 * np.pi / T) * t + b)
|
||
|
||
def planetary_model_16day(t, v0, a, b):
|
||
"""16日行星波模型"""
|
||
T = 384 # 周期为16天（384小时）
|
||
return v0 + a * np.sin((2 * np.pi / T) * t + b)
|
||
|
||
# ----------------------------------------------------------第二大部分————————————————————————————————————————————————————
|
||
# # 根据纬向风(u)、经向风(v),得到武汉左岭镇站/黑龙江漠河站
|
||
# 大气波动（6、8、12、24小时潮汐波，2日、5日、10日、16日行星波）随时空变化热力图
|
||
|
||
# 设置文件路径
|
||
year = 2022
|
||
station = '武汉左岭镇站'
|
||
# station = '黑龙江漠河站'
|
||
|
||
# 文件路径设置
|
||
file_dir = rf'D:\wu\gravity wave\Meteor\流星雷达分析数据\{station}\{year}' # 数据文件路径
|
||
output_dir = rf'D:\wu\gravity wave\Meteor\输出\{station}\{year}\时空' # 参数txt、图片输出路径
|
||
os.makedirs(output_dir, exist_ok=True) # 确保输出路径存在
|
||
|
||
|
||
# 获取目录下所有txt文件
|
||
file_list = [f for f in os.listdir(file_dir) if f.endswith('.txt')]
|
||
|
||
# 初始化空的DataFrame，用于合并所有文件的数据
|
||
final_df = pd.DataFrame()
|
||
|
||
# 批量处理每个txt文件
|
||
for file_name in file_list:
|
||
file_path = os.path.join(file_dir, file_name)
|
||
df = pd.read_csv(file_path, delim_whitespace=True)
|
||
|
||
date_part = os.path.splitext(file_name)[0][-8:] # 假设文件名最后8个字符是日期
|
||
|
||
# 替换异常值为 NaN
|
||
df.loc[df['uwind'] == 1000000, 'uwind'] = np.nan
|
||
df.loc[df['vwind'] == 1000000, 'vwind'] = np.nan
|
||
|
||
|
||
# 筛选出所需的列
|
||
df = df[['height', 'time', 'uwind', 'vwind']]
|
||
df = df.rename(columns={
|
||
'uwind': f'{date_part}_uwind',
|
||
'vwind': f'{date_part}_vwind'
|
||
})
|
||
|
||
if final_df.empty:
|
||
final_df = df
|
||
else:
|
||
final_df = pd.merge(final_df, df, on=['height', 'time'], how='outer')
|
||
|
||
# 生成完整日期范围的列名，补全缺失的天，使得平年365，闰年366
|
||
all_dates = pd.date_range(f'{year}-01-01', f'{year}-12-31').strftime('%Y%m%d')
|
||
expected_columns = [f'{date}_uwind' for date in all_dates] + [f'{date}_vwind' for date in all_dates]
|
||
|
||
# 确保 final_df 包含所有日期的列，缺失列补为 NaN
|
||
for col in expected_columns:
|
||
if col not in final_df.columns:
|
||
final_df[col] = np.nan
|
||
|
||
# 按列名排序（确保日期顺序）
|
||
final_df = final_df[['height', 'time'] + sorted(expected_columns)]
|
||
# 此时，final_df已经是全高度包含完整天的分析数据---------------------------------------------------------
|
||
|
||
|
||
# -----------------------------------------得到潮汐波拟合参数-------------------------------------------
|
||
# 定义单高度的拟合函数
|
||
def fit_wind_speed(data, window=5):
|
||
# 动态调整最小有效数据点数量：参数数量 * 6
|
||
num_params = 9 # 模型参数数量
|
||
|
||
# 转换为 pandas DataFrame 并设置时间为索引
|
||
data["date"] = pd.to_datetime(data["date"])
|
||
data.set_index("date", inplace=True)
|
||
|
||
# 获取唯一日期
|
||
unique_dates = np.unique(data.index.date)
|
||
|
||
# 存储拟合参数
|
||
params = []
|
||
|
||
for date in unique_dates:
|
||
start_date = pd.to_datetime(date) - pd.Timedelta(days=(window - 1) / 2)
|
||
end_date = pd.to_datetime(date) + pd.Timedelta(days=(window - 1) / 2)
|
||
|
||
# 获取这几天的数据
|
||
window_data = data[start_date:end_date]
|
||
|
||
t = np.arange(len(window_data)) + 1 # 时间点
|
||
y = window_data['wind_speed'].values # 风速数据
|
||
|
||
# 将 y 转换为浮点类型，以防有非数值类型数据
|
||
y = pd.to_numeric(y, errors='coerce') # 将非数值转换为 NaN
|
||
|
||
# 去掉包含nan的那一列，筛选有效数据
|
||
valid_mask = ~np.isnan(y) # 创建布尔掩码，标记非 NaN 的位置
|
||
t = t[valid_mask]
|
||
y = y[valid_mask]
|
||
|
||
# 检查筛选后的数据是否足够
|
||
if len(y) < num_params*6: # 确保数据足够
|
||
# print(f"Not enough valid data after filtering for date: {date}")
|
||
params.append([None] * 9) # 如果数据不足，记录None
|
||
continue
|
||
|
||
# 设置初始参数
|
||
initial_guess = [0, 1, 0, 1, 0, 1, 0, 1, 0] # v0, a1, b1, a2, b2, a3, b3, a4, b4
|
||
|
||
# 设置参数界限
|
||
bounds = ([-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]) # 上界
|
||
|
||
|
||
# 拟合模型
|
||
try:
|
||
popt, _ = curve_fit(tidal_model, t, y, p0=initial_guess, bounds=bounds)
|
||
params.append(popt) # 存储拟合参数
|
||
except RuntimeError:
|
||
# print(f"Fitting failed for date: {date}")
|
||
params.append([None] *9) # 如果拟合失败，记录None
|
||
|
||
# 转换拟合参数为 DataFrame
|
||
params_df = pd.DataFrame(
|
||
params,
|
||
columns=["v0", "a1", "b1", "a2", "b2", "a3", "b3", "a4", "b4"],
|
||
index=unique_dates,
|
||
)
|
||
return params_df
|
||
|
||
# 1-拟合纬向风------------------------------------------
|
||
# 提取高度、时间和纬向风数据
|
||
uwind_data = final_df[["height", "time"] + [col for col in final_df.columns if "_uwind" in col]]
|
||
uwind_data.columns = [col[:8] for col in uwind_data.columns]
|
||
|
||
# 遍历所有高度
|
||
heights = uwind_data["height"].unique()
|
||
u_result_dict = {}
|
||
|
||
for height in heights:
|
||
print(f"Processing height: {height} m")
|
||
# 获取当前高度的纬向风数据
|
||
height_data = uwind_data[uwind_data["height"] == height]
|
||
|
||
# 转换为长格式（适配拟合函数）
|
||
long_data = pd.melt(
|
||
height_data,
|
||
id_vars=["time", "height"],
|
||
var_name="date",
|
||
value_name="wind_speed",
|
||
)
|
||
long_data["date"] = pd.to_datetime(long_data["date"], format="%Y%m%d")
|
||
|
||
# 调用拟合函数
|
||
params_df = fit_wind_speed(long_data)
|
||
u_result_dict[height] = params_df
|
||
|
||
# 汇总结果为 DataFrame
|
||
u_final_result = pd.concat(u_result_dict, names=["height", "date"])
|
||
# 将结果保存为 TXT 文档
|
||
file_name1 = f'{year}年_纬向风_潮汐波振幅参数.txt'
|
||
u_final_result.to_csv(
|
||
os.path.join(output_dir, file_name),
|
||
sep='\t',
|
||
index=True,
|
||
float_format='%.2f', # 控制数值列保留 2 位小数
|
||
header=True)
|
||
|
||
# 2-拟合经向风------------------------------------------
|
||
# 提取高度、时间和经向风数据
|
||
vwind_data = final_df[["height", "time"] + [col for col in final_df.columns if "_vwind" in col]]
|
||
vwind_data.columns = [col[:8] for col in vwind_data.columns]
|
||
|
||
# 遍历所有高度
|
||
heights = vwind_data["height"].unique()
|
||
v_result_dict = {}
|
||
|
||
for height in heights:
|
||
print(f"Processing height: {height} m")
|
||
height_data = vwind_data[vwind_data["height"] == height]
|
||
|
||
long_data = pd.melt(
|
||
height_data,
|
||
id_vars=["time", "height"],
|
||
var_name="date",
|
||
value_name="wind_speed",
|
||
)
|
||
long_data["date"] = pd.to_datetime(long_data["date"], format="%Y%m%d")
|
||
|
||
params_df = fit_wind_speed(long_data)
|
||
v_result_dict[height] = params_df
|
||
|
||
# 汇总结果为 DataFrame
|
||
v_final_result = pd.concat(v_result_dict, names=["height", "date"])
|
||
|
||
# 将结果保存为 TXT 文档
|
||
# 保存参数
|
||
file_name2 = f'{year}年_经向风_潮汐波振幅参数.txt'
|
||
v_final_result.to_csv(
|
||
os.path.join(output_dir, file_name),
|
||
sep='\t',
|
||
index=True,
|
||
float_format='%.2f', # 控制数值列保留 2 位小数
|
||
header=True)
|
||
|
||
|
||
# 画热力图-----------------------------------------
|
||
'''
|
||
# 读取数据，设置第一列和第二列为多重索引
|
||
u_final_result = pd.read_csv(os.path.join(output_dir, file_name1), sep='\t')
|
||
u_final_result = u_final_result.set_index([u_final_result.columns[0], u_final_result.columns[1]])
|
||
|
||
v_final_result = pd.read_csv(os.path.join(output_dir, file_name2), sep='\t')
|
||
v_final_result = v_final_result.set_index([v_final_result.columns[0], v_final_result.columns[1]])
|
||
'''
|
||
# 定义振幅参数
|
||
amplitude_param = 'a4' # 代表 24 小时周期振幅 # todo: 选择a4为周日潮,a3半日潮
|
||
|
||
# 提取纬向风数据，重置索引为平坦结构
|
||
u_heatmap_data = u_final_result[amplitude_param].unstack(level=0) # 将高度作为列，日期作为行
|
||
u_heatmap_data.index = pd.to_datetime(u_heatmap_data.index) # 确保日期为时间格式
|
||
u_heatmap_data = u_heatmap_data.iloc[:, ::-1] # 按高度降序排列
|
||
|
||
# 提取经向风数据，重置索引为平坦结构
|
||
v_heatmap_data = v_final_result[amplitude_param].unstack(level=0) # 将高度作为列，日期作为行
|
||
v_heatmap_data.index = pd.to_datetime(v_heatmap_data.index) # 确保日期为时间格式
|
||
v_heatmap_data = v_heatmap_data.iloc[:, ::-1] # 按高度降序排列
|
||
|
||
# 创建画布和子图
|
||
fig, axes = plt.subplots(1, 2, figsize=(18, 8)) # 1行2列的布局
|
||
fig.suptitle('2022年周日潮汐波振幅时空变化', fontsize=20) # 添加总标题
|
||
|
||
# 绘制第一幅热力图（纬向风）
|
||
sns.heatmap(
|
||
u_heatmap_data.T, # 转置，横轴为时间，纵轴为高度
|
||
cmap="viridis", # 配色方案
|
||
vmax=100, # 设置颜色范围的最大值
|
||
cbar_kws={'label': '振幅 (m/s)'}, # 设置颜色条标签
|
||
xticklabels=True, # 保留真实的横轴刻度标签
|
||
#yticklabels=u_heatmap_data.columns.astype(str), # 显示纵轴的高度
|
||
yticklabels=[f'{height / 1000}' for height in u_heatmap_data.columns], # 转换高度为千米并格式化
|
||
ax=axes[0] # 指定子图
|
||
)
|
||
axes[0].set_title('纬向风振幅变化', fontsize=14)
|
||
axes[0].set_xlabel('日期', fontsize=12)
|
||
axes[0].set_ylabel('高度 (km)', fontsize=12)
|
||
|
||
# 绘制第二幅热力图（经向风）
|
||
sns.heatmap(
|
||
v_heatmap_data.T, # 转置，横轴为时间，纵轴为高度
|
||
cmap="viridis", # 配色方案
|
||
vmax=100, # 设置颜色范围的最大值
|
||
cbar_kws={'label': '振幅 (m/s)'}, # 设置颜色条标签
|
||
xticklabels=True, # 保留真实的横轴刻度标签
|
||
#yticklabels=v_heatmap_data.columns.astype(str), # 显示纵轴的高度
|
||
yticklabels=[f'{height / 1000}' for height in v_heatmap_data.columns], # 转换高度为千米并格式化
|
||
ax=axes[1] # 指定子图
|
||
)
|
||
axes[1].set_title('经向风振幅变化', fontsize=14)
|
||
axes[1].set_xlabel('日期', fontsize=12)
|
||
axes[1].set_ylabel('高度 (km)', fontsize=12)
|
||
|
||
# 调整布局
|
||
plt.tight_layout(rect=[0, 0, 1, 0.95]) # 为总标题留出空间
|
||
# 显示图像
|
||
plt.show()
|
||
|
||
# 输出路径和文件名
|
||
picture_file_name = f'{year}年_周日_潮汐波振幅时空变化.png'
|
||
plt.savefig(os.path.join(output_dir, picture_file_name))
|
||
|
||
|
||
# -----------------------------------------得到2日行星波拟合参数-------------------------------------------
|
||
# 定义新的拟合函数
|
||
def fit_wind_speed_v1(data, window=5):
|
||
# 动态调整最小有效数据点数量：参数数量 * 6
|
||
num_params = 3 # 模型参数数量
|
||
# 转换为 pandas DataFrame 并设置时间为索引
|
||
data["date"] = pd.to_datetime(data["date"])
|
||
data.set_index("date", inplace=True)
|
||
|
||
# 获取唯一日期
|
||
unique_dates = np.unique(data.index.date)
|
||
|
||
# 存储拟合参数
|
||
params = []
|
||
|
||
for date in unique_dates:
|
||
start_date = pd.to_datetime(date) - pd.Timedelta(days=(window - 1) / 2)
|
||
end_date = pd.to_datetime(date) + pd.Timedelta(days=(window - 1) / 2)
|
||
|
||
# 获取这几天的数据
|
||
window_data = data[start_date:end_date]
|
||
|
||
t = np.arange(len(window_data)) + 1 # 时间点
|
||
y = window_data['wind_speed'].values # 风速数据
|
||
|
||
# 将 y 转换为浮点类型，以防有非数值类型数据
|
||
y = pd.to_numeric(y, errors='coerce') # 将非数值转换为 NaN
|
||
|
||
# 去掉包含nan的那一列，筛选有效数据
|
||
valid_mask = ~np.isnan(y) # 创建布尔掩码，标记非 NaN 的位置
|
||
t = t[valid_mask]
|
||
y = y[valid_mask]
|
||
|
||
# 检查筛选后的数据是否足够
|
||
if len(y) < num_params * 6: # 确保数据足够
|
||
# print(f"Not enough valid data after filtering for date: {date}")
|
||
params.append([None] * 3) # 如果数据不足，记录None
|
||
continue
|
||
|
||
# 设置初始参数
|
||
initial_guess = [0, 1, 0] # 新模型的初始猜测参数
|
||
bounds = (
|
||
[-np.inf, 0, -np.inf], # 下界
|
||
[np.inf, np.inf, np.inf], # 上界
|
||
)
|
||
|
||
# 拟合模型
|
||
try:
|
||
popt, _ = curve_fit(planetary_model_2day, t, y, p0=initial_guess, bounds=bounds)
|
||
params.append(popt) # 存储拟合参数
|
||
except RuntimeError:
|
||
params.append([None] * 3) # 如果拟合失败，记录None
|
||
|
||
# 转换拟合参数为 DataFrame
|
||
params_df = pd.DataFrame(
|
||
params,
|
||
columns=["v0", "a", "b"],
|
||
index=unique_dates,
|
||
)
|
||
return params_df
|
||
|
||
# 1-拟合纬向风------------------------------------------
|
||
u_result_dict = {}
|
||
|
||
for height in heights:
|
||
print(f"Processing height: {height} m")
|
||
# 获取当前高度的纬向风数据
|
||
height_data = uwind_data[uwind_data["height"] == height]
|
||
|
||
# 转换为长格式（适配拟合函数）
|
||
long_data = pd.melt(
|
||
height_data,
|
||
id_vars=["time", "height"],
|
||
var_name="date",
|
||
value_name="wind_speed",
|
||
)
|
||
long_data["date"] = pd.to_datetime(long_data["date"], format="%Y%m%d")
|
||
|
||
# 调用修改后的拟合函数 fit_wind_speed_v2
|
||
params_df = fit_wind_speed_v1(long_data)
|
||
u_result_dict[height] = params_df
|
||
|
||
# 汇总结果为 DataFrame
|
||
u_final_result = pd.concat(u_result_dict, names=["height", "date"])
|
||
|
||
|
||
# 2-拟合经向风------------------------------------------
|
||
v_result_dict = {}
|
||
|
||
for height in heights:
|
||
print(f"Processing height: {height} m")
|
||
height_data = vwind_data[vwind_data["height"] == height]
|
||
|
||
long_data = pd.melt(
|
||
height_data,
|
||
id_vars=["time", "height"],
|
||
var_name="date",
|
||
value_name="wind_speed",
|
||
)
|
||
long_data["date"] = pd.to_datetime(long_data["date"], format="%Y%m%d")
|
||
|
||
|
||
# 调用修改后的拟合函数 fit_wind_speed_v2
|
||
params_df = fit_wind_speed_v1(long_data)
|
||
v_result_dict[height] = params_df
|
||
|
||
# 汇总结果为 DataFrame
|
||
v_final_result = pd.concat(v_result_dict, names=["height", "date"])
|
||
|
||
|
||
# 定义振幅参数
|
||
amplitude_param = 'a'
|
||
# 提取纬向风数据，重置索引为平坦结构
|
||
u_heatmap_data = u_final_result[amplitude_param].unstack(level=0) # 将高度作为列，日期作为行
|
||
u_heatmap_data.index = pd.to_datetime(u_heatmap_data.index) # 确保日期为时间格式
|
||
u_heatmap_data = u_heatmap_data.iloc[:, ::-1] # 按高度降序排列
|
||
|
||
# 提取经向风数据，重置索引为平坦结构
|
||
v_heatmap_data = v_final_result[amplitude_param].unstack(level=0) # 将高度作为列，日期作为行
|
||
v_heatmap_data.index = pd.to_datetime(v_heatmap_data.index) # 确保日期为时间格式
|
||
v_heatmap_data = v_heatmap_data.iloc[:, ::-1] # 按高度降序排列
|
||
|
||
# 创建画布和子图
|
||
fig, axes = plt.subplots(1, 2, figsize=(18, 8)) # 1行2列的布局
|
||
fig.suptitle('2022年2日_行星波振幅时空变化', fontsize=20) # 添加总标题
|
||
|
||
# 绘制第一幅热力图（纬向风）
|
||
sns.heatmap(
|
||
u_heatmap_data.T, # 转置，横轴为时间，纵轴为高度
|
||
cmap="viridis", # 配色方案
|
||
vmax=100, # 设置颜色范围的最大值
|
||
cbar_kws={'label': '振幅 (m/s)'}, # 设置颜色条标签
|
||
xticklabels=False, # 禁用默认的横轴刻度标签
|
||
#yticklabels=u_heatmap_data.columns.astype(str), # 显示纵轴的高度
|
||
yticklabels=[f'{height // 1000}' for height in u_heatmap_data.columns], # 转换高度为千米并格式化
|
||
ax=axes[0] # 指定子图
|
||
)
|
||
axes[0].set_title('纬向风振幅变化', fontsize=14)
|
||
axes[0].set_xlabel('日期', fontsize=12)
|
||
axes[0].set_ylabel('高度 (km)', fontsize=12)
|
||
|
||
|
||
# 绘制第二幅热力图（经向风）
|
||
sns.heatmap(
|
||
v_heatmap_data.T, # 转置，横轴为时间，纵轴为高度
|
||
cmap="viridis", # 配色方案
|
||
vmax=100, # 设置颜色范围的最大值
|
||
cbar_kws={'label': '振幅 (m/s)'}, # 设置颜色条标签
|
||
xticklabels=False, # 禁用默认的横轴刻度标签
|
||
#yticklabels=v_heatmap_data.columns.astype(str), # 显示纵轴的高度
|
||
yticklabels=[f'{height // 1000}' for height in v_heatmap_data.columns], # 转换高度为千米并格式化
|
||
ax=axes[1] # 指定子图
|
||
)
|
||
axes[1].set_title('经向风振幅变化', fontsize=14)
|
||
axes[1].set_xlabel('日期', fontsize=12)
|
||
axes[1].set_ylabel('高度 (km)', fontsize=12)
|
||
|
||
# 调整布局
|
||
plt.tight_layout(rect=[0, 0, 1, 0.95]) # 为总标题留出空间
|
||
plt.show()
|
||
|
||
# 输出路径和文件名
|
||
picture_file_name = f'{year}年_2日_行星波振幅时空变化.png'
|
||
plt.savefig(os.path.join(output_dir, picture_file_name))
|
||
# -----------------------------------------得到5日行星波拟合参数-------------------------------------------
|
||
# 定义新的拟合函数
|
||
def fit_wind_speed_v2(data, window=15):
|
||
# 动态调整最小有效数据点数量：参数数量 * 6
|
||
num_params = 3 # 模型参数数量
|
||
# 转换为 pandas DataFrame 并设置时间为索引
|
||
data["date"] = pd.to_datetime(data["date"])
|
||
data.set_index("date", inplace=True)
|
||
|
||
# 获取唯一日期
|
||
unique_dates = np.unique(data.index.date)
|
||
|
||
# 存储拟合参数
|
||
params = []
|
||
|
||
for date in unique_dates:
|
||
start_date = pd.to_datetime(date) - pd.Timedelta(days=(window - 1) / 2)
|
||
end_date = pd.to_datetime(date) + pd.Timedelta(days=(window - 1) / 2)
|
||
|
||
# 获取这几天的数据
|
||
window_data = data[start_date:end_date]
|
||
|
||
t = np.arange(len(window_data)) + 1 # 时间点
|
||
y = window_data['wind_speed'].values # 风速数据
|
||
|
||
# 将 y 转换为浮点类型，以防有非数值类型数据
|
||
y = pd.to_numeric(y, errors='coerce') # 将非数值转换为 NaN
|
||
|
||
# 去掉包含nan的那一列，筛选有效数据
|
||
valid_mask = ~np.isnan(y) # 创建布尔掩码，标记非 NaN 的位置
|
||
t = t[valid_mask]
|
||
y = y[valid_mask]
|
||
|
||
# 检查筛选后的数据是否足够
|
||
if len(y) < num_params * 6: # 确保数据足够
|
||
# print(f"Not enough valid data after filtering for date: {date}")
|
||
params.append([None] * 3) # 如果数据不足，记录None
|
||
continue
|
||
|
||
# 设置初始参数
|
||
initial_guess = [0, 1, 0] # 新模型的初始猜测参数
|
||
bounds = (
|
||
[-np.inf, 0, -np.inf], # 下界
|
||
[np.inf, np.inf, np.inf], # 上界
|
||
)
|
||
|
||
# 拟合模型
|
||
try:
|
||
popt, _ = curve_fit(planetary_model_5day, t, y, p0=initial_guess, bounds=bounds)
|
||
params.append(popt) # 存储拟合参数
|
||
except RuntimeError:
|
||
params.append([None] * 3) # 如果拟合失败，记录None
|
||
|
||
# 转换拟合参数为 DataFrame
|
||
params_df = pd.DataFrame(
|
||
params,
|
||
columns=["v0", "a", "b"],
|
||
index=unique_dates,
|
||
)
|
||
return params_df
|
||
|
||
# 1-拟合纬向风------------------------------------------
|
||
u_result_dict = {}
|
||
|
||
for height in heights:
|
||
print(f"Processing height: {height} m")
|
||
# 获取当前高度的纬向风数据
|
||
height_data = uwind_data[uwind_data["height"] == height]
|
||
|
||
# 转换为长格式（适配拟合函数）
|
||
long_data = pd.melt(
|
||
height_data,
|
||
id_vars=["time", "height"],
|
||
var_name="date",
|
||
value_name="wind_speed",
|
||
)
|
||
long_data["date"] = pd.to_datetime(long_data["date"], format="%Y%m%d")
|
||
|
||
# 调用修改后的拟合函数 fit_wind_speed_v2
|
||
params_df = fit_wind_speed_v2(long_data)
|
||
u_result_dict[height] = params_df
|
||
|
||
# 汇总结果为 DataFrame
|
||
u_final_result = pd.concat(u_result_dict, names=["height", "date"])
|
||
'''
|
||
# 将结果保存为 TXT 文档
|
||
file_name = f'{year}年_纬向风_5日行星波振幅参数.txt'
|
||
u_final_result.to_csv(
|
||
os.path.join(output_dir, file_name),
|
||
sep='\t',
|
||
index=True,
|
||
float_format='%.2f', # 控制数值列保留 2 位小数
|
||
header=True)
|
||
'''
|
||
|
||
# 2-拟合经向风------------------------------------------
|
||
v_result_dict = {}
|
||
|
||
for height in heights:
|
||
print(f"Processing height: {height} m")
|
||
height_data = vwind_data[vwind_data["height"] == height]
|
||
|
||
long_data = pd.melt(
|
||
height_data,
|
||
id_vars=["time", "height"],
|
||
var_name="date",
|
||
value_name="wind_speed",
|
||
)
|
||
long_data["date"] = pd.to_datetime(long_data["date"], format="%Y%m%d")
|
||
|
||
|
||
# 调用修改后的拟合函数 fit_wind_speed_v2
|
||
params_df = fit_wind_speed_v2(long_data)
|
||
v_result_dict[height] = params_df
|
||
|
||
# 汇总结果为 DataFrame
|
||
v_final_result = pd.concat(v_result_dict, names=["height", "date"])
|
||
'''
|
||
# 保存参数
|
||
file_name = f'{year}年_经向风_5日行星波振幅参数.txt'
|
||
v_final_result.to_csv(
|
||
os.path.join(output_dir, file_name),
|
||
sep='\t',
|
||
index=True,
|
||
float_format='%.2f', # 控制数值列保留 2 位小数
|
||
header=True)
|
||
'''
|
||
|
||
# 定义振幅参数
|
||
amplitude_param = 'a'
|
||
# 提取纬向风数据，重置索引为平坦结构
|
||
u_heatmap_data = u_final_result[amplitude_param].unstack(level=0) # 将高度作为列，日期作为行
|
||
u_heatmap_data.index = pd.to_datetime(u_heatmap_data.index) # 确保日期为时间格式
|
||
u_heatmap_data = u_heatmap_data.iloc[:, ::-1] # 按高度降序排列
|
||
|
||
# 提取经向风数据，重置索引为平坦结构
|
||
v_heatmap_data = v_final_result[amplitude_param].unstack(level=0) # 将高度作为列，日期作为行
|
||
v_heatmap_data.index = pd.to_datetime(v_heatmap_data.index) # 确保日期为时间格式
|
||
v_heatmap_data = v_heatmap_data.iloc[:, ::-1] # 按高度降序排列
|
||
|
||
# 创建画布和子图
|
||
fig, axes = plt.subplots(1, 2, figsize=(18, 8)) # 1行2列的布局
|
||
fig.suptitle('2022年5日_行星波振幅时空变化', fontsize=20) # 添加总标题
|
||
|
||
# 绘制第一幅热力图（纬向风）
|
||
sns.heatmap(
|
||
u_heatmap_data.T, # 转置，横轴为时间，纵轴为高度
|
||
cmap="viridis", # 配色方案
|
||
vmax=100, # 设置颜色范围的最大值
|
||
cbar_kws={'label': '振幅 (m/s)'}, # 设置颜色条标签
|
||
xticklabels=False, # 禁用默认的横轴刻度标签
|
||
#yticklabels=u_heatmap_data.columns.astype(str), # 显示纵轴的高度
|
||
yticklabels=[f'{height // 1000}' for height in u_heatmap_data.columns], # 转换高度为千米并格式化
|
||
ax=axes[0] # 指定子图
|
||
)
|
||
axes[0].set_title('纬向风振幅变化', fontsize=14)
|
||
axes[0].set_xlabel('日期', fontsize=12)
|
||
axes[0].set_ylabel('高度 (km)', fontsize=12)
|
||
|
||
|
||
# 绘制第二幅热力图（经向风）
|
||
sns.heatmap(
|
||
v_heatmap_data.T, # 转置，横轴为时间，纵轴为高度
|
||
cmap="viridis", # 配色方案
|
||
vmax=100, # 设置颜色范围的最大值
|
||
cbar_kws={'label': '振幅 (m/s)'}, # 设置颜色条标签
|
||
xticklabels=False, # 禁用默认的横轴刻度标签
|
||
#yticklabels=v_heatmap_data.columns.astype(str), # 显示纵轴的高度
|
||
yticklabels=[f'{height // 1000}' for height in v_heatmap_data.columns], # 转换高度为千米并格式化
|
||
ax=axes[1] # 指定子图
|
||
)
|
||
axes[1].set_title('经向风振幅变化', fontsize=14)
|
||
axes[1].set_xlabel('日期', fontsize=12)
|
||
axes[1].set_ylabel('高度 (km)', fontsize=12)
|
||
|
||
# 调整布局
|
||
plt.tight_layout(rect=[0, 0, 1, 0.95]) # 为总标题留出空间
|
||
plt.show()
|
||
|
||
# 输出路径和文件名
|
||
picture_file_name = f'{year}年_5日_行星波振幅时空变化.png'
|
||
plt.savefig(os.path.join(output_dir, picture_file_name))
|
||
|
||
|
||
# -----------------------------------------得到10日行星波拟合参数-------------------------------------------
|
||
|
||
# 定义新的拟合函数
|
||
def fit_wind_speed_v3(data, window=29):
|
||
|
||
# 动态调整最小有效数据点数量：参数数量 * 6
|
||
num_params = 3 # 模型参数数量
|
||
# 转换为 pandas DataFrame 并设置时间为索引
|
||
data["date"] = pd.to_datetime(data["date"])
|
||
data.set_index("date", inplace=True)
|
||
|
||
# 获取唯一日期
|
||
unique_dates = np.unique(data.index.date)
|
||
|
||
# 存储拟合参数
|
||
params = []
|
||
|
||
for date in unique_dates:
|
||
start_date = pd.to_datetime(date) - pd.Timedelta(days=(window - 1) / 2)
|
||
end_date = pd.to_datetime(date) + pd.Timedelta(days=(window - 1) / 2)
|
||
|
||
# 获取这几天的数据
|
||
window_data = data[start_date:end_date]
|
||
|
||
t = np.arange(len(window_data)) + 1 # 时间点
|
||
y = window_data['wind_speed'].values # 风速数据
|
||
|
||
# 将 y 转换为浮点类型，以防有非数值类型数据
|
||
y = pd.to_numeric(y, errors='coerce') # 将非数值转换为 NaN
|
||
|
||
# 去掉包含nan的那一列，筛选有效数据
|
||
valid_mask = ~np.isnan(y) # 创建布尔掩码，标记非 NaN 的位置
|
||
t = t[valid_mask]
|
||
y = y[valid_mask]
|
||
|
||
# 检查筛选后的数据是否足够
|
||
if len(y) < num_params * 6: # 确保数据足够
|
||
# print(f"Not enough valid data after filtering for date: {date}")
|
||
params.append([None] * 3) # 如果数据不足，记录None
|
||
continue
|
||
|
||
# 设置初始参数
|
||
initial_guess = [0, 1, 0] # 新模型的初始猜测参数
|
||
bounds = (
|
||
[-np.inf, 0, -np.inf], # 下界
|
||
[np.inf, np.inf, np.inf], # 上界
|
||
)
|
||
|
||
# 拟合模型
|
||
try:
|
||
popt, _ = curve_fit(planetary_model_10day, t, y, p0=initial_guess, bounds=bounds)
|
||
params.append(popt) # 存储拟合参数
|
||
except RuntimeError:
|
||
params.append([None] * 3) # 如果拟合失败，记录None
|
||
|
||
# 转换拟合参数为 DataFrame
|
||
params_df = pd.DataFrame(
|
||
params,
|
||
columns=["v0", "a", "b"],
|
||
index=unique_dates,
|
||
)
|
||
return params_df
|
||
|
||
# 1-拟合纬向风------------------------------------------
|
||
|
||
u_result_dict = {}
|
||
|
||
for height in heights:
|
||
print(f"Processing height: {height} m")
|
||
# 获取当前高度的纬向风数据
|
||
height_data = uwind_data[uwind_data["height"] == height]
|
||
|
||
# 转换为长格式（适配拟合函数）
|
||
long_data = pd.melt(
|
||
height_data,
|
||
id_vars=["time", "height"],
|
||
var_name="date",
|
||
value_name="wind_speed",
|
||
)
|
||
long_data["date"] = pd.to_datetime(long_data["date"], format="%Y%m%d")
|
||
|
||
# 调用修改后的拟合函数 fit_wind_speed_v3
|
||
params_df = fit_wind_speed_v3(long_data)
|
||
u_result_dict[height] = params_df
|
||
|
||
# 汇总结果为 DataFrame
|
||
u_final_result = pd.concat(u_result_dict, names=["height", "date"])
|
||
|
||
# 2-拟合经向风------------------------------------------
|
||
|
||
v_result_dict = {}
|
||
|
||
for height in heights:
|
||
print(f"Processing height: {height} m")
|
||
height_data = vwind_data[vwind_data["height"] == height]
|
||
|
||
long_data = pd.melt(
|
||
height_data,
|
||
id_vars=["time", "height"],
|
||
var_name="date",
|
||
value_name="wind_speed",
|
||
)
|
||
long_data["date"] = pd.to_datetime(long_data["date"], format="%Y%m%d")
|
||
|
||
|
||
# 调用修改后的拟合函数 fit_wind_speed_v3
|
||
params_df = fit_wind_speed_v3(long_data)
|
||
v_result_dict[height] = params_df
|
||
|
||
# 汇总结果为 DataFrame
|
||
v_final_result = pd.concat(v_result_dict, names=["height", "date"])
|
||
|
||
|
||
# 定义振幅参数
|
||
amplitude_param = 'a'
|
||
|
||
# 提取纬向风数据，重置索引为平坦结构
|
||
u_heatmap_data = u_final_result[amplitude_param].unstack(level=0) # 将高度作为列，日期作为行
|
||
u_heatmap_data.index = pd.to_datetime(u_heatmap_data.index) # 确保日期为时间格式
|
||
u_heatmap_data = u_heatmap_data.iloc[:, ::-1] # 按高度降序排列
|
||
|
||
# 提取经向风数据，重置索引为平坦结构
|
||
v_heatmap_data = v_final_result[amplitude_param].unstack(level=0) # 将高度作为列，日期作为行
|
||
v_heatmap_data.index = pd.to_datetime(v_heatmap_data.index) # 确保日期为时间格式
|
||
v_heatmap_data = v_heatmap_data.iloc[:, ::-1] # 按高度降序排列
|
||
|
||
# 创建画布和子图
|
||
fig, axes = plt.subplots(1, 2, figsize=(18, 8)) # 1行2列的布局
|
||
fig.suptitle('2022年10日_行星波振幅时空变化', fontsize=20) # 添加总标题
|
||
|
||
# 绘制第一幅热力图（纬向风）
|
||
sns.heatmap(
|
||
u_heatmap_data.T, # 转置，横轴为时间，纵轴为高度
|
||
cmap="viridis", # 配色方案
|
||
vmax=100, # 设置颜色范围的最大值
|
||
cbar_kws={'label': '振幅 (m/s)'}, # 设置颜色条标签
|
||
xticklabels=False, # 禁用默认的横轴刻度标签
|
||
#yticklabels=u_heatmap_data.columns.astype(str), # 显示纵轴的高度
|
||
yticklabels=[f'{height // 1000}' for height in u_heatmap_data.columns], # 转换高度为千米并格式化
|
||
ax=axes[0] # 指定子图
|
||
)
|
||
axes[0].set_title('纬向风振幅变化', fontsize=14)
|
||
axes[0].set_xlabel('日期', fontsize=12)
|
||
axes[0].set_ylabel('高度 (m)', fontsize=12)
|
||
|
||
# 绘制第二幅热力图（经向风）
|
||
sns.heatmap(
|
||
v_heatmap_data.T, # 转置，横轴为时间，纵轴为高度
|
||
cmap="viridis", # 配色方案
|
||
vmax=100, # 设置颜色范围的最大值
|
||
cbar_kws={'label': '振幅 (m/s)'}, # 设置颜色条标签
|
||
xticklabels=False, # 禁用默认的横轴刻度标签
|
||
#yticklabels=v_heatmap_data.columns.astype(str), # 显示纵轴的高度
|
||
yticklabels=[f'{height // 1000}' for height in u_heatmap_data.columns], # 转换高度为千米并格式化
|
||
ax=axes[1] # 指定子图
|
||
)
|
||
axes[1].set_title('经向风振幅变化', fontsize=14)
|
||
axes[1].set_xlabel('日期', fontsize=12)
|
||
axes[1].set_ylabel('高度 (m)', fontsize=12)
|
||
|
||
# 调整布局
|
||
plt.tight_layout(rect=[0, 0, 1, 0.95]) # 为总标题留出空间
|
||
# 显示或保存图像
|
||
plt.show()
|
||
|
||
# 输出路径和文件名
|
||
picture_file_name = f'{year}年_10日_行星波振幅时空变化.png'
|
||
plt.savefig(os.path.join(output_dir, picture_file_name))
|
||
|
||
|
||
# -----------------------------------------得到16日行星波拟合参数-------------------------------------------
|
||
|
||
# 定义新的拟合函数
|
||
def fit_wind_speed_v4(data, window=49):
|
||
# 动态调整最小有效数据点数量：参数数量 * 6
|
||
num_params = 3 # 模型参数数量
|
||
# 转换为 pandas DataFrame 并设置时间为索引
|
||
data["date"] = pd.to_datetime(data["date"])
|
||
data.set_index("date", inplace=True)
|
||
|
||
# 获取唯一日期
|
||
unique_dates = np.unique(data.index.date)
|
||
|
||
# 存储拟合参数
|
||
params = []
|
||
|
||
for date in unique_dates:
|
||
start_date = pd.to_datetime(date) - pd.Timedelta(days=(window - 1) / 2)
|
||
end_date = pd.to_datetime(date) + pd.Timedelta(days=(window - 1) / 2)
|
||
|
||
# 获取这几天的数据
|
||
window_data = data[start_date:end_date]
|
||
|
||
t = np.arange(len(window_data)) + 1 # 时间点
|
||
y = window_data['wind_speed'].values # 风速数据
|
||
|
||
# 将 y 转换为浮点类型，以防有非数值类型数据
|
||
y = pd.to_numeric(y, errors='coerce') # 将非数值转换为 NaN
|
||
|
||
# 去掉包含nan的那一列，筛选有效数据
|
||
valid_mask = ~np.isnan(y) # 创建布尔掩码，标记非 NaN 的位置
|
||
t = t[valid_mask]
|
||
y = y[valid_mask]
|
||
|
||
# 检查筛选后的数据是否足够
|
||
if len(y) < num_params * 6: # 确保数据足够
|
||
# print(f"Not enough valid data after filtering for date: {date}")
|
||
params.append([None] * 3) # 如果数据不足，记录None
|
||
continue
|
||
|
||
# 设置初始参数
|
||
initial_guess = [0, 1, 0] # 新模型的初始猜测参数
|
||
bounds = (
|
||
[-np.inf, 0, -np.inf], # 下界
|
||
[np.inf, np.inf, np.inf], # 上界
|
||
)
|
||
|
||
# 拟合模型
|
||
try:
|
||
popt, _ = curve_fit(planetary_model_16day, t, y, p0=initial_guess, bounds=bounds)
|
||
params.append(popt) # 存储拟合参数
|
||
except RuntimeError:
|
||
params.append([None] * 3) # 如果拟合失败，记录None
|
||
|
||
# 转换拟合参数为 DataFrame
|
||
params_df = pd.DataFrame(
|
||
params,
|
||
columns=["v0", "a", "b"],
|
||
index=unique_dates,
|
||
)
|
||
return params_df
|
||
|
||
# 1-拟合纬向风------------------------------------------
|
||
|
||
u_result_dict = {}
|
||
|
||
for height in heights:
|
||
print(f"Processing height: {height} m")
|
||
# 获取当前高度的纬向风数据
|
||
height_data = uwind_data[uwind_data["height"] == height]
|
||
|
||
# 转换为长格式（适配拟合函数）
|
||
long_data = pd.melt(
|
||
height_data,
|
||
id_vars=["time", "height"],
|
||
var_name="date",
|
||
value_name="wind_speed",
|
||
)
|
||
long_data["date"] = pd.to_datetime(long_data["date"], format="%Y%m%d")
|
||
|
||
# 调用修改后的拟合函数 fit_wind_speed_v4
|
||
params_df = fit_wind_speed_v4(long_data)
|
||
u_result_dict[height] = params_df
|
||
|
||
# 汇总结果为 DataFrame
|
||
u_final_result = pd.concat(u_result_dict, names=["height", "date"])
|
||
|
||
# 2-拟合经向风------------------------------------------
|
||
|
||
v_result_dict = {}
|
||
|
||
for height in heights:
|
||
print(f"Processing height: {height} m")
|
||
height_data = vwind_data[vwind_data["height"] == height]
|
||
|
||
long_data = pd.melt(
|
||
height_data,
|
||
id_vars=["time", "height"],
|
||
var_name="date",
|
||
value_name="wind_speed",
|
||
)
|
||
long_data["date"] = pd.to_datetime(long_data["date"], format="%Y%m%d")
|
||
|
||
|
||
# 调用修改后的拟合函数 fit_wind_speed_v4
|
||
params_df = fit_wind_speed_v4(long_data)
|
||
v_result_dict[height] = params_df
|
||
|
||
# 汇总结果为 DataFrame
|
||
v_final_result = pd.concat(v_result_dict, names=["height", "date"])
|
||
|
||
|
||
# 定义振幅参数
|
||
amplitude_param = 'a'
|
||
|
||
# 提取纬向风数据，重置索引为平坦结构
|
||
u_heatmap_data = u_final_result[amplitude_param].unstack(level=0) # 将高度作为列，日期作为行
|
||
u_heatmap_data.index = pd.to_datetime(u_heatmap_data.index) # 确保日期为时间格式
|
||
u_heatmap_data = u_heatmap_data.iloc[:, ::-1] # 按高度降序排列
|
||
|
||
# 提取经向风数据，重置索引为平坦结构
|
||
v_heatmap_data = v_final_result[amplitude_param].unstack(level=0) # 将高度作为列，日期作为行
|
||
v_heatmap_data.index = pd.to_datetime(v_heatmap_data.index) # 确保日期为时间格式
|
||
v_heatmap_data = v_heatmap_data.iloc[:, ::-1] # 按高度降序排列
|
||
|
||
# 创建画布和子图
|
||
fig, axes = plt.subplots(1, 2, figsize=(18, 8)) # 1行2列的布局
|
||
fig.suptitle('2022年16日_行星波振幅时空变化', fontsize=20) # 添加总标题
|
||
|
||
# 绘制第一幅热力图（纬向风）
|
||
sns.heatmap(
|
||
u_heatmap_data.T, # 转置，横轴为时间，纵轴为高度
|
||
cmap="viridis", # 配色方案
|
||
vmax=100, # 设置颜色范围的最大值
|
||
cbar_kws={'label': '振幅 (m/s)'}, # 设置颜色条标签
|
||
xticklabels=False, # 禁用默认的横轴刻度标签
|
||
#yticklabels=u_heatmap_data.columns.astype(str), # 显示纵轴的高度
|
||
yticklabels=[f'{height // 1000}' for height in u_heatmap_data.columns], # 转换高度为千米并格式化
|
||
|
||
ax=axes[0] # 指定子图
|
||
)
|
||
axes[0].set_title('纬向风振幅变化', fontsize=14)
|
||
axes[0].set_xlabel('日期', fontsize=12)
|
||
axes[0].set_ylabel('高度 (m)', fontsize=12)
|
||
|
||
# 绘制第二幅热力图（经向风）
|
||
sns.heatmap(
|
||
v_heatmap_data.T, # 转置，横轴为时间，纵轴为高度
|
||
cmap="viridis", # 配色方案
|
||
vmax=100, # 设置颜色范围的最大值
|
||
cbar_kws={'label': '振幅 (m/s)'}, # 设置颜色条标签
|
||
xticklabels=False, # 禁用默认的横轴刻度标签
|
||
#yticklabels=v_heatmap_data.columns.astype(str), # 显示纵轴的高度
|
||
yticklabels=[f'{height // 1000}' for height in v_heatmap_data.columns], # 转换高度为千米并格式化
|
||
|
||
ax=axes[1] # 指定子图
|
||
)
|
||
axes[1].set_title('经向风振幅变化', fontsize=14)
|
||
axes[1].set_xlabel('日期', fontsize=12)
|
||
axes[1].set_ylabel('高度 (m)', fontsize=12)
|
||
|
||
# 调整布局
|
||
plt.tight_layout(rect=[0, 0, 1, 0.95]) # 为总标题留出空间
|
||
|
||
# 显示或保存图像
|
||
plt.show()
|
||
|
||
# 输出路径和文件名
|
||
picture_file_name = f'{year}年_16日_行星波振幅时空变化.png'
|
||
plt.savefig(os.path.join(output_dir, picture_file_name))
|
||
|
||
|
||
|