93 lines
3.6 KiB
Python
93 lines
3.6 KiB
Python
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 # 正常显示负号
|
|
# 读取一年的数据文件
|
|
# 设置初始参数
|
|
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 cosmic_planetw_plot_perday(
|
|
df: pd.DataFrame,
|
|
T=16,
|
|
start_day=0
|
|
):
|
|
|
|
def u_func(x, *params, t):
|
|
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)
|
|
)
|
|
|
|
# 设置最小数据量的阈值
|
|
min_data_points = 36
|
|
|
|
# 进行多个时间窗口的拟合
|
|
# 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)}")
|
|
plt.figure(figsize=(10, 6))
|
|
# write nothing but a title
|
|
plt.title(f'时间窗口：{start_day} 到 {end_day}')
|
|
plt.grid(True)
|
|
# add a label, says unable to fit
|
|
plt.text(0.5, 0.5, "数据量不足，无法拟合", ha='center', va='center',
|
|
transform=plt.gca().transAxes, fontsize=20) # 跳过当前时间窗口，继续下一个窗口
|
|
|
|
# 提取时间、经度、温度数据
|
|
t = np.array(df_8['Time']) # 时间
|
|
x = np.array(df_8['Longitude_Radians']) # 经度弧度制
|
|
temperature = np.array(df_8['Temperature']) # 温度，因变量
|
|
|
|
# `tempeture` 为因变量，`x` 为自变量，`t` 为参数
|
|
# tempeture is possible to contain NaNs, so we need to drop them
|
|
mask = ~np.isnan(temperature)
|
|
x = x[mask]
|
|
temperature = temperature[mask]
|
|
t = t[mask]
|
|
|
|
# 用T进行拟合
|
|
popt, pcov = curve_fit(lambda x, *params: u_func(x, *params, t=t),
|
|
x, temperature, p0=initial_guess, bounds=bounds, maxfev=50000)
|
|
|
|
# 绘制拟合曲线
|
|
t_fixed = (start_day + end_day) / 2 # 窗口的中间时间
|
|
x_fit = np.linspace(min(x), max(x), 100) # 生成用于绘制拟合曲线的x值
|
|
y_fit = u_func(x_fit, *popt, t=t_fixed) # 计算拟合曲线的y值
|
|
|
|
plt.figure(figsize=(10, 6))
|
|
plt.plot(x_fit, y_fit, label='拟合曲线', color='red', linewidth=2)
|
|
plt.xlabel('经度（弧度制）')
|
|
plt.ylabel('温度')
|
|
plt.title(f'时间窗口：{start_day} 到 {end_day}')
|
|
plt.legend()
|
|
plt.grid(True)
|