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

fix

Browse Source
This commit is contained in:
Dustella 2025-02-21 14:24:05 +08:00
parent 27a24ba6a8
commit ceeed5de92
Signed by: Dustella
GPG Key ID: 35AA0AA3DC402D5C
8 changed files with 94 additions and 29 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

36
modules/balloon/extract_wave.py
View File

@ -4,6 +4,32 @@ from scipy.optimize import curve_fit
from scipy.signal import lombscargle, hilbert
from sklearn.decomposition import PCA
import sympy as sp
from scipy.interpolate import interp1d
from scipy.signal import lombscargle
import os
import glob
from scipy.optimize import curve_fit
import math
from sklearn.decomposition import PCA
import sympy as sp
import matplotlib.pyplot as plt
import seaborn as sns
from matplotlib.patches import Patch
from matplotlib.lines import Line2D
from windrose import WindroseAxes
from scipy.signal import lombscargle, find_peaks
from itertools import permutations
from itertools import product
def get_top_peaks(pgram, ff, num_peaks=3):
peaks, _ = find_peaks(pgram)
peak_values = pgram[peaks]
peak_frequencies = ff[peaks]
sorted_indices = np.argsort(peak_values)[-num_peaks:][::-1] # 按照峰值从大到小排序
top_peaks_frequencies = peak_frequencies[sorted_indices]
top_peaks_values = peak_values[sorted_indices]
return top_peaks_frequencies, top_peaks_values
def calculate_n(Av, Au, Bu, Bv):
@ -118,7 +144,10 @@ def calculate_wave(data: pd.DataFrame, lat: float, g: float):
pgram_ucmp = lombscargle(height, residual_ucmp, ff)
pgram_vcmp = lombscargle(height, residual_vcmp, ff)
# todo:我的频率等于垂直波数（单位10-3rad/m）？
# BEGIN CHANGE
# BACK
# todo: 我的频率等于垂直波数（单位10-3rad/m）？
main_frequency_temp = round(ff[np.argmax(pgram_temp)], 2)
main_frequency_ucmp = round(ff[np.argmax(pgram_ucmp)], 2)
main_frequency_vcmp = round(ff[np.argmax(pgram_vcmp)], 2)
@ -132,6 +161,11 @@ def calculate_wave(data: pd.DataFrame, lat: float, g: float):
rsd_temp = abs(round(((λ1 - mean) / λ1) * 100, 2))
rsd_ucmp = abs(round(((λ2 - mean) / λ2) * 100, 2))
rsd_vcmp = abs(round(((λ3 - mean) / λ3) * 100, 2))
# BACK
# new
# END CHANGE
if not (rsd_temp < 20 and rsd_ucmp < 20 and rsd_vcmp < 20):
return []

27
modules/balloon/read_data.py
View File

@ -2,9 +2,25 @@ import xarray as xr
import pandas as pd
import numpy as np
from scipy.interpolate import interp1d
from scipy.signal import lombscargle
import os
import glob
from scipy.optimize import curve_fit
import math
from sklearn.decomposition import PCA
import sympy as sp
import matplotlib.pyplot as plt
import seaborn as sns
from matplotlib.patches import Patch
from matplotlib.lines import Line2D
from windrose import WindroseAxes
from scipy.signal import lombscargle, find_peaks
from itertools import permutations
from itertools import product
# 定义四舍五入函数
def round_to_nearest_multiple(value, multiple):
return round(value / multiple) * multiple
@ -24,10 +40,12 @@ def read_data(path):
# 移除重复的高度值
extracted_df = extracted_df.drop_duplicates(subset=["alt"])
new_height = np.arange(extracted_df["alt"].min(), extracted_df["alt"].max() + 0.05, 0.05)
new_height = np.arange(extracted_df["alt"].min(
), extracted_df["alt"].max() + 0.05, 0.05)
# 将每个高度值转换为最接近0.05的整数倍，并转化为数组
rounded_heights = [round_to_nearest_multiple(height, 0.05) for height in new_height]
rounded_heights = [round_to_nearest_multiple(
height, 0.05) for height in new_height]
rounded_heights_np = np.array(rounded_heights)
# 初始化一个新的 DataFrame 用于存储插值结果
@ -35,7 +53,8 @@ def read_data(path):
# 对每个因变量进行线性插值
for col in ["press", "temp", "rh", "u", "v", "wspeed"]:
interp_func = interp1d(extracted_df["alt"], extracted_df[col], kind="linear", fill_value="extrapolate")
interp_func = interp1d(
extracted_df["alt"], extracted_df[col], kind="linear", fill_value="extrapolate")
interpolated_data[col] = interp_func(rounded_heights_np)
return interpolated_data

