diff --git a/CONSTANT.py b/CONSTANT.py index 93740cd..8b6a11a 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 3480505..8735f5d 100644 --- a/backend.py +++ b/backend.py @@ -1,10 +1,10 @@ import resource from matplotlib import pyplot as plt -import cosmic -import tidi -import saber -import radar -import balloon +from modules import cosmic +from modules import tidi +from modules import saber +from modules import radar +from modules import balloon from quart import Quart, request from quart_cors import cors from typing import get_args diff --git a/modules/__init__.py b/modules/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/balloon/__init__.py b/modules/balloon/__init__.py similarity index 79% rename from balloon/__init__.py rename to modules/balloon/__init__.py index 688c56d..5af50b5 100644 --- a/balloon/__init__.py +++ b/modules/balloon/__init__.py @@ -2,14 +2,11 @@ import base64 from urllib import response from quart import Blueprint from CONSTANT import DATA_BASEPATH -import balloon.extract_wave -import balloon.read_data -import balloon.plot_once -import balloon.plot_year +import modules.balloon as balloon from quart import jsonify, request, send_file -from balloon.plot_once_backend import render_by_mode_single -from balloon.plot_year_backend import get_all_modes, render_based_on_mode -from balloon.utils import * +from modules.balloon.plot_once_backend import render_by_mode_single +from modules.balloon.plot_year_backend import get_all_modes, render_based_on_mode +from modules.balloon.utils import * BASE_PATH_BALLOON = DATA_BASEPATH.balloon @@ -23,13 +20,6 @@ def get_dataframe_between_year(start_year, end_year): return filtered_res -plot_once = balloon.plot_once.plot_once -plot_year = balloon.plot_year.plot_year -is_terrain_wave = balloon.extract_wave.is_terrain_wave -read_data = balloon.read_data.read_data -extract_wave = balloon.extract_wave.extract_wave - - balloon_module = Blueprint("Balloon", __name__) diff --git a/balloon/extract_wave.py b/modules/balloon/extract_wave.py similarity index 100% rename from balloon/extract_wave.py rename to modules/balloon/extract_wave.py diff --git a/balloon/plot_once.py b/modules/balloon/plot_once.py similarity index 85% rename from balloon/plot_once.py rename to modules/balloon/plot_once.py index e990a21..fb94485 100644 --- a/balloon/plot_once.py +++ b/modules/balloon/plot_once.py @@ -2,17 +2,20 @@ import matplotlib.pyplot as plt import numpy as np import pandas as pd import numpy as np -from balloon.extract_wave import calculate_wave, is_terrain_wave +from modules.balloon.extract_wave import calculate_wave, is_terrain_wave def plot_once(data: pd.DataFrame, path: str, lat: float, g: float): - wave = calculate_wave(data[(data["alt"] >= 15) & (data["alt"] <= 25)], lat, g) + wave = calculate_wave( + data[(data["alt"] >= 15) & (data["alt"] <= 25)], lat, g) if len(wave) == 0: return [] c = is_terrain_wave(data, lat, g) - a, b, omega_upper, w_f, λ_z, λ_h, c_x, c_y, c_z, Ek, E_p, MFu1, MFv1, params_u, params_v, params_T = wave[:16] - ucmp, vcmp, temp, height, poly_ucmp, poly_vcmp, _, poly_temp, residual_ucmp, residual_vcmp, residual_temp, u_fit, v_fit, T_fit, uh, _ = wave[16:] + a, b, omega_upper, w_f, λ_z, λ_h, c_x, c_y, c_z, Ek, E_p, MFu1, MFv1, params_u, params_v, params_T = wave[ + :16] + ucmp, vcmp, temp, height, poly_ucmp, poly_vcmp, _, poly_temp, residual_ucmp, residual_vcmp, residual_temp, u_fit, v_fit, T_fit, uh, _ = wave[ + 16:] plt.figure(figsize=(16, 10)) @@ -81,7 +84,8 @@ def plot_once(data: pd.DataFrame, path: str, lat: float, g: float): plt.plot(u_fit, v_fit) # 绘制U-V曲线 for h, marker in zip(specified_heights, markers): index = np.abs(height - h).argmin() # 找到离指定高度最近的索引 - plt.scatter(u_fit[index], v_fit[index], marker=marker, s=100, label=f"{h} km") + plt.scatter(u_fit[index], v_fit[index], + marker=marker, s=100, label=f"{h} km") # 设置坐标轴范围和比例 plt.xlim(-8, 8) @@ -93,7 +97,8 @@ def plot_once(data: pd.DataFrame, path: str, lat: float, g: float): plt.ylabel("Meridional Wind (m/s)") plt.legend() # 显示图例 plt.grid(True) # 显示网格线 - plt.text(0.05, 1, "(a)", transform=plt.gca().transAxes, fontsize=14, verticalalignment="top") + plt.text(0.05, 1, "(a)", transform=plt.gca().transAxes, + fontsize=14, verticalalignment="top") plt.title("径向风-纬向风矢量图") # b. 水平风-温度的矢端曲线 @@ -101,7 +106,8 @@ def plot_once(data: pd.DataFrame, path: str, lat: float, g: float): plt.plot(uh, T_fit) # 绘制水平风-温度曲线 for h, marker in zip(specified_heights, markers): index = np.abs(height - h).argmin() # 找到离指定高度最近的索引 - plt.scatter(uh[index], T_fit[index], marker=marker, s=100, label=f"{h} km") + plt.scatter(uh[index], T_fit[index], + marker=marker, s=100, label=f"{h} km") # 设置坐标轴范围和比例 plt.xlim(-4, 4) @@ -113,7 +119,8 @@ def plot_once(data: pd.DataFrame, path: str, lat: float, g: float): plt.ylabel("Temp (K)") plt.legend() plt.grid(True) - plt.text(0.05, 1, "(b)", transform=plt.gca().transAxes, fontsize=14, verticalalignment="top") + plt.text(0.05, 1, "(b)", transform=plt.gca().transAxes, + fontsize=14, verticalalignment="top") plt.title("温度-水平风矢量图") # 设置中文显示和负号正常显示 diff --git a/balloon/plot_once_backend.py b/modules/balloon/plot_once_backend.py similarity index 95% rename from balloon/plot_once_backend.py rename to modules/balloon/plot_once_backend.py index baa8e27..131559b 100644 --- a/balloon/plot_once_backend.py +++ b/modules/balloon/plot_once_backend.py @@ -3,7 +3,7 @@ from io import BytesIO import matplotlib.pyplot as plt import numpy as np -from balloon.extract_wave import calculate_wave, is_terrain_wave +from modules.balloon.extract_wave import calculate_wave, is_terrain_wave lat = 52.21 g = 9.76 diff --git a/balloon/plot_year.py b/modules/balloon/plot_year.py similarity index 100% rename from balloon/plot_year.py rename to modules/balloon/plot_year.py diff --git a/balloon/plot_year_backend.py b/modules/balloon/plot_year_backend.py similarity index 100% rename from balloon/plot_year_backend.py rename to modules/balloon/plot_year_backend.py diff --git a/balloon/read_data.py b/modules/balloon/read_data.py similarity index 100% rename from balloon/read_data.py rename to modules/balloon/read_data.py diff --git a/balloon/utils.py b/modules/balloon/utils.py similarity index 91% rename from balloon/utils.py rename to modules/balloon/utils.py index d60a71f..d4db2a5 100644 --- a/balloon/utils.py +++ b/modules/balloon/utils.py @@ -3,10 +3,12 @@ import re import pandas as pd -import balloon import time import logging +from modules.balloon.extract_wave import extract_wave, is_terrain_wave +from modules.balloon.read_data import read_data + filter_columns = [ "file_name", "c", @@ -107,14 +109,14 @@ def get_ballon_full_df_by_year(start_year, end_year): logging.debug(f"Processing file {idx}/{len(paths)}: {file}") # Read data - data = balloon.read_data(file) + data = read_data(file) read_time = time.time() logging.debug( f"Read data in {read_time - file_start_time:.2f} seconds") # Extract wave try: - wave = balloon.extract_wave(data, lat, g) + wave = extract_wave(data, lat, g) extract_time = time.time() logging.debug( f"Extracted wave in {extract_time - read_time:.2f} seconds") @@ -128,7 +130,7 @@ def get_ballon_full_df_by_year(start_year, end_year): continue # Determine terrain wave - c = balloon.is_terrain_wave(data, lat, g) + c = is_terrain_wave(data, lat, g) terrain_time = time.time() year_pattern = r"products-RS92-GDP.2-LIN-(\d{4})" diff --git a/cosmic/__init__.py b/modules/cosmic/__init__.py similarity index 81% rename from cosmic/__init__.py rename to modules/cosmic/__init__.py index 232b376..e82ddb3 100644 --- a/cosmic/__init__.py +++ b/modules/cosmic/__init__.py @@ -2,43 +2,46 @@ from io import BytesIO from quart import Blueprint, request, send_file from matplotlib import pyplot as plt -from cosmic.single import SingleCosmicWavePlot -from cosmic.temp_render import temp_render +from modules.cosmic.single import SingleCosmicWavePlot +from modules.cosmic.temp_render import temp_render from quart.utils import run_sync cosmic_module = Blueprint("Cosmic", __name__) + @cosmic_module.route('/metadata') def get_meta(): return [] + @cosmic_module.route('/temp_render') async def render(): T = request.args.get("T", 16) await run_sync(temp_render)(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/single') async def single_render(): year = request.args.get("year", 2008) day = request.args.get("day", 1) mode = request.args.get("mode", "mean_ktemp_Nz") - - p :SingleCosmicWavePlot = await run_sync(SingleCosmicWavePlot)(year=int(year), day=int(day)) + + p: SingleCosmicWavePlot = await run_sync(SingleCosmicWavePlot)(year=int(year), day=int(day)) if mode == "mean_ktemp_Nz": await run_sync(p.plot_results_mean_ktemp_Nz)() elif mode == "mean_ktemp_Ptz": await run_sync(p.plot_results_mean_ktemp_Ptz)() - else : + else: raise ValueError("Invalid mode") - + buf = BytesIO() plt.savefig(buf, format="png") buf.seek(0) - - return await send_file(buf, mimetype="image/png") \ No newline at end of file + + return await send_file(buf, mimetype="image/png") diff --git a/cosmic/multiple.py b/modules/cosmic/multiple.py similarity index 100% rename from cosmic/multiple.py rename to modules/cosmic/multiple.py diff --git a/cosmic/single.py b/modules/cosmic/single.py similarity index 100% rename from cosmic/single.py rename to modules/cosmic/single.py diff --git a/cosmic/temp_render.py b/modules/cosmic/temp_render.py similarity index 100% rename from cosmic/temp_render.py rename to modules/cosmic/temp_render.py diff --git a/radar/__init__.py b/modules/radar/__init__.py similarity index 85% rename from radar/__init__.py rename to modules/radar/__init__.py index f090d59..b432276 100644 --- a/radar/__init__.py +++ b/modules/radar/__init__.py @@ -4,8 +4,8 @@ from io import BytesIO from quart import Blueprint, request, send_file from matplotlib import pyplot as plt from CONSTANT import DATA_BASEPATH -from radar.plot_original import final_render_v2 -from radar.plot_prod import final_plot_v2 +from modules.radar.plot_original import final_render_v2 +from modules.radar.plot_prod import final_plot_v2 from quart.utils import run_sync BASE_PATH_RADAR = DATA_BASEPATH.radar @@ -37,18 +37,18 @@ async def render_v1(): year = request.args.get('year') H = request.args.get('H') station = request.args.get('station') - + model_name = request.args.get('model_name') mode = request.args.get('mode') - + renderer = await run_sync(final_render_v2)(int(H), int(year), station, wind_type) - if mode=="day": + if mode == "day": day = request.args.get('day') - renderer.render_day(day,model_name) - elif mode=="month": + renderer.render_day(day, model_name) + elif mode == "month": month = request.args.get('month') - renderer.render_month(int(month),model_name) - elif mode=="year": + renderer.render_month(int(month), model_name) + elif mode == "year": renderer.render_year(model_name) else: raise ValueError("mode not supported") @@ -78,7 +78,6 @@ async def render_v2(): month_range = (start_month, end_month) else: month_range = (1, 12) - buffer = await run_sync(final_plot_v2)(int(year), station, model_name, month_range) buffer.seek(0) diff --git a/radar/plot_original.py b/modules/radar/plot_original.py similarity index 100% rename from radar/plot_original.py rename to modules/radar/plot_original.py diff --git a/radar/plot_prod.py b/modules/radar/plot_prod.py similarity index 100% rename from radar/plot_prod.py rename to modules/radar/plot_prod.py diff --git a/saber/__init__.py b/modules/saber/__init__.py similarity index 91% rename from saber/__init__.py rename to modules/saber/__init__.py index c6e6cce..789230a 100644 --- a/saber/__init__.py +++ b/modules/saber/__init__.py @@ -3,11 +3,11 @@ 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 * +from modules.saber.archive.gravity_plot import do_gravity_plot +from modules.saber.archive.gravity_process import process_gravity_data +from modules.saber.process import DataProcessor +from modules.saber.render import Renderer +from modules.saber.utils import * from quart.utils import run_sync diff --git a/saber/archive/gravity_plot.py b/modules/saber/archive/gravity_plot.py similarity index 100% rename from saber/archive/gravity_plot.py rename to modules/saber/archive/gravity_plot.py diff --git a/saber/archive/gravity_process.py b/modules/saber/archive/gravity_process.py similarity index 100% rename from saber/archive/gravity_process.py rename to modules/saber/archive/gravity_process.py diff --git a/saber/archive/legacy.py b/modules/saber/archive/legacy.py similarity index 100% rename from saber/archive/legacy.py rename to modules/saber/archive/legacy.py diff --git a/saber/archive/saber_render.py b/modules/saber/archive/saber_render.py similarity index 97% rename from saber/archive/saber_render.py rename to modules/saber/archive/saber_render.py index a8e3776..3c78b03 100644 --- a/saber/archive/saber_render.py +++ b/modules/saber/archive/saber_render.py @@ -3,7 +3,7 @@ import matplotlib.pyplot as plt import matplotlib.dates as mdates from CONSTANT import DATA_BASEPATH -from saber.utils import * +from modules.saber.utils import * # from matplotlib.colors import LinearSegmentedColormap # 设置字体为支持中文的字体 plt.rcParams['font.family'] = 'SimHei' # 设置为黑体(需要你的环境中有该字体) diff --git a/saber/archive/zjy_wave_fit_plot.py b/modules/saber/archive/zjy_wave_fit_plot.py similarity index 100% rename from saber/archive/zjy_wave_fit_plot.py rename to modules/saber/archive/zjy_wave_fit_plot.py diff --git a/saber/process.py b/modules/saber/process.py similarity index 96% rename from saber/process.py rename to modules/saber/process.py index 7e5f233..4136732 100644 --- a/saber/process.py +++ b/modules/saber/process.py @@ -3,7 +3,7 @@ from datetime import datetime from typing import Dict, List, Optional, Tuple import numpy as np from CONSTANT import DATA_BASEPATH -from saber.utils import * +from modules.saber.utils import * # lat_range=(latitude_min, latitude_max), # alt_range=(altitude_min, altitude_max), diff --git a/saber/render.py b/modules/saber/render.py similarity index 97% rename from saber/render.py rename to modules/saber/render.py index a9eba4a..830a3c8 100644 --- a/saber/render.py +++ b/modules/saber/render.py @@ -6,7 +6,7 @@ from matplotlib.figure import Figure from matplotlib.axes import Axes import matplotlib.dates as mdates -from saber.process import DataProcessor, WaveData, YearlyData +from modules.saber.process import DataProcessor, WaveData, YearlyData @dataclass diff --git a/saber/utils.py b/modules/saber/utils.py similarity index 100% rename from saber/utils.py rename to modules/saber/utils.py diff --git a/tidi/__init__.py b/modules/tidi/__init__.py similarity index 91% rename from tidi/__init__.py rename to modules/tidi/__init__.py index 22a24a9..f2a2fed 100644 --- a/tidi/__init__.py +++ b/modules/tidi/__init__.py @@ -5,8 +5,8 @@ from quart.utils import run_sync from matplotlib import pyplot as plt from CONSTANT import DATA_BASEPATH -from tidi.plot import TidiPlotv2 -from tidi.staged.plot import tidi_render +from modules.tidi.plot import TidiPlotv2 +from modules.tidi.staged.plot import tidi_render tidi_module = Blueprint("Tidi", __name__) diff --git a/tidi/plot.py b/modules/tidi/plot.py similarity index 100% rename from tidi/plot.py rename to modules/tidi/plot.py diff --git a/tidi/staged/plot.py b/modules/tidi/staged/plot.py similarity index 96% rename from tidi/staged/plot.py rename to modules/tidi/staged/plot.py index 07ff5a5..0b494df 100644 --- a/tidi/staged/plot.py +++ b/modules/tidi/staged/plot.py @@ -13,7 +13,7 @@ import matplotlib.pyplot as plt import matplotlib from CONSTANT import DATA_BASEPATH -from tidi.staged.process import tidi_process_idk +from modules.tidi.staged.process import tidi_process_idk # 设置中文显示和负号正常显示 matplotlib.rcParams['font.sans-serif'] = ['SimHei'] # 显示中文 matplotlib.rcParams['axes.unicode_minus'] = False # 正常显示负号 diff --git a/tidi/staged/plot_old.py b/modules/tidi/staged/plot_old.py similarity index 100% rename from tidi/staged/plot_old.py rename to modules/tidi/staged/plot_old.py diff --git a/tidi/staged/process.py b/modules/tidi/staged/process.py similarity index 100% rename from tidi/staged/process.py rename to modules/tidi/staged/process.py diff --git a/staged/cosmic行星波参数全年逐日绘图.py b/staged/cosmic行星波参数全年逐日绘图.py deleted file mode 100644 index 3cd3e6f..0000000 --- a/staged/cosmic行星波参数全年逐日绘图.py +++ /dev/null @@ -1,115 +0,0 @@ -#此代码是对数据处理后的txt数据进行行星波参数提取绘图 - -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 # 正常显示负号 -# 读取一年的数据文件 - -df = pd.read_csv(r'./cosmic.txt', sep='\s+') -# 删除有 NaN 的行 -df = df.dropna(subset=['Temperature']) -# 设置初始参数 -# 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 参数 - -# 设置参数界限 -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 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 = [] - -# 设置最小数据量的阈值 -min_data_points = 36 - -# 进行多个时间窗口的拟合 -#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 - - # 选择当前窗口的数据 - 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 # 跳过当前时间窗口,继续下一个窗口 - - # 提取时间、经度、温度数据 - 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) - - # 将拟合结果添加到列表中 - 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即为拟合的参数汇总 - -# -------------------------------画图---------------------------- -#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)} - -# 获取索引并转换为 numpy 数组 -x_values = fit_df.index.to_numpy() - -# 对每一列生成独立的图 -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('振幅') - - # 设置横坐标的动态调整 - 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) - - plt.show() # 显示图形 \ No newline at end of file diff --git a/staged/cosmic重力波多天.py b/staged/cosmic重力波多天.py deleted file mode 100644 index f829d0a..0000000 --- a/staged/cosmic重力波多天.py +++ /dev/null @@ -1,547 +0,0 @@ -import os -import numpy as np -import pandas as pd -from scipy.interpolate import interp1d -from scipy.optimize import curve_fit -import netCDF4 as nc -import matplotlib.pyplot as plt -import seaborn as sns -# 设置支持中文的字体 -plt.rcParams['font.sans-serif'] = ['Microsoft YaHei'] # 设置字体为微软雅黑 -plt.rcParams['axes.unicode_minus'] = False # 正常显示负号 - -# 定义处理单个文件的函数 -def process_single_file(base_folder_path, 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): - dfs = [] - # 遍历文件夹中的文件 - 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'][:] - # 创建DataFrame - df = pd.DataFrame({ - 'Longitude': lon, - 'Latitude': lat, - 'Altitude': altitude, - 'Temperature': temp - }) - dataset.close() - # 剔除高度大于60的行 - df = df[df['Altitude'] <= 60] - # 对每个文件的数据进行插值 - alt_interp = np.linspace(df['Altitude'].min(), df['Altitude'].max(), 3000) - f_alt = interp1d(df['Altitude'], df['Altitude'], kind='linear', fill_value="extrapolate") - f_lon = interp1d(df['Altitude'], df['Longitude'], kind='linear', fill_value="extrapolate") - f_lat = interp1d(df['Altitude'], df['Latitude'], kind='linear', fill_value="extrapolate") - f_temp = interp1d(df['Altitude'], df['Temperature'], kind='linear', fill_value="extrapolate") - # 计算插值结果 - interpolated_alt = f_alt(alt_interp) - interpolated_lon = f_lon(alt_interp) - interpolated_lat = f_lat(alt_interp) - interpolated_temp = f_temp(alt_interp) - # 创建插值后的DataFrame - interpolated_df = pd.DataFrame({ - 'Altitude': interpolated_alt, - 'Longitude': interpolated_lon, - 'Latitude': interpolated_lat, - 'Temperature': interpolated_temp - }) - # 将插值后的DataFrame添加到列表中 - dfs.append(interpolated_df) - except Exception as e: - print(f"处理文件 {finfo} 时出错: {e}") - # 按行拼接所有插值后的DataFrame - final_df = pd.concat(dfs, axis=0, ignore_index=True) - # 获取 DataFrame 的长度 - num_rows = len(final_df) - # 生成一个每3000个数从0到2999的序列并重复 - altitude_values = np.tile(np.arange(3000), num_rows // 3000 + 1)[:num_rows] - # 将生成的值赋给 DataFrame 的 'Altitude' 列 - final_df['Altitude'] = altitude_values - # 摄氏度换算开尔文 - final_df['Temperature'] = final_df['Temperature'] + 273.15 - - # 筛选出纬度在30到40度之间的数据 - latitude_filtered_df = final_df[(final_df['Latitude'] >= 30) & (final_df['Latitude'] <= 40)] - - # 划分经度网格,20°的网格 - lon_min, lon_max = latitude_filtered_df['Longitude'].min(), latitude_filtered_df['Longitude'].max() - lon_bins = np.arange(lon_min, lon_max + 20, 20) # 创建经度网格边界 - # 将数据分配到网格中 - latitude_filtered_df['Longitude_Grid'] = np.digitize(latitude_filtered_df['Longitude'], lon_bins) - 1 - - # 对相同高度的温度取均值,忽略NaN - altitude_temperature_mean = latitude_filtered_df.groupby('Altitude')['Temperature'].mean().reset_index() - # 重命名列,使其更具可读性 - altitude_temperature_mean.columns = ['Altitude', 'Mean_Temperature'] - - # 定义高度的范围(这里从0到最短段) - altitude_range = range(0, 3000) - all_heights_mean_temperature = [] # 用于存储所有高度下的温度均值结果 - for altitude in altitude_range: - # 筛选出当前高度的所有数据 - altitude_df = latitude_filtered_df[latitude_filtered_df['Altitude'] == altitude] - # 对Longitude_Grid同一区间的温度取均值 - temperature_mean_by_grid = altitude_df.groupby('Longitude_Grid')['Temperature'].mean().reset_index() - # 重命名列,使其更具可读性 - temperature_mean_by_grid.columns = ['Longitude_Grid', 'Mean_Temperature'] - # 添加高度信息列,方便后续区分不同高度的结果 - temperature_mean_by_grid['Altitude'] = altitude - # 将当前高度的结果添加到列表中 - all_heights_mean_temperature.append(temperature_mean_by_grid) - # 将所有高度的结果合并为一个DataFrame - combined_mean_temperature_df = pd.concat(all_heights_mean_temperature, ignore_index=True) - - # 基于Altitude列合并两个DataFrame,只保留能匹配上的行 - merged_df = pd.merge(combined_mean_temperature_df, altitude_temperature_mean, on='Altitude', how='inner') - # 计算差值(减去wn0的扰动) - merged_df['Temperature_Difference'] = merged_df['Mean_Temperature_x'] - merged_df['Mean_Temperature_y'] - - # 按Altitude分组 - grouped = merged_df.groupby('Altitude') - - def single_harmonic(x, A, phi): - return A * np.sin(2 * np.pi / (18 / k) * x + phi) - # 初始化存储每个高度的最佳拟合参数、拟合曲线、残差值以及背景温度的字典 - fit_results = {} - fitted_curves = {} - residuals = {} - background_temperatures = {} - for altitude, group in grouped: - y_data = group['Temperature_Difference'].values - x_data = np.arange(len(y_data)) - wn0_data = group['Mean_Temperature_y'].values # 获取同一高度下的wn0数据 - # 检查Temperature_Difference列是否全部为NaN - if np.all(np.isnan(y_data)): - fit_results[altitude] = {'A': [np.nan] * 5, 'phi': [np.nan] * 5} - fitted_curves[altitude] = [np.nan * x_data] * 5 - residuals[altitude] = np.nan * x_data - background_temperatures[altitude] = np.nan * x_data - else: - # 替换NaN值为非NaN值的均值 - y_data = np.where(np.isnan(y_data), np.nanmean(y_data), y_data) - # 初始化存储WN参数和曲线的列表 - wn_params = [] - wn_curves = [] - # 计算wn0(使用Mean_Temperature_y列数据) - wn0 = wn0_data - - # 对WN1至WN5进行拟合 - for k in range(1, 6): - # 更新单谐波函数中的k值 - harmonic_func = lambda x, A, phi: single_harmonic(x, A, phi) - # 使用curve_fit进行拟合 - popt, pcov = curve_fit(harmonic_func, x_data, y_data, p0=[np.nanmax(y_data) - np.nanmin(y_data), 0]) - A_fit, phi_fit = popt - # 存储拟合结果 - wn_params.append({'A': A_fit, 'phi': phi_fit}) - # 使用拟合参数生成拟合曲线 - WN = harmonic_func(x_data, A_fit, phi_fit) - wn_curves.append(WN) - # 计算残差值 - y_data = y_data - WN # 使用残差值作为下一次拟合的y_data - # 存储结果 - fit_results[altitude] = wn_params - fitted_curves[altitude] = wn_curves - residuals[altitude] = y_data - # 计算同一高度下的背景温度(wn0 + wn1 + wn2 + wn3 + wn4 + wn5) - wn_sum = np.sum([wn0] + wn_curves, axis=0) - background_temperatures[altitude] = wn_sum - - # 将每个字典转换成一个 DataFrame - df = pd.DataFrame(residuals) - # 使用前向填充(用上一个有效值填充 NaN) - df.ffill(axis=1, inplace=True) - # 初始化一个新的字典来保存处理结果 - result = {} - # 定义滤波范围 - lambda_low = 2 # 2 km - lambda_high = 15 # 15 km - f_low = 2 * np.pi / lambda_high - f_high = 2 * np.pi / lambda_low - # 循环处理df的每一行(每个高度) - for idx, residuals_array in df.iterrows(): - # 提取有效值 - valid_values = np.ma.masked_array(residuals_array, np.isnan(residuals_array)) - compressed_values = valid_values.compressed() # 去除NaN值后的数组 - N = len(compressed_values) # 有效值的数量 - # 如果有效值为空(即所有值都是NaN),则将结果设置为NaN - if N == 0: - result[idx] = np.full_like(residuals_array, np.nan) - else: - # 时间序列和频率 - dt = 0.02 # 假设的时间间隔 - n = np.arange(N) - f = n / (N * dt) - # 傅里叶变换 - y = np.fft.fft(compressed_values) - # 频率滤波 - yy = y.copy() - freq_filter = (f >= f_low) & (f <= f_high) - yy[~freq_filter] = 0 - # 逆傅里叶变换 - perturbation_after = np.real(np.fft.ifft(yy)) - # 将处理结果插回到result字典中 - result[idx] = perturbation_after - - # 处理背景温度和扰动温度数据格式 - heights = list(background_temperatures.keys()) - data_length = len(next(iter(background_temperatures.values()))) - background_matrix = np.zeros((data_length, len(heights))) - for idx, height in enumerate(heights): - background_matrix[:, idx] = background_temperatures[height] - - heights = list(result.keys()) - data_length = len(next(iter(result.values()))) - perturbation_matrix = np.zeros((data_length, len(heights))) - for idx, height in enumerate(heights): - perturbation_matrix[:, idx] = result[height] - perturbation_matrix = perturbation_matrix.T - - # 计算 Brunt-Väisälä 频率和势能 - heights_for_calc = np.linspace(0, 60, 3000) * 1000 - - def brunt_vaisala_frequency(g, BT_z, c_p, heights): - # 计算位温随高度的变化率 - dBT_z_dz = np.gradient(BT_z, heights) - # 计算 Brunt-Väisälä 频率,根号内取绝对值 - frequency_squared = (g / BT_z) * ((g / c_p) + dBT_z_dz) - frequency = np.sqrt(np.abs(frequency_squared)) - return frequency - - # 计算势能 - def calculate_gravitational_potential_energy(g, BT_z, N_z, PT_z): - # 计算势能 - return 0.5 * ((g / N_z) ** 2) * ((PT_z / BT_z) ** 2) - g = 9.81 # 重力加速度 - c_p = 1004.5 # 比热容 - N_z_matrix = [] - PT_z_matrix = [] - for i in range(background_matrix.shape[0]): - BT_z = np.array(background_matrix[i]) - PT_z = np.array(perturbation_matrix[i]) - N_z = brunt_vaisala_frequency(g, BT_z, c_p, heights_for_calc) - PW = calculate_gravitational_potential_energy(g, BT_z, N_z, PT_z) - N_z_matrix.append(N_z) - PT_z_matrix.append(PW) - ktemp_Nz = np.vstack(N_z_matrix) - ktemp_Ptz = np.vstack(PT_z_matrix) - mean_ktemp_Nz = np.mean(ktemp_Nz, axis=0) - mean_ktemp_Ptz = np.mean(ktemp_Ptz, axis=0) - - return mean_ktemp_Nz, mean_ktemp_Ptz - else: - print(f"文件夹 {folder_path} 不存在。") - return None, None - -# 主循环,处理1到365个文件 -base_folder_path = r"E:\COSMIC\2008" -all_mean_ktemp_Nz = [] -all_mean_ktemp_Ptz = [] -for file_index in range(101, 104): - try: - mean_ktemp_Nz, mean_ktemp_Ptz = process_single_file(base_folder_path, file_index) - 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) - except ValueError as e: - print(f"Error processing file index {file_index}: {e}, skipping this file.") - continue - -# 转换每个数组为二维形状 -final_mean_ktemp_Nz = np.vstack([arr.reshape(1, -1) for arr in all_mean_ktemp_Nz]) -final_mean_ktemp_Ptz = np.vstack([arr.reshape(1, -1) for arr in all_mean_ktemp_Ptz]) -# 使用条件索引替换大于50的值为NaN -final_mean_ktemp_Ptz[final_mean_ktemp_Ptz > 50] = np.nan -# heights 为每个高度的值 -heights = np.linspace(0, 60, 3000) -df_final_mean_ktemp_Ptz = pd.DataFrame(final_mean_ktemp_Ptz) -df_final_mean_ktemp_Nz = pd.DataFrame(final_mean_ktemp_Nz) -#-------------------------------------------------绘制年统计图------------------------------------ -#-----------绘制浮力频率年统计图----------------------- -data=df_final_mean_ktemp_Nz.T -# 对每个元素进行平方(计算N2) -data= data ** 2 -data=data*10000#(绘图好看) -#将大于 10 的值替换为 NaN(个别异常值) -data[data > 10] = np.nan -# 绘制热力图的函数 -def plot_heatmap(data, heights, title): - plt.figure(figsize=(10, 8)) - # 绘制热力图,数据中的行代表高度,列代表天数 - sns.heatmap(data, cmap='coolwarm', xticklabels=1, yticklabels=50, cbar_kws={'label': 'Value'}) - plt.xlabel('Day') - plt.ylabel('Height (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) - plt.xticks(x_tick_positions, x_tick_labels) - # 设置y轴的刻度,使其显示对应的高度 - plt.yticks(np.linspace(0, data.shape[0] - 1, 6), np.round(np.linspace(heights[0], heights[-1], 6), 2)) - # 反转 y 轴,使 0 在底部 - plt.gca().invert_yaxis() - plt.show() -# 调用函数绘制热力图 -plot_heatmap(data, heights, 'Heatmap of final_mean_ktemp_Nz(10^(-4))') -#----------------------------------------------------------------------------- -#-------------绘制重力势能年统计图------------------------------------------------ -data1=df_final_mean_ktemp_Ptz.T -# 绘制热力图的函数 -def plot_heatmap(data, heights, title): - plt.figure(figsize=(10, 8)) - # 绘制热力图,数据中的行代表高度,列代表天数 - sns.heatmap(data, cmap='coolwarm', xticklabels=1, yticklabels=50, cbar_kws={'label': 'Value'}) - plt.xlabel('Day') - plt.ylabel('Height (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) - plt.xticks(x_tick_positions, x_tick_labels) - # 设置y轴的刻度,使其显示对应的高度 - plt.yticks(np.linspace(0, data.shape[0] - 1, 6), np.round(np.linspace(heights[0], heights[-1], 6), 2)) - # 反转 y 轴,使 0 在底部 - plt.gca().invert_yaxis() - plt.show() -# 调用函数绘制热力图 -plot_heatmap(data1, heights, 'Heatmap of final_mean_ktemp_Ptz(J/kg)') -#------------------------绘制月统计图--------------------------------------------------------------------------------- -#----------绘制浮力频率月统计图------------------------------------------------- -# 获取总列数 -num_columns = data.shape[1] -# 按30列分组计算均值 -averaged_df = [] -# 逐步处理每30列 -for i in range(0, num_columns, 30): - # 获取当前范围内的列,并计算均值 - subset = data.iloc[:, i:i+30] # 获取第i到i+29列 - mean_values = subset.mean(axis=1) # 对每行计算均值 - averaged_df.append(mean_values) # 将均值添加到列表 -# 将结果转化为一个新的 DataFrame -averaged_df = pd.DataFrame(averaged_df).T -# 1. 每3000行取一个均值 -# 获取总行数 -num_rows = averaged_df.shape[0] -# 创建一个新的列表来存储每3000行的均值 -averaged_by_rows_df = [] -# 逐步处理每3000行 -for i in range(0, num_rows, 3000): - # 获取当前范围内的行 - subset = averaged_df.iloc[i:i+3000, :] # 获取第i到i+99行 - mean_values = subset.mean(axis=0) # 对每列计算均值 - averaged_by_rows_df.append(mean_values) # 将均值添加到列表 -# 将结果转化为一个新的 DataFrame -averaged_by_rows_df = pd.DataFrame(averaged_by_rows_df) -# 绘制折线图 -plt.figure(figsize=(10, 6)) # 设置图形的大小 -plt.plot(averaged_by_rows_df.columns, averaged_by_rows_df.mean(axis=0), marker='o', color='b', label='平均值') -# 添加标题和标签 -plt.title('每月平均N^2的折线图') -plt.xlabel('月份') -plt.ylabel('N^2(10^-4)') -plt.legend() -# 显示图形 -plt.grid(True) -plt.xticks(rotation=45) # 让x轴标签(月份)倾斜,以便更清晰显示 -plt.tight_layout() -plt.show() -#------------重力势能的月统计----------------------------------- -# 获取总列数 -num_columns = data1.shape[1] -# 按30列分组计算均值 -averaged_df = [] -# 逐步处理每30列 -for i in range(0, num_columns, 30): - # 获取当前范围内的列,并计算均值 - subset = data1.iloc[:, i:i+30] # 获取第i到i+29列 - mean_values = subset.mean(axis=1) # 对每行计算均值 - averaged_df.append(mean_values) # 将均值添加到列表 -# 将结果转化为一个新的 DataFrame -averaged_df = pd.DataFrame(averaged_df).T -# 1. 每3000行取一个均值 -# 获取总行数 -num_rows = averaged_df.shape[0] -# 创建一个新的列表来存储每3000行的均值 -averaged_by_rows_df = [] -# 逐步处理每3000行 -for i in range(0, num_rows, 3000): - # 获取当前范围内的行 - subset = averaged_df.iloc[i:i+3000, :] # 获取第i到i+99行 - mean_values = subset.mean(axis=0) # 对每列计算均值 - averaged_by_rows_df.append(mean_values) # 将均值添加到列表 -# 将结果转化为一个新的 DataFrame -averaged_by_rows_df = pd.DataFrame(averaged_by_rows_df) -# 绘制折线图 -plt.figure(figsize=(10, 6)) # 设置图形的大小 -plt.plot(averaged_by_rows_df.columns, averaged_by_rows_df.mean(axis=0), marker='o', color='b', label='平均值') -# 添加标题和标签 -plt.title('每月平均重力势能的折线图') -plt.xlabel('月份') -plt.ylabel('重力势能(J/Kg)') -plt.legend() -# 显示图形 -plt.grid(True) -plt.xticks(rotation=45) # 让x轴标签(月份)倾斜,以便更清晰显示 -plt.tight_layout() -plt.show() - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -# 获取总列数 -total_columns = data.shape[1] -# 用于存储每一组30列计算得到的均值列数据(最终会构成新的DataFrame) -mean_columns = [] -# 分组序号,用于生成列名时区分不同的均值列,从1开始 -group_index = 1 -# 按照每30列一组进行划分(不滑动) -for start_col in range(0, total_columns, 30): - end_col = start_col + 30 - if end_col > total_columns: - end_col = total_columns - # 选取当前组的30列(如果不足30列,按实际剩余列数选取) - group_data = data.iloc[:, start_col:end_col] - # 按行对当前组的列数据求和 - sum_per_row = group_data.sum(axis=1) - # 计算平均(每一组的平均,每行都有一个平均结果) - mean_per_row = sum_per_row / (end_col - start_col) - # 生成新的列名,格式为'Mean_分组序号',例如'Mean_1'、'Mean_2'等 - new_column_name = f'Mean_{group_index}' - group_index += 1 - # 将当前组计算得到的均值列添加到列表中 - mean_columns.append(mean_per_row) -# 将所有的均值列合并为一个新的DataFrame(列方向合并) -new_mean_df = pd.concat(mean_columns, axis=1) -# 按行对new_mean_df所有列的数据进行求和,得到一个Series,索引与new_mean_df的索引一致,每个元素是每行的总和 -row_sums = new_mean_df.sum(axis=1) -# 计算所有行总和的均值 -mean_value = row_sums.mean() -# 设置中文字体为黑体,解决中文显示问题(Windows系统下),如果是其他系统或者有其他字体需求可适当调整 -plt.rcParams['font.sans-serif'] = ['SimHei'] -# 解决负号显示问题 -plt.rcParams['axes.unicode_minus'] = False -# 提取月份作为x轴标签(假设mean_value的索引就是月份信息) -months = mean_value.index.tolist() -# 提取均值数据作为y轴数据 -energy_values = mean_value.tolist() -# 创建图形和坐标轴对象 -fig, ax = plt.subplots(figsize=(10, 6)) -# 绘制折线图 -ax.plot(months, energy_values, marker='o', linestyle='-', color='b') -# 设置坐标轴标签和标题 -ax.set_xlabel('月份') -ax.set_ylabel('平均浮力频率') -ax.set_title('每月浮力频率变化趋势') -# 设置x轴刻度,让其旋转一定角度以便更好地显示所有月份标签,避免重叠 -plt.xticks(rotation=45) -# 显示网格线,增强图表可读性 -ax.grid(True) -# 显示图形 -plt.show() -#--------------------------------绘制重力势能月统计图------------------------------ -# 获取总列数 -total_columns = data1.shape[1] -# 用于存储每一组30列计算得到的均值列数据(最终会构成新的DataFrame) -mean_columns = [] -# 分组序号,用于生成列名时区分不同的均值列,从1开始 -group_index = 1 -# 按照每30列一组进行划分(不滑动) -for start_col in range(0, total_columns, 30): - end_col = start_col + 30 - if end_col > total_columns: - end_col = total_columns - # 选取当前组的30列(如果不足30列,按实际剩余列数选取) - group_data = data1.iloc[:, start_col:end_col] - # 按行对当前组的列数据求和 - sum_per_row = group_data.sum(axis=1) - # 计算平均(每一组的平均,每行都有一个平均结果) - mean_per_row = sum_per_row / (end_col - start_col) - # 生成新的列名,格式为'Mean_分组序号',例如'Mean_1'、'Mean_2'等 - new_column_name = f'Mean_{group_index}' - group_index += 1 - # 将当前组计算得到的均值列添加到列表中 - mean_columns.append(mean_per_row) -# 将所有的均值列合并为一个新的DataFrame(列方向合并) -new_mean_df = pd.concat(mean_columns, axis=1) -# 按行对new_mean_df所有列的数据进行求和,得到一个Series,索引与new_mean_df的索引一致,每个元素是每行的总和 -row_sums = new_mean_df.sum(axis=1) -# 计算所有行总和的均值 -mean_value = row_sums.mean() -# 设置中文字体为黑体,解决中文显示问题(Windows系统下),如果是其他系统或者有其他字体需求可适当调整 -plt.rcParams['font.sans-serif'] = ['SimHei'] -# 解决负号显示问题 -plt.rcParams['axes.unicode_minus'] = False -# 提取月份作为x轴标签(假设mean_value的索引就是月份信息) -months = mean_value.index.tolist() -# 提取均值数据作为y轴数据 -energy_values = mean_value.tolist() -# 创建图形和坐标轴对象 -fig, ax = plt.subplots(figsize=(10, 6)) -# 绘制折线图 -ax.plot(months, energy_values, marker='o', linestyle='-', color='b') -# 设置坐标轴标签和标题 -ax.set_xlabel('月份') -ax.set_ylabel('平均浮力频率') -ax.set_title('每月浮力频率变化趋势') -# 设置x轴刻度,让其旋转一定角度以便更好地显示所有月份标签,避免重叠 -plt.xticks(rotation=45) -# 显示网格线,增强图表可读性 -ax.grid(True) -# 显示图形 -plt.show() \ No newline at end of file