From 05a2ea6c166a2ed2abd6a432c0e2db2d8c6be351 Mon Sep 17 00:00:00 2001 From: Dustella Date: Wed, 5 Mar 2025 11:40:19 +0800 Subject: [PATCH] sync --- .gitignore | 5 +- CONSTANT.py | 10 +-- backend.py | 56 +++++++++---- modules/balloon/__init__.py | 12 +-- modules/balloon/extract_wave.py | 71 ++++++++++++---- modules/balloon/gravityw_plot_perday.py | 32 ++++---- modules/balloon/utils.py | 80 +++++++++++++----- modules/cosmic/__init__.py | 26 ++++++ modules/cosmic/gravityw_multiday.py | 9 +- modules/cosmic/gravityw_perday.py | 8 +- modules/cosmic/planetw_daily.py | 2 +- modules/cosmic/planetw_daily_process.py | 105 ++++++++++++++---------- modules/cosmic/planetw_perday.py | 79 ++++++++++++++++++ modules/saber/__init__.py | 22 +++++ modules/saber/gravityw_render.py | 4 +- modules/saber/planetw_perday.py | 86 +++++++++++++++++++ modules/saber/planetw_plot.py | 2 +- modules/saber/planetw_process.py | 2 +- modules/tidi/gravity_wave_year.py | 10 +-- modules/tidi/planet_wave_daily.py | 2 +- ping.py | 2 +- 21 files changed, 478 insertions(+), 147 deletions(-) create mode 100644 modules/cosmic/planetw_perday.py create mode 100644 modules/saber/planetw_perday.py diff --git a/.gitignore b/.gitignore index d503dba..a690175 100644 --- a/.gitignore +++ b/.gitignore @@ -18,4 +18,7 @@ dist *.spec notebooks passcode -data \ No newline at end of file +data + +res +*.txt diff --git a/CONSTANT.py b/CONSTANT.py index 8b6a11a..9e3eec9 100644 --- a/CONSTANT.py +++ b/CONSTANT.py @@ -3,8 +3,8 @@ from dataclasses import dataclass @dataclass class DATA_BASEPATH: - balloon = "../data/balloon" - cosmic = "../data/cosmic" - radar = "../data/radar" - saber = "../data/saber" - tidi = "../data/tidi" + balloon = "./data/balloon" + cosmic = "./data/cosmic" + radar = "./data/radar" + saber = "./data/saber" + tidi = "./data/tidi" diff --git a/backend.py b/backend.py index 7842928..a98161d 100644 --- a/backend.py +++ b/backend.py @@ -1,3 +1,4 @@ +import os import resource from matplotlib import pyplot as plt from modules import cosmic @@ -12,7 +13,7 @@ import sys import matplotlib.font_manager as fm import ping -app = Quart(__name__) +app = Quart(__name__, static_folder="./res") app = cors(app, send_origin_wildcard=True, allow_origin="*") # limit the whole app not to use over 8G ram @@ -22,29 +23,47 @@ app = cors(app, send_origin_wildcard=True, allow_origin="*") plt.switch_backend('agg') fm.fontManager.addfont("./SimHei.ttf") +with open("./passcode.txt", "r") as f: + code = f.read() + if code: + code = code.strip() + else: + code = "0101" + @app.before_request def auth(): # check for method # if it is OPTIONS, do not check for auth - if request.method == "OPTIONS": + if request.method == "OPTIONS" or not request.path.startswith("/api"): return - code = request.headers.get("Authorization") - print(code) - if code != "0101": + if request.path.startswith("/api/ping"): + return + _code = request.headers.get("Authorization") + if _code != code: return "Unauthorized", 401 -@app.route("/ping") -def dping(): - return "pong" +@app.route('/') +async def return_index_html(): + return await app.send_static_file("index.html") -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") -app.register_blueprint(tidi.tidi_module, url_prefix="/tidi") -app.register_blueprint(cosmic.cosmic_module, url_prefix="/cosmic") +@app.route("/") +async def index(path): + # check if there is path in static folder + # if not, return index.html + if app.static_folder: + # check if the file exists + if os.path.exists(os.path.join(app.static_folder, path)): + return await app.send_static_file(path) + return await app.send_static_file("index.html") + +app.register_blueprint(balloon.balloon_module, url_prefix="/api/balloon") +app.register_blueprint(radar.radar_module, url_prefix="/api/radar") +app.register_blueprint(saber.saber_module, url_prefix="/api/saber") +app.register_blueprint(tidi.tidi_module, url_prefix="/api/tidi") +app.register_blueprint(cosmic.cosmic_module, url_prefix="/api/cosmic") # allow cors ping.setup_websocket(app) @@ -52,10 +71,11 @@ ping.setup_websocket(app) if __name__ == '__main__': # get args '--prod' args = sys.argv - if 'prod' in args: - app.run("0.0.0.0", port=5000, debug=False) - pass - elif 'debug' in args: + print(os.getcwd()) + # read from ENV Z_PORT + env_port = os.getenv("Z_PORT") + port = 5000 if env_port is None else int(env_port) + if 'debug' in args: app.run("0.0.0.0", port=18200, debug=True) else: - raise Exception("Invalied Mode") + app.run("0.0.0.0", port=port, debug=False) diff --git a/modules/balloon/__init__.py b/modules/balloon/__init__.py index 171d098..80a51df 100644 --- a/modules/balloon/__init__.py +++ b/modules/balloon/__init__.py @@ -10,19 +10,13 @@ from quart import jsonify, request, send_file from modules.balloon.gravityw_plot_perday import render_by_mode_single from modules.balloon.gravityw_plot_year import get_all_modes, render_based_on_mode from modules.balloon.utils import * +from quart.utils import run_sync BASE_PATH_BALLOON = DATA_BASEPATH.balloon all_year_data = pd.read_parquet(f"{BASE_PATH_BALLOON}/ballon_data_lin.parquet") -def get_dataframe_between_year(start_year, end_year): - res = all_year_data - filtered_res = res[(res['file_name'].str.extract(r'LIN-(\d{4})')[0].astype(int) >= start_year) & - (res['file_name'].str.extract(r'LIN-(\d{4})')[0].astype(int) <= end_year)] - return filtered_res - - balloon_module = Blueprint("Balloon", __name__) @@ -60,8 +54,10 @@ async def render_full_year(): end_year = request.args.get('end_year') mode = request.args.get('mode') season = request.args.get('season') + station = request.args.get('station') + print(start_year, end_year) - df = get_dataframe_between_year(int(start_year), int(end_year)) + df = await run_sync(get_ballon_full_df_by_year)(int(start_year), int(end_year), station, False) buff = render_based_on_mode(df, mode, season) return await send_file(buff, mimetype='image/png') diff --git a/modules/balloon/extract_wave.py b/modules/balloon/extract_wave.py index 6ab7e40..c48f6a8 100644 --- a/modules/balloon/extract_wave.py +++ b/modules/balloon/extract_wave.py @@ -30,6 +30,27 @@ def get_top_peaks(pgram, ff, num_peaks=3): top_peaks_frequencies = peak_frequencies[sorted_indices] top_peaks_values = peak_values[sorted_indices] return top_peaks_frequencies, top_peaks_values +# 定义一个函数来计算 RSD + + +def calculate_rsd(x, y, z): + + # todo: 由波数计算波长 + # lamde 垂直波长，只不过单位由10e-3 rad/m变为/km，值差不多 + λ1 = round(1 / ((x / 2.5) * 0.4), 2) + λ2 = round(1 / ((y / 2.5) * 0.4), 2) + λ3 = round(1 / ((z / 2.5) * 0.4), 2) + + # 计算均值 + mean = round(sum([λ1, λ2, λ3]) / 3, 2) + + # 计算相对标准差 + rsd1 = abs(round(((λ1 - mean) / λ1) * 100, 2)) + rsd2 = abs(round(((λ2 - mean) / λ2) * 100, 2)) + rsd3 = abs(round(((λ3 - mean) / λ3) * 100, 2)) + + # 判断是否满足 RSD 小于 20%并返回结果和均值 + return rsd1 < 20 and rsd2 < 20 and rsd3 < 20, mean def calculate_n(Av, Au, Bu, Bv): @@ -148,29 +169,49 @@ def calculate_wave(data: pd.DataFrame, lat: float, g: float): # 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) + # 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) - # 由波数计算波长 - λ1 = round(1 / ((main_frequency_temp / 2.5) * 0.4), 2) - λ2 = round(1 / ((main_frequency_ucmp / 2.5) * 0.4), 2) - λ3 = round(1 / ((main_frequency_vcmp / 2.5) * 0.4), 2) + # # 由波数计算波长 + # λ1 = round(1 / ((main_frequency_temp / 2.5) * 0.4), 2) + # λ2 = round(1 / ((main_frequency_ucmp / 2.5) * 0.4), 2) + # λ3 = round(1 / ((main_frequency_vcmp / 2.5) * 0.4), 2) - mean = round(sum([λ1, λ2, λ3]) / 3, 2) - 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)) + # mean = round(sum([λ1, λ2, λ3]) / 3, 2) + # 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 + # if not (rsd_temp < 20 and rsd_ucmp < 20 and rsd_vcmp < 20): + # return [] - # new + # begin new + + peaks_ucmp_frequencies, _ = get_top_peaks(pgram_ucmp, ff) + peaks_vcmp_frequencies, _ = get_top_peaks(pgram_vcmp, ff) + peaks_temp_frequencies, _ = get_top_peaks(pgram_temp, ff) + + # 获取所有可能的频率组合（9 种组合） + combinations = list(product(peaks_ucmp_frequencies, + peaks_vcmp_frequencies, + peaks_temp_frequencies)) + result = "0" # 初始结果为 "0" + for combo in combinations: + x, y, z = combo # 获取每个组合的三个频率 + valid_rsd, mean = calculate_rsd(x, y, z) # 计算 RSD，并获取均值 + if valid_rsd: # 如果有任何组合满足 RSD 小于 20%，设置 result 为 "1" + result = "1" + break # 找到一个符合条件的组合后，退出循环，避免继续检查其他组合 + if result != "1": + return [] + + # end new # END CHANGE - if not (rsd_temp < 20 and rsd_ucmp < 20 and rsd_vcmp < 20): - return [] - # 定义正弦波模型函数 + def sine_wave(height, A, B, omega): omega = 2 * np.pi / mean return A * np.sin(omega * height + B) diff --git a/modules/balloon/gravityw_plot_perday.py b/modules/balloon/gravityw_plot_perday.py index 7d38326..d4a39dd 100644 --- a/modules/balloon/gravityw_plot_perday.py +++ b/modules/balloon/gravityw_plot_perday.py @@ -18,16 +18,16 @@ def plot_polynomial_fit(ucmp, vcmp, temp, height, poly_ucmp, poly_vcmp, poly_tem plt.subplot(1, 3, 1) plt.plot(ucmp, height, linestyle="-") plt.plot(poly_ucmp(height), height) - plt.ylabel("Height(km)") - plt.xlabel("Zonal wind (m/s)") + plt.ylabel("高度(km)") + plt.xlabel("纬向风(m/s)") plt.grid(True) # v 风扰动量 plt.subplot(1, 3, 2) plt.plot(vcmp, height, linestyle="-") plt.plot(poly_vcmp(height), height) - plt.ylabel("Height(km)") - plt.xlabel("Meridional wind (m/s)") + plt.ylabel("高度(km)") + plt.xlabel("经向风 (m/s)") plt.grid(True) plt.title("观测的二阶多项式拟合", fontsize=16) @@ -36,8 +36,8 @@ def plot_polynomial_fit(ucmp, vcmp, temp, height, poly_ucmp, poly_vcmp, poly_tem plt.subplot(1, 3, 3) plt.plot(temp, height, linestyle="-") plt.plot(poly_temp(height), height) - plt.ylabel("Height(km)") - plt.xlabel("Temperature(K)") + plt.ylabel("高度(km)") + plt.xlabel("温度(K)") plt.grid(True) plt.tight_layout() @@ -50,16 +50,16 @@ def plot_sine_fit(residual_ucmp, residual_vcmp, residual_temp, height, u_fit, v_ plt.subplot(1, 3, 1) plt.plot(residual_ucmp, height, marker="o", linestyle="None") plt.plot(u_fit, height) - plt.ylabel("Height(km)") - plt.xlabel("Zonal wind (m/s)") + plt.ylabel("高度(km)") + plt.xlabel("纬向风 (m/s)") plt.grid(True) # v 风扰动量 plt.subplot(1, 3, 2) plt.plot(residual_vcmp, height, marker="o", linestyle="None") plt.plot(v_fit, height) - plt.ylabel("Height(km)") - plt.xlabel("Meridional wind (m/s)") + plt.ylabel("高度(km)") + plt.xlabel("经向风(m/s)") plt.grid(True) plt.title("扰动分量的正弦波拟合", fontsize=16) @@ -68,8 +68,8 @@ def plot_sine_fit(residual_ucmp, residual_vcmp, residual_temp, height, u_fit, v_ plt.subplot(1, 3, 3) plt.plot(residual_temp, height, marker="o", linestyle="None") plt.plot(T_fit, height) - plt.ylabel("Height(km)") - plt.xlabel("Temperature(K)") + plt.ylabel("高度(km)") + plt.xlabel("温度(K)") plt.grid(True) plt.tight_layout() @@ -87,8 +87,8 @@ def plot_uv_vector(u_fit, v_fit, height, specified_heights, markers): # plt.ylim(-4, 4) plt.axvline(0, color="gray", linestyle="--") plt.axhline(0, color="gray", linestyle="--") - plt.xlabel("Zonal Wind (m/s)") - plt.ylabel("Meridional Wind (m/s)") + plt.xlabel("经向风 (m/s)") + plt.ylabel("纬向风(m/s)") plt.legend() plt.grid(True) plt.title("径向风-纬向风矢量图") @@ -102,8 +102,8 @@ def plot_temp_horizontal_wind(uh, T_fit, height, specified_heights, markers): plt.scatter(uh[index], T_fit[index], marker=marker, s=100, label=f"{h} km") - plt.xlim(-4, 4) - plt.ylim(-2, 2) + # plt.xlim(-4, 4) + # plt.ylim(-2, 2) plt.axvline(0, color="gray", linestyle="--") plt.axhline(0, color="gray", linestyle="--") plt.xlabel("Horizontal wind (m/s)") diff --git a/modules/balloon/utils.py b/modules/balloon/utils.py index d4db2a5..15602d7 100644 --- a/modules/balloon/utils.py +++ b/modules/balloon/utils.py @@ -1,4 +1,5 @@ import glob +import os import re import pandas as pd @@ -6,6 +7,7 @@ import pandas as pd import time import logging +from CONSTANT import DATA_BASEPATH from modules.balloon.extract_wave import extract_wave, is_terrain_wave from modules.balloon.read_data import read_data @@ -60,51 +62,81 @@ comboDate = [ ] +def get_dataframe_between_year(all_year_data, start_year, end_year, station): + res = all_year_data + filtered_res = res[(res['file_name'].str.extract(rf'{station}-(\d{{4}})')[0].astype(int) >= start_year) & + (res['file_name'].str.extract(rf'{station}-(\d{{4}})')[0].astype(int) <= end_year)] + return filtered_res + + def get_ballon_files(): try: - data = glob.glob("data/探空气球/**/*.nc", recursive=True) + data = glob.glob(f"{DATA_BASEPATH.balloon}/**/*.nc", recursive=True) except FileNotFoundError: return [] return data -all_ballon_files = get_ballon_files() - - -def get_ballon_path_by_year(start_year, end_year): +def get_ballon_path_by_year(start_year, end_year, station=None): + all_ballon_files = get_ballon_files() return list(filter( - lambda x: any(f"LIN-{year}" in x for year in range( + lambda x: any(f"{station}-{year}" in x for year in range( start_year, end_year + 1)), all_ballon_files )) -def get_ballon_full_df_by_year(start_year, end_year): +def get_ballon_full_df_by_year(start_year, end_year, station=None, ignore_cache=False): # Set up logging logging.basicConfig(level=logging.DEBUG, format='%(asctime)s %(levelname)s: %(message)s') + if (station is None): + raise ValueError("Station is required") + + cache_dir = f"{DATA_BASEPATH.balloon}/cache/{station}-b{start_year}-e{end_year}.parquet" + cache_base_dir = f"{DATA_BASEPATH.balloon}/cache/" + if os.path.exists(cache_dir) and not ignore_cache: + logging.debug(f"Reading cache from {cache_dir}") + return pd.read_parquet(cache_dir) + # check if there is any file begin with cache base dir, or say ,filename begin with station + # if so, read the file and return + if os.path.exists(cache_base_dir): + for file in os.listdir(cache_base_dir): + # should match {station}-bxxxx-exxxx.parquet, xxxx is year + if re.match(f"{station}-b\d{{4}}-e\d{{4}}.parquet", file): + # extract begin year and end year from filename + b_year = int(file.split("-")[1][1:]) + e_year = int(file.split("-")[2][1:].split(".")[0]) + logging.debug( + f"Found cache file {file}, b_year={b_year}, e_year={e_year}") + if b_year <= start_year and e_year >= end_year: + wider_df = pd.read_parquet(f"{cache_base_dir}/{file}") + return get_dataframe_between_year(wider_df, start_year, end_year, station) + start_time = time.time() logging.debug( f"Starting get_ballon_full_df_by_year with start_year={start_year}, end_year={end_year}") # Timing the path retrieval t0 = time.time() - paths = get_ballon_path_by_year(start_year, end_year) + paths = get_ballon_path_by_year(start_year, end_year, station) + # if no path found, raise + if len(paths) == 0: + raise ValueError( + f"No balloon data found for {station} between {start_year} and {end_year}") t1 = time.time() logging.debug(f"Retrieved {len(paths)} paths in {t1 - t0:.2f} seconds") # optimization: add cache. only select need to be reprocessed - with open("./cache/ballon_lin_has_wave", "r") as f: - cache_has_waves = f.readlines() - cache_has_waves = [x.strip() for x in cache_has_waves] + # with open("./cache/ballon_lin_has_wave", "r") as f: + # cache_has_waves = f.readlines() + # cache_has_waves = [x.strip() for x in cache_has_waves] year_df = pd.DataFrame() for idx, file in enumerate(paths, 1): - if len(cache_has_waves) > 0 and file not in cache_has_waves: - logging.debug(f"Skipping {file} as it has no wave data") - continue + file_start_time = time.time() logging.debug(f"Processing file {idx}/{len(paths)}: {file}") @@ -130,11 +162,12 @@ def get_ballon_full_df_by_year(start_year, end_year): continue # Determine terrain wave - c = is_terrain_wave(data, lat, g) - terrain_time = time.time() - - year_pattern = r"products-RS92-GDP.2-LIN-(\d{4})" - year = int(re.search(year_pattern, file).group(1)) + try: + c = is_terrain_wave(data, lat, g) + terrain_time = time.time() + except Exception as e: + logging.error(f"Error determining terrain wave from {file}: {e}") + continue logging.debug( f"Determined terrain wave in {terrain_time - extract_time:.2f} seconds") @@ -157,6 +190,15 @@ def get_ballon_full_df_by_year(start_year, end_year): total_time = time.time() - start_time logging.debug( f"Completed get_ballon_full_df_by_year in {total_time:.2f} seconds") + year_df['hori_wave_len'] = year_df['hori_wave_len'].apply( + lambda x: float(x) if x != ' ' else None) + + # save to cache + # if parent folder does not exist, create it + if not os.path.exists(os.path.dirname(cache_dir)): + os.makedirs(os.path.dirname(cache_dir)) + year_df.to_parquet(cache_dir) + return year_df diff --git a/modules/cosmic/__init__.py b/modules/cosmic/__init__.py index 9f46409..eb9d910 100644 --- a/modules/cosmic/__init__.py +++ b/modules/cosmic/__init__.py @@ -10,6 +10,7 @@ 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 +from modules.cosmic.planetw_perday import cosmic_planetw_plot_perday cosmic_module = Blueprint("Cosmic", __name__) @@ -48,6 +49,31 @@ async def render(): return await send_file(buf, mimetype="image/png") +@cosmic_module.route('/render/planet_wave/perday') +async def render_by_mode_single(): + year = request.args.get("year", 2008) + T = request.args.get("T", 16) + 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_plot_perday)( + temp_df, + T=int(T), + ) + buf = BytesIO() + + plt.savefig(buf, format="png") + buf.seek(0) + return await send_file(buf, mimetype="image/png") + + @cosmic_module.route('/render/gravity_wave/perday') @cosmic_module.route('/render/single') async def single_render(): diff --git a/modules/cosmic/gravityw_multiday.py b/modules/cosmic/gravityw_multiday.py index bbb9f23..95ff830 100644 --- a/modules/cosmic/gravityw_multiday.py +++ b/modules/cosmic/gravityw_multiday.py @@ -66,10 +66,9 @@ def plot_heatmap(data, heights, title): # 绘制热力图，数据中的行代表高度，列代表天数 sns.heatmap(data, cmap='coolwarm', xticklabels=1, yticklabels=50, cbar_kws={'label': 'Value'}) - plt.xlabel('Day') - plt.ylabel('Height (km)') + plt.xlabel('天') + plt.ylabel('高度 (km)') plt.title(title) - # 设置x轴的刻度，使其每30天显示一个标签 num_days = data.shape[1] x_tick_positions = np.arange(0, num_days, 30) x_tick_labels = np.arange(0, num_days, 30) @@ -98,13 +97,13 @@ class GravityMultidayPlot: # 调用函数绘制热力图 data = self.data plot_heatmap(data, self.heights, - 'Heatmap of final_mean_ktemp_Nz（10^(-4)）') + '"最终平均温度场的热图 (单位：10^(-4))，随浮力频率Nz变化') def plot_heatmap_tempPtz(self): # -------------绘制重力势能年统计图------------------------------------------------ data = self.df_final_mean_ktemp_Ptz.T plot_heatmap(data, self.heights, - 'Heatmap of final_mean_ktemp_Ptz(J/kg)') + '最终平均动能温度-位势温度的热图 (单位：焦耳/千克)') def plot_monthly_tempNz(self): # ------------------------绘制月统计图--------------------------------------------------------------------------------- diff --git a/modules/cosmic/gravityw_perday.py b/modules/cosmic/gravityw_perday.py index 94c8cd1..7989e51 100644 --- a/modules/cosmic/gravityw_perday.py +++ b/modules/cosmic/gravityw_perday.py @@ -359,8 +359,8 @@ def plot_results_mean_ktemp_Nz(mean_ktemp_Nz, heights): plt.clf() plt.figure(figsize=(10, 6)) plt.plot(mean_ktemp_Nz, heights / 1000) # 高度单位换算为km，方便展示 - plt.xlabel('Average (N_z)') - plt.ylabel('H(km)') + plt.xlabel('平均浮力频率 (N_z)') + plt.ylabel('高度(km)') # plt.gca().invert_yaxis() # 使高度坐标轴从上到下递增，符合常规习惯 # plt.show() @@ -374,8 +374,8 @@ def plot_results_mean_ktemp_Ptz(mean_ktemp_Ptz, heights): plt.clf() plt.figure(figsize=(10, 6)) plt.plot(mean_ktemp_Ptz, heights / 1000) # 高度单位换算为km，方便展示 - plt.xlabel('Average (PT_z)') - plt.ylabel('H (km)') + plt.xlabel('平均势能 (PT_z)') + plt.ylabel('高度(km)') # plt.gca().invert_yaxis() # 使高度坐标轴从上到下递增，符合常规习惯 # plt.show() diff --git a/modules/cosmic/planetw_daily.py b/modules/cosmic/planetw_daily.py index f751bf9..5795992 100644 --- a/modules/cosmic/planetw_daily.py +++ b/modules/cosmic/planetw_daily.py @@ -108,7 +108,7 @@ def cosmic_planetw_daily_plot( plt.figure(figsize=(8, 6)) # 创建新的图形 plt.plot(x_values, fit_df[col].values) plt.title(f'k={k} 振幅图') - plt.xlabel('Day') + plt.xlabel('天') plt.ylabel('振幅') # 设置横坐标的动态调整 diff --git a/modules/cosmic/planetw_daily_process.py b/modules/cosmic/planetw_daily_process.py index 220535c..3587dd1 100644 --- a/modules/cosmic/planetw_daily_process.py +++ b/modules/cosmic/planetw_daily_process.py @@ -24,7 +24,7 @@ def cosmic_planet_daily_process( return pd.read_parquet(cache_path) # 遍历文件夹序号1到365 - for i in range(1, 365): + for i in range(1, 165): # 根据i的值调整文件夹名称 if i < 10: folder_name = f"atmPrf_repro2021_2008_00{i}" # 一位数，前面加两个0 @@ -39,49 +39,63 @@ def cosmic_planet_daily_process( # 检查文件夹是否存在 if os.path.exists(folder_path): # 遍历文件夹中的文件 + folder_level_cache_path = f"{folder_path}/planet_{target_h}.parquet" + if os.path.exists(folder_level_cache_path): + result_folder_level = pd.read_parquet(folder_level_cache_path) + dfs.append(result_folder_level) + continue + result_folder_level = [] 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}") + if not file_name.endswith('.0390_nc'): + continue + + 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) + result_folder_level.append(df_filtered) + except Exception as e: + print(f"处理文件 {finfo} 时出错: {e}") + # save to parquet + result_folder_level = pd.concat( + result_folder_level, axis=0, ignore_index=True) + result_folder_level.to_parquet(folder_level_cache_path) else: print(f"文件夹 {folder_path} 不存在。") @@ -125,8 +139,11 @@ def cosmic_planet_daily_process( return total_minutes / 1440 # 应用函数转换 - closest_per_day['Time'] = closest_per_day['Datetime_str'].apply( - convert_to_days) + try: + closest_per_day['Time'] = closest_per_day['Datetime_str'].apply( + convert_to_days) + except Exception as e: + print(e) # 时间t day_data = closest_per_day['Time'] # 计算背景温度，时间Day和弧度经度的函数 diff --git a/modules/cosmic/planetw_perday.py b/modules/cosmic/planetw_perday.py new file mode 100644 index 0000000..74a140f --- /dev/null +++ b/modules/cosmic/planetw_perday.py @@ -0,0 +1,79 @@ +import pandas as pd +import numpy as np +from scipy.optimize import curve_fit +import matplotlib.pyplot as plt +# 解决绘图中中文不能显示的问题 +import matplotlib +matplotlib.rcParams['font.sans-serif'] = ['SimHei'] # 显示中文 +matplotlib.rcParams['axes.unicode_minus'] = False # 正常显示负号 +# 读取一年的数据文件 +# 设置初始参数 +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] # 上界 +) +# 定义拟合函数 + + +def cosmic_planetw_plot_perday( + df: pd.DataFrame, + T=16, + +): + + def u_func(x, *params, t): + 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) + ) + + # 设置最小数据量的阈值 + min_data_points = 36 + + # 进行多个时间窗口的拟合 + 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)] + + # 检查当前窗口的数据量 + if len(df_8) < min_data_points: + print(f"数据量不足，无法拟合：{start_day} 到 {end_day}，数据点数量：{len(df_8)}") + continue # 跳过当前时间窗口，继续下一个窗口 + + # 提取时间、经度、温度数据 + t = np.array(df_8['Time']) # 时间 + x = np.array(df_8['Longitude_Radians']) # 经度弧度制 + temperature = np.array(df_8['Temperature']) # 温度，因变量 + + # 用T进行拟合 + popt, pcov = curve_fit(lambda x, *params: u_func(x, *params, t=t), + x, temperature, p0=initial_guess, bounds=bounds, maxfev=50000) + + # 绘制拟合曲线 + t_fixed = (start_day + end_day) / 2 # 窗口的中间时间 + x_fit = np.linspace(min(x), max(x), 100) # 生成用于绘制拟合曲线的x值 + y_fit = u_func(x_fit, *popt, t=t_fixed) # 计算拟合曲线的y值 + + plt.figure(figsize=(10, 6)) + plt.plot(x_fit, y_fit, label='拟合曲线', color='red', linewidth=2) + plt.xlabel('经度（弧度制）') + plt.ylabel('温度') + plt.title(f'时间窗口：{start_day} 到 {end_day}') + plt.legend() + plt.grid(True) diff --git a/modules/saber/__init__.py b/modules/saber/__init__.py index 8971971..8463bcc 100644 --- a/modules/saber/__init__.py +++ b/modules/saber/__init__.py @@ -3,6 +3,7 @@ from io import BytesIO from quart import Blueprint, request, send_file from matplotlib import pyplot as plt from CONSTANT import DATA_BASEPATH +from modules.saber.planetw_perday import saber_planetw_plot_perday from modules.saber.planetw_plot import saber_planetw_plot from modules.saber.planetw_process import saber_planetw_process from modules.saber.gravityw_process import SaberGravitywProcessor @@ -164,3 +165,24 @@ async def do_gravity_year(): ) await run_sync(saber_planetw_plot)(data, T, k=int(k)) return await extract_payload() + + +@saber_module.route("/render/planet_wave/per_day") +async def do_pl_year(): + year = request.args.get("year") + T = request.args.get("T") + # k = request.args.get("k") + H = request.args.get("H") + lat_target = request.args.get("lat_target") + start_day = request.args.get("start_day", 0) + 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, + target_h=int(H), + target_lat=int(lat_target) + ) + await run_sync(saber_planetw_plot_perday)(data, T, int(start_day)) + return await extract_payload() diff --git a/modules/saber/gravityw_render.py b/modules/saber/gravityw_render.py index 66569f8..52ccc95 100644 --- a/modules/saber/gravityw_render.py +++ b/modules/saber/gravityw_render.py @@ -157,14 +157,14 @@ class SaberGravitywRenderer: plt.subplot(1, 2, 1) plt.plot(x1[::-1], y, label='原始信号') plt.title('（a）Nz') - plt.xlabel('Nz', labelpad=10) # 增加标签间距 + plt.xlabel('势能', labelpad=10) # 增加标签间距 plt.ylabel('高度 (km)', labelpad=10) # 增加标签间距 # 原始信号的时间序列 plt.subplot(1, 2, 2) plt.plot(x2[::-1], y, label='原始信号') plt.title('（b）Ptz') - plt.xlabel('Ptz', labelpad=10) # 增加标签间距 + plt.xlabel('浮力频率', labelpad=10) # 增加标签间距 plt.ylabel('高度 (km)', labelpad=10) # 增加标签间距 # 调整子图之间的边距 diff --git a/modules/saber/planetw_perday.py b/modules/saber/planetw_perday.py new file mode 100644 index 0000000..151edad --- /dev/null +++ b/modules/saber/planetw_perday.py @@ -0,0 +1,86 @@ +import pandas as pd +import numpy as np +from scipy.optimize import curve_fit +import matplotlib.pyplot as plt +# 解决绘图中中文不能显示的问题 +import matplotlib +matplotlib.rcParams['font.sans-serif'] = ['SimHei'] # 显示中文 +matplotlib.rcParams['axes.unicode_minus'] = False # 正常显示负号 +# 读取一年的数据文件 +# 设置初始参数 +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] # 上界 +) +# 定义拟合函数 + + +def saber_planetw_plot_perday( + df: pd.DataFrame, + T=16, + start_day=0 + +): + + def u_func(x, *params, t): + 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) + ) + + # 设置最小数据量的阈值 + min_data_points = 36 + + # 进行多个时间窗口的拟合 + # 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)] + + # 检查当前窗口的数据量 + if len(df_8) < min_data_points: + print(f"数据量不足，无法拟合：{start_day} 到 {end_day}，数据点数量：{len(df_8)}") + plt.figure(figsize=(10, 6)) + # write nothing but a title + plt.title(f'时间窗口：{start_day} 到 {end_day}') + plt.grid(True) + # add a label, says unable to fit + plt.text(0.5, 0.5, "数据量不足，无法拟合", ha='center', va='center', + transform=plt.gca().transAxes, fontsize=20) + return + # 提取时间、经度、温度数据 + t = np.array(df_8['Time']) # 时间 + x = np.array(df_8['Longitude_Radians']) # 经度弧度制 + temperature = np.array(df_8['Temperature']) # 温度，因变量 + + # 用T进行拟合 + popt, pcov = curve_fit(lambda x, *params: u_func(x, *params, t=t), + x, temperature, p0=initial_guess, bounds=bounds, maxfev=50000) + + # 绘制拟合曲线 + t_fixed = (start_day + end_day) / 2 # 窗口的中间时间 + x_fit = np.linspace(min(x), max(x), 100) # 生成用于绘制拟合曲线的x值 + y_fit = u_func(x_fit, *popt, t=t_fixed) # 计算拟合曲线的y值 + + plt.figure(figsize=(10, 6)) + plt.plot(x_fit, y_fit, label='拟合曲线', color='red', linewidth=2) + plt.xlabel('经度（弧度制）') + plt.ylabel('温度') + plt.title(f'时间窗口：{start_day} 到 {end_day}') + plt.legend() + plt.grid(True) diff --git a/modules/saber/planetw_plot.py b/modules/saber/planetw_plot.py index 11690f2..5cfd395 100644 --- a/modules/saber/planetw_plot.py +++ b/modules/saber/planetw_plot.py @@ -105,7 +105,7 @@ def saber_planetw_plot( col = k_to_a[f'k={k}'] plt.plot(x_values, fit_df[col].values) plt.suptitle(f'k = {k} 振幅图') - plt.xlabel('Day') + plt.xlabel('天') plt.ylabel('振幅') # 设置横坐标的动态调整 diff --git a/modules/saber/planetw_process.py b/modules/saber/planetw_process.py index 65eae00..558fc4c 100644 --- a/modules/saber/planetw_process.py +++ b/modules/saber/planetw_process.py @@ -23,7 +23,7 @@ months = [ def saber_planetw_process( year, target_h=70, - target_lat=60 + target_lat=60 ): # 定义文件路径模板和年份 diff --git a/modules/tidi/gravity_wave_year.py b/modules/tidi/gravity_wave_year.py index b1953fd..a7ff005 100644 --- a/modules/tidi/gravity_wave_year.py +++ b/modules/tidi/gravity_wave_year.py @@ -869,8 +869,8 @@ class TidiGravityWPlotMonthly: labels=months[::interval], rotation=45) # rotation旋转可不加 # 添加轴标签 - plt.xlabel('Month') # X轴标签 - plt.ylabel('Height') # Y轴标签 + plt.xlabel('月') # X轴标签 + plt.ylabel('高度') # Y轴标签 # 显示图形 # plt.show() @@ -920,9 +920,9 @@ class TidiGravityWPlotMonthly: plt.plot(months, monthly_average, marker='o', linestyle='-', color='b') # 添加标题和标签 - plt.title("Monthly Average ENERGY(log)") - plt.xlabel("Month") - plt.ylabel("Average Energy") + plt.title("月平均能量(结果取log)") + plt.xlabel("月") + plt.ylabel("平均能量") # 显示图表 plt.xticks(rotation=45) # 让月份标签更清晰可读 plt.grid(True) diff --git a/modules/tidi/planet_wave_daily.py b/modules/tidi/planet_wave_daily.py index c57caf0..981b23a 100644 --- a/modules/tidi/planet_wave_daily.py +++ b/modules/tidi/planet_wave_daily.py @@ -127,7 +127,7 @@ def TidiPlotPlanetWDaily( plt.figure(figsize=(8, 6)) # 创建新的图形 plt.plot(x_values, fit_df[col].values) plt.title(f'k={k} 振幅图') - plt.xlabel('Day') + plt.xlabel('天') plt.ylabel('振幅') # 设置横坐标的动态调整 diff --git a/ping.py b/ping.py index 44a5689..572acaf 100644 --- a/ping.py +++ b/ping.py @@ -19,7 +19,7 @@ class WebSocketConfig: def setup_websocket( app: Quart, - url_prefix: str = '/ping/ws', + url_prefix: str = '/api/ping/ws', config: Optional[WebSocketConfig] = None ) -> None: """