diff --git a/backend.py b/backend.py index 24fa054..3480505 100644 --- a/backend.py +++ b/backend.py @@ -1,3 +1,4 @@ +import resource from matplotlib import pyplot as plt import cosmic import tidi @@ -13,9 +14,14 @@ import matplotlib.font_manager as fm app = Quart(__name__) app = cors(app, send_origin_wildcard=True, allow_origin="*") +# limit the whole app not to use over 8G ram +# use resource module to do this +# resource.setrlimit(resource.RLIMIT_AS, (8 * 1024 * 1024 * 1024, -1)) plt.switch_backend('agg') fm.fontManager.addfont("./SimHei.ttf") + + @app.before_request def auth(): # check for method @@ -27,10 +33,12 @@ def auth(): 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") @@ -46,6 +54,6 @@ if __name__ == '__main__': app.run("0.0.0.0", port=5000, debug=False) pass elif 'debug' in args: - app.run("0.0.0.0",port=18200,debug=True) + app.run("0.0.0.0", port=18200, debug=True) else: raise Exception("Invalied Mode") diff --git a/saber/__init__.py b/saber/__init__.py index 0da10db..c6e6cce 100644 --- a/saber/__init__.py +++ b/saber/__init__.py @@ -3,6 +3,8 @@ from io import BytesIO from quart import Blueprint, request, send_file from matplotlib import pyplot as plt from CONSTANT import DATA_BASEPATH +from saber.archive.gravity_plot import do_gravity_plot +from saber.archive.gravity_process import process_gravity_data from saber.process import DataProcessor from saber.render import Renderer from saber.utils import * @@ -19,12 +21,13 @@ saber_module = Blueprint('saber', __name__) # lat_range=lat_range, alt_range=alt_range, lambda_range=lambda_range, lvboin=lvboin) renderer = Renderer() + def get_processer(): lat_str = request.args.get("lat_range") if lat_str is None: lat_str = "30.0,40.0" lat_range = tuple(map(float, lat_str.split(","))) - + alt_str = request.args.get("alt_range") if alt_str is None: alt_str = "20.0,105.0" @@ -36,7 +39,7 @@ def get_processer(): lvboin = request.args.get("lvboin") == "true" if lvboin is None: lvboin = True - + _p = DataProcessor( lat_range=lat_range, alt_range=alt_range, lambda_range=lambda_range, lvboin=lvboin) return _p @@ -46,7 +49,7 @@ async def extract_payload(): buffer = BytesIO() plt.savefig(buffer, format="png") buffer.seek(0) - return await send_file(buffer, mimetype="image/png") + return await send_file(buffer, mimetype="image/png") _all_saber_files = glob.glob(f"{DATA_BASEPATH.saber}/**/**.nc", recursive=True) @@ -78,7 +81,7 @@ async def do_plot_wave_day_fitting(): processor = await run_sync(get_processer)() data = await run_sync(processor.process_day)(ncfile, int(day)) await run_sync(renderer.plot_wave_fitting)(data, int(height)) - return await extract_payload() + return await extract_payload() @saber_module.route("/render/day_fft_ifft_plot") @@ -92,7 +95,7 @@ async def do_day_fft_ifft_plot(): data = await run_sync(processor.process_day)(ncfile, int(day)) # 高度滤波处理 await run_sync(renderer.day_fft_ifft_plot)(wave_data=data, cycle_no=int(cycle_no)) - return await extract_payload() + return await extract_payload() @saber_module.route("/render/day_cycle_power_wave_plot") @@ -115,9 +118,9 @@ async def do_month_power_wave_plot(): ncfile = await run_sync(data_nc_load)(path) processor = await run_sync(get_processer)() - data =await run_sync(processor.process_month)(ncfile) + data = await run_sync(processor.process_month)(ncfile) await run_sync(renderer.month_power_wave_plot)(wave_data=data, date_time=ncfile.date_time) - return await extract_payload() + return await extract_payload() @saber_module.route("/render/year_power_wave_plot") @@ -129,3 +132,16 @@ async def do_year_power_wave_plot(): ]) await run_sync(renderer.year_power_wave_plot)(year_wave=data) return await extract_payload() + + +@saber_module.route("/render/gravity_year") +async def do_gravity_year(): + year = request.args.get("year") + T = request.args.get("T") + year = int(year) + T = int(T) + if T not in [5, 10, 16]: + raise ValueError("T must be in [5, 10, 16]") + data = await run_sync(process_gravity_data)(year) + await run_sync(do_gravity_plot)(data, T) + return await extract_payload() diff --git a/saber/archive/gravity_plot.py b/saber/archive/gravity_plot.py index 318fb78..87ccec3 100644 --- a/saber/archive/gravity_plot.py +++ b/saber/archive/gravity_plot.py @@ -11,7 +11,7 @@ import matplotlib matplotlib.rcParams['font.sans-serif'] = ['SimHei'] # 显示中文 matplotlib.rcParams['axes.unicode_minus'] = False # 正常显示负号 # 读取一年的数据文件 -df = pd.read_csv(r'C:\pythonProject3\combined_data.txt', sep='\s+') +# df = pd.read_csv(r'C:\pythonProject3\combined_data.txt', sep='\s+') # 设置初始参数 # 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 @@ -26,92 +26,98 @@ bounds = ( 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) - ) +def do_gravity_plot(df: pd.DataFrame, 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) + ) + # 用于存储拟合参数结果的列表 + all_fit_results = [] -# 用于存储拟合参数结果的列表 -all_fit_results = [] + # 设置最小数据量的阈值 + min_data_points = 36 -# 设置最小数据量的阈值 -min_data_points = 36 + # 进行多个时间窗口的拟合 + # T应该取5、10、16 -# 进行多个时间窗口的拟合 -# 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 -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)] - # 选择当前窗口的数据 - 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 # 跳过当前时间窗口，继续下一个窗口 - # 检查当前窗口的数据量 - 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 = 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) - # 用T进行拟合 - popt, pcov = curve_fit(u_func, x, temperature, - p0=initial_guess, bounds=bounds, maxfev=50000) + # 将拟合结果添加到列表中 + all_fit_results.append(popt) - # 将拟合结果添加到列表中 - 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即为拟合的参数汇总 -# 将结果转换为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 -# -------------------------------画图---------------------------- -# 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)} -# 创建一个字典映射 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() -# 获取索引并转换为 numpy 数组 -x_values = fit_df.index.to_numpy() + # 创建一个包含多个子图的图形 + fig, axs = plt.subplots(len(k_to_a), 1, figsize=(10, 2 * len(k_to_a))) -# 对每一列生成独立的图 -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('振幅') + # 对每一列生成独立的子图 + for ax, (k, col) in zip(axs, k_to_a.items()): + ax.plot(x_values, fit_df[col].values) + ax.set_title(f'{k} 振幅图') + ax.set_xlabel('Day') + ax.set_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] + # 设置横坐标的动态调整 + 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) + ax.set_xticks(tick_positions) + ax.set_xticklabels(tick_labels) - plt.show() # 显示图形 + plt.tight_layout() + # plt.show() + + # plt.show() # 显示图形 diff --git a/saber/archive/gravity_process.py b/saber/archive/gravity_process.py index 22f6e4a..50a8528 100644 --- a/saber/archive/gravity_process.py +++ b/saber/archive/gravity_process.py @@ -1,11 +1,13 @@ # 高度、纬度可以指定，得到指定高度纬度下的数据，高度在101行，一般输入70、90、110，纬度在115行，纬度一般输入-20、30、60 +import os import netCDF4 as nc import numpy as np import netCDF4 as nc import numpy as np import matplotlib.pyplot as plt from datetime import datetime, timedelta +import pandas as pd from scipy.optimize import curve_fit from scipy.optimize import least_squares @@ -23,6 +25,13 @@ def process_gravity_data(year): base_path = DATA_BASEPATH.saber + \ "/{year}/SABER_Temp_O3_{month}{year}_v2.0.nc" + cache_path = DATA_BASEPATH.saber + \ + f"/{year}/gravity_result.parquet" + + if os.path.exists(cache_path): + # load from cache + return pd.read_parquet(cache_path) + # 初始化空数组来存储拼接后的数据 tplongitude = None tplatitude = None @@ -165,7 +174,12 @@ def process_gravity_data(year): combined_data = np.column_stack( (closest_longitudes_radians, closest_times, closest_temperatures)) - return combined_data + df = pd.DataFrame(combined_data, columns=[ + "Longitude_Radians", "Time", "Temperature"]) + + # dump to cache + df.to_parquet(cache_path) + return df # 将数据保存为txt文件 diff --git a/tidi/__init__.py b/tidi/__init__.py index df69a97..22a24a9 100644 --- a/tidi/__init__.py +++ b/tidi/__init__.py @@ -11,6 +11,7 @@ from tidi.staged.plot import tidi_render tidi_module = Blueprint("Tidi", __name__) + @tidi_module.route('/metadata') def get_all_years(): res = glob.glob(f"{DATA_BASEPATH.tidi}/**/**.txt", recursive=True) @@ -20,28 +21,35 @@ def get_all_years(): "path": res } + @tidi_module.route('/render/wave') async def render_wave(): mode = request.args.get('mode') year = request.args.get('year') k = request.args.get('k') T = request.args.get('T') + height = request.args.get('height') + lat_range = request.args.get('lat_range') + # like `0 ~ 5` year = int(year) k = int(k) T = int(T) - - await run_sync(tidi_render)(mode, year, k, T) + height = float(height) + lat_range = tuple(map(int, lat_range.split('~'))) + + await run_sync(tidi_render)(mode, year, k, height, lat_range, T) buffer = BytesIO() plt.savefig(buffer, format="png") buffer.seek(0) - - return await send_file(buffer, mimetype="image/png") + + return await send_file(buffer, mimetype="image/png") + @tidi_module.route('/render/month_stats_v1') async def render_stats_v1(): year = request.args.get('year') year = int(year) - + plotter = await run_sync(TidiPlotv2)(year) await run_sync(plotter.plot_v1)() buffer = BytesIO() @@ -49,11 +57,12 @@ async def render_stats_v1(): buffer.seek(0) return await send_file(buffer, mimetype="image/png") + @tidi_module.route('/render/month_stats_v2') async def render_stats_v2(): year = request.args.get('year') year = int(year) - + plotter = await run_sync(TidiPlotv2)(year) await run_sync(plotter.plot_month)() buffer = BytesIO() diff --git a/tidi/staged/plot.py b/tidi/staged/plot.py index f557c01..9b36069 100644 --- a/tidi/staged/plot.py +++ b/tidi/staged/plot.py @@ -1,8 +1,9 @@ # 行星波 -#此代码是对数据处理后的txt数据进行行星波参数提取绘图 -#2006_TIDI_V_Meridional_data.txt也可以做同样处理，一个是经向风提取的行星波，一个是纬向风的 +# 此代码是对数据处理后的txt数据进行行星波参数提取绘图 +# 2006_TIDI_V_Meridional_data.txt也可以做同样处理，一个是经向风提取的行星波，一个是纬向风的 +import os import pandas as pd import numpy as np from scipy.optimize import curve_fit @@ -12,25 +13,34 @@ import matplotlib.pyplot as plt import matplotlib from CONSTANT import DATA_BASEPATH +from tidi.staged.process import tidi_process_idk # 设置中文显示和负号正常显示 matplotlib.rcParams['font.sans-serif'] = ['SimHei'] # 显示中文 matplotlib.rcParams['axes.unicode_minus'] = False # 正常显示负号 # 读取一年的数据文件 # 设置初始参数 # 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]) # 上界 # 定义拟合函数 -def tidi_render(mode, year, k, T=10): +def tidi_render( + mode, + year, + k, + input_height, + lat_range=(0, 5), + T=10 +): # 用于存储拟合参数结果的列表 all_fit_results = [] @@ -38,27 +48,37 @@ def tidi_render(mode, year, k, T=10): min_data_points = 36 if mode != 'V_Zonal' and mode != 'V_Meridional': raise ValueError('mode must be V_Zonal or V_Meridional') - if k not in list(range(-4, 5)): + if k not in list(range(-4, 5)): raise ValueError('k must be in range(-4, 5)') - df = pd.read_csv(f'{DATA_BASEPATH.tidi}/{year}/TIDI_{mode}_data.txt', sep='\s+') + data_cache_path = f'{DATA_BASEPATH.tidi}/{year}/TIDI_{mode}h{input_height}r{lat_range[0]}-{lat_range[1]}_data.parquet' + # if cache exist, read from cache + if os.path.exists(data_cache_path): + df = pd.read_parquet(data_cache_path) + else: + df = tidi_process_idk( + year, input_height, lat_range + )[mode] + df.to_parquet(data_cache_path) + # 删除有 NaN 的行 # V_Meridional df = df.dropna(subset=[mode]) # 进行多个时间窗口的拟合 - #T应该取5、10、16 + # T应该取5、10、16 # todo:对于5日波，滑动窗口选择15天 + 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) ) for start_day in range(0, 365-3*T): # 最后一个窗口为[351, 366] end_day = start_day + 3 * T # 每个窗口的结束时间为 start_day + 3*T @@ -80,17 +100,19 @@ def tidi_render(mode, year, k, T=10): temperature = np.array(df_8[mode]) # 经向风速，因变量 # 用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 @@ -121,7 +143,8 @@ def tidi_render(mode, year, k, T=10): plt.xticks(ticks=tick_positions, labels=tick_labels) # plt.show() # 显示图形 - + + if __name__ == '__main__': tidi_render('V_Zonal', 2015, 1) tidi_render('V_Meridional', 2015, 1) @@ -203,7 +226,3 @@ if __name__ == '__main__': # # 显示图形 # plt.tight_layout() # plt.show() - - - - diff --git a/tidi/staged/process.py b/tidi/staged/process.py index 0a10ef4..86aa104 100644 --- a/tidi/staged/process.py +++ b/tidi/staged/process.py @@ -5,184 +5,201 @@ import pandas as pd from CONSTANT import DATA_BASEPATH -#第156行代码，纬度可以自己选择，但是都是每5°做一个选择 -#70km-120km每2.5km一个，共21个，第六个就是82.5km高度的情况，高度也可以自己选择 +# 第156行代码，纬度可以自己选择，但是都是每5°做一个选择 +# 70km-120km每2.5km一个，共21个，第六个就是82.5km高度的情况，高度也可以自己选择 -fixed_height_index = 6 -# 初始化空列表来存储每天的数据 -height_list = [] -lat_list = [] -lon_list = [] -vmeridional_list = [] -vzonal_list = [] - -# 文件路径 -base_path = f"{DATA_BASEPATH.tidi}/2022/" - -# 循环读取从第1天到第365天的数据 -for day in range(1, 366): # 365天数据，确保循环次数为365 - # 构建文件名（使用三位数格式） - day_str = f"{day:03d}" # 格式化数字为三位数，前面补0 - height_file = f"{day_str}_Height.mat" - lat_file = f"{day_str}_Lat.mat" - lon_file = f"{day_str}_Lon.mat" - vmeridional_file = f"{day_str}_VMerdional.mat" - vzonal_file = f"{day_str}_Vzonal.mat" - - # 检查文件是否存在 - try: - if not os.path.exists(os.path.join(base_path, height_file)): - raise FileNotFoundError(f"{height_file} not found") - if not os.path.exists(os.path.join(base_path, lat_file)): - raise FileNotFoundError(f"{lat_file} not found") - if not os.path.exists(os.path.join(base_path, lon_file)): - raise FileNotFoundError(f"{lon_file} not found") - if not os.path.exists(os.path.join(base_path, vmeridional_file)): - raise FileNotFoundError(f"{vmeridional_file} not found") - if not os.path.exists(os.path.join(base_path, vzonal_file)): - raise FileNotFoundError(f"{vzonal_file} not found") - - # 读取.mat文件 - height_data = loadmat(os.path.join(base_path, height_file)) - lat_data = loadmat(os.path.join(base_path, lat_file)) - lon_data = loadmat(os.path.join(base_path, lon_file)) - vmeridional_data = loadmat(os.path.join(base_path, vmeridional_file)) - vzonal_data = loadmat(os.path.join(base_path, vzonal_file)) - - # 将数据转换为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']) - - # 将每天的数据添加到列表中 - height_list.append(height_df) - lat_list.append(lat_df) - lon_list.append(lon_df) - vmeridional_list.append(vmeridional_df) - vzonal_list.append(vzonal_df) - except FileNotFoundError as e: - print(f"Error: {e}") - continue # 跳过当前文件，继续处理其他文件 - -# 合并所有天的数据 -height_df = pd.concat(height_list, ignore_index=True) -lat_df = pd.concat(lat_list, ignore_index=True) -lon_df = pd.concat(lon_list, ignore_index=True) -vmeridional_df = pd.concat(vmeridional_list, axis=1) # 按列合并 -vzonal_df = pd.concat(vzonal_list, axis=1) # 按列合并 - -# 初始化空列表来存储每天的数据 -UTime_list = [] -# 初始化累加的时间偏移量 -time_offset = 0 - -# 循环读取从第1天到第365天的数据 -for day in range(1, 366): # 365天数据 - # 构建文件名（使用三位数格式） - day_str = f"{day:03d}" # 格式化数字为三位数，前面补0 - UTime_file = f"{day_str}_UTime.mat" - - # 检查UTime文件是否存在 - try: - if not os.path.exists(os.path.join(base_path, UTime_file)): - raise FileNotFoundError(f"{UTime_file} not found") - - # 读取.mat文件 - UTime_data = loadmat(os.path.join(base_path, UTime_file)) - - # 将数据转换为DataFrame - UTime_df = pd.DataFrame(UTime_data['UTime']) - - # 对UTime进行累加处理 - UTime_df += time_offset - - # 将每天的数据添加到列表中 - UTime_list.append(UTime_df) - - # 更新时间偏移量，为下一天增加24小时 - time_offset += 24 - except FileNotFoundError as e: - print(f"Error: {e}") - continue # 跳过当前文件，继续处理其他文件 - -# 合并所有天的数据 -UTime_df = pd.concat(UTime_list, ignore_index=True) -# 将UTime_df的单位从小时转换为天（除以24） -UTime_df = UTime_df / 24 +HEIGHTS_CONSTANT = [ + 70, 72.5, 75, 77.5, 80, 82.5, 85, 87.5, 90, 92.5, 95, 97.5, 100, 102.5, 105, 107.5, 110, 112.5, 115, 117.5, 120 +] -# 获取固定高度下的纬度数据 -fixed_height_lat = lat_df.iloc[:, 0].values -# 获取固定高度下的经度数据 -fixed_height_lon = lon_df.iloc[:, 0].values -# 获取固定高度下的经向风数据 -fixed_height_vmeridional = vmeridional_df.iloc[fixed_height_index, :].values -# 获取固定高度下的纬向风数据 -fixed_height_vzonal = vzonal_df.iloc[fixed_height_index, :].values -# 获取对应的世界时数据（这里直接使用整个UTime数据，你可以根据实际情况判断是否需要处理下格式等） -fixed_height_utime = UTime_df.iloc[:, 0].values -# 将这些数据整合到一个新的DataFrame（可选操作，方便后续查看等） -combined_data = pd.DataFrame({ - 'Latitude': fixed_height_lat, - 'Longitude': fixed_height_lon, - 'V_Meridional': fixed_height_vmeridional, - 'V_Zonal': fixed_height_vzonal, - 'UTime': fixed_height_utime -}) +def tidi_process_idk( + year, + input_height: float = 82.5, + lat_range: tuple = (0, 5), +): + # 初始化空列表来存储每天的数据 + if lat_range[1] - lat_range[0] != 5: + raise ValueError("lat_range must be 5°") + if input_height not in HEIGHTS_CONSTANT: + raise ValueError(f"input_height must be in {HEIGHTS_CONSTANT}") + fixed_height_index = HEIGHTS_CONSTANT.index(input_height) + height_list = [] + lat_list = [] + lon_list = [] + vmeridional_list = [] + vzonal_list = [] -# 简单打印下整合后的数据前几行查看下内容 -print(combined_data) + # 文件路径 + base_path = f"{DATA_BASEPATH.tidi}/{year}/" + # if dir not exist, raise error + if not os.path.exists(base_path): + raise FileNotFoundError(f"{base_path} not found") -# 确定纬度和经度的范围 -lat_min, lat_max = np.min(fixed_height_lat), np.max(fixed_height_lat) -lon_min, lon_max = np.min(fixed_height_lon), np.max(fixed_height_lon) + # 循环读取从第1天到第365天的数据 + for day in range(1, 366): # 365天数据，确保循环次数为365 + # 构建文件名（使用三位数格式） + day_str = f"{day:03d}" # 格式化数字为三位数，前面补0 + height_file = f"{day_str}_Height.mat" + lat_file = f"{day_str}_Lat.mat" + lon_file = f"{day_str}_Lon.mat" + vmeridional_file = f"{day_str}_VMerdional.mat" + vzonal_file = f"{day_str}_Vzonal.mat" -# 计算网格的边界 -lat_bins = np.linspace(lat_min, lat_max, 10) # 纬度分为6个网格，需要7个边界 -lon_bins = np.linspace(lon_min, lon_max, 13) # 经度分为20个网格，需要21个边界 + # 检查文件是否存在 + try: + if not os.path.exists(os.path.join(base_path, height_file)): + raise FileNotFoundError(f"{height_file} not found") + if not os.path.exists(os.path.join(base_path, lat_file)): + raise FileNotFoundError(f"{lat_file} not found") + if not os.path.exists(os.path.join(base_path, lon_file)): + raise FileNotFoundError(f"{lon_file} not found") + if not os.path.exists(os.path.join(base_path, vmeridional_file)): + raise FileNotFoundError(f"{vmeridional_file} not found") + if not os.path.exists(os.path.join(base_path, vzonal_file)): + raise FileNotFoundError(f"{vzonal_file} not found") -# 将数据分配到网格中 -grid_data = [] -for i in range(len(fixed_height_lat)): - lat_idx = np.digitize(fixed_height_lat[i], lat_bins) - 1 # 找到纬度所在的网格索引 - lon_idx = np.digitize(fixed_height_lon[i], lon_bins) - 1 # 找到经度所在的网格索引 - grid_idx = lat_idx * 12 + lon_idx # 计算网格的唯一索引（9x12=108） - grid_data.append([fixed_height_lat[i], fixed_height_lon[i], fixed_height_vmeridional[i], fixed_height_vzonal[i], fixed_height_utime[i], grid_idx]) -# 将网格数据转换为DataFrame -grid_df = pd.DataFrame(grid_data, columns=['Latitude', 'Longitude', 'V_Meridional', 'V_Zonal', 'UTime', 'Grid_Index']) -# 打印网格数据的前几行查看 -print(grid_df.head()) -# 筛选纬度在 0 到 5° 范围内的数据 -filtered_df = grid_df[(grid_df['Latitude'] >= 0) & (grid_df['Latitude'] <= 5)] -# 将经度从度转换为弧度 -filtered_df['Longitude_Radians'] = np.radians(filtered_df['Longitude']) -# 选择需要输出的列 -output_columns = ['Longitude_Radians', 'UTime', 'V_Zonal'] -# 将数据输出为txt文件 -output_file = base_path+'TIDI_V_Zonal_data.txt'#经向风速分量 -filtered_df[output_columns].to_csv(output_file, sep='\t', index=False) -# 打印结果确认 -print(f"数据已保存到 {output_file}") -# 选择需要输出的列 -output_columns1 = ['Longitude_Radians', 'UTime', 'V_Meridional'] -# 将数据输出为txt文件 -output_file1 = base_path+'TIDI_V_Meridional_data.txt'#纬向风速分量 -filtered_df[output_columns1].to_csv(output_file1, sep='\t', index=False) -# 打印结果确认 -print(f"数据已保存到 {output_file1}") + # 读取.mat文件 + height_data = loadmat(os.path.join(base_path, height_file)) + lat_data = loadmat(os.path.join(base_path, lat_file)) + lon_data = loadmat(os.path.join(base_path, lon_file)) + vmeridional_data = loadmat( + os.path.join(base_path, vmeridional_file)) + vzonal_data = loadmat(os.path.join(base_path, vzonal_file)) + # 将数据转换为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']) + # 将每天的数据添加到列表中 + height_list.append(height_df) + lat_list.append(lat_df) + lon_list.append(lon_df) + vmeridional_list.append(vmeridional_df) + vzonal_list.append(vzonal_df) + except FileNotFoundError as e: + print(f"Error: {e}") + continue # 跳过当前文件，继续处理其他文件 + # 合并所有天的数据 + height_df = pd.concat(height_list, ignore_index=True) + lat_df = pd.concat(lat_list, ignore_index=True) + lon_df = pd.concat(lon_list, ignore_index=True) + vmeridional_df = pd.concat(vmeridional_list, axis=1) # 按列合并 + vzonal_df = pd.concat(vzonal_list, axis=1) # 按列合并 + # 初始化空列表来存储每天的数据 + UTime_list = [] + # 初始化累加的时间偏移量 + time_offset = 0 + # 循环读取从第1天到第365天的数据 + for day in range(1, 366): # 365天数据 + # 构建文件名（使用三位数格式） + day_str = f"{day:03d}" # 格式化数字为三位数，前面补0 + UTime_file = f"{day_str}_UTime.mat" + # 检查UTime文件是否存在 + try: + if not os.path.exists(os.path.join(base_path, UTime_file)): + raise FileNotFoundError(f"{UTime_file} not found") + # 读取.mat文件 + UTime_data = loadmat(os.path.join(base_path, UTime_file)) + # 将数据转换为DataFrame + UTime_df = pd.DataFrame(UTime_data['UTime']) + # 对UTime进行累加处理 + UTime_df += time_offset + # 将每天的数据添加到列表中 + UTime_list.append(UTime_df) + # 更新时间偏移量，为下一天增加24小时 + time_offset += 24 + except FileNotFoundError as e: + print(f"Error: {e}") + continue # 跳过当前文件，继续处理其他文件 + + # 合并所有天的数据 + UTime_df = pd.concat(UTime_list, ignore_index=True) + # 将UTime_df的单位从小时转换为天（除以24） + UTime_df = UTime_df / 24 + + # 获取固定高度下的纬度数据 + fixed_height_lat = lat_df.iloc[:, 0].values + # 获取固定高度下的经度数据 + fixed_height_lon = lon_df.iloc[:, 0].values + # 获取固定高度下的经向风数据 + fixed_height_vmeridional = vmeridional_df.iloc[fixed_height_index, :].values + # 获取固定高度下的纬向风数据 + fixed_height_vzonal = vzonal_df.iloc[fixed_height_index, :].values + # 获取对应的世界时数据（这里直接使用整个UTime数据，你可以根据实际情况判断是否需要处理下格式等） + fixed_height_utime = UTime_df.iloc[:, 0].values + # 将这些数据整合到一个新的DataFrame（可选操作，方便后续查看等） + combined_data = pd.DataFrame({ + 'Latitude': fixed_height_lat, + 'Longitude': fixed_height_lon, + 'V_Meridional': fixed_height_vmeridional, + 'V_Zonal': fixed_height_vzonal, + 'UTime': fixed_height_utime + }) + + # 简单打印下整合后的数据前几行查看下内容 + print(combined_data) + + # 确定纬度和经度的范围 + lat_min, lat_max = np.min(fixed_height_lat), np.max(fixed_height_lat) + lon_min, lon_max = np.min(fixed_height_lon), np.max(fixed_height_lon) + + # 计算网格的边界 + lat_bins = np.linspace(lat_min, lat_max, 10) # 纬度分为6个网格，需要7个边界 + lon_bins = np.linspace(lon_min, lon_max, 13) # 经度分为20个网格，需要21个边界 + + # 将数据分配到网格中 + grid_data = [] + for i in range(len(fixed_height_lat)): + lat_idx = np.digitize(fixed_height_lat[i], lat_bins) - 1 # 找到纬度所在的网格索引 + lon_idx = np.digitize(fixed_height_lon[i], lon_bins) - 1 # 找到经度所在的网格索引 + grid_idx = lat_idx * 12 + lon_idx # 计算网格的唯一索引（9x12=108） + grid_data.append([fixed_height_lat[i], fixed_height_lon[i], fixed_height_vmeridional[i], + fixed_height_vzonal[i], fixed_height_utime[i], grid_idx]) + # 将网格数据转换为DataFrame + grid_df = pd.DataFrame(grid_data, columns=[ + 'Latitude', 'Longitude', 'V_Meridional', 'V_Zonal', 'UTime', 'Grid_Index']) + # 打印网格数据的前几行查看 + print(grid_df.head()) + # 筛选纬度在 0 到 5° 范围内的数据 + lat_range_down, lat_range_upper = lat_range + + filtered_df = grid_df[(grid_df['Latitude'] >= lat_range_down) & + (grid_df['Latitude'] <= lat_range_upper)] + # 将经度从度转换为弧度 + filtered_df['Longitude_Radians'] = np.radians(filtered_df['Longitude']) + # 选择需要输出的列 + output_columns = ['Longitude_Radians', 'UTime', 'V_Zonal'] + # 将数据输出为txt文件 + output_file = base_path+'TIDI_V_Zonal_data.txt' # 经向风速分量 + filtered_df[output_columns].to_csv(output_file, sep='\t', index=False) + + zonal_data = filtered_df[output_columns] + # 选择需要输出的列 + output_columns1 = ['Longitude_Radians', 'UTime', 'V_Meridional'] + # 将数据输出为txt文件 + merdi_data = filtered_df[output_columns1] + + return { + "V_Zonal": zonal_data, + "V_Meridional": merdi_data + } + # output_file1 = base_path+'TIDI_V_Meridional_data.txt' # 纬向风速分量 + # final_result.to_csv(output_file1, sep='\t', index=False) + # # 打印结果确认 + # print(f"数据已保存到 {output_file1}") # # 对网格数据按照Grid_Index进行分组，并计算每个网格的统计数据 @@ -209,4 +226,4 @@ print(f"数据已保存到 {output_file1}") # print("无法拟合行星波") # else: # # 打印筛选后的网格数据的前几行查看 -# print(filtered_grid_stats.head()) \ No newline at end of file +# print(filtered_grid_stats.head())