fix: cosmic filter

This commit is contained in:
Dustella 2025-04-23 15:40:09 +08:00
parent 5408317ee5
commit a6e598a26f
Signed by: Dustella
GPG Key ID: 35AA0AA3DC402D5C
8 changed files with 71 additions and 13 deletions

View File

@ -34,12 +34,22 @@ def read_data(path):
columns_to_extract = ["alt", "press", "temp", "rh", "u", "v", "wspeed"] columns_to_extract = ["alt", "press", "temp", "rh", "u", "v", "wspeed"]
extracted_df = df[columns_to_extract].copy() extracted_df = df[columns_to_extract].copy()
# 过渡所有 NaN 的行
# 进行单位转换 # 进行单位转换
extracted_df["alt"] = extracted_df["alt"] / 1000 # km extracted_df["alt"] = extracted_df["alt"] / 1000 # km
extracted_df["rh"] = extracted_df["rh"] * 100 # % extracted_df["rh"] = extracted_df["rh"] * 100 # %
# 移除重复的高度值 # 移除重复的高度值
extracted_df = extracted_df.drop_duplicates(subset=["alt"]) extracted_df = extracted_df.drop_duplicates(subset=["alt"])
# Replace infinity values and then drop NaN values
extracted_df = extracted_df.replace([np.inf, -np.inf], np.nan)
extracted_df = extracted_df.replace(math.nan, np.nan)
extracted_df = extracted_df.dropna(how="any")
# if data is empty, return None
if extracted_df.empty:
return None
new_height = np.arange(extracted_df["alt"].min( new_height = np.arange(extracted_df["alt"].min(
), extracted_df["alt"].max() + 0.05, 0.05) ), extracted_df["alt"].max() + 0.05, 0.05)

View File

@ -64,8 +64,9 @@ comboDate = [
def get_dataframe_between_year(all_year_data, start_year, end_year, station): def get_dataframe_between_year(all_year_data, start_year, end_year, station):
res = all_year_data res = all_year_data
filtered_res = res[(res['file_name'].str.extract(rf'{station}-(\d{{4}})')[0].astype(int) >= start_year) & # print(all_year_data["file_name"])
(res['file_name'].str.extract(rf'{station}-(\d{{4}})')[0].astype(int) <= end_year)] filtered_res = res[(res['file_name'].str.extract(rf'-(\d{{4}})\/')[0].astype(int) >= start_year) &
(res['file_name'].str.extract(rf'-(\d{{4}})\/')[0].astype(int) <= end_year)]
return filtered_res return filtered_res
@ -145,6 +146,8 @@ def get_ballon_full_df_by_year(start_year, end_year, station=None, ignore_cache=
# Read data # Read data
data = read_data(file) data = read_data(file)
if data is None:
continue
read_time = time.time() read_time = time.time()
logging.debug( logging.debug(
f"Read data in {read_time - file_start_time:.2f} seconds") f"Read data in {read_time - file_start_time:.2f} seconds")

View File

@ -39,6 +39,16 @@ def process_single_file(base_folder_path, i, lat_range=(30, 40)):
# print(f"正在处理文件: {finfo}") # print(f"正在处理文件: {finfo}")
try: try:
dataset = nc.Dataset(finfo, 'r') dataset = nc.Dataset(finfo, 'r')
try:
bad = int(dataset.getncattr('bad'))
except AttributeError:
bad = None
if bad == 1:
print(f"文件 {finfo} 被标记为坏文件,跳过。")
continue
# 提取变量数据 # 提取变量数据
temp = dataset.variables['Temp'][:] temp = dataset.variables['Temp'][:]
altitude = dataset.variables['MSL_alt'][:] altitude = dataset.variables['MSL_alt'][:]
@ -54,6 +64,10 @@ def process_single_file(base_folder_path, i, lat_range=(30, 40)):
dataset.close() dataset.close()
# 剔除高度大于60的行 # 剔除高度大于60的行
df = df[df['Altitude'] <= 60] df = df[df['Altitude'] <= 60]
df = df[(df["Temperature"] <= 49.85)
& (df["Temperature" >= -120.15])]
# 对每个文件的数据进行插值 # 对每个文件的数据进行插值
alt_interp = np.linspace( alt_interp = np.linspace(
df['Altitude'].min(), df['Altitude'].max(), 3000) df['Altitude'].min(), df['Altitude'].max(), 3000)
@ -77,6 +91,7 @@ def process_single_file(base_folder_path, i, lat_range=(30, 40)):
'Latitude': interpolated_lat, 'Latitude': interpolated_lat,
'Temperature': interpolated_temp 'Temperature': interpolated_temp
}) })
# 将插值后的DataFrame添加到列表中 # 将插值后的DataFrame添加到列表中
dfs.append(interpolated_df) dfs.append(interpolated_df)
except Exception as e: except Exception as e:

