search

OpenDrift / opendrift

Open source framework for ocean trajectory modelling
https://opendrift.github.io
GNU General Public License v2.0
246 stars 120 forks source link

Problems with svelocities x, y in the oceandrift module. They are scrambled #438

Closed GuiCruz closed 3 years ago

GuiCruz commented 3 years ago

Hello. I'm trying to use the output data from the Delft3D modeling software. However, as the grid that delft3d uses is curvilinear, I need to extract the data and interpolate in a regular grid so that Opendrift can read the data correctly.

I did a little routine to extract data from Delft3d and create a netcdf file that oceandrift can read correctly, but genereic_reader shuffles my data towards east and west. Since actually the velocity field is always constant north to the south.

If anyone can help me I would appreciate it very much. I think there are other people trying to couple delft3d with opendrift

PS:if necessary I can attach the data and the grid of the delft3d model

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
fn  = 'D:/Dropbox/Doutorado/LC08_RS/Ideal_model/trim-Model_Ideal_river_25x100km_Wind10ms.nc'
fn_p = 'D:/Dropbox/Doutorado/LC08_RS/Ideal_model/trim-Model_Ideal_river_25x100km_Wind10ms_NA.nc'
fn_shp = 'D:/Dropbox/Doutorado/LC08_RS/Ideal_model/Coastline_25x100km.shp'

ds = nc.Dataset(fn)
x = ds.variables['XCOR'][10:81,1:51]
y = ds.variables['YCOR'][10:81,1:51]
k = 1 
#ds.variables['KMAXOUT_RESTR'][:]
t = ds.variables['time'][:]
#the data from delft3d is time, K-layer,Y,X
U = ds.variables['U1'][:,1,10:81,1:51]
V = ds.variables['V1'][:,1,10:81,1:51]

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

[gx_Deg,gy_Deg] = utm.to_latlon(x, y, 22, 'H')
[x_utm,y_utm,n,z] = utm.from_latlon(gx_Deg,gy_Deg)
x_f = np.ravel(x_utm)
y_f = np.ravel(y_utm)

#lats_gy = []
#lons_gx = []

utm_gy = np.arange(np.min(y_utm),np.max(y_utm), 500)
utm_gx = np.arange(np.min(x_utm), np.max(x_utm), 500)

#lats_gy = np.arange(np.min(gy_Deg),np.max(gy_Deg), 0.001)
#lons_gx = np.arange(np.min(gx_Deg), np.max(gx_Deg), 0.005)

xi, yi = np.meshgrid(utm_gx, utm_gy, sparse=False, indexing='xy')
r_u = []
r_v = []
i = 0
while i < dim0:

    u_fp = np.ravel(U[i,:,:])
    v_fp = np.ravel(V[i,:,:])
    ui = griddata((y_f, x_f), u_fp, (yi, xi),method='cubic',fill_value = 0)
    vi = griddata((y_f, x_f), v_fp, (yi, xi),method='cubic',fill_value = 0)
    r_u.append(ui)
    r_v.append(vi)
    i+=1
    print(i)

gx_1 = xi[0,:]
gy_1 = yi[:,0]    
r_u = np.expand_dims(r_u, axis=1)
r_v = np.expand_dims(r_v, axis=1)
np.shape(r_u)

checking my data

fig = plt.figure()
ax = plt.axes(projection='3d')
ax.plot_trisurf(np.ravel(xi), np.ravel(yi), np.ravel(ui), cmap=cm.jet)

