230 lines
10 KiB
Python
230 lines
10 KiB
Python
import os
|
||
import numpy as np
|
||
from scipy.io import loadmat
|
||
import pandas as pd
|
||
|
||
from CONSTANT import DATA_BASEPATH
|
||
|
||
# 第156行代码，纬度可以自己选择，但是都是每5°做一个选择
|
||
# 70km-120km每2.5km一个，共21个，第六个就是82.5km高度的情况，高度也可以自己选择
|
||
|
||
HEIGHTS_CONSTANT = [
|
||
70, 72.5, 75, 77.5, 80, 82.5, 85, 87.5, 90, 92.5, 95, 97.5, 100, 102.5, 105, 107.5, 110, 112.5, 115, 117.5, 120
|
||
]
|
||
|
||
|
||
def tidi_process_idk(
|
||
year,
|
||
input_height: float = 82.5,
|
||
lat_range: tuple = (0, 5),
|
||
):
|
||
# 初始化空列表来存储每天的数据
|
||
if lat_range[1] - lat_range[0] != 5:
|
||
raise ValueError("lat_range must be 5°")
|
||
if input_height not in HEIGHTS_CONSTANT:
|
||
raise ValueError(f"input_height must be in {HEIGHTS_CONSTANT}")
|
||
fixed_height_index = HEIGHTS_CONSTANT.index(input_height)
|
||
height_list = []
|
||
lat_list = []
|
||
lon_list = []
|
||
vmeridional_list = []
|
||
vzonal_list = []
|
||
|
||
# 文件路径
|
||
base_path = f"{DATA_BASEPATH.tidi}/{year}/"
|
||
# if dir not exist, raise error
|
||
if not os.path.exists(base_path):
|
||
raise FileNotFoundError(f"{base_path} not found")
|
||
|
||
# 循环读取从第1天到第365天的数据
|
||
for day in range(1, 366): # 365天数据，确保循环次数为365
|
||
# 构建文件名（使用三位数格式）
|
||
day_str = f"{day:03d}" # 格式化数字为三位数，前面补0
|
||
height_file = f"{day_str}_Height.mat"
|
||
lat_file = f"{day_str}_Lat.mat"
|
||
lon_file = f"{day_str}_Lon.mat"
|
||
vmeridional_file = f"{day_str}_VMerdional.mat"
|
||
vzonal_file = f"{day_str}_Vzonal.mat"
|
||
|
||
# 检查文件是否存在
|
||
try:
|
||
if not os.path.exists(os.path.join(base_path, height_file)):
|
||
raise FileNotFoundError(f"{height_file} not found")
|
||
if not os.path.exists(os.path.join(base_path, lat_file)):
|
||
raise FileNotFoundError(f"{lat_file} not found")
|
||
if not os.path.exists(os.path.join(base_path, lon_file)):
|
||
raise FileNotFoundError(f"{lon_file} not found")
|
||
if not os.path.exists(os.path.join(base_path, vmeridional_file)):
|
||
raise FileNotFoundError(f"{vmeridional_file} not found")
|
||
if not os.path.exists(os.path.join(base_path, vzonal_file)):
|
||
raise FileNotFoundError(f"{vzonal_file} not found")
|
||
|
||
# 读取.mat文件
|
||
height_data = loadmat(os.path.join(base_path, height_file))
|
||
lat_data = loadmat(os.path.join(base_path, lat_file))
|
||
lon_data = loadmat(os.path.join(base_path, lon_file))
|
||
vmeridional_data = loadmat(
|
||
os.path.join(base_path, vmeridional_file))
|
||
vzonal_data = loadmat(os.path.join(base_path, vzonal_file))
|
||
|
||
# 将数据转换为DataFrame
|
||
height_df = pd.DataFrame(height_data['Height'])
|
||
lat_df = pd.DataFrame(lat_data['Lat'])
|
||
lon_df = pd.DataFrame(lon_data['Lon'])
|
||
vmeridional_df = pd.DataFrame(vmeridional_data['VMerdional'])
|
||
vzonal_df = pd.DataFrame(vzonal_data['Vzonal'])
|
||
|
||
# 将每天的数据添加到列表中
|
||
height_list.append(height_df)
|
||
lat_list.append(lat_df)
|
||
lon_list.append(lon_df)
|
||
vmeridional_list.append(vmeridional_df)
|
||
vzonal_list.append(vzonal_df)
|
||
except FileNotFoundError as e:
|
||
print(f"Error: {e}")
|
||
continue # 跳过当前文件，继续处理其他文件
|
||
|
||
# 合并所有天的数据
|
||
height_df = pd.concat(height_list, ignore_index=True)
|
||
lat_df = pd.concat(lat_list, ignore_index=True)
|
||
lon_df = pd.concat(lon_list, ignore_index=True)
|
||
vmeridional_df = pd.concat(vmeridional_list, axis=1) # 按列合并
|
||
vzonal_df = pd.concat(vzonal_list, axis=1) # 按列合并
|
||
|
||
# 初始化空列表来存储每天的数据
|
||
UTime_list = []
|
||
# 初始化累加的时间偏移量
|
||
time_offset = 0
|
||
|
||
# 循环读取从第1天到第365天的数据
|
||
for day in range(1, 366): # 365天数据
|
||
# 构建文件名（使用三位数格式）
|
||
day_str = f"{day:03d}" # 格式化数字为三位数，前面补0
|
||
UTime_file = f"{day_str}_UTime.mat"
|
||
|
||
# 检查UTime文件是否存在
|
||
try:
|
||
if not os.path.exists(os.path.join(base_path, UTime_file)):
|
||
raise FileNotFoundError(f"{UTime_file} not found")
|
||
|
||
# 读取.mat文件
|
||
UTime_data = loadmat(os.path.join(base_path, UTime_file))
|
||
|
||
# 将数据转换为DataFrame
|
||
UTime_df = pd.DataFrame(UTime_data['UTime'])
|
||
|
||
# 对UTime进行累加处理
|
||
UTime_df += time_offset
|
||
|
||
# 将每天的数据添加到列表中
|
||
UTime_list.append(UTime_df)
|
||
|
||
# 更新时间偏移量，为下一天增加24小时
|
||
time_offset += 24
|
||
except FileNotFoundError as e:
|
||
print(f"Error: {e}")
|
||
continue # 跳过当前文件，继续处理其他文件
|
||
|
||
# 合并所有天的数据
|
||
UTime_df = pd.concat(UTime_list, ignore_index=True)
|
||
# 将UTime_df的单位从小时转换为天（除以24）
|
||
UTime_df = UTime_df / 24
|
||
|
||
# 获取固定高度下的纬度数据
|
||
fixed_height_lat = lat_df.iloc[:, 0].values
|
||
# 获取固定高度下的经度数据
|
||
fixed_height_lon = lon_df.iloc[:, 0].values
|
||
# 获取固定高度下的经向风数据
|
||
fixed_height_vmeridional = vmeridional_df.iloc[fixed_height_index, :].values
|
||
# 获取固定高度下的纬向风数据
|
||
fixed_height_vzonal = vzonal_df.iloc[fixed_height_index, :].values
|
||
# 获取对应的世界时数据（这里直接使用整个UTime数据，你可以根据实际情况判断是否需要处理下格式等）
|
||
fixed_height_utime = UTime_df.iloc[:, 0].values
|
||
# 将这些数据整合到一个新的DataFrame（可选操作，方便后续查看等）
|
||
combined_data = pd.DataFrame({
|
||
'Latitude': fixed_height_lat,
|
||
'Longitude': fixed_height_lon,
|
||
'V_Meridional': fixed_height_vmeridional,
|
||
'V_Zonal': fixed_height_vzonal,
|
||
'UTime': fixed_height_utime
|
||
})
|
||
|
||
# 简单打印下整合后的数据前几行查看下内容
|
||
print(combined_data)
|
||
|
||
# 确定纬度和经度的范围
|
||
lat_min, lat_max = np.min(fixed_height_lat), np.max(fixed_height_lat)
|
||
lon_min, lon_max = np.min(fixed_height_lon), np.max(fixed_height_lon)
|
||
|
||
# 计算网格的边界
|
||
lat_bins = np.linspace(lat_min, lat_max, 10) # 纬度分为6个网格，需要7个边界
|
||
lon_bins = np.linspace(lon_min, lon_max, 13) # 经度分为20个网格，需要21个边界
|
||
|
||
# 将数据分配到网格中
|
||
grid_data = []
|
||
for i in range(len(fixed_height_lat)):
|
||
lat_idx = np.digitize(fixed_height_lat[i], lat_bins) - 1 # 找到纬度所在的网格索引
|
||
lon_idx = np.digitize(fixed_height_lon[i], lon_bins) - 1 # 找到经度所在的网格索引
|
||
grid_idx = lat_idx * 12 + lon_idx # 计算网格的唯一索引（9x12=108）
|
||
grid_data.append([fixed_height_lat[i], fixed_height_lon[i], fixed_height_vmeridional[i],
|
||
fixed_height_vzonal[i], fixed_height_utime[i], grid_idx])
|
||
# 将网格数据转换为DataFrame
|
||
grid_df = pd.DataFrame(grid_data, columns=[
|
||
'Latitude', 'Longitude', 'V_Meridional', 'V_Zonal', 'UTime', 'Grid_Index'])
|
||
# 打印网格数据的前几行查看
|
||
print(grid_df.head())
|
||
# 筛选纬度在 0 到 5° 范围内的数据
|
||
lat_range_down, lat_range_upper = lat_range
|
||
|
||
filtered_df = grid_df[(grid_df['Latitude'] >= lat_range_down) &
|
||
(grid_df['Latitude'] <= lat_range_upper)]
|
||
# 将经度从度转换为弧度
|
||
filtered_df['Longitude_Radians'] = np.radians(filtered_df['Longitude'])
|
||
# 选择需要输出的列
|
||
output_columns = ['Longitude_Radians', 'UTime', 'V_Zonal']
|
||
# 将数据输出为txt文件
|
||
output_file = base_path+'TIDI_V_Zonal_data.txt' # 经向风速分量
|
||
filtered_df[output_columns].to_csv(output_file, sep='\t', index=False)
|
||
|
||
zonal_data = filtered_df[output_columns]
|
||
# 选择需要输出的列
|
||
output_columns1 = ['Longitude_Radians', 'UTime', 'V_Meridional']
|
||
# 将数据输出为txt文件
|
||
merdi_data = filtered_df[output_columns1]
|
||
|
||
return {
|
||
"V_Zonal": zonal_data,
|
||
"V_Meridional": merdi_data
|
||
}
|
||
# output_file1 = base_path+'TIDI_V_Meridional_data.txt' # 纬向风速分量
|
||
# final_result.to_csv(output_file1, sep='\t', index=False)
|
||
# # 打印结果确认
|
||
# print(f"数据已保存到 {output_file1}")
|
||
|
||
|
||
# # 对网格数据按照Grid_Index进行分组，并计算每个网格的统计数据
|
||
# grid_stats = grid_df.groupby('Grid_Index').agg({
|
||
# 'Latitude': 'mean', # 计算纬度的平均值
|
||
# 'Longitude': 'mean', # 计算经度的平均值
|
||
# 'V_Meridional': 'mean', # 计算经向风的平均值、最大值和最小值
|
||
# 'V_Zonal': 'mean', # 计算纬向风的平均值、最大值和最小值
|
||
# 'UTime': 'mean' # 计算UTime的平均值
|
||
# }).reset_index()
|
||
#
|
||
# # 重命名列，使其更具可读性
|
||
# grid_stats.columns = ['Grid_Index', 'Mean_Latitude', 'Mean_Longitude', 'Mean_VMeridional', 'Mean_VZonal', 'Mean_UTime']
|
||
# # 将经度转换为2π范围
|
||
# grid_stats['Mean_Longitude_2pi'] = (grid_stats['Mean_Longitude'] / 360) * (2 * np.pi)
|
||
# # 定义纬度的最小值和最大值（起码要40°区间，否则数据太少无法拟合）
|
||
# latmin =0
|
||
# latmax=40
|
||
# # 筛选出Mean_Latitude在lat_min到lat_max之间的数据
|
||
# filtered_grid_stats = grid_stats[(grid_stats['Mean_Latitude'] >= latmin) & (grid_stats['Mean_Latitude'] <= latmax)]
|
||
#
|
||
# # 检查筛选后的数据行数是否小于18
|
||
# if filtered_grid_stats.shape[0] < 18:
|
||
# print("无法拟合行星波")
|
||
# else:
|
||
# # 打印筛选后的网格数据的前几行查看
|
||
# print(filtered_grid_stats.head())
|