View File

@ -38,6 +38,16 @@ def read_and_preprocess_data(base_folder_path, day_num):
print(f"正在处理文件: {finfo}") print(f"正在处理文件: {finfo}")
try: try:
dataset = nc.Dataset(finfo, 'r') dataset = nc.Dataset(finfo, 'r')
try:
bad = int(dataset.getncattr('bad'))
except AttributeError:
bad = None
if bad == 1:
print(f"文件 {finfo} 被标记为坏文件,跳过。")
continue
# 提取变量数据 # 提取变量数据
temp = dataset.variables['Temp'][:] temp = dataset.variables['Temp'][:]
altitude = dataset.variables['MSL_alt'][:] altitude = dataset.variables['MSL_alt'][:]
@ -53,6 +63,11 @@ def read_and_preprocess_data(base_folder_path, day_num):
dataset.close() dataset.close()
# 剔除高度大于60的行 # 剔除高度大于60的行
df = df[df['Altitude'] <= 60] df = df[df['Altitude'] <= 60]
# filter Temperature > df = df[(df["Temperature"] <= 49.85)
# & (df["Temperature" >= -120.15])]
df = df[(df["Temperature"] <= 49.85) &
(df["Temperature"] >= -120.15)]
# 对每个文件的数据进行插值 # 对每个文件的数据进行插值
alt_interp = np.linspace( alt_interp = np.linspace(

View File

@ -44,6 +44,16 @@ def _process_per_folder(base_folder_path, i, target_h):
finfo = os.path.join(folder_path, file_name) finfo = os.path.join(folder_path, file_name)
try: try:
dataset = nc.Dataset(finfo, 'r') dataset = nc.Dataset(finfo, 'r')
try:
bad = int(dataset.getncattr('bad'))
except AttributeError:
bad = None
if bad == 1:
print(f"文件 {finfo} 被标记为坏文件,跳过。")
continue
# 提取变量数据 # 提取变量数据
temp = dataset.variables['Temp'][:] temp = dataset.variables['Temp'][:]
altitude = dataset.variables['MSL_alt'][:] altitude = dataset.variables['MSL_alt'][:]
@ -75,9 +85,12 @@ def _process_per_folder(base_folder_path, i, target_h):
}) })
dataset.close() dataset.close()
# 仅筛选高度 # 仅筛选高度
df_filtered = df[(df["Temperature"] <= 49.85)
& (df["Temperature" >= -120.15])]
df_filtered = df[(df['Altitude'] >= target_h - 0.5) df_filtered = df[(df['Altitude'] >= target_h - 0.5)
& (df['Altitude'] <= target_h + 0.5)] & (df['Altitude'] <= target_h + 0.5)]
result_folder_level.append(df_filtered) result_folder_level.append(df_filtered)
except Exception as e: except Exception as e:
print(f"处理文件 {finfo} 时出错: {e}") print(f"处理文件 {finfo} 时出错: {e}")

View File

@ -342,7 +342,7 @@ def plot_sine_wave_for_day(params_df, date, wind_type, H, model_name):
# 添加图例和标签 # 添加图例和标签
plt.title( plt.title(
f'{date.strftime("%Y-%m-%d")} {wind_type}风速拟合潮汐波', fontsize=16) f'{date.strftime("%Y-%m-%d")} {wind_type}风速, 高度 {H/1000} km 拟合潮汐波', fontsize=16)
plt.xlabel('时间 (小时)', fontsize=14) plt.xlabel('时间 (小时)', fontsize=14)
plt.ylabel('幅度 (m/s)', fontsize=14) plt.ylabel('幅度 (m/s)', fontsize=14)
plt.legend() plt.legend()
@ -386,7 +386,7 @@ def plot_sine_wave_for_day(params_df, date, wind_type, H, model_name):
# 添加图例和标签 # 添加图例和标签
plt.title( plt.title(
f'{date.strftime("%Y-%m-%d")} {wind_type}风速拟合 {model_name}', fontsize=16) f'{date.strftime("%Y-%m-%d")} {wind_type}风速, 高度 {H/1000} km 拟合 {model_name}', fontsize=16)
plt.xlabel('时间 (小时)', fontsize=14) plt.xlabel('时间 (小时)', fontsize=14)
plt.ylabel('幅度 (m/s)', fontsize=14) plt.ylabel('幅度 (m/s)', fontsize=14)
plt.legend() plt.legend()

View File

