352 lines
12 KiB
Python
352 lines
12 KiB
Python
import xarray as xr
|
||
import numpy as np
|
||
import sys
|
||
import ctypes
|
||
from datetime import datetime, timedelta
|
||
# from WindDataHandler import DataFileInfo,WindDataHandler
|
||
from PySide6.QtWidgets import (QApplication, QMainWindow, QPushButton,
|
||
QVBoxLayout, QWidget, QProgressDialog)
|
||
from PySide6.QtCore import QTimer, Qt
|
||
from PySide6.QtWidgets import QMessageBox
|
||
import argparse
|
||
|
||
|
||
import ctypes
|
||
|
||
# 定义与 C++ 中 DataFileInfo 对应的结构体
|
||
class DataFileInfo(ctypes.Structure):
|
||
_fields_ = [
|
||
("Height", ctypes.c_int64),
|
||
("Width", ctypes.c_int64),
|
||
("Num", ctypes.c_int64),
|
||
("firstTimestamp", ctypes.c_int64),
|
||
("lastTimestamp", ctypes.c_int64),
|
||
("minLon", ctypes.c_double),
|
||
("maxLon", ctypes.c_double),
|
||
("minLat", ctypes.c_double),
|
||
("maxLat", ctypes.c_double),
|
||
("ESPGCODE", ctypes.c_int64),
|
||
("T11", ctypes.c_double),
|
||
("T12", ctypes.c_double),
|
||
("T13", ctypes.c_double),
|
||
("T21", ctypes.c_double),
|
||
("T22", ctypes.c_double),
|
||
("T23", ctypes.c_double),
|
||
("fileSize", ctypes.c_int64)
|
||
]
|
||
|
||
|
||
class WindDataHandler:
|
||
def __init__(self, dll_path='./LAMPWindData.dll'):
|
||
self.dll_path = dll_path
|
||
self.wind_lib = ctypes.CDLL(dll_path)
|
||
self._setup_prototypes()
|
||
self.wind_lib.getDataFileInfo.argtypes = [ctypes.c_char_p] # 参数是 const char*
|
||
self.wind_lib.getDataFileInfo.restype = DataFileInfo # 返回类型是 DataFileInfo 结构体
|
||
self.current_file_info = None
|
||
|
||
|
||
def _setup_prototypes(self):
|
||
"""设置 DLL 函数的参数和返回类型"""
|
||
# ... 函数原型设置代码见上一节，此处省略 ...
|
||
|
||
def get_data_file_info(self, filepath):
|
||
"""获取风场数据文件信息"""
|
||
# 将Python字符串转换为字节字符串
|
||
filepath_bytes = filepath.encode('utf-8')
|
||
info = self.wind_lib.getDataFileInfo(filepath_bytes)
|
||
self.current_file_info = info
|
||
return info
|
||
|
||
def get_wind_data_file_time_arr(self, filepath, info):
|
||
"""获取时间数组"""
|
||
filepath_bytes = filepath.encode('utf-8')
|
||
if info.Num <= 0:
|
||
return []
|
||
|
||
# 创建足够大的数组来存储时间戳
|
||
time_arr = (ctypes.c_int64 * info.Num)()
|
||
result = self.wind_lib.get_WindDataFileTimeArr(filepath_bytes, info, time_arr)
|
||
|
||
if result == 0: # 假设返回0表示成功
|
||
return list(time_arr)
|
||
else:
|
||
print(f"Error getting time array: {result}")
|
||
return []
|
||
|
||
def read_wind_data_file(self, filepath, info, time_id):
|
||
"""读取指定时间ID的风场数据 (U/V分量)"""
|
||
filepath_bytes = filepath.encode('utf-8')
|
||
layer_size = info.Height * info.Width
|
||
|
||
# 创建数组来存储U和V分量数据
|
||
u_arr = (ctypes.c_double * layer_size)()
|
||
v_arr = (ctypes.c_double * layer_size)()
|
||
|
||
result = self.wind_lib.Read_WindDataFile(filepath_bytes, info, time_id, u_arr, v_arr)
|
||
|
||
if result >= 0:
|
||
# 将ctypes数组转换为Python列表或numpy数组（如果可用）
|
||
return list(u_arr), list(v_arr)
|
||
else:
|
||
print(f"Error reading wind data: {result}")
|
||
return None, None
|
||
|
||
def write_wind_data_file(self, filepath, info, timeInt64,time_id, u_arr, v_arr):
|
||
"""将风场数据写入文件"""
|
||
filepath_bytes = filepath.encode('utf-8')
|
||
layer_size = info.Height * info.Width
|
||
|
||
# 确保输入数据长度正确
|
||
if len(u_arr) != layer_size or len(v_arr) != layer_size:
|
||
raise ValueError(f"UArr and VArr must have length {layer_size}")
|
||
|
||
# 将Python列表转换为ctypes数组
|
||
u_arr_ctypes = (ctypes.c_double * layer_size)(*u_arr)
|
||
v_arr_ctypes = (ctypes.c_double * layer_size)(*v_arr)
|
||
|
||
result = self.wind_lib.Write_WindDataFile(filepath_bytes, info, timeInt64,time_id, u_arr_ctypes, v_arr_ctypes)
|
||
return result
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
class ERA5ToWindDataConverter:
|
||
def __init__(self, dll_path='./LAMPWindData.dll'):
|
||
self.wind_handler = WindDataHandler(dll_path)
|
||
|
||
def read_era5_data(self, nc_file_path):
|
||
"""读取ERA5数据文件"""
|
||
try:
|
||
# 使用xarray读取NetCDF文件[6,8](@ref)
|
||
ds = xr.open_dataset(nc_file_path)
|
||
print("ERA5数据集变量:", list(ds.variables.keys()))
|
||
|
||
# 提取u10和v10分量（10米高度风场）[1,5](@ref)
|
||
u10 = ds['u10'].values # 东西风分量
|
||
v10 = ds['v10'].values # 南北风分量
|
||
lon = ds['longitude'].values # 经度
|
||
lat = ds['latitude'].values # 纬度
|
||
time = ds['time'].values # 时间
|
||
|
||
ds.close()
|
||
|
||
return u10, v10, lon, lat, time
|
||
except Exception as e:
|
||
print(f"读取ERA5数据失败: {e}")
|
||
return None, None, None, None, None
|
||
|
||
|
||
def create_data_file_info(self, lon, lat, time):
|
||
"""创建DataFileInfo结构体"""
|
||
info = DataFileInfo()
|
||
|
||
# 设置网格尺寸[1](@ref)
|
||
info.Height = len(lat) # 纬度维度
|
||
info.Width = len(lon) # 经度维度
|
||
info.Num = len(time) # 时间步数
|
||
|
||
# 设置地理范围[1](@ref)
|
||
info.minLon = float(lon.min())
|
||
info.maxLon = float(lon.max())
|
||
info.minLat = float(lat.min())
|
||
info.maxLat = float(lat.max())
|
||
|
||
|
||
time_deltas = time
|
||
info.firstTimestamp = int(time_deltas[0])
|
||
info.lastTimestamp = int(time_deltas[-1])
|
||
|
||
# 设置坐标系（WGS84的EPSG代码为4326）
|
||
info.ESPGCODE = 4326
|
||
|
||
# 变换参数（单位矩阵表示无变换）
|
||
info.T11, info.T12, info.T13 = lon[0],lon[1]-lon[0],0.0 # 计算 lon X_geo = GT(0) + X_col * GT(1) + Y_row * GT(2)
|
||
info.T21, info.T22, info.T23 = lat[0], 0.0, lat[1]-lat[0] # 计算 lat Y_geo = GT(3) + X_col * GT(4) + Y_row * GT(5)
|
||
|
||
return info
|
||
|
||
def convert_era5_time(self, time_values):
|
||
"""转换ERA5时间格式为datetime对象[1,5](@ref)"""
|
||
# ERA5时间通常是从1900-01-01开始的小时数
|
||
return [time_val.astype("int64") for time_val in time_values]
|
||
# base_time = datetime(1900, 1, 1, 0, 0, 0)
|
||
#
|
||
# converted_times = []
|
||
# for time_val in time_values:
|
||
# if isinstance(time_val, np.datetime64):
|
||
# # 如果是numpy datetime64格式
|
||
# dt = time_val.astype(datetime)
|
||
# else:
|
||
# # 如果是数值格式（从1900-01-01开始的小时数）
|
||
# hours = float(time_val)
|
||
# dt = base_time + timedelta(hours=hours)
|
||
# converted_times.append(dt)
|
||
#
|
||
# return converted_times
|
||
|
||
def convert_era5_to_winddata(self, era5_file, output_file):
|
||
"""将ERA5数据转换为WindData格式"""
|
||
|
||
progress_dialog = QProgressDialog(
|
||
"文件写入进度，请稍候...", # 标签文本
|
||
"取消", # 取消按钮文本
|
||
0, 100, # 最小值和最大值
|
||
None # 父窗口
|
||
)
|
||
# 1. 读取ERA5数据
|
||
u10, v10, lon, lat, time = self.read_era5_data(era5_file)
|
||
if u10 is None:
|
||
return False
|
||
|
||
# 2. 转换时间格式
|
||
time_datetime = self.convert_era5_time(time)
|
||
|
||
|
||
# 4. 创建文件信息头
|
||
file_info = self.create_data_file_info(lon, lat, time_datetime)
|
||
|
||
# 5. 写入WindData文件
|
||
try:
|
||
|
||
progress_dialog.setWindowTitle("文件写入进度，请稍候...")
|
||
progress_dialog.setFixedSize(500, 50)
|
||
progress_dialog.setCancelButton(None)
|
||
progress_dialog.setWindowModality(Qt.NonModal)
|
||
progress_dialog.setMinimumDuration(0)
|
||
progress_dialog.setRange(0, len(time_datetime))
|
||
progress_dialog.show()
|
||
# 逐个时间步写入数据
|
||
for time_idx in range(len(time_datetime)):
|
||
if progress_dialog is None or not progress_dialog.isVisible() or progress_dialog.wasCanceled():
|
||
print("窗口已关闭")
|
||
continue
|
||
# 获取当前时间步的u、v分量
|
||
u_slice = u10[time_idx, :, :].flatten()
|
||
v_slice = v10[time_idx, :, :].flatten()
|
||
|
||
timeInt=time_datetime[time_idx]
|
||
print("时间：",timeInt)
|
||
# 写入数据
|
||
result = self.wind_handler.write_wind_data_file(
|
||
output_file, file_info,ctypes.c_int64(timeInt), time_idx, list(u_slice), list(v_slice)
|
||
)
|
||
progress_dialog.setValue(time_idx)
|
||
QApplication.processEvents()
|
||
if result != 0:
|
||
print(f"写入时间步 {time_idx} 失败: {result}")
|
||
return False
|
||
|
||
print(f"转换成功！输出文件: {output_file}")
|
||
if progress_dialog is None or not progress_dialog.isVisible() or progress_dialog.wasCanceled():
|
||
return True
|
||
else:
|
||
progress_dialog.close()
|
||
return True
|
||
|
||
except Exception as e:
|
||
print(f"转换过程出错: {e}")
|
||
QMessageBox.critical(None, "错误", f"文件保存错误: {e}\n{output_file}")
|
||
return False
|
||
|
||
def compareResult(self,era5_file, bin_file):
|
||
u10, v10, lon, lat, time = self.read_era5_data(era5_file)
|
||
time_datetime = self.convert_era5_time(time)
|
||
info = converter.wind_handler.get_data_file_info(bin_file)
|
||
print(f"文件信息: {info.Height}x{info.Width} 网格, {info.Num} 个时间步")
|
||
# 获取时间数组
|
||
time_arr = converter.wind_handler.get_wind_data_file_time_arr(bin_file, info)
|
||
print(f"时间范围: {time_arr[0]} 到 {time_arr[-1]}")
|
||
|
||
print("比对时间")
|
||
for tid in range(info.Num):
|
||
if time_arr[tid]-time_datetime[tid]==0:
|
||
continue
|
||
else:
|
||
print("时间错误，",time_arr[tid],time_datetime[tid])
|
||
|
||
print("比对实际数据值")
|
||
|
||
for tid in range(info.Num):
|
||
u_slice = u10[tid, :, :].flatten()
|
||
v_slice = v10[tid, :, :].flatten()
|
||
|
||
u_wind,v_wind=self.wind_handler.read_wind_data_file(bin_file,info,tid)
|
||
|
||
u_wind=np.array(u_wind)
|
||
v_wind=np.array(v_wind)
|
||
|
||
u_d=u_slice-u_wind
|
||
v_d=v_slice-v_wind
|
||
|
||
u_d_min=np.max(np.abs(u_d))
|
||
v_d_min=np.max(np.abs(v_d))
|
||
|
||
if u_d_min>0 or v_d_min>0:
|
||
print("tid:(u_d,v_d)",tid,u_d_min,v_d_min)
|
||
|
||
|
||
pass
|
||
|
||
|
||
def getParams():
|
||
parser = argparse.ArgumentParser()
|
||
parser.add_argument('-i', '--infile', type=str,default=r"D:\Programs\SpacetyLabelAIPlante\Manual-Labeling-Tool\WindFileData\data\total_precipitation202001.nc",help='输入ERA5 的netcdf文件')
|
||
parser.add_argument('-o', '--outfile', type=str, default=r"D:\Programs\SpacetyLabelAIPlante\Manual-Labeling-Tool\WindFileData\output\converted_wind_data.lampwindbin",help='输出自定义风场文件')
|
||
parser.add_argument('-d', '--dllfile', type=str, default=r"D:\Programs\SpacetyLabelAIPlante\Manual-Labeling-Tool\WindFileData\dllFolder\WindDataOperator.dll",help='输出自定义风场文件')
|
||
|
||
args = parser.parse_args()
|
||
|
||
return args
|
||
|
||
|
||
# 使用示例
|
||
if __name__ == "__main__":
|
||
|
||
app = QApplication(sys.argv)
|
||
|
||
parser = getParams()
|
||
dllpath=parser.dllfile
|
||
ncpath=parser.infile
|
||
binpath=parser.outfile
|
||
# 创建转换器实例
|
||
converter = ERA5ToWindDataConverter(dllpath)
|
||
|
||
# 执行转换
|
||
era5_file = ncpath # ERA5输入文件
|
||
output_file = binpath # 输出文件
|
||
|
||
success = converter.convert_era5_to_winddata(era5_file, output_file)
|
||
print('写入结果',success)
|
||
if success:
|
||
QMessageBox.information(None, "操作成功", "文件保存成功！\n{}".format(output_file))
|
||
# 验证转换结果
|
||
info = converter.wind_handler.get_data_file_info(output_file)
|
||
print(f"文件信息: {info.Height}x{info.Width} 网格, {info.Num} 个时间步")
|
||
else:
|
||
# 错误框
|
||
pass
|
||
|
||
sys.exit(app.exec())
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|