search

OpenDrift / opendrift

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

Telemac 3D + OpenDrift (Readers) #213

Open josefilho77 opened 4 years ago

josefilho77 commented 4 years ago

Hello everyone,

Has anyone already worked with Telemac3D alongwith OpenDrift? I am working on my PhD thesis and recently discovered OpenDrift. I managed to install it and it works fine when I am using the examples.

However, I can not figure out which way would be most straight to use slf as input for OpenDrift...

Do you think that is better to convert slf to NETcdf first and use the Unstructured Reader or is there another way to do that?

I apologize for any incovenience.

Cordially,

José

knutfrode commented 4 years ago

Hi José,

I am not familiar with Telemac nor .slf-files. However, OpenDrift includes one reader for regular CF-compliant netCDF files and one for corresponding unstructured netCDF files. The regular reader is much more mature and robust, whereas the unstructured reader is less developed. I believe it is even returns surface currents only, and not 3D currents.

So there are 3 alternatives:

josefilho77 commented 4 years ago

Dear Knut Frode, Thank you for your kind answer. I am trying to follow your suggestions and write a new .slf reader in 3D dimension using Ugrid convention. Could you explain me which are the steps to achieve that? I am just a new user in Python and started to write a new script however I am stucked in the part of loading data. I mean I can load the data, declarate variables and make a 2D triangulation plot. Which are supposed to be the next steps? Interpolate in all vertical plans, write variables to the new netCDF file?

Thanks in advance,

José

Here is my code until now:

import sys from os import path import numpy as np import csv import datetime from scipy.interpolate import griddata import shapefile from matplotlib.patches import Polygon import matplotlib.pyplot as plt

from utils.files import getFileContent

Activating PYTEL

sys.path.append('/Users/Jose_Filho/scripts/python27')#pytel path

from parsers.parserSELAFIN import SELAFIN,getValueHistorySLF,\ getValuePolylineSLF,subsetVariablesSLF from parsers.parserStrings import parseArrayPaires from parsers.parserSortie import Sortie