18
modules/cosmic/__init__.py
View File

@ -35,8 +35,11 @@ async def single_render():
year = request.args.get("year", 2008)
day = request.args.get("day", 1)
mode = request.args.get("mode", "mean_ktemp_Nz")
lat_range = request.args.get("lat_range", "30.0 ~ 40.0")
lat_range = tuple(map(float, lat_range.split("~")))
p: CosmicGravitywPlot = await run_sync(CosmicGravitywPlot)(year=int(year), day=int(day))
p: CosmicGravitywPlot = await run_sync(CosmicGravitywPlot)(
year=int(year), day=int(day), lat_range=lat_range)
if mode == "mean_ktemp_Nz":
await run_sync(p.plot_results_mean_ktemp_Nz)()
elif mode == "mean_ktemp_Ptz":
@ -57,10 +60,15 @@ async def multiday_render():
start_day = request.args.get("start_day", 1)
end_day = request.args.get("end_day", 204)
mode = request.args.get("mode", "位温分布")
lat_range = request.args.get("lat_range", "30.0 ~ 40.0")
lat_range = tuple(map(float, lat_range.split("~")))
p: GravityMultidayPlot = await run_sync(GravityMultidayPlot)(year=int(year),
day_range=(start_day, end_day))
if mode == "布伦特-维萨拉频率分布":
p: GravityMultidayPlot = await run_sync(GravityMultidayPlot)(
year=int(year),
day_range=(start_day, end_day),
lat_range=lat_range
)
if mode == "浮力频率均值":
await run_sync(p.plot_heatmap_tempNz)()
elif mode == "不同高度下的逐日统计分析":
await run_sync(p.plot_heatmap_tempPtz)()
@ -68,6 +76,8 @@ async def multiday_render():
await run_sync(p.plot_floatage_trend)()
elif mode == "每月平均重力势能的折线图":
await run_sync(p.plot_monthly_energy)()
elif mode == "每月平均N^2的折线图":
await run_sync(p.plot_monthly_tempNz)()
else:
raise ValueError("Invalid mode")

12
modules/cosmic/gravityw_multiday.py
View File

@ -21,7 +21,7 @@ plt.rcParams['font.sans-serif'] = ['Simhei'] # 设置字体为微软雅黑
# 主循环，处理1到365个文件
def get_multiday_data(year=2008, day_range=DAY_RANGE):
def get_multiday_data(year=2008, day_range=DAY_RANGE, lat_range=(30, 40)):
base_folder_path = f"{DATA_BASEPATH.cosmic}/{year}"
all_mean_ktemp_Nz = []
all_mean_ktemp_Ptz = []
@ -29,7 +29,7 @@ def get_multiday_data(year=2008, day_range=DAY_RANGE):
for file_index in range(day_begin, day_end):
try:
mean_ktemp_Nz, mean_ktemp_Ptz = process_single_file(
base_folder_path, file_index)
base_folder_path, file_index, lat_range)
if mean_ktemp_Nz is not None and mean_ktemp_Ptz is not None:
all_mean_ktemp_Nz.append(mean_ktemp_Nz)
all_mean_ktemp_Ptz.append(mean_ktemp_Ptz)
@ -83,11 +83,11 @@ def plot_heatmap(data, heights, title):
class GravityMultidayPlot:
def __init__(self, year, day_range):
def __init__(self, year, day_range, lat_range=None):
self.year = year
df_final_mean_ktemp_Nz, df_final_mean_ktemp_Ptz, data, heights = get_multiday_data(
year, day_range)
year, day_range, lat_range)
self.df_final_mean_ktemp_Nz = df_final_mean_ktemp_Nz
self.df_final_mean_ktemp_Ptz = df_final_mean_ktemp_Ptz
self.data = data
@ -148,7 +148,6 @@ class GravityMultidayPlot:
plt.grid(True)
plt.xticks(rotation=45) # 让x轴标签（月份）倾斜，以便更清晰显示
plt.tight_layout()
plt.show()
def plot_monthly_energy(self):
data1 = self.data1
@ -192,7 +191,6 @@ class GravityMultidayPlot:
plt.grid(True)
plt.xticks(rotation=45) # 让x轴标签（月份）倾斜，以便更清晰显示
plt.tight_layout()
plt.show()
def plot_floatage_trend(self):
data = self.data
@ -246,7 +244,6 @@ class GravityMultidayPlot:
# 显示网格线，增强图表可读性
ax.grid(True)
# 显示图形
plt.show()
def plot_floatage_trend(self):
data1 = self.data1
@ -300,4 +297,3 @@ class GravityMultidayPlot:
# 显示网格线，增强图表可读性
ax.grid(True)
# 显示图形
plt.show()

