diff --git a/modules/cosmic/__init__.py b/modules/cosmic/__init__.py index 8a9136e..9f46409 100644 --- a/modules/cosmic/__init__.py +++ b/modules/cosmic/__init__.py @@ -1,12 +1,16 @@ from io import BytesIO +import pandas as pd from quart import Blueprint, request, send_file from matplotlib import pyplot as plt +from CONSTANT import DATA_BASEPATH from modules.cosmic.gravityw_multiday import GravityMultidayPlot from modules.cosmic.gravityw_perday import CosmicGravitywPlot from modules.cosmic.planetw_daily import cosmic_planetw_daily_plot from quart.utils import run_sync +from modules.cosmic.planetw_daily_process import cosmic_planet_daily_process + cosmic_module = Blueprint("Cosmic", __name__) @@ -19,9 +23,24 @@ def get_meta(): @cosmic_module.route('/render/planet_wave/daily') @cosmic_module.route('/temp_render') async def render(): + year = request.args.get("year", 2008) T = request.args.get("T", 16) k = request.args.get("k", 1) - await run_sync(cosmic_planetw_daily_plot)(T=int(T), k=int(k)) + target_h = request.args.get("target_h", 40) + target_latitude = request.args.get("target_lat", 30) + + # path: str = f"{DATA_BASEPATH.cosmic}/cosmic.txt" + temp_df = await run_sync(cosmic_planet_daily_process)( + year=int(year), + target_h=int(target_h), + target_latitude=int(target_latitude) + ) + + await run_sync(cosmic_planetw_daily_plot)( + temp_df, + T=int(T), + k=int(k) + ) buf = BytesIO() plt.savefig(buf, format="png") diff --git a/modules/cosmic/planetw_daily.py b/modules/cosmic/planetw_daily.py index 476a9b1..f751bf9 100644 --- a/modules/cosmic/planetw_daily.py +++ b/modules/cosmic/planetw_daily.py @@ -18,7 +18,7 @@ matplotlib.rcParams['axes.unicode_minus'] = False # 正常显示负号 def cosmic_planetw_daily_plot( - path: str = f"{DATA_BASEPATH.cosmic}/cosmic.txt", + df: pd.DataFrame, T=16, k=0 ): @@ -36,7 +36,6 @@ def cosmic_planetw_daily_plot( + a9 * np.sin((2 * np.pi / T) * t + 4 * x + b9) ) - df = pd.read_csv(path, sep='\s+') # 删除有 NaN 的行 df = df.dropna(subset=['Temperature']) # 设置初始参数 diff --git a/modules/cosmic/planetw_daily_process.py b/modules/cosmic/planetw_daily_process.py new file mode 100644 index 0000000..220535c --- /dev/null +++ b/modules/cosmic/planetw_daily_process.py @@ -0,0 +1,150 @@ +# 高度和target_latitude可以自己选择,纬度可以选+-30,+-60,高度可以每10km +import netCDF4 as nc +import pandas as pd +import numpy as np +import os +from scipy.interpolate import interp1d +import matplotlib.pyplot as plt + +from CONSTANT import DATA_BASEPATH + + +def cosmic_planet_daily_process( + year=2008, + target_latitude=30, + target_h=40 +): + # 用于存储所有DataFrame的列表 + dfs = [] + # 设置文件夹的路径模板 + base_folder_path = f"{DATA_BASEPATH.cosmic}/{year}" + + cache_path = f"{base_folder_path}/planet_{target_h}_{target_latitude}.parquet" + if os.path.exists(cache_path): + return pd.read_parquet(cache_path) + + # 遍历文件夹序号1到365 + for i in range(1, 365): + # 根据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): + # 遍历文件夹中的文件 + 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'][:] + # 读取month, hour, minute + month = dataset.month + hour = dataset.hour + minute = dataset.minute + # 将i的值作为day的值 + day = i + # 将day, hour, minute拼接成一个字符串 + datetime_str = f"{day:03d}{hour:02d}{minute:02d}" + # 确保所有变量都是一维数组 + temp = np.squeeze(temp) + altitude = np.squeeze(altitude) + lat = np.squeeze(lat) + lon = np.squeeze(lon) + # 检查所有数组的长度是否相同 + assert len(temp) == len(altitude) == len(lat) == len( + lon), "Arrays must be the same length" + # 创建DataFrame,并将datetime_str作为第一列添加 + df = pd.DataFrame({ + 'Datetime_str': datetime_str, + 'Longitude': lon, + 'Latitude': lat, + 'Altitude': altitude, + 'Temperature': temp + }) + dataset.close() + # 仅筛选高度 + + df_filtered = df[(df['Altitude'] >= target_h - 0.5) + & (df['Altitude'] <= target_h + 0.5)] + dfs.append(df_filtered) + except Exception as e: + print(f"处理文件 {finfo} 时出错: {e}") + else: + print(f"文件夹 {folder_path} 不存在。") + + # 按行拼接所有DataFrame + final_df = pd.concat(dfs, axis=0, ignore_index=True) + # 假设 final_df 是你的大 DataFrame + final_df.to_csv("output365.txt", index=False, sep=",") # 逗号分隔 + print("数据已保存为 output365.txt") + # 计算差值找到最接近的索引来筛选纬度为30的数据 + # 目标值 + + # 计算 'Latitude' 列与目标值的差异 + final_df['Latitude_diff'] = np.abs(final_df['Latitude'] - target_latitude) + # 按天分组并找到每一天最接近的行 + + def find_closest_row(group): + return group.loc[group['Latitude_diff'].idxmin()] + + # 使用groupby按Day分组并找到每个分组的最接近的记录 + closest_per_day = final_df.groupby('Datetime_str', as_index=False).apply( + find_closest_row, include_groups=False) + # 删除不需要的列 + closest_per_day = closest_per_day.drop(columns=['Latitude_diff']) + closest_per_day = closest_per_day[(closest_per_day['Latitude'] >= ( + target_latitude - 2)) & (closest_per_day['Latitude'] <= (target_latitude + 2))] + # 将经度转换为弧度制单位 + closest_per_day['Longitude_rad'] = np.radians(closest_per_day['Longitude']) + print(closest_per_day) + # 函数:将 Datetime_str 转换为以天为单位的时间 + + def convert_to_days(datetime_str): + # 获取日期、小时和分钟 + dd = int(datetime_str[:3]) # 日期部分 (DDD) + hh = int(datetime_str[3:5]) # 小时部分 (HH) + mm = int(datetime_str[5:]) # 分钟部分 (MM) + + # 计算总分钟数 + total_minutes = dd * 1440 + hh * 60 + mm + + # 转换为天数 + return total_minutes / 1440 + + # 应用函数转换 + closest_per_day['Time'] = closest_per_day['Datetime_str'].apply( + convert_to_days) + # 时间t + day_data = closest_per_day['Time'] + # 计算背景温度,时间Day和弧度经度的函数 + background_temperature = closest_per_day['Temperature'].mean() + # 计算温度扰动 + closest_per_day['Temperature_perturbation'] = closest_per_day['Temperature'] - \ + background_temperature + # 选择需要的列并重命名 + selected_columns = closest_per_day[[ + 'Longitude_rad', 'Time', 'Temperature']] + selected_columns.columns = ['Longitude_Radians', 'Time', 'Temperature'] + # 输出到 txt 文件 + # save to parquet + selected_columns.to_parquet(cache_path) + + return selected_columns + + # output_file_path = 'cosmic.txt' # 你可以修改这个路径来指定文件的保存位置 +# selected_columns.to_csv(output_file_path, sep='\t', index=False) +# # 输出成功提示 +# print("输出成功,文件已保存为 'cosmic.txt'") diff --git a/modules/saber/__init__.py b/modules/saber/__init__.py index a771fe5..8971971 100644 --- a/modules/saber/__init__.py +++ b/modules/saber/__init__.py @@ -152,11 +152,15 @@ async def do_gravity_year(): T = request.args.get("T") k = request.args.get("k") H = request.args.get("H") - lat_range = request.args.get("lat_range") + lat_target = request.args.get("lat_target") 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(saber_planetw_process)(year) + data = await run_sync(saber_planetw_process)( + year, + target_h=int(H), + target_lat=int(lat_target) + ) await run_sync(saber_planetw_plot)(data, T, k=int(k)) return await extract_payload() diff --git a/modules/saber/planetw_process.py b/modules/saber/planetw_process.py index 771c8f5..65eae00 100644 --- a/modules/saber/planetw_process.py +++ b/modules/saber/planetw_process.py @@ -20,7 +20,11 @@ months = [ ] -def saber_planetw_process(year): +def saber_planetw_process( + year, + target_h=70, + target_lat=60 +): # 定义文件路径模板和年份 base_path = DATA_BASEPATH.saber + \ @@ -116,7 +120,7 @@ def saber_planetw_process(year): adjusted_time_day = adjusted_time_hour / 24.0 # 转换为天 # 找到第一个事件中与70最接近的高度的索引 # 输入高度 - target_h = 70 + first_event_altitudes = tpaltitude[0] diff = np.abs(first_event_altitudes - target_h) closest_height_index = np.argmin(diff) @@ -130,7 +134,7 @@ def saber_planetw_process(year): # 1. 创建经度分组(0到360°,每20°一个区间) longitude_bins = np.arange(0, 361, 20) # 2. 创建纬度目标值 - target_lat = 60 # 假设目标纬度为60° + # 假设目标纬度为60° # 3. 用于存储结果 closest_event_indices = []