""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" INPUT / OUTPUT Folder and TELEMAC RESULT FILE """""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""

Directory

local='//Users//Jose_Filho//opendrift_telemacV3//'

Input file

arquivo=local+'PERNAMBUCO3D_0206.slf'

arquivocosta=local+'Lc_BTS_agua'

Loading selafin

slf = SELAFIN(arquivo)

xmin = np.min(slf.MESHX) xmax = np.max(slf.MESHX) ymin = np.min(slf.MESHY) ymax = np.max(slf.MESHY) x = slf.MESHX y = slf.MESHY mesh = np.array(slf.IKLE2)

elements = np.dstack((slf.MESHX[slf.IKLE2],slf.MESHY[slf.IKLE2])) nomes=slf.VARNAMES# variables names times=slf.tags["times"]

data0slf=datetime.datetime(slf.DATETIME[0],slf.DATETIME[1],slf.DATETIME[2],slf.DATETIME[3],slf.DATETIME[4])#data original do selafin

Creating list of dates

for i,name in enumerate(slf.VARNAMES): print "%i %s" %(i, name)

for i,time in enumerate(slf.tags["times"]): print "%i %s" %(i, time)

values = slf.getVALUES(i) ztri = values[0,0:27470]

fig = plt.figure() ax=fig.add_subplot(111)

zmintri = ztri.min() zmaxtri = ztri.max()

levels = np.arange(zmintri,zmaxtri,(zmaxtri - zmintri)/slf.NPLAN) triplotcontourf = plt.tricontourf(x,y,ztri,levels) plt.colorbar() u = values[:,1] v = values[:,2] cs = plt.quiver(x,y,u,v)

ax.set_xlim(xmin,xmax) ax.set_ylim(ymin,ymax) ax.axis('equal') plt.show() plt.close()

and here is the dict keys of my slf file: 'TITLE', 'fole', 'file', 'MESHX', 'MESHY', 'CLDUNITS', 'CLDNAMES', 'NPOIN2', 'IKLE3', 'alterZnames', 'IPOB3', 'IPOB2', 'VARUNITS', 'neighbours', 'NDP2', 'NDP3', 'NBV2', 'NPLAN', 'NBV1', 'VARINDEX', 'tags', 'NVAR', 'IPARAM', 'DATETIME', 'edges', 'IKLE2', 'NPOIN3', 'NELEM3', 'NELEM2', 'tree', 'VARNAMES'

gauteh commented 3 years ago

See #243. There is now an implementation of a reader for FVCOM unstructured files that might be of use for implementing a reader for Telemac.

Boorhin commented 3 years ago

Just bouncing on this thread rather than the other one. The pytel script here is outdated as the functions have changed (Telemac v8p2) FVCOM is a finite volume model and Telemac a finite element one so you would not be able to have natively the same grids. Of course one could interpolate the data from nodes to volumes etc. I am interested in developing that aspect

Boorhin commented 3 years ago

This is just a little sneak peek of what I have been writing to convert a 3D Telemac slf file (unstructured grid) into a netcdf trying to follow convention 1.8 and ugrid. I have no certainty it is the way it should but here is where I am at. I need to check if it is correct regarding the standards and if it will work with OpenDrift. I will put it on my Github when I have a minute. There are lots of comments because there are lots of things that I don't think are useful (up to you to check)

###
# Exporting a Telemac 3D selafin file to an unstructured Netcdf file
# Following CF 1-8 and ugrid conventions
import numpy as np
import netCDF4 as nC
from datetime import datetime
import sys
from os import path, environ, sep
from data_manip.formats.selafin import Selafin

# Selafin variables
local='/home/julien/DATA/Marine plastics/PlasticsATBay/Simulations/oceano/NW_Scot/'#Directory
f=local+'r3d_tide-ES_real_gen_new_7.slf' #Input file
slf = Selafin(f)

steps= len(slf.tags['times'])
timestep= slf.tags['times'][1]
starttime=datetime(slf.datetime[0],slf.datetime[1],slf.datetime[2],slf.datetime[3],slf.datetime[4])

d = nC.Dataset('somefile.nc', 'w')

d.setncattr('Conventions','CF-1.8') #https://cfconventions.org/Data/cf-conventions/cf-conventions-1.8/cf-conventions.pdf
d.setncattr('title',slf.title)
d.setncattr('history',str(datetime.now())[:16])
d.setncattr('institution',"Plastic@bay - MTS-CFD") #Specifies where the original data was produced
d.setncattr('source',"Telemac 3D v8p2") # The  method  of  production  of  the  original  data.
d.setncattr('CoordinateSystem', 'Cartesian')
d.setncattr('CoordinateProjection', 'UTM30')

###make topology
node= d.createDimension('nMesh2_node',slf.npoin2)
nele= d.createDimension('nMesh2_face',slf.nelem2)
maxnode=d.createDimension('nMaxMesh2_face_nodes',slf.ndp2)
siglay= d.createDimension('Mesh2_layers',slf.nplan)
time= d.createDimension('time',None)#len(slf.tags['times'])
Two = d.createDimension('Two',2)

# time constrain
Mesh2_time= d.createVariable('time',np.float,('time'))
Mesh2_time.long_name='time'
Mesh2_time.standard_name='time'
Mesh2_time.units='seconds since '+str(starttime)
Mesh2_time.base_date=str(starttime)
Mesh2_time[:] = slf.tags['times']

## variables.
# Mesh topology
Mesh2= d.createVariable('Mesh2', np.int)
Mesh2.cf_role = "mesh_topology"
Mesh2.long_name = "Topology data of 2D unstructured mesh"
Mesh2.topology_dimension = 2
Mesh2.node_coordinates = "Mesh2_node_x Mesh2_node_y"
Mesh2.face_node_connectivity = "Mesh2_face_nodes"
Mesh2.face_dimension = "nele"

Mesh2_face_nodes =  d.createVariable('Mesh2_face_nodes',np.int,('nMesh2_face','nMaxMesh2_face_nodes'), fill_value=999999)
Mesh2_face_nodes.cf_role = "face_node_connectivity"
Mesh2_face_nodes.long_name = "Maps every face to its corner nodes."
Mesh2_face_nodes.start_index = 0
Mesh2_face_nodes[:]=slf.ikle2

#creating a var of projection
Mesh2_proj = d.createVariable('crs',np.int)
Mesh2_proj.grid_mapping_name = 'transverse_mercator'
Mesh2_proj.scale_factor_at_central_meridian =0.9996
Mesh2_proj.longitude_of_central_meridian = 3
Mesh2_proj.latitude_of_projection_origin = 0
Mesh2_proj.false_easting =500000.

# Mesh node coordinates (needs changed from crs)
Mesh2_node_x = d.createVariable('Mesh2_node_x', np.float, ('nMesh2_node'))
Mesh2_node_x.axis="X"
Mesh2_node_x.long_name = "x-coordinate in Cartesian system"
Mesh2_node_x.units = "m"
Mesh2_node_x[:]=slf.meshx

Mesh2_node_y = d.createVariable('Mesh2_node_y', np.float, ('nMesh2_node'))
Mesh2_node_y.axis="Y"
Mesh2_node_y.long_name = "y-coordinate in Cartesian system"
Mesh2_node_y.units = "m"
Mesh2_node_y[:]=slf.meshy
# Vertical coordinate
Mesh2_layers= d.createVariable('Mesh2_layers', np.float,('Mesh2_layers'))
Mesh2_layers.standard_name = "ocean_sigma_coordinate"
Mesh2_layers.long_name = "sigma at layer midpoints"
Mesh2_layers.positive = "up"
Mesh2_layers.formula_terms = "sigma. Mesh2_layers eta. Mesh2_surface depth. Mesh2_depth"
Mesh2_layers[:]= np.linspace(0,1, slf.nplan) ###<= needs to be checked for FE

Mesh2_depth= d.createVariable('Mesh2_depth', np.float,('nMesh2_node'))
Mesh2_depth.standard_name = "sea_floor_depth_below_geoid"
Mesh2_depth.units = "m"
Mesh2_depth.positive = "up"
Mesh2_depth.mesh = "Mesh2"
Mesh2_depth.location = "node"
Mesh2_depth.coordinates = "time Mesh2_node_x Mesh2_node_y"
for i in range(steps):
    Mesh2_depth[i]=slf.get_values(i)[0][:slf.npoin2]

Mesh2_surface= d.createVariable('Mesh2_surface', np.float,('nMesh2_node'))
Mesh2_surface.standard_name = "sea_surface_height_above_geoid"
Mesh2_surface.units = "m"
Mesh2_surface.mesh = "Mesh2"
Mesh2_surface.location = "node"
Mesh2_surface.coordinates = "time Mesh2_node_x Mesh2_node_y"
for i in range(steps):
    Mesh2_surface[i]=slf.get_values(i)[0][-slf.npoin2:]

#creating variable currents X
Mesh2_u = d.createVariable('u','f4',('time','Mesh2_layers','nMesh2_node'))
Mesh2_u.standard_name='eastward_sea_water_velocity'
Mesh2_u.units='m s-1'
Mesh2_u.coordinates = 'time siglay Mesh2_node_x Mesh2_node_y'
for i in range(steps):
    Mesh2_u[i]=slf.get_values(i)[1]

#creating variable currents Y
Mesh2_v = d.createVariable('v','f4',('time','Mesh2_layers','nMesh2_node'))
Mesh2_v.standard_name='northward_sea_water_velocity'
Mesh2_v.units='m s-1'
Mesh2_v.coordinates = 'time siglay Mesh2_node_x Mesh2_node_y'
for i in range(steps):
    Mesh2_v[i]=slf.get_values(i)[2]

#creating variable currents Z
Mesh2_w = d.createVariable('w','f4',('time','Mesh2_layers','nMesh2_node'))
Mesh2_w.standard_name='upward_sea_water_velocity'
Mesh2_w.units='m s-1'
Mesh2_w.coordinates = 'time siglay Mesh2_node_x Mesh2_node_y'
for i in range(steps):
    Mesh2_w[i]=slf.get_values(i)[3]

d.close()
### Optional stuff not used
#Mesh2.edge_node_connectivity = "Mesh2_edge_nodes"  // attribute required if variables will be defined on edges
#Mesh2.edge_dimension = "nMesh2_edge"
#Mesh2.edge_coordinates = "Mesh2_edge_x Mesh2_edge_y"  // optional attribute (requires edge_node_connectivity)
#Mesh2.face_coordinates = "Mesh2_face_x Mesh2_face_y"  // optional attribute
#Mesh2.face_edge_connectivity = "Mesh2_face_edges"  // optional attribute (requires edge_node_connectivity)
#Mesh2.face_face_connectivity = "Mesh2_face_links"  // optional attribute
#Mesh2.edge_face_connectivity = "Mesh2_edge_face_links"  // optional attribute (requires edge_node_connectivity)
# integer Mesh2_edge_nodes(nMesh2_edge, Two)
# Mesh2_edge_nodes.cf_role = "edge_node_connectivity"
# Mesh2_edge_nodes.long_name = "Maps every edge to the two nodes that it connects."
# Mesh2_edge_nodes.start_index = 1

# // Optional mesh topology variables
# integer Mesh2_face_edges(nMesh2_face, nMaxMesh2_face_nodes)
# Mesh2_face_edges.cf_role = "face_edge_connectivity"
# Mesh2_face_edges.long_name = "Maps every face to its edges."
# Mesh2_face_edges._FillValue = 999999
# Mesh2_face_edges.start_index = 1
# integer Mesh2_face_links(nMesh2_face, nMaxMesh2_face_nodes)
# Mesh2_face_links.cf_role = "face_face_connectivity"
# Mesh2_face_links.long_name = "neighbor faces for faces"
# Mesh2_face_links.start_index = 1
# Mesh2_face_links._FillValue = -999
# Mesh2_face_links.comment = "missing edges as well as missing neighbor faces are indicated using _FillValue"
# integer Mesh2_edge_face_links(nMesh2_edge, Two)
# Mesh2_edge_face_links.cf_role = "edge_face_connectivity"
# Mesh2_edge_face_links.long_name = "neighbor faces for edges"
# Mesh2_edge_face_links.start_index = 1
# Mesh2_edge_face_links._FillValue = -999
# Mesh2_edge_face_links.comment = "missing neighbor faces are indicated using _FillValue"
# // Optional mesh face and edge coordinate variables
# double Mesh2_face_x(nMesh2_face)
# Mesh2_face_x.standard_name = "longitude"
# Mesh2_face_x.long_name = "Characteristics longitude of 2D mesh face."
# Mesh2_face_x.units = "degrees_east"
# Mesh2_face_x.bounds = "Mesh2_face_xbnds"
# double Mesh2_face_y(nMesh2_face)
# Mesh2_face_y.standard_name = "latitude"
# Mesh2_face_y.long_name = "Characteristics latitude of 2D mesh face."
# Mesh2_face_y.units = "degrees_north"
# Mesh2_face_y.bounds = "Mesh2_face_ybnds"
# double Mesh2_face_xbnds(nMesh2_face,nMaxMesh2_face_nodes)
# Mesh2_face_xbnds.standard_name = "longitude"
# Mesh2_face_xbnds.long_name = "Longitude bounds of 2D mesh face (i.e. corner coordinates)."
# Mesh2_face_xbnds.units = "degrees_east"
# Mesh2_face_xbnds._FillValue = 9.9692099683868690E36
# double Mesh2_face_ybnds(nMesh2_face,nMaxMesh2_face_nodes)
# Mesh2_face_ybnds.standard_name = "latitude"
# Mesh2_face_ybnds.long_name = "Latitude bounds of 2D mesh face (i.e. corner coordinates)."
# Mesh2_face_ybnds.units = "degrees_north"
# Mesh2_face_ybnds._FillValue = 9.9692099683868690E36
# double Mesh2_edge_x(nMesh2_edge)
# Mesh2_edge_x.standard_name = "longitude"
# Mesh2_edge_x.long_name = "Characteristic longitude of 2D mesh edge (e.g. midpoint of the edge)."
# Mesh2_edge_x.units = "degrees_east"
# double Mesh2_edge_y(nMesh2_edge)
# Mesh2_edge_y.standard_name = "latitude"
# Mesh2_edge_y.long_name = "Characteristic latitude of 2D mesh edge (e.g. midpoint of the edge)."
# Mesh2_edge_y.units = "degrees_north"
# // bounds variables for edges skipped
josefilho77 commented 3 years ago

Hi, Boorhin. Thank you for your contribution. I am also trying to do this task for my T3D results...I managed to convert and interpolate the unstructured data from t3d into a regular grid and write it into a Ncdf file which works nicely with opendrift regular netcdf reader . Let me know if do you need a notebook of that.

Meanwhile, I suggest you to use the most recent python api from Telemac, which includes Telemac_File.py (extraction scripts). This script gives you information you need to write the unstructured data, including the triangulation aspects.

Let me know if do you want to develop it along with me.

Best regards

Boorhin commented 3 years ago

Thanks I will have a look at Telemac_file.py. v8p2 is the latest stable version as far as I know I am not interested in structured grids as it is not manageable in my study. I would advise you to read the netcdf and ugrid conventions as requested by the dev. Maybe share your script on github so we could improve it from both sides of Opendrift and Telemac?