内容管理网站广州网站优化工具
该文参考博主@Anywayyyyy.的博文《如何使用python处理nc数据制作Mike风场文件–以ERA5数据为例》
链接: link
代码结构主体参考上文,但由于@Anywayyyyy.博主使用的mikeio库是1.0以下版本,1.0以上版本进行了大量改动导致无法直接使用,而安装老版本的mikeio库也会出现各种问题,因此笔者查找了源码及注释,给出了新版本mikeio库制作风场dfs2文件的方法,代码如下。
"""
把era5下载的风场改写为mike的输入文件dfs2格式
"""
import mikeio
from mikeio import Dfs2
from mikeio.eum import ItemInfo, EUMUnit
from mikeio.eum import EUMType
import numpy as np
import netCDF4 as nc
import datetime as dtfile = nc.Dataset('D:/MIKE Zero Projects/waigaoqiao/era5_wind_and_pressure.nc', 'r') # 文件路径及名称
time = np.array(file.variables['time'][:]) # 读取时间信息
lon = np.array(file.variables['longitude'][:]) # 读取经度信息
lat = np.flipud(np.array(file.variables['latitude'][:])) # 读取纬度信息
u = np.flip(np.array(file.variables['u10'][:]), axis=1) # 读取东西风速u
v = np.flip(np.array(file.variables['v10'][:]), axis=1) # 读取南北风速v
p = np.flip(np.array(file.variables['sp'][:]), axis=1) # 读取表面压强
由于era5风场数据的纬度是从大到小排列的,与平时的习惯不同,为了防止后面出现各种问题,在导入数据后将其以及其他变量进行翻转以符合正常习惯。
u = u * 1.1 # 将风速放大1.1倍
v = v * 1.1time_dt = [] # 指定一个新的列表
tstart = dt.datetime(1900, 1, 1, 0) # nc文件的起始时刻
for i in time:time_dt.append(tstart+dt.timedelta(hours=int(i))) # 将nc文件中的时刻转换为正常的时间格式dfsfilename = 'D:/MIKE Zero Projects/waigaoqiao/2021uvp_for_fm.dfs2' # 给定文件名
x0 = lon[0]
y0 = lat[0]
dx = 0.25 # 设置网格间距,ERA5默认网格间距均为0.25度
dy = 0.25
time0 = time_dt[0] # 设置其实时间,用到了上一步的时间处理
delta_t = 3600 # 设置时间间隔,ERA5数据时间间隔为1h,3600s# 创建变量, EUMType.Wind_Velocity说明变量类型为风速, EUMUnit.meter_per_sec说明变量单位为m/s
items = [ItemInfo("u", EUMType.Wind_Velocity, EUMUnit.meter_per_sec),ItemInfo("v", EUMType.Wind_Velocity, EUMUnit.meter_per_sec),ItemInfo("p", EUMType.Pressure, EUMUnit.pascal)]geometry = mikeio.Grid2D(x0=x0, nx=len(lon), dx=dx, y0=y0, ny=len(lat), dy=dy, projection="LONG/LAT")dfs = Dfs2() # 实例化一个dfs变量d = [np.array(u), np.array(v), np.array(p)] # 将要输出至dfs2的数组打包
aa = mikeio.Dataset(data=d, time=time_dt, items=items, geometry=geometry)# dfs.write(filename=dfsfilename, data=d, start_time=time0, dt=delta_t, items=items, coordinate=coordinate, dx=dx, dy=dy) # 输出变量
dfs.write(data=aa, filename=dfsfilename, dt=delta_t, dx=dx, dy=dy, title='wind_uvp') # 输出变量
写的比较粗糙,有各种问题欢迎大家评论区留言,笔者会根据问题再做修改。