Delft 3D根据netCDF创建空间风场文件

使用数据为ECMWF的ERA5数据，仅选择10m u component of wind和10m u component of windERA5 hourly time-series data on single levels from 1940 to presenthttps://cds.climate.copernicus.eu/datasets/reanalysis-era5-single-levels-timeseries?tab=overview

 下文是官方示例，详见Delft3D-FLOW_User_Manual文档

FileVersion = 1.03 
filetype = meteo_on_equidistant_grid 
NODATA_value = -999.000 
n_cols = 3 
n_rows = 4 
grid_unit = degree 
x_llcenter = -12.000 
y_llcenter = 48.000 
dx = 0.12500 
dy = 0.083333333 
n_quantity = 1 
quantity1 = x_wind 
unit1 = m s-1 
TIME = 0.0 hours since 2008-01-15 04:35:00 +00:00 
2 3.0 3.6 
3 4.5 2 
2.2 1 2.3 
1.2 0.7 -0.4 
TIME = 6.0 hours since 2008-01-15 04:35:00 +00:00 
-1.1 -2.3 -3.6 
-3.2 0.8 1.1 
2.2 -1 -1.6 
1.2 -0.7 -0.4

以下是python代码 ，nc文件和python文件放一起即可

import sys
import numpy as np
from netCDF4 import Dataset, num2date
from datetime import datetime, timedelta
import os

def total_seconds_cal(start, end):
    """
    计算两个时间点之间的秒数差，支持字符串或 datetime 格式输入
    """
    if isinstance(start, str):
        start = datetime.strptime(start, "%Y-%m-%d %H:%M:%S")
    if isinstance(end, str):
        end = datetime.strptime(end, "%Y-%m-%d %H:%M:%S")
    return int((start - end).total_seconds())

# 输入的文件名（自己的文件）
nc_fileID = "2016_1-10.nc"

# 检查文件是否存在
if not os.path.exists(nc_fileID):
    sys.exit(f'\nERROR! File "{nc_fileID}" not found')

# 读取 NetCDF 文件并打印可用的变量
ncfile = Dataset(nc_fileID)
print('\nAvailable variables in selected .NC file:\n')

vars = ncfile.variables
for n, var_name in enumerate(vars):
    print(f'# {n}    {var_name} {vars[var_name].shape}')

# 检查是否存在 'valid_time' 变量
if 'valid_time' not in vars:
    sys.exit("\nERROR! 'valid_time' variable not found in the dataset")

# 设置读取的变量
variables_to_read = ['u10', 'v10']

# 检查变量是否存在
for var in variables_to_read:
    if var not in vars:
        sys.exit(f'\nERROR! Variable "{var}" not found in the dataset')

# 读取变量（u10 和 v10）
nc_data = {var: vars[var][:].astype(np.float32) for var in variables_to_read}

print('\nData read complete\n')

# 获取时间参考
time_ref = num2date(0, units=ncfile.variables['valid_time'].units)
print(f"Time reference: {time_ref}")

# 计算文件valid_time的起止时间
start_time_inFile = time_ref + timedelta(seconds=int(ncfile.variables['valid_time'][0]))

# 起止时间（根据需要更改）
extracted_start_time = "2016-01-01 00:00:00"
extracted_end_time = "2016-06-30 23:00:00"

# 计算总秒数（相对于 time_ref）
extracted_start_seconds = total_seconds_cal(extracted_start_time, time_ref)
extracted_end_seconds = total_seconds_cal(extracted_end_time, time_ref)

# 获取时间索引
valid_times = vars['valid_time'][:].astype(int)
try:
    start_index = np.where(valid_times == extracted_start_seconds)[0][0]
    end_index = np.where(valid_times == extracted_end_seconds)[0][0]
except IndexError:
    sys.exit("\nERROR! Extracted time not found in dataset")

# 计算 dx 和 dy (经纬度差)
dx = np.abs(vars['longitude'][1] - vars['longitude'][0])
dy = np.abs(vars['latitude'][1] - vars['latitude'][0])

# 处理变量数据
for var_name in variables_to_read:
    variable_mapping = {'u10': ('.amu', 'x_wind'), 'v10': ('.amv', 'y_wind')}
    file_extension, quantity = variable_mapping[var_name]
    outfile_name = f'{var_name}{file_extension}'
    
    with open(outfile_name, 'w') as outfile:
        # 写入header
        outfile.write('FileVersion = 1.03\n')
        outfile.write('filetype = meteo_on_equidistant_grid\n')
        outfile.write('NODATA_value = -999.0\n')
        
        n_cols, n_rows = vars[var_name].shape[2], vars[var_name].shape[1]
        outfile.write(f'n_cols = {n_cols}\n')
        outfile.write(f'n_rows = {n_rows}\n')
        
        outfile.write('grid_unit = degree\n')
        outfile.write(f'x_llcenter = {ncfile.variables["longitude"][0]}\n')
        outfile.write(f'y_llcenter = {ncfile.variables["latitude"][-1]}\n')
        
        outfile.write(f'dx = {dx}\n')
        outfile.write(f'dy = {dy}\n')
        
        outfile.write(f'n_quantity = 1\n')
        outfile.write(f'quantity1 = {quantity}\n')
        outfile.write('unit1 = m s-1\n')
        
        # 写入每个时间步长的数据
        for t in range(start_index, end_index):
            time_value = np.ma.filled(vars['valid_time'][t], 0)
            if isinstance(time_value, np.ndarray):
                time_value = time_value.item()
            current_time = time_ref + timedelta(seconds=time_value)
            hours = total_seconds_cal(current_time, start_time_inFile) // 3600
            outfile.write(f'TIME = {hours}.0 hours since {start_time_inFile} +00:00\n')
            np.savetxt(outfile, nc_data[var_name][t, :, :], fmt='%.6f')
            
print('Procession finished.')