# -*- coding: utf-8 -*-
"""
Functions for extracting velocity, velocity gradient, and vorticity 2D fields
from Open Telemac-Mascret model output.
Chris Old
IES, School of Engineering, University of Edinburgh
Sep 2023
"""
# ----------------------------------------------------------------------------
# IMPORTS
# ----------------------------------------------------------------------------
# Standard Python Dependencies
import os
import sys
import time
# Non-Standard Python Dependencies
import numpy as np
# Local Module Dependencies
from ocmw.core.fileTools import generateGlobalAttributes
from ocmw.core.similitude import courantNumber, froudeNumber, ekmanNumber, rossbyNumber
from ocmw.core.physics import circulation2D
from ocmw.core.in_poly_funcs import cn_PnPoly
from ocmw.core.meshTools import getNodesAroundNode, gradientAtNode, weightsMatrix
from .otm_extraction import OTMDate2Num, OTMDate2Str, genOTMTimeStamps, default_epoch
from .otm_extraction import getSlfNodesWeights, extractMesh, enhanceMesh, slfObj
from .ocmw_extract import save_ocmw_file, load_ocmw_file, ocmwExtn
# Other Dependencies
#--------------------------------------------------------------------------
# GLOBAL CONSTANTS
#--------------------------------------------------------------------------
# ----------------------------------------------------------------------------
# CLASS DEFINITIONS
# ----------------------------------------------------------------------------
#--------------------------------------------------------------------------
# FUNCTION DEFINITIONS
#--------------------------------------------------------------------------
[docs]
def initialiseVectorField(slf):
"""
Get Selafin structure data to set parameters for the vector field data
Parameters
----------
slf : Selafin obj
Selafin file handle.
Returns
-------
times : numpy.ndarray, float
Timestamp as an offset in seconds from idate.
idate : numpy.ndarray, int32
Start time as an array of integers [year,month,day,hour,minute,second].
nsteps : int
Number of time steps in file.
nnodes : int
Number of mesh nodes.
meshx : numpy.ndarray, float
X coordinate of mesh node locations.
meshy : numpy.ndarray, float
Y coordinate of mesh node locations.
elems : numpy.ndarray, int32
Indices of nodes that define each mesh triangular element.
edges : numpy.ndarray, int32
Indices of node that define every element edge.
"""
slf.set_mpl_tri(True)
slf.set_kd_tree(True)
times = slf.tags['times']
idate = slf.datetime
nsteps = len(times)
nnodes = slf.npoin2
meshx = slf.meshx
meshy = slf.meshy
elems = slf.ikle3
edges = slf.edges
return times,idate,nsteps,nnodes,meshx,meshy,elems,edges
[docs]
def getVariableAtNodes(slf,nodes,varIndex):
"""
Get timeseries of data for the requested nodes and field variable
Parameters
----------
slf : Selfin obj
Handle to the Selafin file being processed.
nodes : numpy.ndarray, int
Indices of nodes to be processed.
varIndex : int
Index of OTM variable in Selafin data structure.
Returns
-------
var : numpy.ndarray
Array containing variable timeseries for all requested nodes.
"""
var = np.squeeze(slf.get_series(nodes+1,[varIndex],False))
return var
[docs]
def getVariablesAtTime(slf,frame,varIndex):
"""
Get variable data for all mesh nodes at the requested time step
Parameters
----------
slf : Selfin obj
Handle to the Selafin file being processed.
frame : int
Index of the required model time step.
varIndex : int
Index of OTM variable in Selafin data structure.
Returns
-------
vars : numpy.ndarray, float
Array of variable data for requested time step.
"""
vars = slf.getVariablesAt(frame,varIndex)
return vars
[docs]
def vectorFieldAtLocation(slf,Loc,d_off=0.0):
"""
Get the 2-D vector field data at a specifc geospatial location
Parameters
----------
slf : Selfin obj
Handle to the Selafin file being processed.
Loc : list, float
Geospatial coordinates of location [x,y].
d_off : float, optional
Mean sea level correction if used in model configuration. The default is 0.0.
Returns
-------
vect : dict
Dictionary strcutre containing the vector field data at the location.
"""
NNghb, NWghts, Locs = getSlfNodesWeights(slf,Loc[0],Loc[1])
npts = len(NNghb)
times,idate,nsteps,nnodes,meshx,meshy,elems,edges = initialiseVectorField(slf)
elev = np.zeros([nsteps,npts])
velx = np.zeros([nsteps,npts])
vely = np.zeros([nsteps,npts])
dpth = np.zeros([nsteps,npts])
dhdx = np.zeros([nsteps,npts])
dhdy = np.zeros([nsteps,npts])
dudx = np.zeros([nsteps,npts])
dudy = np.zeros([nsteps,npts])
dvdx = np.zeros([nsteps,npts])
dvdy = np.zeros([nsteps,npts])
for inode in np.arange(npts):
node = NNghb[inode]
nodes = getNodesAroundNode(node,edges)
nodesX = meshx[nodes]
nodesY = meshy[nodes]
indx, = np.where(nodes == node)
u = getVariableAtNodes(slf,nodes,0)
v = getVariableAtNodes(slf,nodes,1)
h = getVariableAtNodes(slf,nodes,3)
b = getVariableAtNodes(slf,nodes,4)
dpth[:,inode] = d_off - b[indx,:]
elev[:,inode] = h[indx,:]
velx[:,inode] = u[indx,:]
vely[:,inode] = v[indx,:]
for istep in np.arange(nsteps):
val = h[:,istep]
grad,res,rnk,sng = gradientAtNode(node,nodes,nodesX,nodesY,val)
dhdx[istep,inode] = grad[0]
dhdy[istep,inode] = grad[1]
val = u[:,istep]
grad,res,rnk,sng = gradientAtNode(node,nodes,nodesX,nodesY,val)
dudx[istep,inode] = grad[0]
dudy[istep,inode] = grad[1]
val = v[:,istep]
grad,res,rnk,sng = gradientAtNode(node,nodes,nodesX,nodesY,val)
dvdx[istep,inode] = grad[0]
dvdy[istep,inode] = grad[1]
del node, nodes, u, v, h
vort = dvdx - dudy
# Interpolate onto location of interest
wgts = weightsMatrix(NWghts,nsteps)
epochStr = OTMDate2Str(default_epoch)
vect = {}
vect['loc'] = Loc
vect['nearestNodes'] = NNghb
vect['nodeWeights'] = NWghts
vect['nodeLocs'] = Locs
vect['timeUnits'] = 'days since '+epochStr
vect['epoch'] = epochStr
vect['times'] = genOTMTimeStamps(idate, times)
vect['elev'] = elev
vect['dpth'] = dpth
vect['U'] = velx
vect['V'] = vely
vect['dhdx'] = dhdx
vect['dhdy'] = dhdy
vect['dudx'] = dudx
vect['dudy'] = dudy
vect['dvdx'] = dvdx
vect['dvdy'] = dvdy
vect['vort'] = vort
vect['elev_loc'] = np.sum(elev.transpose()*wgts,0)
vect['dpth_loc'] = np.sum(dpth.transpose()*wgts,0)
vect['U_loc'] = np.sum(velx.transpose()*wgts,0)
vect['V_loc'] = np.sum(vely.transpose()*wgts,0)
vect['dhdx_loc'] = np.sum(dhdx.transpose()*wgts,0)
vect['dhdy_loc'] = np.sum(dhdy.transpose()*wgts,0)
vect['dudx_loc'] = np.sum(dudx.transpose()*wgts,0)
vect['dudy_loc'] = np.sum(dudy.transpose()*wgts,0)
vect['dvdx_loc'] = np.sum(dvdx.transpose()*wgts,0)
vect['dvdy_loc'] = np.sum(dvdy.transpose()*wgts,0)
vect['vort_loc'] = np.sum(vort.transpose()*wgts,0)
return vect
[docs]
def characteriseAtNode(slf,node):
"""
Calculate a timeseries of characteristic parameters for a given model node
Parameters
----------
slf : Selfin obj
Handle to the Selafin file being processed.
node : int
Node index.
Returns
-------
vect : dict
Dictionary structure containing timeseries of characteristic parameter
at the node requested.
"""
times,idate,nsteps,nnodes,meshx,meshy,elems,edges = initialiseVectorField(slf)
Cr = np.zeros([nsteps,])
Fr = np.zeros([nsteps,])
Ek = np.zeros([nsteps,])
Ro = np.zeros([nsteps,])
nodes = getNodesAroundNode(node,edges)
nodesX = meshx[nodes]
nodesY = meshy[nodes]
indx, = np.where(nodes == node)
u = getVariableAtNodes(slf,nodes,0)
v = getVariableAtNodes(slf,nodes,1)
d = getVariableAtNodes(slf,nodes,2)
dpth = np.squeeze(d[indx,:])
velx = np.squeeze(u[indx,:])
vely = np.squeeze(v[indx,:])
vel = np.sqrt(np.square(velx)+np.square(vely))
for istep in np.arange(nsteps):
val = u[:,istep]
grad,res,rnk,sng = gradientAtNode(node,nodes,nodesX,nodesY,val)
dudy = grad[1]
val = v[:,istep]
grad,res,rnk,sng = gradientAtNode(node,nodes,nodesX,nodesY,val)
dvdx = grad[0]
omega = 0.5 * (dvdx - dudy)
Cr[istep] = courantNumber(node, nodes, nodesX, nodesY, velx[istep], vely[istep], 10.0)
Fr[istep] = froudeNumber(dpth[istep], vel[istep])
Ek[istep] = ekmanNumber(dpth[istep], omega)
Ro[istep] = rossbyNumber(dpth[istep], omega)
# Store results in dictionary structure
vect = {}
vect['node'] = node
vect['nearestNodes'] = nodes
vect['nodesX'] = nodesX
vect['nodesY'] = nodesY
vect['times'] = np.reshape(OTMDate2Num(idate) + (times / (24.*60.0*60.)),[nsteps,1])
vect['Cr'] = Cr
vect['Fr'] = Fr
vect['Ek'] = Ek
vect['Ro'] = Ro
return vect
[docs]
def getVelocityField(slfHnd, mesh, dispFreq=5000):
"""
Generate a series of single time 2D velocity field data for animations
Parameters
----------
slf : Selfin obj
Handle to the Selafin file being processed.
mesh : dict
Dictionary structure containing enhanced mesh
dispFreq : float, optional
How often to display cell processing progress. The default is 5000.
Returns
-------
fname : str
Name of model output file processed.
"""
# Construct file name parts
fparts = slfHnd.file['name'].replace('\\','/').split('/')
pidx = fparts.index('RESULTS')
fpath = "/".join(fparts[0:pidx])+'/VELOCITY'
if not os.path.isdir(fpath):
os.mkdir(fpath)
fname = fparts[-1].split('.')[0]
# Get data dimensions
nNode = slfHnd.npoin2
nElem = slfHnd.nelem3
elems = slfHnd.ikle3.copy()
nStep = len(slfHnd.tags['times'])
# Initialise data structures
print('Initializing variables...')
velx_node = np.zeros([nNode,nStep],dtype=np.single)
vely_node = np.zeros([nNode,nStep],dtype=np.single)
vmag_node = np.zeros([nNode,nStep],dtype=np.single)
# Load velocity data
print('Loading velocity data...')
for istp in np.arange(nStep, dtype=np.int64):
velx_node[:,istp] = slfHnd.get_variables_at(istp,[0])[0]
vely_node[:,istp] = slfHnd.get_variables_at(istp,[1])[0]
vmag_node[:,istp] = np.sqrt(np.square(velx_node[:,istp])+
np.square(vely_node[:,istp]))
print('Writing node data to individual frame files')
for istp in np.arange(nStep, dtype=np.int64):
ocmwFile = fname+"_vels_nodes_step"+str(istp).zfill(4)+".mat"
vel_nodes = {}
vel_nodes['velx'] = velx_node[:,istp]
vel_nodes['vely'] = vely_node[:,istp]
vel_nodes['vmag'] = vmag_node[:,istp]
save_ocmw_file(fpath,ocmwFile,vel_nodes)
# Interploate velocity data onto element centres
print('Interpolate onto cell centres...')
# Get centroid B weights from mesh
if 'cWghts' not in mesh.keys():
mesh = enhanceMesh(mesh)
cWghts = mesh['cWghts'].copy()
velx_cntr = np.zeros([nElem,nStep],dtype=np.single)
vely_cntr = np.zeros([nElem,nStep],dtype=np.single)
vmag_cntr = np.zeros([nElem,nStep],dtype=np.single)
print('Number of cells to process = ',nElem)
for ielm in np.arange(nElem, dtype=np.int64):
cw = cWghts[ielm,:]
v = velx_node[elems[ielm,:],:]
velx_cntr[ielm,:] = np.sum((v.transpose()*cw),1)
v = vely_node[elems[ielm,:],:]
vely_cntr[ielm,:] = np.sum((v.transpose()*cw),1)
v = vmag_node[elems[ielm,:],:]
vmag_cntr[ielm,:] = np.sum((v.transpose()*cw),1)
if not(np.mod(ielm+1,dispFreq)):
sys.stdout.write("\r"+"cells processed = "+str(ielm+1))
sys.stdout.flush()
print("\r"+"cells processed = "+str(nElem), flush=True)
print('Writing element data to individual frame files')
for istp in np.arange(nStep, dtype=np.int64):
ocmwFile = fname+"_vels_elems_step"+str(istp).zfill(4)+".mat"
vel_elems = {}
vel_elems['velx'] = velx_cntr[:,istp]
vel_elems['vely'] = vely_cntr[:,istp]
vel_elems['vmag'] = vmag_cntr[:,istp]
save_ocmw_file(fpath,ocmwFile,vel_elems)
print('Process complete.')
return fname
[docs]
def frameSetIndices(nFrames,nCores,setNum):
"""
Get subset of time steps for batch process generation of animation frames
from 2D model fields.
Parameters
----------
nFrames : int
Number of frames in 2D Selafin file.
nCores : int
Number of cores to spread processing across.
setNum : int
Sube set to get incides for.
Returns
-------
frameSet : list
Time stpe indices for frames to be generated.
"""
if nCores > 1:
nfrm = nFrames//nCores
nrem = np.mod(nFrames,nCores)
if setNum < nrem:
frameSet = (np.arange(nfrm+1)+setNum*(nfrm+1)).tolist()
else:
frameSet = (np.arange(nfrm)+nrem*(nfrm+1)+(setNum-nrem)*nfrm).tolist()
else:
frameSet = (np.arange(nFrames)).tolist()
return frameSet
[docs]
def getCirculationField(slfHnd, mesh, tsteps=[], dispFreq=5000):
"""
Generate a series of single time 2D circulation data files for animations.
Parameters
----------
slf : Selfin obj
Handle to the Selafin file being processed.
mesh : dict
Dictionary structure containing enhanced mesh
dispFreq : float, optional
How often to display cell processing progress. The default is 5000.
Returns
-------
fname : str
Name of model output file processed.
"""
# Construct file name parts
fparts = slfHnd.file['name'].replace('\\','/').split('/')
pidx = fparts.index('RESULTS')
fpath = "/".join(fparts[0:pidx])+'/VELOCITY'
if not os.path.isdir(fpath):
os.mkdir(fpath)
fname = fparts[-1].split('.')[0]
# Get data dimensions
times = slfHnd.tags['times']
idate = slfHnd.datetime
nStep = len(times)
meshx = slfHnd.meshx
meshy = slfHnd.meshy
nElem = slfHnd.nelem3
elems = slfHnd.ikle3.copy()
print('Number of frames to process = '+str(len(tsteps))+'\n\n')
# Calculate element circulation
print('Calculate cell circulation...')
if len(tsteps) == 0:
tsteps = np.arange(nStep, dtype=np.int64).tolist()
circ = np.zeros([nElem,],dtype=np.single)
#for istp in np.arange(nStep, dtype=np.int64):
for istp in tsteps:
tic = time.time()
print('timestep '+str(istp)+':')
t = genOTMTimeStamps(idate, times[istp])
print('Loading velocity data...')
velx_node = slfHnd.get_variables_at(istp,[0])[0]
vely_node = slfHnd.get_variables_at(istp,[1])[0]
print('Number of cells to process = ',nElem)
for ielm in np.arange(nElem, dtype=np.int64):
x = meshx[elems[ielm,:]].tolist()
y = meshy[elems[ielm,:]].tolist()
vx = velx_node[elems[ielm,:]].tolist()
vy = vely_node[elems[ielm,:]].tolist()
circ[ielm] = circulation2D(x, y, vx, vy)
if not(np.mod(ielm+1,dispFreq)):
sys.stdout.write("\r"+"cells processed = "+str(ielm+1))
sys.stdout.flush()
print("\r"+"cells processed = "+str(nElem), flush=True)
ocmwFile = fname+"_circ_elems_step"+str(istp).zfill(4)+ocmwExtn
epochStr = OTMDate2Str(default_epoch)
circ_elems = {}
circ_elems['time'] = t
circ_elems['timeUnits'] = 'days since '+epochStr
circ_elems['epoch'] = epochStr
circ_elems['mesh'] = mesh
circ_elems['circ'] = circ
circ_elems['circ_dens'] = circ/mesh['area']
save_ocmw_file(fpath,ocmwFile,circ_elems)
toc = time.time() - tic
print('Time Elasped = ',toc,' s')
print('Process complete.')
return fname
[docs]
def mergeCirculationFrames(dataPath,fileList):
"""
Merge a set of circulation data frames into a single data strcuture.
Parameters
----------
dataPath : str
Path to location of files to be merged.
fileList : list, str
List of file names for the files to be merged.
Returns
-------
circ_frames : dict
Dictionary data strcutre containing the merged data frames.
"""
nframes = len(fileList)
epochStr = OTMDate2Str(default_epoch)
for indx,f in enumerate(fileList):
data = load_ocmw_file(dataPath,f)
if indx == 0:
t = np.zeros((nframes,))
circ = np.zeros((data['circ'].shape[0],nframes))
circd = np.zeros((data['circ_dens'].shape[0],nframes))
t[indx] = data['time']
circ[:,indx] = data['circ']
circd[:,indx] = data['circ_dens']
circ_frames = {}
circ_frames['time'] = t
if 'timeUnits' in data.keys():
circ_frames['timeUnits'] = data['timeUnits']
else:
circ_frames['timeUnits'] = 'days since '+epochStr
if 'epoch' in data.keys():
circ_frames['epoch'] = data['epoch']
else:
circ_frames['epoch'] = epochStr
circ_frames['mesh'] = data['mesh']
circ_frames['time'] = t
circ_frames['circ'] = circ
circ_frames['circ_dens'] = circd
return circ_frames
[docs]
def processNodeSet(fPath,fName,nCores=1,setNum=0):
"""
Calculate field gradients and vorticty for a set of nodes in a batch
process and save to a set file for merging.
Parameters
----------
fPath : str
DESCRIPTION.
fName : str
Name of Selafin file containing 2D field data.
nCores : int, optional
Number of cores to split processing across. The default is 1.
setNum : int, optional
Index of batch processing set. The default is 0.
Returns
-------
None.
"""
slf = slfObj(os.path.join(fPath,fName).replace('\\','/'))
times,idate,nsteps,nnodes,meshx,meshy,elems,edges = initialiseVectorField(slf)
t = np.reshape(OTMDate2Num(idate) + (times / (24.*60.0*60.)),[nsteps,1])
npts = np.int64(np.floor(slf.npoin3/nCores))
node0 = setNum*npts
# Account for extra points to last core
if np.mod(slf.npoin3,nCores) != 0:
if (setNum+1) == nCores:
npts = npts + np.mod(slf.npoin3,nCores)
print(' Number of Nodes to Process: '+str(npts))
print(' First Node in Set: '+str(node0))
if nCores == 1:
ocmwFile = slf.file['name'].split('/')[-1].split('.')[0]+'_gradients.mat'
else:
ocmwFile = slf.file['name'].split('/')[-1].split('.')[0]+'_gradients_'+'SET'+str(setNum).zfill(4)+'.mat'
elev = np.zeros([nsteps,npts])
velx = np.zeros([nsteps,npts])
vely = np.zeros([nsteps,npts])
dpth = np.zeros([nsteps,npts])
dhdx = np.zeros([nsteps,npts])
dhdy = np.zeros([nsteps,npts])
dudx = np.zeros([nsteps,npts])
dudy = np.zeros([nsteps,npts])
dvdx = np.zeros([nsteps,npts])
dvdy = np.zeros([nsteps,npts])
tic = time.time()
for inode in np.arange(npts):
node = node0 + inode
nodes = getNodesAroundNode(node,edges)
nodesX = meshx[nodes]
nodesY = meshy[nodes]
indx = np.where(nodes == node)
u = getVariableAtNodes(slf,nodes,0)
v = getVariableAtNodes(slf,nodes,1)
h = getVariableAtNodes(slf,nodes,3)
b = getVariableAtNodes(slf,nodes,4)
dpth[:,inode] = b[indx,:]
elev[:,inode] = h[indx,:]
velx[:,inode] = u[indx,:]
vely[:,inode] = v[indx,:]
for istep in np.arange(nsteps):
grad,res,rnk,sng = gradientAtNode(node,nodes,nodesX,nodesY,h[:,istep])
dhdx[istep,inode] = grad[0]
dhdy[istep,inode] = grad[1]
grad,res,rnk,sng = gradientAtNode(node,nodes,nodesX,nodesY,u[:,istep])
dudx[istep,inode] = grad[0]
dudy[istep,inode] = grad[1]
grad,res,rnk,sng = gradientAtNode(node,nodes,nodesX,nodesY,v[:,istep])
dvdx[istep,inode] = grad[0]
dvdy[istep,inode] = grad[1]
if (node > 0) & (np.mod(node,1000)==0):
print('Processed '+str(node+1)+' nodes : '+str(time.time()-tic)+' s')
del node, nodes, u, v, h
vort = dvdx - dudy
toc = time.time() - tic
print('Time Elasped = ',toc,' s')
# Close Selafin file
slf.__del__()
# Store and save results
dataDict = {}
dataDict['times'] = t
dataDict['nsteps'] = nsteps
dataDict['nodes'] = node0+np.arange(npts)
dataDict['meshx'] = meshx[node0+np.arange(npts)].copy()
dataDict['meshy'] = meshy[node0+np.arange(npts)].copy()
dataDict['elev'] = elev
dataDict['dpth'] = dpth
dataDict['U'] = velx
dataDict['V'] = vely
dataDict['dhdx'] = dhdx
dataDict['dhdy'] = dhdy
dataDict['dudx'] = dudx
dataDict['dudy'] = dudy
dataDict['dvdx'] = dvdx
dataDict['dvdy'] = dvdy
dataDict['vort'] = vort
extractPath = fPath.replace('RESULTS','EXTRACT')
save_ocmw_file(extractPath,ocmwFile,dataDict)
return
[docs]
def processNodesInPolygon(fPath,fName,polygon):
"""
Calculate field gradients and vorticty at nodes inside a polygon and save
to an OCMW format output file.
Parameters
----------
fPath : str
DESCRIPTION.
fName : str
Name of Selafin file containing 2D field data.
polygon : list, float
List of coordinates for vertices of the polygon [[x,y],[x,y],...].
Returns
-------
None.
"""
slf = slfObj(os.path.join(fPath,fName).replace('\\','/'))
times,idate,nsteps,nnodes,meshx,meshy,elems,edges = initialiseVectorField(slf)
mesh = extractMesh(slf, saveMesh=False)
points = np.zeros((nnodes,2),dtype=float)
points[:,0] = meshx
points[:,1] = meshy
inside = [cn_PnPoly(point,polygon) for point in points]
nodes = np.where(inside)[0]
npts = np.int64(len(nodes))
print(' Number of Nodes to Process: '+str(npts))
ocmwFile = slf.file['name'].split('/')[-1].split('.')[0].replace('_2D','_poly_dyn2d.mat')
elev = np.zeros([nsteps,npts])
velx = np.zeros([nsteps,npts])
vely = np.zeros([nsteps,npts])
dpth = np.zeros([nsteps,npts])
dhdx = np.zeros([nsteps,npts])
dhdy = np.zeros([nsteps,npts])
dudx = np.zeros([nsteps,npts])
dudy = np.zeros([nsteps,npts])
dvdx = np.zeros([nsteps,npts])
dvdy = np.zeros([nsteps,npts])
tic = time.time()
for inode, node in enumerate(nodes):
nnodes = getNodesAroundNode(node,edges)
nodesX = meshx[nnodes]
nodesY = meshy[nnodes]
indx = np.where(nnodes == node)
u = getVariableAtNodes(slf,nnodes,0)
v = getVariableAtNodes(slf,nnodes,1)
h = getVariableAtNodes(slf,nnodes,3)
b = getVariableAtNodes(slf,nnodes,4)
dpth[:,inode] = b[indx,:]
elev[:,inode] = h[indx,:]
velx[:,inode] = u[indx,:]
vely[:,inode] = v[indx,:]
for istep in np.arange(nsteps):
grad,res,rnk,sng = gradientAtNode(node,nnodes,nodesX,nodesY,h[:,istep])
dhdx[istep,inode] = grad[0]
dhdy[istep,inode] = grad[1]
grad,res,rnk,sng = gradientAtNode(node,nnodes,nodesX,nodesY,u[:,istep])
dudx[istep,inode] = grad[0]
dudy[istep,inode] = grad[1]
grad,res,rnk,sng = gradientAtNode(node,nnodes,nodesX,nodesY,v[:,istep])
dvdx[istep,inode] = grad[0]
dvdy[istep,inode] = grad[1]
del node, nnodes, u, v, h
vort = dvdx - dudy
toc = time.time() - tic
print('Time Elasped = ',toc,' s')
# Close Selafin file
slf.__del__()
# Store and save results
epochStr = OTMDate2Str(default_epoch)
dyn2d = {}
dyn2d['mesh'] = mesh
dyn2d['times'] = genOTMTimeStamps(idate, times)
dyn2d['timeUnits'] = 'days since '+epochStr
dyn2d['epoch'] = epochStr
dyn2d['nsteps'] = nsteps
dyn2d['polygon'] = polygon
dyn2d['nodes'] = nodes
dyn2d['locX'] = meshx[nodes].copy()
dyn2d['locY'] = meshy[nodes].copy()
dyn2d['elev'] = elev
dyn2d['dpth'] = dpth
dyn2d['U'] = velx
dyn2d['V'] = vely
dyn2d['dhdx'] = dhdx
dyn2d['dhdy'] = dhdy
dyn2d['dudx'] = dudx
dyn2d['dudy'] = dudy
dyn2d['dvdx'] = dvdx
dyn2d['dvdy'] = dvdy
dyn2d['vort'] = vort
glbattr = generateGlobalAttributes('Telemac - Extraction of 2D dynamics at a location')
dyn2d['globalAttributes'] = glbattr
extractPath = fPath.replace('RESULTS','EXTRACT')
save_ocmw_file(extractPath,ocmwFile,dyn2d)
return
[docs]
def processLocation(filePath,runStr,suffix,Loc,depthOffset:float=0.0,modelVer=None,saveData=False):
"""
Calculate dynamic parameters for a specifc geospatial location
Parameters
----------
filePath : str
File path to Selfin file containing 2D fields.
runStr : str
Selfin file with out "_2D.slf".
suffix : str
Identifier string for location.
Loc : list, float
Coordinates of the geospatial location [x,y].
depthOffset : float, optional
Mean sea level correction used in model configuration. The default is 0.0.
modelVer : int, optional
Model version number. The default is None.
saveData : bool, optional
Toggle for activating saving data to a file. The default is False.
Returns
-------
vect : dict
Dictionary structure containing the parameters at the requested location.
"""
fName = runStr+'_2D.slf'
slf = slfObj(os.path.join(filePath,fName).replace('\\','/'))
vect = vectorFieldAtLocation(slf,Loc,d_off=depthOffset)
glbattr = generateGlobalAttributes('Telemac - Extraction of 2D dynamics at a location')
vect['globalAttributes'] = glbattr
if saveData:
if modelVer is None:
mvStr = 'v'+str(modelVer).zfill(2)
ocmwFile = runStr+'_'+suffix+'_DYN_'+mvStr+ocmwExtn
else:
ocmwFile = runStr+'_'+suffix+'_DYN'+ocmwExtn
extractPath = filePath.replace('RESULTS','EXTRACT')
save_ocmw_file(extractPath,ocmwFile,vect)
return vect