6
modules/cosmic/gravityw_multiday_process.py
View File

@ -10,7 +10,7 @@ import seaborn as sns
import matplotlib.font_manager as fm
def process_single_file(base_folder_path, i):
def process_single_file(base_folder_path, i, lat_range=(30, 40)):
# 构建当前文件夹的路径
if i < 10:
folder_name = f"atmPrf_repro2021_2008_00{i}" # 一位数，前面加两个0
@ -93,9 +93,11 @@ def process_single_file(base_folder_path, i):
# 摄氏度换算开尔文
final_df['Temperature'] = final_df['Temperature'] + 273.15
lat_begin, lat_end = lat_range
# 筛选出纬度在30到40度之间的数据
latitude_filtered_df = final_df[(
final_df['Latitude'] >= 30) & (final_df['Latitude'] <= 40)]
final_df['Latitude'] >= lat_begin) & (final_df['Latitude'] <= lat_end)]
# 划分经度网格，20°的网格
lon_min, lon_max = latitude_filtered_df['Longitude'].min(

9
modules/cosmic/gravityw_perday.py
View File

@ -100,7 +100,7 @@ def read_and_preprocess_data(base_folder_path, day_num):
# 模块2: 数据筛选与网格划分
def filter_and_grid_data(final_df):
def filter_and_grid_data(final_df, lat_range=(30, 40)):
"""
对输入的数据框进行纬度筛选以及经度网格划分等操作
@ -112,8 +112,9 @@ def filter_and_grid_data(final_df):
altitude_temperature_mean (pd.DataFrame): 按高度分组求平均温度后的数据框
"""
# 筛选出纬度在30到40度之间的数据
lat_begin, lat_end = lat_range
latitude_filtered_df = final_df[(
final_df['Latitude'] >= 30) & (final_df['Latitude'] <= 40)]
final_df['Latitude'] >= lat_begin) & (final_df['Latitude'] <= lat_end)]
# 划分经度网格，20°的网格
lon_min, lon_max = latitude_filtered_df['Longitude'].min(
), latitude_filtered_df['Longitude'].max()
@ -381,7 +382,7 @@ def plot_results_mean_ktemp_Ptz(mean_ktemp_Ptz, heights):
class CosmicGravitywPlot:
def __init__(self, year, day):
def __init__(self, year, day, lat_range):
BASE_PATH_COSMIC = DATA_BASEPATH.cosmic
base_folder_path = f"./{BASE_PATH_COSMIC}/{year}"
@ -390,7 +391,7 @@ class CosmicGravitywPlot:
final_df = read_and_preprocess_data(base_folder_path, day_num)
# 模块2调用
latitude_filtered_df, altitude_temperature_mean = filter_and_grid_data(
final_df)
final_df, lat_range)
# 模块3调用
merged_df = calculate_wn0_and_perturbation(
latitude_filtered_df, altitude_temperature_mean)

6
modules/saber/__init__.py
View File

@ -24,10 +24,14 @@ renderer = SaberGravitywRenderer()
def get_gravityw_processer():
lat_str = request.args.get("lat_range")
# if found ~ in str, replace it with ,
if "~" in lat_str:
lat_str = lat_str.replace("~", ",")
if lat_str is None:
lat_str = "30.0,40.0"
print(lat_str)
lat_range = tuple(map(float, lat_str.split(",")))
print(lat_range)
alt_str = request.args.get("alt_range")
if alt_str is None:
alt_str = "20.0,105.0"

9
modules/saber/planetw_plot.py
View File

@ -104,9 +104,9 @@ def saber_planetw_plot(
# 对每一列生成独立的子图
col = k_to_a[f'k={k}']
plt.plot(x_values, fit_df[col].values)
plt.set_title(f'k = {k} 振幅图')
plt.set_xlabel('Day')
plt.set_ylabel('振幅')
plt.suptitle(f'k = {k} 振幅图')
plt.xlabel('Day')
plt.ylabel('振幅')
# 设置横坐标的动态调整
adjusted_x_values = x_values + (3 * T + 1) / 2
@ -118,8 +118,7 @@ def saber_planetw_plot(
tick_positions = adjusted_x_values
tick_labels = [f'{int(val)}' for val in tick_positions]
plt.set_xticks(tick_positions)
plt.set_xticklabels(tick_labels)
plt.xticks(tick_positions, tick_labels)
plt.tight_layout()
# plt.show()