@ -142,7 +142,7 @@ class SaberGravitywProcessor:
def __init__(self, lat_range: Tuple[float, float], def __init__(self, lat_range: Tuple[float, float],
alt_range: Tuple[float, float], alt_range: Tuple[float, float],
lambda_range: Tuple[float, float], lambda_range: Tuple[float, float],
lvboin: float): lvboin: float = None):
self.lat_min, self.lat_max = lat_range self.lat_min, self.lat_max = lat_range
self.alt_min, self.alt_max = alt_range self.alt_min, self.alt_max = alt_range
self.lambda_low, self.lambda_high = lambda_range self.lambda_low, self.lambda_high = lambda_range

View File

@ -1,4 +1,5 @@
from dataclasses import dataclass from dataclasses import dataclass
from multiprocessing import process
from typing import Optional, Tuple, List from typing import Optional, Tuple, List
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import numpy as np import numpy as np
@ -6,6 +7,7 @@ from matplotlib.figure import Figure
from matplotlib.axes import Axes from matplotlib.axes import Axes
import matplotlib.dates as mdates import matplotlib.dates as mdates
from modules.cosmic import render
from modules.saber.gravityw_process import SaberGravitywProcessor, WaveData, YearlyData from modules.saber.gravityw_process import SaberGravitywProcessor, WaveData, YearlyData
@ -205,8 +207,8 @@ class SaberGravitywRenderer:
axs[0].set_title('(a) NZ') axs[0].set_title('(a) NZ')
axs[0].set_xlabel('Time') axs[0].set_xlabel('Time')
axs[0].set_ylabel('Height') axs[0].set_ylabel('Height')
axs[0].set_yticks(np.linspace(30, 100, 8)) axs[0].set_yticks(np.linspace(20, 100, 9))
axs[0].set_yticklabels(np.round(np.linspace(30, 100, 8), 1)) axs[0].set_yticklabels(np.round(np.linspace(20, 100, 9), 1))
axs[0].set_xticks(x) axs[0].set_xticks(x)
axs[0].set_xticklabels(x) axs[0].set_xticklabels(x)
@ -217,8 +219,8 @@ class SaberGravitywRenderer:
axs[1].set_title('(b) PTZ') axs[1].set_title('(b) PTZ')
axs[1].set_xlabel('Time') axs[1].set_xlabel('Time')
axs[1].set_ylabel('Height') axs[1].set_ylabel('Height')
axs[1].set_yticks(np.linspace(30, 100, 8)) axs[1].set_yticks(np.linspace(20, 100, 9))
axs[1].set_yticklabels(np.round(np.linspace(30, 100, 8), 1)) axs[1].set_yticklabels(np.round(np.linspace(20, 100, 9), 1))
axs[1].set_xticks(x) axs[1].set_xticks(x)
axs[1].set_xticklabels(x) axs[1].set_xticklabels(x)
@ -273,8 +275,8 @@ class SaberGravitywRenderer:
axs[0].set_title('(a) NZ') axs[0].set_title('(a) NZ')
axs[0].set_xlabel('Time') axs[0].set_xlabel('Time')
axs[0].set_ylabel('Height') axs[0].set_ylabel('Height')
axs[0].set_yticks(np.linspace(30, 100, 8)) axs[0].set_yticks(np.linspace(20, 100, 9))
axs[0].set_yticklabels(np.round(np.linspace(30, 100, 8), 1)) axs[0].set_yticklabels(np.round(np.linspace(20, 100, 9), 1))
axs[0].set_xticks(x_ticks) # 设置新的横坐标刻度 axs[0].set_xticks(x_ticks) # 设置新的横坐标刻度
axs[0].set_xticklabels(x_labels, rotation=45) axs[0].set_xticklabels(x_labels, rotation=45)
@ -285,8 +287,8 @@ class SaberGravitywRenderer:
axs[1].set_title('(b) PTZ') axs[1].set_title('(b) PTZ')
axs[1].set_xlabel('Time') axs[1].set_xlabel('Time')
axs[1].set_ylabel('Height') axs[1].set_ylabel('Height')
axs[1].set_yticks(np.linspace(30, 100, 8)) axs[1].set_yticks(np.linspace(20, 100, 9))
axs[1].set_yticklabels(np.round(np.linspace(30, 100, 8), 1)) axs[1].set_yticklabels(np.round(np.linspace(20, 100, 9), 1))
axs[1].set_xticks(x_ticks) # 设置新的横坐标刻度 axs[1].set_xticks(x_ticks) # 设置新的横坐标刻度
axs[1].set_xticklabels(x_labels, rotation=45) axs[1].set_xticklabels(x_labels, rotation=45)