zephyr/zephyr-backend
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

feat: cosmic thing

Browse Source
This commit is contained in:
Dustella 2025-02-21 15:05:38 +08:00
parent ceeed5de92
commit b7c522004b
Signed by: Dustella
GPG Key ID: 35AA0AA3DC402D5C
5 changed files with 184 additions and 8 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

21
modules/cosmic/__init__.py
View File

@ -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")

3
modules/cosmic/planetw_daily.py
View File

@ -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'])
# 设置初始参数

150
modules/cosmic/planetw_daily_process.py Normal file
View File

@ -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'")

8
modules/saber/__init__.py
View File

@ -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()

10
modules/saber/planetw_process.py
View File

@ -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 = []