Converting NetCDF file from Delft3D to Netcdf for Opendrift

[GuiCruz]

Hello,

I am trying to convert my .netcdf file from Delft3D, but the file format and variable names are not compatible with the netcdf format that OpenDrift reads. To address this, I am working on converting this file to a new netcdf file that OpenDrift can read. The file is loading correctly, but for some reason, the point’s trajectory is stationary, as if there were no velocity field.

Could anyone help me figure out what I might be doing wrong? I've attached the file with the original data. I am only using the depth-integrated mean current velocity.

`import netCDF4 as nc
import numpy as np
import utm
from opendrift.models.oceandrift import OceanDrift
from opendrift.readers import reader_netCDF_CF_unstructured
from opendrift.readers import reader_netCDF_CF_generic
from scipy.interpolate import griddata
import matplotlib.pyplot as plt
from pyproj import Proj, transform
from metpy.interpolate import (interpolate_to_grid, remove_nan_observations,
remove_repeat_coordinates)

fn = 'D:/Dropbox/Doutorado/LC08_RS/Model_artigo2/trim-Modelo_Dir207_Mag5.nc'
ds = nc.Dataset(fn)

x = ds.variables['XCOR'][9:91, 2:72]
y = ds.variables['YCOR'][9:91, 2:72]

U = np.mean(ds.variables['U1'][-1, :, 9:91, 2:72], axis=0)
V = np.mean(ds.variables['V1'][-1, :, 9:91, 2:72], axis=0)

plt.figure(figsize=(6, 6))
plt.quiver(x, y, U, V, scale_units='xy', angles='xy')
plt.title('Campo de Vetores Interpolado - Tempo 5')
plt.xlabel('Coordenadas UTM (X)')
plt.ylabel('Coordenadas UTM (Y)')

resolucao_metros = 100 # Defina a resolução desejada em metros

x_flat = np.ravel(x)
y_flat = np.ravel(y)
U_flat = np.ravel(U)
V_flat = np.ravel(V)

xmin, xmax = np.min(x_flat), np.max(x_flat)
ymin, ymax = np.min(y_flat), np.max(y_flat)

num_pontos_x = int((xmax - xmin) / resolucao_metros) + 1
num_pontos_y = int((ymax - ymin) / resolucao_metros) + 1

xi = np.linspace(xmin, xmax, resolucao_metros)
yi = np.linspace(ymin, ymax, resolucao_metros)
xi, yi = np.meshgrid(xi, yi)

U_grid = griddata((y_flat, x_flat), U_flat, (yi, xi), method='linear')
V_grid = griddata((y_flat, x_flat), V_flat, (yi, xi), method='linear')

gx_1 = xi[0,:]
gy_1 = yi[:,0]

M = np.sqrt(U_grid2 + V_grid2)

plt.figure(figsize=(6, 6))
plt.quiver(gx_1, gy_1 , U_grid, V_grid, angles='xy')
plt.title('Campo de Vetores Interpolado - Tempo 5')
plt.xlabel('Coordenadas UTM (X)')
plt.ylabel('Coordenadas UTM (Y)')
#plt.axis('equal')

t = np.arange(0, 18000, 600)
t

r_u = np.expand_dims(U_grid, axis=0)
r_v = np.expand_dims(V_grid, axis=0)

r_u = np.expand_dims(r_u, axis=0)
r_v = np.expand_dims(r_v, axis=0)

r_u = np.concatenate([r_u, r_u, r_u, r_u, r_u, r_u, r_u, r_u, r_u, r_u, r_u, r_u, r_u, r_u, r_u, r_u, r_u, r_u, r_u, r_u, r_u, r_u, r_u, r_u, r_u, r_u, r_u, r_u, r_u, r_u])
r_v = np.concatenate([r_v, r_v, r_v, r_v, r_v, r_v, r_v, r_v, r_v, r_v, r_v, r_v, r_v, r_v, r_v, r_v, r_v, r_v, r_v, r_v, r_v, r_v, r_v, r_v, r_v, r_v, r_v, r_v, r_v, r_v])

np.shape(r_u)

#DS_Open.close()
DS_Open = nc.Dataset('./teste.nc', 'w', format='NETCDF4_CLASSIC', diskless=True, persist=True)
DS_Open.Conventions = 'CF-1.7'
DS_Open.title = 'teste'
DS_Open.institution = 'Unidata'
DS_Open.source = 'WRF-1.5'
DS_Open.history = ''
DS_Open.references = ''
DS_Open.comment = ''
DS_Open.createDimension('time', None)
DS_Open.createDimension('depth', 1)
DS_Open.createDimension('X', gx_1.size)
DS_Open.createDimension('Y', gy_1.size)

UTM WGS84 Zona 22S)

crs = DS_Open.createVariable('crs', np.int32, ())
crs.grid_mapping_name = 'transverse_mercator'
crs.false_easting = 500000.0
crs.false_northing = 10000000.0
crs.longitude_of_central_meridian = -51.0
crs.latitude_of_projection_origin = 0.0
crs.scale_factor_at_central_meridian = 0.9996
crs.semi_major_axis = 6378137.0
crs.inverse_flattening = 298.257223563
crs.proj4 = "+proj=utm +zone=22 +south +datum=WGS84 +units=m +no_defs"
crs.spatial_ref = "EPSG:32722" # WKT para EPSG:32722 (UTM 22S)

D_var = DS_Open.createVariable('depth', np.float32, ('depth',))
D_var[:] = 0
D_var.units = 'm'
D_var.positive = 'down'
D_var.axis = 'Z' # Optional
D_var.standard_name = 'depth'

time_var = DS_Open.createVariable('time', np.int32, ('time',))
time_var[:] = t
time_var.long_name = 'time since initialization'
time_var.units = 'seconds since 2020-07-01 00:00:00'
time_var.calendar = 'proleptic_gregorian'
time_var.standard_name = 'time' # Optional

U_var = DS_Open.createVariable('u', datatype=np.float32,
dimensions=('time', 'depth', 'X', 'Y'),
zlib=True)
U_var[:] = r_u
U_var.units = 'meter second-1'
U_var.standard_name = 'y_sea_water_velocity'
U_var.long_name = 'Sea water y velocity'

V_var = DS_Open.createVariable('v', datatype=np.float32,
dimensions=('time', 'depth', 'X', 'Y'),
zlib=True)
V_var[:] = r_v
V_var.units = 'meter second-1'
V_var.standard_name = 'x_sea_water_velocity'
V_var.long_name = 'Sea water x velocity'

y_var = DS_Open.createVariable('y', np.float32, ('X',))
y_var[:] = gy_1
y_var.units = 'm'
y_var.standard_name = 'projection_y_coordinate'
y_var.long_name = 'y-coordinate in projected coordinate system'

x_var = DS_Open.createVariable('x', np.float32, ('Y',))
x_var[:] = gx_1
x_var.units = 'm'
x_var.standard_name = 'projection_x_coordinate'
x_var.long_name = 'x-coordinate in projected coordinate system'

DS_Open.variables

`