plt.show()
#DS_Open.close()
DS_Open = nc.Dataset('./OD_25x100_Wind05ms_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('UTM', np.float32, ())#
crs.proj4 = '+proj=utm +zone=22 +south +datum=WGS84 +units=m +no_defs '

y_var = DS_Open.createVariable('lat', np.float64, ('Y',))
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('lon', np.float64, ('X',))
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'

D_var = DS_Open.createVariable('depth', np.int32, ('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 = 'gregorian'
time_var.field = 'time, scalar, series'
time_var.axis = 'T'  # Optional
time_var.standard_name = 'time'  # Optional

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

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

DS_Open.variables

{'UTM': <class 'netCDF4._netCDF4.Variable'>
 float32 UTM()
     proj4: +proj=utm +zone=22 +south +datum=WGS84 +units=m +no_defs 
 unlimited dimensions: 
 current shape = ()
 filling on, default _FillValue of 9.969209968386869e+36 used,
 'lat': <class 'netCDF4._netCDF4.Variable'>
 float64 lat(Y)
     units: m
     standard_name: projection_y_coordinate
     long_name: y-coordinate in projected coordinate system
 unlimited dimensions: 
 current shape = (163,)
 filling on, default _FillValue of 9.969209968386869e+36 used,
 'lon': <class 'netCDF4._netCDF4.Variable'>
 float64 lon(X)
     units: m
     standard_name: projection_x_coordinate
     long_name: x-coordinate in projected coordinate system
 unlimited dimensions: 
 current shape = (40,)
 filling on, default _FillValue of 9.969209968386869e+36 used,
 'depth': <class 'netCDF4._netCDF4.Variable'>
 int32 depth(depth)
     units: m
     positive: down
     axis: Z
     standard_name: depth
 unlimited dimensions: 
 current shape = (1,)
 filling on, default _FillValue of -2147483647 used,
 'time': <class 'netCDF4._netCDF4.Variable'>
 int32 time(time)
     long_name: time since initialization
     units: seconds since 2020-07-01 00:00:00
     calendar: gregorian
     field: time, scalar, series
     axis: T
     standard_name: time
 unlimited dimensions: time
 current shape = (97,)
 filling on, default _FillValue of -2147483647 used,
 'u': <class 'netCDF4._netCDF4.Variable'>
 float64 u(time, depth, Y, X)
     units: meter second-1
     standard_name: x_sea_water_velocity
     long_name: Sea water x velocity
 unlimited dimensions: time
 current shape = (97, 1, 163, 40)
 filling on, default _FillValue of 9.969209968386869e+36 used,
 'v': <class 'netCDF4._netCDF4.Variable'>
 float64 v(time, depth, Y, X)
     units: meter second-1
     standard_name: y_sea_water_velocity
     long_name: Sea water y velocity
 unlimited dimensions: time
 current shape = (97, 1, 163, 40)
 filling on, default _FillValue of 9.969209968386869e+36 used}
o = OceanDrift(loglevel=0)
reader = reader_netCDF_CF_generic.Reader('./OD_25x100_Wind05ms_opendrift.nc')
#reader_shp = reader_shape.Reader.from_shpfiles(fn_shp)
#reader_landmask = reader_global_landmask.Reader(extent=[np.min(gx_Deg),np.min(gy_Deg),np.max(gx_Deg),np.max(gy_Deg)])  # lonmin, latmin, lonmax, latmax
o.add_reader([reader])
lats = np.linspace(-32.32541009, -32.32541009, 100)
lons = np.linspace(-52.18633433, -51.97683977, 100)

o.seed_elements(lon=lons, lat=lats,time=datetime.datetime(2020, 7, 3, 0, 0))
#o.set_config('general:use_auto_landmask', False)
#o.set_config('general:coastline_action', 'stranding')

o.run()
o.plot(background=['x_sea_water_velocity', 'y_sea_water_velocity'])

link of a printscreen with the results: [https://paste.pics/AQAP8]

knutfrode commented 3 years ago

Are the original current components oriented along the UTM x- and y-directions, or are they east- and north-components provided on a UTM grid?

It could help if you could post ncdump of your Delft-file, and the output from >>> print(reader)

GuiCruz commented 3 years ago

The components are oriented in Y and X direction, but I specified when interpolating with griddata as YX )
ui = griddata((y_f, x_f), u_fp, (yi, xi),method='cubic',fill_value = 0) vi = griddata((y_f, x_f), v_fp, (yi, xi),method='cubic',fill_value = 0) )

and the new netcdf file the dimensions are also specified as Y, X: ( U_var = DS_Open.createVariable('u', datatype=np.float64, dimensions=('time', 'depth', 'Y', 'X'), )

I checked my data with this command and they are normal just like the original model

fig = plt.figure()
ax = plt.axes(projection='3d')
ax.plot_trisurf(np.ravel(xi), np.ravel(yi), np.ravel(ui), cmap=cm.jet)
plt.show()

Ouput from delft3d netcdf file:

ds.variables

{'XCOR': <class 'netCDF4._netCDF4.Variable'>
 float32 XCOR(MC, NC)
     standard_name: projection_x_coordinate
     long_name: X-coordinate of grid points
     units: m
 unlimited dimensions: 
 current shape = (92, 112)
 filling on, default _FillValue of 9.969209968386869e+36 used,
 'YCOR': <class 'netCDF4._netCDF4.Variable'>
 float32 YCOR(MC, NC)
     standard_name: projection_y_coordinate
     long_name: Y-coordinate of grid points
     units: m
 unlimited dimensions: 
 current shape = (92, 112)
 filling on, default _FillValue of 9.969209968386869e+36 used,
 'XZ': <class 'netCDF4._netCDF4.Variable'>
 float32 XZ(M, N)
     standard_name: projection_x_coordinate
     long_name: X-coordinate of cell centres
     units: m
 unlimited dimensions: 
 current shape = (92, 112)
 filling on, default _FillValue of 9.969209968386869e+36 used,
 'YZ': <class 'netCDF4._netCDF4.Variable'>
 float32 YZ(M, N)
     standard_name: projection_y_coordinate
     long_name: Y-coordinate of cell centres
     units: m
 unlimited dimensions: 
 current shape = (92, 112)
 filling on, default _FillValue of 9.969209968386869e+36 used,
 'ALFAS': <class 'netCDF4._netCDF4.Variable'>
 float32 ALFAS(M, N)
     long_name: Orientation ksi-axis w.r.t. pos.x-axis at water level point
     units: arc_degrees
     coordinates: XZ YZ
     grid: grid
     location: face
 unlimited dimensions: 
 current shape = (92, 112)
 filling on, default _FillValue of 9.969209968386869e+36 used,
 'KCU': <class 'netCDF4._netCDF4.Variable'>
 int32 KCU(MC, N)
     long_name: Mask array for U-velocity points
     units: 1
     grid: grid
     location: edge1
 unlimited dimensions: 
 current shape = (92, 112)
 filling on, default _FillValue of -2147483647 used,
 'KCV': <class 'netCDF4._netCDF4.Variable'>
 int32 KCV(M, NC)
     long_name: Mask array for V-velocity points
     units: 1
     grid: grid
     location: edge2
 unlimited dimensions: 
 current shape = (92, 112)
 filling on, default _FillValue of -2147483647 used,
 'KCS': <class 'netCDF4._netCDF4.Variable'>
 int32 KCS(M, N)
     long_name: Non-active/active water-level point
     units: 1
     coordinates: XZ YZ
     grid: grid
     location: face
 unlimited dimensions: 
 current shape = (92, 112)
 filling on, default _FillValue of -2147483647 used,
 'DP0': <class 'netCDF4._netCDF4.Variable'>
 float32 DP0(MC, NC)
     long_name: Initial bottom depth (positive down)
     units: m
     coordinates: XCOR YCOR
     grid: grid
     location: node
 unlimited dimensions: 
 current shape = (92, 112)
 filling on, default _FillValue of 9.969209968386869e+36 used,
 'DPS0': <class 'netCDF4._netCDF4.Variable'>
 float32 DPS0(M, N)
     long_name: Initial bottom depth at zeta points (positive down)
     units: m
     coordinates: XZ YZ
     grid: grid
     location: face
 unlimited dimensions: 
 current shape = (92, 112)
 filling on, default _FillValue of 9.969209968386869e+36 used,
 'DPU0': <class 'netCDF4._netCDF4.Variable'>
 float32 DPU0(MC, N)
     long_name: Initial bottom depth at u points (positive down)
     units: m
     grid: grid
     location: edge1
 unlimited dimensions: 
 current shape = (92, 112)
 filling on, default _FillValue of 9.969209968386869e+36 used,
 'DPV0': <class 'netCDF4._netCDF4.Variable'>
 float32 DPV0(M, NC)
     long_name: Initial bottom depth at v points (positive down)
     units: m
     grid: grid
     location: edge2
 unlimited dimensions: 
 current shape = (92, 112)
 filling on, default _FillValue of 9.969209968386869e+36 used,
 'NAMCON': <class 'netCDF4._netCDF4.Variable'>
 |S1 NAMCON(LSTSCI, strlen20)
     long_name: Name of constituent quantity
 unlimited dimensions: 
 current shape = (1, 20)
 filling on, default _FillValue of  used,
 'NAMTUR': <class 'netCDF4._netCDF4.Variable'>
 |S1 NAMTUR(LTUR, strlen20)
     long_name: Name of turbulent quantity
 unlimited dimensions: 
 current shape = (2, 20)
 filling on, default _FillValue of  used,
 'SIG_LYR': <class 'netCDF4._netCDF4.Variable'>
 float32 SIG_LYR(SIG_LYR)
     standard_name: ocean_sigma_coordinate
     long_name: Sigma coordinates of layer centres
     units: 1
     formula_terms: sigma: SIG_LYR eta: S1 depth: DPS0
 unlimited dimensions: 
 current shape = (5,)
 filling on, default _FillValue of 9.969209968386869e+36 used,
 'SIG_INTF': <class 'netCDF4._netCDF4.Variable'>
 float32 SIG_INTF(SIG_INTF)
     standard_name: ocean_sigma_coordinate
     long_name: Sigma coordinates of layer interfaces
     units: 1
     formula_terms: sigma: SIG_INTF eta: S1 depth: DPS0
 unlimited dimensions: 
 current shape = (6,)
 filling on, default _FillValue of 9.969209968386869e+36 used,
 'GSQS': <class 'netCDF4._netCDF4.Variable'>
 float32 GSQS(M, N)
     long_name: Horizontal area of computational cell
     units: m2
     coordinates: XZ YZ
     grid: grid
     location: face
 unlimited dimensions: 
 current shape = (92, 112)
 filling on, default _FillValue of 9.969209968386869e+36 used,
 'PPARTITION': <class 'netCDF4._netCDF4.Variable'>
 int32 PPARTITION(M, N)
     long_name: Partition
     units: 1
     coordinates: XZ YZ
     grid: grid
     location: face
 unlimited dimensions: 
 current shape = (92, 112)
 filling on, default _FillValue of -2147483647 used,
 'KMAXOUT': <class 'netCDF4._netCDF4.Variable'>
 int32 KMAXOUT(KMAXOUT)
     long_name: User selected output layer interfaces
     units: 1
     compress: SIG_INTF
 unlimited dimensions: 
 current shape = (6,)
 filling on, default _FillValue of -2147483647 used,
 'KMAXOUT_RESTR': <class 'netCDF4._netCDF4.Variable'>
 int32 KMAXOUT_RESTR(KMAXOUT_RESTR)
     long_name: User selected output layer centres
     units: 1
     compress: SIG_LYR
 unlimited dimensions: 
 current shape = (5,)
 filling on, default _FillValue of -2147483647 used,
 'RHOCONST': <class 'netCDF4._netCDF4.Variable'>
 float32 RHOCONST()
     long_name: User specified constant density
     units: kg/m3
 unlimited dimensions: 
 current shape = ()
 filling on, default _FillValue of 9.969209968386869e+36 used,
 'GRAVITY': <class 'netCDF4._netCDF4.Variable'>
 float32 GRAVITY()
     long_name: Gravitational acceleration
     units: m/s2
 unlimited dimensions: 
 current shape = ()
 filling on, default _FillValue of 9.969209968386869e+36 used,
 'grid': <class 'netCDF4._netCDF4.Variable'>
 int32 grid()
     units: 1
     cf_role: grid_topology
     topology_dimension: 2
     node_dimensions: MC NC
     face_dimensions: M:MC (padding: low) N:NC (padding: low)
     face_coordinates: XZ YZ
     node_coordinates: XCOR YCOR
     vertical_dimensions: SIG_LYR:SIG_INTF (padding: none)
 unlimited dimensions: 
 current shape = ()
 filling on, default _FillValue of -2147483647 used,
 'time': <class 'netCDF4._netCDF4.Variable'>
 float32 time(time)
     standard_name: time
     long_name: time
     units: seconds since 2020-07-01 00:00:00
     calendar: proleptic_gregorian
 unlimited dimensions: time
 current shape = (97,)
 filling on, default _FillValue of 9.969209968386869e+36 used,
 'S1': <class 'netCDF4._netCDF4.Variable'>
 float32 S1(time, M, N)
     long_name: Water-level in zeta point
     units: m
     coordinates: XZ YZ
     grid: grid
     location: face
 unlimited dimensions: time
 current shape = (97, 92, 112)
 filling on, default _FillValue of 9.969209968386869e+36 used,
 'KFU': <class 'netCDF4._netCDF4.Variable'>
 int32 KFU(time, MC, N)
     long_name: Non-active/active in U-point
     units: 1
     grid: grid
     location: edge1
 unlimited dimensions: time
 current shape = (97, 92, 112)
 filling on, default _FillValue of -2147483647 used,
 'KFV': <class 'netCDF4._netCDF4.Variable'>
 int32 KFV(time, M, NC)
     long_name: Non-active/active in V-point
     units: 1
     grid: grid
     location: edge2
 unlimited dimensions: time
 current shape = (97, 92, 112)
 filling on, default _FillValue of -2147483647 used,
 'U1': <class 'netCDF4._netCDF4.Variable'>
 float32 U1(time, KMAXOUT_RESTR, MC, N)
     long_name: U-velocity per layer in U-point (Eulerian)
     units: m/s
     grid: grid
     location: edge1
 unlimited dimensions: time
 current shape = (97, 5, 92, 112)
 filling on, default _FillValue of 9.969209968386869e+36 used,
 'V1': <class 'netCDF4._netCDF4.Variable'>
 float32 V1(time, KMAXOUT_RESTR, M, NC)
     long_name: V-velocity per layer in V-point (Eulerian)
     units: m/s
     grid: grid
     location: edge2
 unlimited dimensions: time
 current shape = (97, 5, 92, 112)
 filling on, default _FillValue of 9.969209968386869e+36 used,
 'W': <class 'netCDF4._netCDF4.Variable'>
 float32 W(time, KMAXOUT, M, N)
     long_name: W-omega per layer in zeta point
     units: m/s
     coordinates: XZ YZ
     grid: grid
     location: face
 unlimited dimensions: time
 current shape = (97, 6, 92, 112)
 filling on, default _FillValue of 9.969209968386869e+36 used,
 'WPHY': <class 'netCDF4._netCDF4.Variable'>
 float32 WPHY(time, KMAXOUT_RESTR, M, N)
     long_name: W-velocity per layer in zeta point
     units: m/s
     coordinates: XZ YZ
     grid: grid
     location: face
 unlimited dimensions: time
 current shape = (97, 5, 92, 112)
 filling on, default _FillValue of 9.969209968386869e+36 used,
 'R1': <class 'netCDF4._netCDF4.Variable'>
 float32 R1(time, LSTSCI, KMAXOUT_RESTR, M, N)
     long_name: Concentrations per layer in zeta point
     units: 1
     coordinates: XZ YZ
     grid: grid
     location: face
 unlimited dimensions: time
 current shape = (97, 1, 5, 92, 112)
 filling on, default _FillValue of 9.969209968386869e+36 used,
 'RTUR1': <class 'netCDF4._netCDF4.Variable'>
 float32 RTUR1(time, LTUR, KMAXOUT, M, N)
     long_name: Turbulent quantity per layer in zeta point
     units: 1
     coordinates: XZ YZ
     grid: grid
     location: face
 unlimited dimensions: time
 current shape = (97, 2, 6, 92, 112)
 filling on, default _FillValue of 9.969209968386869e+36 used,
 'TAUKSI': <class 'netCDF4._netCDF4.Variable'>
 float32 TAUKSI(time, MC, N)
     long_name: Bottom stress in U-point
     units: N/m2
     grid: grid
     location: edge1
 unlimited dimensions: time
 current shape = (97, 92, 112)
 filling on, default _FillValue of 9.969209968386869e+36 used,
 'TAUETA': <class 'netCDF4._netCDF4.Variable'>
 float32 TAUETA(time, M, NC)
     long_name: Bottom stress in V-point
     units: N/m2
     grid: grid
     location: edge2
 unlimited dimensions: time
 current shape = (97, 92, 112)
 filling on, default _FillValue of 9.969209968386869e+36 used,
 'TAUMAX': <class 'netCDF4._netCDF4.Variable'>
 float32 TAUMAX(time, M, N)
     long_name: Tau_max in zeta points (scalar)
     units: N/m2
     coordinates: XZ YZ
     grid: grid
     location: face
 unlimited dimensions: time
 current shape = (97, 92, 112)
 filling on, default _FillValue of 9.969209968386869e+36 used,
 'VICWW': <class 'netCDF4._netCDF4.Variable'>
 float32 VICWW(time, KMAXOUT, M, N)
     long_name: Vertical eddy viscosity-3D in zeta point
     units: m2/s
     coordinates: XZ YZ
     grid: grid
     location: face
 unlimited dimensions: time
 current shape = (97, 6, 92, 112)
 filling on, default _FillValue of 9.969209968386869e+36 used,
 'DICWW': <class 'netCDF4._netCDF4.Variable'>
 float32 DICWW(time, KMAXOUT, M, N)
     long_name: Vertical eddy diffusivity-3D in zeta point
     units: m2/s
     coordinates: XZ YZ
     grid: grid
     location: face
 unlimited dimensions: time
 current shape = (97, 6, 92, 112)
 filling on, default _FillValue of 9.969209968386869e+36 used,
 'RICH': <class 'netCDF4._netCDF4.Variable'>
 float32 RICH(time, KMAXOUT, M, N)
     long_name: Richardson number
     units: 1
     coordinates: XZ YZ
     grid: grid
     location: face
 unlimited dimensions: time
 current shape = (97, 6, 92, 112)
 filling on, default _FillValue of 9.969209968386869e+36 used,
 'RHO': <class 'netCDF4._netCDF4.Variable'>
 float32 RHO(time, KMAXOUT_RESTR, M, N)
     long_name: Density per layer in zeta point
     units: kg/m3
     coordinates: XZ YZ
     grid: grid
     location: face
 unlimited dimensions: time
 current shape = (97, 5, 92, 112)
 filling on, default _FillValue of 9.969209968386869e+36 used,
 'UMNLDF': <class 'netCDF4._netCDF4.Variable'>
 float32 UMNLDF(time, MC, N)
     long_name: Filtered U-velocity
     units: m/s
     grid: grid
     location: edge1
 unlimited dimensions: time
 current shape = (97, 92, 112)
 filling on, default _FillValue of 9.969209968386869e+36 used,
 'VMNLDF': <class 'netCDF4._netCDF4.Variable'>
 float32 VMNLDF(time, M, NC)
     long_name: Filtered V-velocity
     units: m/s
     grid: grid
     location: edge2
 unlimited dimensions: time
 current shape = (97, 92, 112)
 filling on, default _FillValue of 9.969209968386869e+36 used,
 'VICUV': <class 'netCDF4._netCDF4.Variable'>
 float32 VICUV(time, KMAXOUT_RESTR, M, N)
     long_name: Horizontal eddy viscosity in zeta point
     units: m2/s
     coordinates: XZ YZ
     grid: grid
     location: face
 unlimited dimensions: time
 current shape = (97, 5, 92, 112)
 filling on, default _FillValue of 9.969209968386869e+36 used}

output from reader:

===========================
Reader: ./OD_25x100_Wind05ms_opendrift.nc
Projection: 
  +proj=utm +zone=22 +south +datum=WGS84 +units=m +no_defs 
Coverage: [m]
  xmin: 388208.622339   xmax: 407708.622339   step: 500   numx: 40
  ymin: 6411333.593972   ymax: 6492333.593972   step: 500   numy: 163
  Corners (lon, lat):
    (-52.18, -31.70)  (-51.97, -31.70)
    (-52.19, -32.43)  (-51.98, -32.43)
Vertical levels [m]: 
  [0]
Available time range:
  start: 2020-07-01 00:00:00   end: 2020-07-05 00:00:00   step: 1:00:00
    97 times (0 missing)
Variables:
  time
  x_sea_water_velocity
  y_sea_water_velocity
===========================
 0:00:01.6  total
 0:00:00.2  preparing
 0:00:00.6  reading
 0:00:00.4  interpolation
 0:00:00.0  interpolation_time
 0:00:00.0  masking
knutfrode commented 3 years ago

I cannot see any immediate errors here. Could you double check that you are not swapping x- and y axes somewhere? E.g. utm.to_latlon returns LAT, LON, which you attribute to resp gx_Deg, gy_Deg, but a few lines below you refer to gx_Deg in relation to longitute (although commented out) and not as latitude.

Is there any other software you could use to check that your regridded file is correct?

GuiCruz commented 3 years ago

I think I found the problem.

The problem occurs when I will interpolate the data with the griddata command.

Do you know of any way to interpolate this data correctly?

knutfrode commented 3 years ago

Perhaps you can use LinearNDInterpolator like in reader_netCDF_CF_unstructured.py https://github.com/OpenDrift/opendrift/blob/master/opendrift/readers/reader_netCDF_CF_unstructured.py#L238 An alternative is KDTree.

Note: the unstructured reader is presently beeing improved, so that hopefully in not too long you may use this directly on the Delft3d output.

I am closing now this issue, please reopen if further questions.

GuiCruz commented 3 years ago

Hello again @knutfrode . I'm having some trouble adding a shapefile file to the model. I followed the tutorial, but for some reason it does not appear when plotting the graph and the velocities do not appear in the background too..

o = OceanDrift(loglevel=0)
reader = reader_netCDF_CF_generic.Reader('./OD_25x100_Wind10ms_opendrift.nc')
reader_shp = reader_shape.Reader.from_shpfiles(fn_shp)
reader_landmask = reader_global_landmask.Reader(extent=[np.min(gx_Deg),np.min(gy_Deg),np.max(gx_Deg),np.max(gy_Deg)])  # lonmin, latmin, lonmax, latmax
o.add_reader([reader_landmask, reader_shp, reader])
#lons = np.linspace(130.919571462843436, 130.579213297218644, 100)
#lats = np.linspace(-57.968343361095428, -57.968343361095428, 100)
lats = np.linspace(-31.82541009, -31.82541009, 100)
lons = np.linspace(-52.18633433, -51.97683977, 100)

o.seed_elements(lon=lats, lat=lons,time=datetime.datetime(2020, 7, 3, 0, 0))
o.set_config('general:use_auto_landmask', False)
o.set_config('general:coastline_action', 'stranding')

o.run()

03:24:56 DEBUG: Calculation SRS set to: None
03:24:56 DEBUG: Adding 17 config items from basemodel
03:24:56 DEBUG: Adding 4 config items from basemodel
03:24:56 DEBUG: Adding 34 config items from basemodel
03:24:56 INFO: OpenDriftSimulation initialised (version 1.4.2)
03:24:56 DEBUG: Adding 13 config items from oceandrift
03:24:56 DEBUG:   Overwriting config item seed:z
03:24:56 INFO: Opening dataset: ./OD_25x100_Wind10ms_opendrift.nc
03:24:56 INFO: Opening file with Dataset
03:24:56 DEBUG: Finding coordinate variables.
03:24:56 DEBUG: Parsing variable: WGS84
03:24:56 DEBUG: Parsing CF grid mapping dictionary: {'grid_mapping_name': 'WGS84'}
03:24:56 INFO: Could not parse CF grid_mapping
03:24:56 DEBUG: Parsing variable: depth
03:24:56 DEBUG: Parsing variable: time
03:24:56 DEBUG: Parsing variable: u
03:24:56 DEBUG: Parsing variable: v
03:24:56 DEBUG: Parsing variable: lon
03:24:56 DEBUG: Parsing variable: lat
03:24:56 INFO: No projection found, using lon/lat arrays
03:24:56 DEBUG: Reading lon lat 2D arrays, since projection is not given
03:24:56 INFO: Making Splines for lon,lat to x,y conversion...
03:24:57 DEBUG: Setting buffer size 33 for reader ./OD_25x100_Wind10ms_opendrift.nc, assuming a maximum average speed of 5 m/s.
03:24:57 DEBUG: Reading shapefile: D:/Dropbox/Doutorado/LC08_RS/Ideal_model/Coastline_25x100km.shp
03:24:57 INFO: Pre-processing 2 geometries
ERROR:root:could not verify read permissions for group and others on landmask.
Traceback (most recent call last):
  File "C:\Users\Cliente\anaconda3\envs\opendrift\lib\site-packages\opendrift_landmask_data-0.6-py3.8.egg\opendrift_landmask_data\mask.py", line 77, in __check_permissions__
    if not os.stat(self.lockf).st_mode & 0o777 == 0o777:
FileNotFoundError: [WinError 2] O sistema não pode encontrar o arquivo especificado: 'C:\\Users\\Public\\Documents\\Wondershare\\CreatorTemp\\landmask\\.mask.dat.lock'
03:25:01 DEBUG: Calculation SRS set to: +proj=longlat +ellps=WGS84 +no_defs
03:25:01 DEBUG: Using srs for common grid: +proj=lonlat +ellps=WGS84
03:25:01 DEBUG: Added reader global_landmask
03:25:01 DEBUG: Added reader shape
03:25:01 DEBUG: Setting buffer size 9 for reader ./OD_25x100_Wind10ms_opendrift.nc, assuming a maximum average speed of 1 m/s.
03:25:01 DEBUG: Added reader ./OD_25x100_Wind10ms_opendrift.nc
03:25:01 DEBUG: 
------------------------------------------------------
Software and hardware:
  OpenDrift version 1.4.2
  Platform: Windows, 10
  15.881317138671875 GB memory
  4 processors (Intel64 Family 6 Model 142 Stepping 9, GenuineIntel)
  NumPy version 1.19.4
  SciPy version 1.5.3
  Matplotlib version 3.1.3
  NetCDF4 version 1.5.4
  Python version 3.7.8 | packaged by conda-forge | (default, Nov 17 2020, 23:02:53) [MSC v.1916 64 bit (AMD64)]
------------------------------------------------------

03:25:01 DEBUG: No output file is specified, neglecting export_buffer_length
03:25:01 INFO: Fallback values will be used for the following variables which have no readers: 
03:25:01 INFO:  x_wind: 0.000000
03:25:01 INFO:  y_wind: 0.000000
03:25:01 INFO:  upward_sea_water_velocity: 0.000000
03:25:01 INFO:  ocean_vertical_diffusivity: 0.000000
03:25:01 INFO:  sea_surface_wave_significant_height: 0.000000
03:25:01 INFO:  sea_surface_wave_stokes_drift_x_velocity: 0.000000
03:25:01 INFO:  sea_surface_wave_stokes_drift_y_velocity: 0.000000
03:25:01 INFO:  sea_surface_wave_period_at_variance_spectral_density_maximum: 0.000000
03:25:01 INFO:  sea_surface_wave_mean_period_from_variance_spectral_density_second_frequency_moment: 0.000000
03:25:01 INFO:  surface_downward_x_stress: 0.000000
03:25:01 INFO:  surface_downward_y_stress: 0.000000
03:25:01 INFO:  turbulent_kinetic_energy: 0.000000
03:25:01 INFO:  turbulent_generic_length_scale: 0.000000
03:25:01 INFO:  sea_floor_depth_below_sea_level: 10000.000000
03:25:01 INFO: Duration, steps or end time not specified, running until end of first reader: None
03:25:01 INFO: Duration, steps or end time not specified, running until end of first reader: None
03:25:01 INFO: Duration, steps or end time not specified, running until end of first reader: 2020-07-05 00:00:00
03:25:01 DEBUG: Preparing readers for simulation coverage ([-34.358622637059504, -53.743089459393474, -29.29219904873151, -50.420084216143636]) and time (2020-07-03 00:00:00 to 2020-07-05 00:00:00)
03:25:01 DEBUG:     Preparing global_landmask
03:25:01 DEBUG:     Preparing shape
03:25:01 DEBUG:     Preparing ./OD_25x100_Wind10ms_opendrift.nc
03:25:01 INFO: Using existing reader for land_binary_mask
03:25:01 DEBUG: ----------------------------------------
03:25:01 DEBUG: Variable group ['land_binary_mask']
03:25:01 DEBUG: ----------------------------------------
03:25:01 DEBUG: Calling reader global_landmask
03:25:01 DEBUG: ----------------------------------------
03:25:01 DEBUG: Data needed for 100 elements
03:25:01 DEBUG: Reader time:
        None (before)
        None (after)
03:25:01 DEBUG: Fetching variables from global_landmask
03:25:01 DEBUG: Fetched env-before
03:25:01 DEBUG: No time interpolation needed - right on time.
03:25:01 DEBUG: Reader SRS is the same as calculation SRS - rotation of vectors is not needed.
03:25:01 DEBUG: Obtained data for all elements.
03:25:01 DEBUG: ---------------------------------------
03:25:01 DEBUG: Finished processing all variable groups
03:25:01 DEBUG: ------------ SUMMARY -------------
03:25:01 DEBUG:     land_binary_mask: 0 (min) 0 (max)
03:25:01 DEBUG: ---------------------------------
03:25:01 DEBUG:         0 active elements
03:25:01 INFO: All points are in ocean
03:25:01 DEBUG: to be seeded: 100, already seeded 0
03:25:01 DEBUG: Released 100 new elements.
03:25:01 DEBUG: ======================================================================
03:25:01 INFO: 2020-07-03 00:00:00 - step 1 of 48 - 100 active elements (0 deactivated)
03:25:01 DEBUG: 0 elements scheduled.
03:25:01 DEBUG: ======================================================================

o.plot(background=['x_sea_water_velocity', 'y_sea_water_velocity'])

03:25:13 DEBUG: Adding GSHHS shapes..
03:25:16 DEBUG: Calculating lonlat->xy sequentially
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-55-a503b5efec9d> in <module>()
----> 1 o.plot(background=['x_sea_water_velocity', 'y_sea_water_velocity'])

D:\Dropbox\Doutorado\Jupyter\opendrift\opendrift\models\basemodel.py in plot(self, background, buffer, corners, linecolor, filename, show, vmin, vmax, compare, cmap, lvmin, lvmax, skip, scale, show_scalar, contourlines, trajectory_dict, colorbar, linewidth, lcs, show_particles, show_initial, density_pixelsize_m, surface_color, submerged_color, markersize, title, legend, legend_loc, lscale, fast, hide_landmask, **kwargs)
   3477             else:
   3478                 map_x, map_y, scalar, u_component, v_component = \
-> 3479                     self.get_map_background(ax, background, time=time)
   3480                                         #self.time_step_output)
   3481 

D:\Dropbox\Doutorado\Jupyter\opendrift\opendrift\models\basemodel.py in get_map_background(self, ax, background, time)
   3618                 reader_x, reader_y, u_component, v_component,
   3619                 reader.proj, ccrs.PlateCarree(
-> 3620                 globe=ccrs.Globe(datum='WGS84', ellipse='WGS84')
   3621                 ).proj4_init)
   3622         else:

D:\Dropbox\Doutorado\Jupyter\opendrift\opendrift\readers\basereader.py in rotate_vectors(self, reader_x, reader_y, u_component, v_component, proj_from, proj_to)
    749         if type(proj_from) is not pyproj.Proj:
    750             proj_from = pyproj.Proj('+proj=latlong +R=6370997.0 +ellps=WGS84')
--> 751             reader_x, reader_y = self.xy2lonlat(reader_x, reader_y)
    752         if type(proj_to) is str:
    753             proj_to = pyproj.Proj(proj_to)

D:\Dropbox\Doutorado\Jupyter\opendrift\opendrift\readers\basereader.py in xy2lonlat(self, x, y)
    800 
    801             # NB: mask coordinates outside domain
--> 802             x[x < self.xmin] = np.nan
    803             x[x > self.xmax] = np.nan
    804             y[y < self.ymin] = np.nan

ValueError: cannot convert float NaN to integer