From ed22e7d3a338093c04a231e0b5278b062588e200 Mon Sep 17 00:00:00 2001 From: Dustella Date: Wed, 19 Feb 2025 12:28:03 +0800 Subject: [PATCH] noted orignal filename --- cosmic/multiple.py | 5 ++- cosmic/single.py | 13 +++++--- cosmic/temp_render.py | 50 ++++++++++++++++-------------- saber/archive/gravity_plot.py | 1 + saber/archive/gravity_process.py | 1 + saber/archive/zjy_wave_fit_plot.py | 2 ++ tidi/plot.py | 49 ++++++++++++++++------------- tidi/staged/plot.py | 2 +- tidi/staged/process.py | 2 ++ 9 files changed, 73 insertions(+), 52 deletions(-) diff --git a/cosmic/multiple.py b/cosmic/multiple.py index a4306be..9a56d00 100644 --- a/cosmic/multiple.py +++ b/cosmic/multiple.py @@ -1,3 +1,6 @@ +# 这里全是重力波相关的 +# 原始文件名 ：cosmic重力波多天.py + import os import numpy as np import pandas as pd @@ -269,7 +272,7 @@ def process_single_file(base_folder_path, i): # 主循环，处理1到3个文件 -base_folder_path = r"./cosmic/data/2008" +base_folder_path = r"./data/cosmic/2008" all_mean_ktemp_Nz = [] all_mean_ktemp_Ptz = [] for file_index in range(1, 365): diff --git a/cosmic/single.py b/cosmic/single.py index 4f04767..fa54cd8 100644 --- a/cosmic/single.py +++ b/cosmic/single.py @@ -1,3 +1,5 @@ +# 重力波 一天cosmic重力波 + import os import time import numpy as np @@ -360,7 +362,8 @@ def plot_results_mean_ktemp_Nz(mean_ktemp_Nz, heights): plt.ylabel('H(km)') # plt.gca().invert_yaxis() # 使高度坐标轴从上到下递增，符合常规习惯 # plt.show() - + + def plot_results_mean_ktemp_Ptz(mean_ktemp_Ptz, heights): """ @@ -377,10 +380,10 @@ def plot_results_mean_ktemp_Ptz(mean_ktemp_Ptz, heights): class SingleCosmicWavePlot: - + def __init__(self, year, day): BASE_PATH_COSMIC = DATA_BASEPATH.cosmic - + base_folder_path = f"./{BASE_PATH_COSMIC}/{year}" day_num = day # 模块1调用 @@ -407,9 +410,9 @@ class SingleCosmicWavePlot: self.mean_ktemp_Ptz = mean_ktemp_Ptz self.heights = heights # 调用绘图模块函数进行绘图 - + def plot_results_mean_ktemp_Nz(self): plot_results_mean_ktemp_Nz(self.mean_ktemp_Nz, self.heights) - + def plot_results_mean_ktemp_Ptz(self): plot_results_mean_ktemp_Ptz(self.mean_ktemp_Ptz, self.heights) diff --git a/cosmic/temp_render.py b/cosmic/temp_render.py index 7d1e242..bc4a75a 100644 --- a/cosmic/temp_render.py +++ b/cosmic/temp_render.py @@ -1,4 +1,5 @@ -#此代码是对数据处理后的txt数据进行行星波参数提取绘图 +# 原始文件 cosmic行星波参数全年逐日绘图 +# 此代码是对数据处理后的txt数据进行行星波参数提取绘图 import pandas as pd import numpy as np @@ -16,36 +17,35 @@ matplotlib.rcParams['axes.unicode_minus'] = False # 正常显示负号 # 读取一年的数据文件 -def temp_render(path:str = f"{DATA_BASEPATH.cosmic}/cosmic.txt", T = 16): +def temp_render(path: str = f"{DATA_BASEPATH.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) + 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 参数 + 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]) # 上界 - - - + [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 = [] @@ -54,7 +54,7 @@ def temp_render(path:str = f"{DATA_BASEPATH.cosmic}/cosmic.txt", T = 16): min_data_points = 36 # 进行多个时间窗口的拟合 - #T应该取5、10、16 + # 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] @@ -77,17 +77,19 @@ def temp_render(path:str = f"{DATA_BASEPATH.cosmic}/cosmic.txt", T = 16): temperature = np.array(df_8['Temperature']) # 温度，因变量 # 用T进行拟合 - popt, pcov = curve_fit(u_func, x, temperature, p0=initial_guess, bounds=bounds, maxfev=50000) + 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'] + 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的行星波振幅 + # 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 @@ -117,4 +119,4 @@ def temp_render(path:str = f"{DATA_BASEPATH.cosmic}/cosmic.txt", T = 16): plt.xticks(ticks=tick_positions, labels=tick_labels) - plt.show() # 显示图形 \ No newline at end of file + plt.show() # 显示图形 diff --git a/saber/archive/gravity_plot.py b/saber/archive/gravity_plot.py index 87ccec3..26cfbee 100644 --- a/saber/archive/gravity_plot.py +++ b/saber/archive/gravity_plot.py @@ -1,3 +1,4 @@ +# saber 行星波参数一年逐日图 # 此代码是对数据处理后的txt数据进行行星波参数提取绘图 import pandas as pd diff --git a/saber/archive/gravity_process.py b/saber/archive/gravity_process.py index 50a8528..8116d0b 100644 --- a/saber/archive/gravity_process.py +++ b/saber/archive/gravity_process.py @@ -1,3 +1,4 @@ +# 原始文件 saber行星波处理.py # 高度、纬度可以指定，得到指定高度纬度下的数据，高度在101行，一般输入70、90、110，纬度在115行，纬度一般输入-20、30、60 import os diff --git a/saber/archive/zjy_wave_fit_plot.py b/saber/archive/zjy_wave_fit_plot.py index 112c750..5942ee7 100644 --- a/saber/archive/zjy_wave_fit_plot.py +++ b/saber/archive/zjy_wave_fit_plot.py @@ -1,3 +1,5 @@ +# saber 重力波 + from io import BytesIO import netCDF4 as nc import numpy as np diff --git a/tidi/plot.py b/tidi/plot.py index 1aa7588..336eb5c 100644 --- a/tidi/plot.py +++ b/tidi/plot.py @@ -1,3 +1,5 @@ +# 原文件：TIDI循环重力波 + import os import pandas as pd import numpy as np @@ -377,14 +379,13 @@ def process_vzonal_day(day, year=2015): # 初始化一个空的DataFrame来存储所有天的结果 - # ------------------------------------------------------------------------------------------- # --------meridional------------------------------------------------------------------------- def process_vmeridional_day(day, year=2015): try: # 读取数据 base_path = DATA_BASEPATH.tidi - + height_data = loadmat(rf"{base_path}/{year}/{day:03d}_Height.mat") lat_data = loadmat(rf"{base_path}/{year}/{day:03d}_Lat.mat") lon_data = loadmat(rf"{base_path}/{year}/{day:03d}_Lon.mat") @@ -757,16 +758,19 @@ def day_to_month(day): 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 cache_path = f"{DATA_BASEPATH.tidi}/cache" if os.path.exists(f"{cache_path}/{year}/all_days_vzonal_results.parquet") \ - and os.path.exists(f"{cache_path}/{year}/all_days_vmeridional_results.parquet"): - all_days_vzonal_results = pd.read_parquet(f"{cache_path}/{year}/all_days_vzonal_results.parquet") - all_days_vmeridional_results = pd.read_parquet(f"{cache_path}/{year}/all_days_vmeridional_results.parquet") + and os.path.exists(f"{cache_path}/{year}/all_days_vmeridional_results.parquet"): + all_days_vzonal_results = pd.read_parquet( + f"{cache_path}/{year}/all_days_vzonal_results.parquet") + all_days_vmeridional_results = pd.read_parquet( + f"{cache_path}/{year}/all_days_vmeridional_results.parquet") else: - + all_days_vzonal_results = pd.DataFrame() # 循环处理每一天的数据 @@ -775,7 +779,8 @@ class TidiPlotv2: all_days_vzonal_results[rf"{day:02d}"] = u2 # 将结果按列拼接 - all_days_vzonal_results.columns = [f"{day:02d}" for day in range(1, 365)] + all_days_vzonal_results.columns = [ + f"{day:02d}" for day in range(1, 365)] # 初始化一个空的DataFrame来存储所有天的结果 all_days_vmeridional_results = pd.DataFrame() @@ -785,33 +790,35 @@ class TidiPlotv2: all_days_vmeridional_results[rf"{day:02d}"] = v2 # 将结果按列拼接 - all_days_vmeridional_results.columns = [f"{day:02d}" for day in range(1, 365)] - + all_days_vmeridional_results.columns = [ + f"{day:02d}" for day in range(1, 365)] + # cache the results # if dir not exists, create it if not os.path.exists(f"{cache_path}/{year}"): os.makedirs(f"{cache_path}/{year}") - - all_days_vzonal_results.to_parquet(f"{cache_path}/{year}/all_days_vzonal_results.parquet") - all_days_vmeridional_results.to_parquet(f"{cache_path}/{year}/all_days_vmeridional_results.parquet") - + + all_days_vzonal_results.to_parquet( + f"{cache_path}/{year}/all_days_vzonal_results.parquet") + all_days_vmeridional_results.to_parquet( + f"{cache_path}/{year}/all_days_vmeridional_results.parquet") + 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), + 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] + 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') @@ -830,7 +837,7 @@ class TidiPlotv2: 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] @@ -859,7 +866,7 @@ class TidiPlotv2: # 横坐标过长，设置等间隔展示 interval = 34 # 横坐标显示间隔 plt.xticks(ticks=range(0, len(dates), interval), - labels=months[::interval], rotation=45) # rotation旋转可不加 + labels=months[::interval], rotation=45) # rotation旋转可不加 # 添加轴标签 plt.xlabel('Month') # X轴标签 @@ -867,7 +874,7 @@ class TidiPlotv2: # 显示图形 # plt.show() - + def plot_month(self): HP = self.HP # --------------绘制月统计图------------------------------------------------------------------- @@ -905,7 +912,7 @@ class TidiPlotv2: # 生成x轴的月份标签 months = ["Jan", "Feb", "Mar", "Apr", "May", "Jun", - "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"] + "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"] # clear the last plot plt.clf() diff --git a/tidi/staged/plot.py b/tidi/staged/plot.py index 9b36069..07ff5a5 100644 --- a/tidi/staged/plot.py +++ b/tidi/staged/plot.py @@ -1,4 +1,4 @@ -# 行星波 +# 原始文件 TIDI行星波一年逐日图 # 此代码是对数据处理后的txt数据进行行星波参数提取绘图 # 2006_TIDI_V_Meridional_data.txt也可以做同样处理，一个是经向风提取的行星波，一个是纬向风的 diff --git a/tidi/staged/process.py b/tidi/staged/process.py index 86aa104..415402b 100644 --- a/tidi/staged/process.py +++ b/tidi/staged/process.py @@ -1,3 +1,5 @@ +# 原始文件 TIDI行星波数据处理.py + import os import numpy as np from scipy.io import loadmat