trim-Modelo_Dir207_Mag5.zip

[GuiCruz]

I managed to convert the Delft3D NetCDF file for OpenDrift, but I had to make some modifications to the code. The issue seemed to be with the TIME variable I created and the time_step in o.run(time_step=600). I didn't fully understand why this happened, but here is the code that worked:

`
import datetime
import utm
from pyproj import CRS
from opendrift.readers import reader_netCDF_CF_generic
from opendrift.readers import reader_netCDF_CF_unstructured
from opendrift.readers import reader_shape
from opendrift.readers import reader_global_landmask
from opendrift.models.oceandrift import OceanDrift
from scipy.interpolate import griddata
import numpy as np
import netCDF4 as nc
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D #
from matplotlib import cm
from metpy.interpolate import (interpolate_to_grid, remove_nan_observations,
remove_repeat_coordinates)

fn = 'C:/Users/guioc/Documents/Art3_OD/trim-Modelo_Dir207_Mag5.nc'
ds = nc.Dataset(fn)
x = ds.variables['XCOR'][9:92,1:51]
y = ds.variables['YCOR'][9:92,1:51]

k = 1
U = np.mean(ds.variables['U1'][-11:-1,:,9:92,1:51],1)
V = np.mean(ds.variables['V1'][-11:-1,:,9:92,1:51],1)

dim0 = len(U)
print(dim0)
dim2 = len(U[0])
print(dim2)
dim3 = len(U[0][0])
print(dim3)

x_f = np.ravel(x)
y_f = np.ravel(y)

r_u = []
r_v = []
r_m = []
i = 0
while i < dim0:

u_fp = np.ravel(U[0,:,:])
v_fp = np.ravel(V[0,:,:])

for t in range (len(u_fp)):
    if u_fp[t]==-999.999:
        u_fp[t]=0

for t in range (len(v_fp)):
    if v_fp[t]==-999.999:
        v_fp[t]=0
    
gx, gy, u_tp = interpolate_to_grid(x_f, y_f, u_fp, interp_type='linear', hres=50)
gx, gy, v_tp = interpolate_to_grid(x_f, y_f, v_fp, interp_type='linear', hres=50)
v_tp = v_tp*-1
r_u.append(u_tp)
r_v.append(v_tp)
i+=1
print(i)

[gy_Deg,gx_Deg] = utm.to_latlon(gx, gy, 22, 'H')
[x_utm,y_utm,n,z] = utm.from_latlon(gx_Deg,gy_Deg)

gx_1 = x_utm[0,:]
gy_1 = y_utm[:,0]
r_u = np.expand_dims(r_u, axis=1)
r_v = np.expand_dims(r_v, axis=1)

gx_Deg = np.nan_to_num(gx_Deg)
gy_Deg = np.nan_to_num(gy_Deg)
r_u = np.nan_to_num(r_u)
r_v = np.nan_to_num(r_v)

t = np.arange(0, 6000, 600)
np.shape(t)

fig, ax = plt.subplots(dpi=600)
M = np.sqrt(u_tp2+v_tp2)
q = ax.quiver(gx_Deg, gy_Deg,u_tp,v_tp,M,cmap=plt.cm.jet, scale_units='xy')#zoom: [200:250,0:80]
colorbar = plt.colorbar(q)
colorbar.ax.tick_params(labelsize=3)
plt.axis('off')
plt.margins(0,0)
ax.set_aspect('equal')
plt.scatter(-52.08, -32.192)

DS_Open = nc.Dataset('./OD_25x100_Wind10ms_opendrift.nc', 'w', format='NETCDF4_CLASSIC', diskless=True, persist=True)
DS_Open.Conventions = 'CF-1.7'
DS_Open.title = '25x100_Wind05ms_opendrift'
DS_Open.institution = 'Unidata'
DS_Open.source = 'WRF-1.5'
DS_Open.history = ''
DS_Open.references = ''
DS_Open.comment = ''
DS_Open.createDimension('time', None)
DS_Open.createDimension('depth', 1)
DS_Open.createDimension('X', gx_1.size)
DS_Open.createDimension('Y', gy_1.size)

crs = DS_Open.createVariable('latitude_longitude', np.int32, ())
crs.grid_mapping_name = 'latitude_longitude'

D_var = DS_Open.createVariable('depth', np.float32, ('depth',))
D_var[:] = 0
D_var.units = 'm'
D_var.positive = 'down'
D_var.axis = 'Z' # Optional
D_var.standard_name = 'depth'

time_var = DS_Open.createVariable('time', np.int32, ('time',))
time_var[:] = t
time_var.long_name = 'time since initialization'
time_var.units = 'seconds since 2020-07-01 00:00:00'
time_var.calendar = 'proleptic_gregorian'
time_var.standard_name = 'time' # Optional

U_var = DS_Open.createVariable('u', datatype=np.float32,
dimensions=('time', 'depth', 'Y', 'X'),
zlib=True)
U_var[:] = r_u
U_var.units = 'meter second-1'
U_var.standard_name = 'y_sea_water_velocity'
U_var.long_name = 'Sea water y velocity'

V_var = DS_Open.createVariable('v', datatype=np.float32,
dimensions=('time', 'depth', 'Y', 'X'),
zlib=True)
V_var[:] = r_v
V_var.units = 'meter second-1'
V_var.standard_name = 'x_sea_water_velocity'
V_var.long_name = 'Sea water x velocity'

lon_var = DS_Open.createVariable('lon', np.float32, ('Y', 'X'))
lon_var[:] = gx_Deg
lon_var.units = ''
lon_var.standard_name = 'longitude' # Optional
lon_var.long_name = 'longitude'

lat_var = DS_Open.createVariable('lat', np.float32, ('Y', 'X'))
lat_var[:] = gy_Deg
lat_var.units = ''
lat_var.standard_name = 'latitude' # Optional
lat_var.long_name = 'latitude'

DS_Open.variables

DS_Open.close()
o = OceanDrift()

reader = reader_netCDF_CF_generic.Reader('./OD_25x100_Wind10ms_opendrift.nc')
o.add_reader([reader])

o.seed_elements(lon=-52.08, lat=-32.192, number=100, radius =1000, z= 0, time=reader.start_time)
o.run(time_step = 600, end_time=reader.end_time)

o.plot(markersize=5)
`