# -*- coding: utf-8 -*-
"""
Open Telemac-Mascaret output data extraction built on the HR Wallingford parser tools supplied with OTM.
Chris Old
IES, School of Engineering, University of Edinburgh
Feb 2020
"""
# ----------------------------------------------------------------------------
# IMPORTS
# ----------------------------------------------------------------------------
# Standard Python Dependencies
import os
import time
import inspect
# Non-Standard Python Dependencies
import numpy as np
# Local Module Dependencies
from ocmw.core.timeFuncs import datetime, dateNumFromDatetime
from ocmw.core.fileTools import getListOfFiles, generateGlobalAttributes
from ocmw.core.geometry import rotateVectorField
from ocmw.core.in_poly_funcs import cn_PnPoly
from ocmw.core.meshTools import isInsideTriangle, getTriangleArea
from ocmw.core.meshTools import getTriangleCentroid, getCentroidWeights
from ocmw.core.meshTools import getBarycentricWeights, weightsMatrix
from ocmw.core.meshTools import getNearestNode
from ocmw.otm_gpl.selafin import Selafin as slfObj
from .ocmw_extract import ocmwExtn, checkOCMWExtn, save_ocmw_file
# Other Dependencies
#--------------------------------------------------------------------------
# GLOBAL CONSTANTS
#--------------------------------------------------------------------------
default_epoch = np.array([1900, 1, 1, 0, 0, 0])
default_epoch_str = '1900-01-01 00:00:00'
OTM2DVarNames = ['VELOCITY U ',
'VELOCITY V ',
'CELERITY ',
'WATER DEPTH ',
'FREE SURFACE ',
'BOTTOM ',
'FROUDE NUMBER ',
'SCALAR FLOWRATE ',
'TRACER ',
'TURBULENT ENERG.',
'DISSIPATION ',
'VISCOSITY ',
'FLOWRATE ALONG X',
'FLOWRATE ALONG Y',
'SCALAR VELOCITY ',
'WIND ALONG X ',
'WIND ALONG Y ',
'AIR PRESSURE ',
'BOTTOM FRICTION ',
'RIGID BED ',
'BED THICKNESS ',
'EROSION FLUX ',
'DEPOSITION FLUX ',
'PRIVE 1 ',
'PRIVE 2 ',
'PRIVE 3 ',
'PRIVE 4 ',
'FRICTION VELOCIT',
'SOLID DISCHARGE ',
'SOLID DIS IN X ',
'SOLID DIS IN Y ',
'HIGH WATER MARK ',
'HIGH WATER TIME ',
'AIR TEMPERATURE ',
'U SURFACE ',
'V SURFACE ',
'W SURFACE ',
'M SURFACE ',
'DRIFT ALONG X ',
'DRIFT ALONG Y ',
'COURANT NUMBER ',
'HIGHEST VELOCITY',
'TIME OF HIGH VEL',
'FRICTION VEL. ',
'TAU_S ',
'1/R ',
'WALLDIST ',
'REFERENCE LEVEL ']
OCMW2DVarNames = ['U','V','C','H','SELV','BATHY','FROUDE','Q','TRAC',
'TKE','EPS','NU','QX','QY','MAG','WINDX','WINDY',
'PAIR','BFRIC','BED','BEDTHICK','EROD','DEPO',
'P2D1','P2D2','P2D3','P2D4',
'USTAR','QSED','QSEDX','QSEDY','MAXZ','TMXZ',
'TAIR','USURF','VSURF','WSURF','MSURF',
'DRFTX','DRFTY','COURNUM','MAXV','TMXV',
'USTAR','TAUS','REVRADCURV','WDIST','ZRL']
OTM3DVarNames = ['ELEVATION Z ',
'VELOCITY U ',
'VELOCITY V ',
'VELOCITY W ',
'NUX FOR VELOCITY',
'NUY FOR VELOCITY',
'NUZ FOR VELOCITY',
'TURBULENT ENERGY',
'DISSIPATION ',
'RICHARDSON NUMB ',
'RELATIVE DENSITY',
'DYNAMIC PRESSURE',
'HYDROSTATIC PRES',
'PRIVE 1 ',
'PRIVE 2 ',
'PRIVE 3 ',
'PRIVE 4 ',
'USTOKES ',
'VSTOKES ',
'WSTOKES ']
OCMW3DVarNames = ['Z','U','V','W','NUX','NUY','NUZ',
'TKE','EPS','RI','RHO','PDYN','PHYD'
'P3D1','P3D2','P3D3','P3D4',
'USTOK','VSTOK','WSTOK']
# ----------------------------------------------------------------------------
# CLASS DEFINITIONS
# ----------------------------------------------------------------------------
[docs]
class slfFile(slfObj):
"""
Class structure for reading OTM Selafin files.
Inhereting properties from the Selafin object in the OTM python scripts.
"""
def __init__(self, fileName):
slfObj.__init__(self, fileName)
self.set_mpl_tri()
self.set_kd_tree()
self.nsteps = len(self.tags['times'])
self.__slfGetTimes()
self.__slfGetMesh()
def __slfGetTimeOrigin(self):
timestamp = self.datetime
dt = datetime(timestamp[0],timestamp[1],timestamp[2],timestamp[3],timestamp[4],timestamp[5])
datenum = dateNumFromDatetime(dt)
return datenum
def __slfGetTimes(self):
t0 = self.__slfGetTimeOrigin()
times = t0 + (self.tags['times'] / (24.*60.0*60.))
self.times = times
return
def __slfGetMesh(self):
mesh = {}
mesh['nNodes'] = self.npoin2
mesh['meshx'] = self.meshx
mesh['meshy'] = self.meshy
mesh['nElems'] = self.nelem2
mesh['elems'] = self.ikle2
mesh['nlays'] = self.nplan
X = mesh['meshx'].copy()
Y = mesh['meshy'].copy()
Elems = mesh['elems'].copy()
nElems = np.shape(Elems)[0]
cWghts = np.zeros([nElems,3])
centroid = np.zeros([nElems,2])
area = np.zeros([nElems,])
for ielem in np.arange(nElems):
elem = mesh['elems'][ielem,:].copy()
p1 = [X[elem[0]],Y[elem[0]]]
p2 = [X[elem[1]],Y[elem[1]]]
p3 = [X[elem[2]],Y[elem[2]]]
centroid[ielem,:] = getTriangleCentroid(p1,p2,p3)
area[ielem] = getTriangleArea(p1,p2,p3)
cWghts[ielem,:] = getCentroidWeights(p1,p2,p3)
mesh['area'] = np.float32(area)
mesh['cLoc'] = np.float32(centroid)
mesh['cWghts'] = np.float32(cWghts)
self.mesh = mesh
return
[docs]
def properties(self):
"""
Print a list of teh class properties to std output
Returns
-------
None.
"""
props = []
for i in inspect.getmembers(self):
if not i[0].startswith('_'):
if not inspect.ismethod(i[1]):
props.append(i[0])
print(props)
return
[docs]
def methods(self):
"""
Print a list of class methods to std output
Returns
-------
None.
"""
meths = []
for i in inspect.getmembers(self):
if not i[0].startswith('_'):
if inspect.ismethod(i[1]):
meths.append(i[0])
print(meths)
return
#--------------------------------------------------------------------------
# FUNCTION DEFINITIONS
#--------------------------------------------------------------------------
[docs]
def InitialiseFile(FileName):
"""
Initialise Selafin file for reading, inlcuding setting up mesh search tree.
Parameters
----------
FileName : Str
Selafin filename including path to location
Returns
-------
SLF : Selafin obj
OTM object based on the Selafin class
"""
SLF = slfObj(FileName)
SLF.set_mpl_tri()
SLF.set_kd_tree()
return SLF
[docs]
def OTMDate2Num(timestamp):
"""
Convert OTM timestamp into a python date number
Parameters
----------
timestamp : numpy.array, int32
Array of date time fields [year,month,day,hour,minute,seccond].
Returns
-------
datenum : float
Date number.
"""
dt = datetime(timestamp[0],timestamp[1],timestamp[2],timestamp[3],timestamp[4],timestamp[5])
datenum = dateNumFromDatetime(dt)
return datenum
[docs]
def OTMDate2Str(timestamp):
"""
Convert OTM timestamp to date strings
Parameters
----------
timestamp : numpy.array, int32
Array of date time fields [year,month,day,hour,minute,seccond].
Returns
-------
datestr : str
Date string.
"""
dt = datetime(timestamp[0],timestamp[1],timestamp[2],timestamp[3],timestamp[4],timestamp[5])
datestr = dt.strftime('%Y/%m/%d %H:%M:%S')
return datestr
[docs]
def genOTMTimeStamps(idate,times,epoch=default_epoch):
"""
Convert OTM timestamp information into an array of python date numbers.
Parameters
----------
idate : numpy.array, int32
Array of date time fields [year,month,day,hour,minute,seccond].
times : numpy.array, float
Timestep offset from idate in seconds
epoch : numpy.array, int32, optional
Time epoch to use when generating timestampes.
The default is default_epoch.
Returns
-------
timeStamp : numpy.array, float
Timestamp as date number in days from epoch.
"""
epochOffset = OTMDate2Num(epoch)
startDate = OTMDate2Num(idate)-epochOffset
nstps = np.size(times)
timeStamp = np.reshape(startDate + (times / (24.*60.0*60.)),[nstps,1])
return timeStamp
[docs]
def getSlfTriangle(slfHnd,LocX,LocY):
"""
Get the triangular element and node the contain the point (LocX,LocY)
from SLF file handle
Parameters
----------
slfHnd : Selafin obj
Handle to the Selafin obj for the file being processed
LocX : float
x coordinate of the point of interest
LocY : float
y coordinate of the point of interest
Returns
-------
NNodes : list[int64]
Indices for the nodes surrounding the location
Elem : int64
Index for the mesh element corresponding to the nodes
"""
X = np.transpose(slfHnd.meshx.copy())
Y = np.transpose(slfHnd.meshy.copy())
NNode = getNearestNode(X,Y,LocX,LocY,slfHnd.nplan)
elems = np.where(slfHnd.ikle2[:,:] == NNode[0])[0]
nelems = np.size(elems)
nodes = slfHnd.ikle2[elems,:]
p = [LocX,LocY]
for Elem in np.arange(nelems):
p1 = [X[nodes[Elem,0]],Y[nodes[Elem,0]]]
p2 = [X[nodes[Elem,1]],Y[nodes[Elem,1]]]
p3 = [X[nodes[Elem,2]],Y[nodes[Elem,2]]]
if isInsideTriangle(p1,p2,p3,p):
break
else:
continue
if Elem == nelems:
NNodes = []
Elem = None
else:
NNodes = nodes[Elem,:].copy()
if len(NNode) > 1:
dn = np.diff(NNode)[0]
nl = len(NNode)
NNodes = np.resize(NNodes,3*nl)
for lay in np.arange(nl-1):
indx = np.arange(3)+(lay+1)*3
NNodes[indx] = NNodes[0:3]+(lay+1)*dn
return NNodes, elems[Elem]
[docs]
def getSlfNodesWeights(slfHnd,LocX,LocY):
"""
Get the barycentric weigthing factors to interpolate from element nodes
onto the point (LocX,LocY) from Selafin file
Parameters
----------
slfHnd : Selafin obj
Handle to the Selafin obj for the file being processed
LocX : float
x coordinate of the point of interest
LocY : float
y coordinate of the point of interest
Returns
-------
Elem : TYPE
Index of the element containing the point
NNodes : list, int64
Indices of the nodes for the element surrounding the point (LocX,LocY)
NWghts : list, float
Barycebtric weights from element nodes to point (LocX,LocY)
NLocs : list, list[float,float]
Location of nodes surroundng the point (LocX,LocY)
"""
NNodes, Elem = getSlfTriangle(slfHnd,LocX,LocY)
meshx = np.transpose(slfHnd.meshx.copy())
meshy = np.transpose(slfHnd.meshy.copy())
p1 = [meshx[NNodes[0]],meshy[NNodes[0]]]
p2 = [meshx[NNodes[1]],meshy[NNodes[1]]]
p3 = [meshx[NNodes[2]],meshy[NNodes[2]]]
p = [LocX,LocY]
NLocs = [p1,p2,p3]
NWghts = getBarycentricWeights(p1,p2,p3,p)
return Elem,NNodes,NWghts,NLocs
[docs]
def getMeshTriangle(mesh,LocX,LocY):
"""
Get the triangular element and node the contain the point (LocX,LocY)
from mesh dictionary
Parameters
----------
slfHnd : Selafin obj
Handle to the Selafin obj for the file being processed
LocX : float
x coordinate of the point of interest
LocY : float
y coordinate of the point of interest
Returns
-------
NNodes : list, int64
Indices for the nodes surrounding the location
Elem : int64
Index for the mesh element corresponding to the nodes
"""
X = np.transpose(mesh['meshx'].copy())
Y = np.transpose(mesh['meshy'].copy())
nLays = mesh['nlays']
mElems = mesh['elems']
NNode = getNearestNode(X,Y,LocX,LocY,nLays)
elems = np.where(mElems[:,:] == NNode[0])[0]
nelems = np.size(elems)
nodes = mElems[elems,:]
p = [LocX,LocY]
for Elem in np.arange(nelems):
p1 = [X[nodes[Elem,0]],Y[nodes[Elem,0]]]
p2 = [X[nodes[Elem,1]],Y[nodes[Elem,1]]]
p3 = [X[nodes[Elem,2]],Y[nodes[Elem,2]]]
if isInsideTriangle(p1,p2,p3,p):
break
else:
continue
if Elem == nelems:
NNodes = []
Elem = None
else:
NNodes = nodes[Elem,:].copy()
if len(NNode) > 1:
dn = np.diff(NNode)[0]
nl = len(NNode)
NNodes = np.resize(NNodes,3*nl)
for lay in np.arange(nl-1):
indx = np.arange(3)+(lay+1)*3
NNodes[indx] = NNodes[0:3]+(lay+1)*dn
return NNodes, elems[Elem]
[docs]
def getMeshNodesWeights(mesh,LocX,LocY):
"""
Get the barycentric weigthing factors to interpolate from element nodes
onto the point (LocX,LocY) from mesh dictionary data
Parameters
----------
slfHnd : Selafin obj
Handle to the Selafin obj for the file being processed
LocX : float
x coordinate of the point of interest
LocY : float
y coordinate of the point of interest
Returns
-------
Elem : TYPE
Index of teh element containing the point
NNodes : list, int64
Indices of the nodes for the element surrounding the point (LocX,LocY)
NWghts : list, float
Barycebtric weights from element nodes to point (LocX,LocY)
NLocs : list, list[float,float]
Location of nodes surroundng the point (LocX,LocY)
"""
NNodes, Elem = getMeshTriangle(mesh,LocX,LocY)
meshx = np.transpose(mesh['meshx'])
meshy = np.transpose(mesh['meshy'])
p1 = [meshx[NNodes[0]],meshy[NNodes[0]]]
p2 = [meshx[NNodes[1]],meshy[NNodes[1]]]
p3 = [meshx[NNodes[2]],meshy[NNodes[2]]]
p = [LocX,LocY]
NLocs = [p1,p2,p3]
NWghts = getBarycentricWeights(p1,p2,p3,p)
return Elem,NNodes,NWghts,NLocs
# ------------------------------------------------
# Functions for saving data to OCMW format
# ------------------------------------------------
[docs]
def generateOutputFileStr(filePath,runStr):
"""
Generate output path + filename string
Parameters
----------
filePath : str
Path to Selfin files
runStr : str
Prefix identifying file to be processes
i.e. filename excluding "_2D.slf" or "_3D.slf"
Returns
-------
outFName : str
Full path and filename for outputing extracted data.
"""
outPath = filePath.replace('RESULTS','EXTRACT')
outName = runStr+ocmwExtn
outFName = os.path.join(outPath,outName)
return outFName
[docs]
def save2ocmw(outfname,glbAttr,Loc,nodes,Data2D,Data3D):
"""
Create data dictionary of extracted node and location data and save to
an OCMW format file.
Parameters
----------
outfname : str
Full path and filename for data output.
glbAttr : str
Global attribute dictionary generated using "generateGlobalAttributes".
Loc : list, float
DESCRIPTION.
nodes : list, [float,float]
Element corner nodes surrounding location.
Data2D : dict
Data dictionary containg the node and interpolated model 2D fields.
Data3D : dict
Data dictionary containg the node and interpolated model 3D fields.
Returns
-------
None.
"""
store = {}
store['globalAttributes'] = glbAttr
store['dataLoc'] = Loc
store['nearestNodes'] = nodes[0]
store['nodeWeights'] = nodes[1]
store['nodeLocs'] = nodes[2]
store['timeUnits'] = 'days since '+Data2D['TS_epoch']
store['epoch'] = Data2D['TS_epoch']
store['times'] = Data2D['TS_t']
store['depth_nodes'] = Data2D['TS_NODES_d']
store['height_nodes'] = Data2D['TS_NODES_h']
store['depAvg_U_nodes'] = Data2D['TS_NODES_U']
store['depAvg_V_nodes'] = Data2D['TS_NODES_V']
store['z_nodes'] = Data3D['TS_NODES_z']
store['vel_east_nodes'] = Data3D['TS_NODES_U']
store['vel_north_nodes'] = Data3D['TS_NODES_V']
store['vel_up_nodes'] = Data3D['TS_NODES_W']
store['depth'] = Data2D['TS_d']
store['height'] = Data2D['TS_h']
store['z'] = Data3D['TS_z']
store['depAvg_U'] = Data2D['TS_U']
store['depAvg_V'] = Data2D['TS_V']
store['vel_east'] = Data3D['TS_U']
store['vel_north'] = Data3D['TS_V']
store['vel_up'] = Data3D['TS_W']
#Check output file extension
fpath,fname = checkOCMWExtn(outfname)
save_ocmw_file(fpath,fname,store)
return
'''
def saveLoc2ocmw(outfname,glbAttr,loc,crs,Data2D,Data3D):
"""
Create data dictionary of extracted location data and save to an OCMW file
Parameters
----------
outfname : str
Full path and filename for data output.
glbAttr : str
Global attribute dictionary generated using "generateGlobalAttributes".
Loc : list, [float,float]
Location coordinates.
crs : str
Coordinate reference system used for location data.
Data2D : dict
Data dictionary containg the interpolated model 2D fields.
Data3D : dict
Data dictionary containg the interpolated model 3D fields.
Returns
-------
None.
"""
store = {}
store['globalAttributes'] = glbAttr
store['dataLoc'] = loc
store['crs'] = crs
store['meshElem'] = Data2D['TS_elem']
store['timeUnits'] = 'days since '+Data2D['TS_epoch']
store['epoch'] = Data2D['TS_epoch']
store['times'] = Data2D['TS_t']
store['depth'] = Data2D['TS_d']
store['height'] = Data2D['TS_h']
store['depAvg_U'] = Data2D['TS_U']
store['depAvg_V'] = Data2D['TS_V']
if Data3D is not None:
store['z'] = Data3D['TS_z']
store['vel_east'] = Data3D['TS_U']
store['vel_north'] = Data3D['TS_V']
store['vel_up'] = Data3D['TS_W']
#Check output file extension
fpath,fname = checkOCMWExtn(outfname)
save_ocmw_file(fpath,fname,store)
return
def saveNodes2ocmw(outfname,glbAttr,nodes,Data2D,Data3D):
"""
Create data dictionary of extracted node data and save to an OCMW file
Parameters
----------
outfname : str
Full path and filename for data output.
glbAttr : str
Global attribute dictionary generated using "generateGlobalAttributes".
nodes : list, float
Indices of nodes where data are extracted.
Data2D : dict
Data dictionary containg the nodes model 2D fields.
Data3D : dict
Data dictionary containg the nodes model 3D fields.
Returns
-------
None.
"""
store = {}
store['globalAttributes'] = glbAttr
store['nodes'] = nodes
store['nodesLoc'] = Data2D['nodesLoc']
store['timeUnits'] = 'days since '+Data2D['TS_epoch']
store['epoch'] = Data2D['TS_epoch']
store['times'] = Data2D['TS_t']
store['depth'] = Data2D['TS_d']
store['height'] = Data2D['TS_h']
store['depAvg_U'] = Data2D['TS_U']
store['depAvg_V'] = Data2D['TS_V']
if Data3D is not None:
store['z'] = Data3D['TS_z']
store['vel_east'] = Data3D['TS_U']
store['vel_north'] = Data3D['TS_V']
store['vel_up'] = Data3D['TS_W']
#Check output file extension
fpath,fname = checkOCMWExtn(outfname)
save_ocmw_file(fpath,fname,store)
return
'''
[docs]
def saveLoc2ocmw(outfname,glbAttr,loc,crs,Data2D,Data3D):
"""
Create data dictionary of extracted location data and save to an OCMW file
Parameters
----------
outfname : str
Full path and filename for data output.
glbAttr : str
Global attribute dictionary generated using "generateGlobalAttributes".
Loc : list, [float,float]
Location coordinates.
crs : str
Coordinate reference system used for location data.
Data2D : dict
Data dictionary containg the interpolated model 2D fields.
Data3D : dict
Data dictionary containg the interpolated model 3D fields.
Returns
-------
None.
"""
store = {}
store['globalAttributes'] = glbAttr
store['dataLoc'] = loc
store['crs'] = crs
store['meshElem'] = Data2D['locsElem']
store['timeUnits'] = 'days since '+Data2D['TS_epoch']
store['epoch'] = Data2D['TS_epoch']
store['times'] = Data2D['TS_t']
store['depth'] = Data2D['H']
store['height'] = Data2D['SELV']
store['depAvg_U'] = Data2D['U']
store['depAvg_V'] = Data2D['V']
for akey in Data2D.keys():
if akey not in ['U','V','H','SELV','TS_t','TS_epoch','locsElem']:
store[akey] = Data2D[akey]
if Data3D is not None:
store['z'] = Data3D['Z']
store['vel_east'] = Data3D['U']
store['vel_north'] = Data3D['V']
store['vel_up'] = Data3D['W']
for akey in Data3D.keys():
if akey not in ['Z','U','V','W']:
store[akey] = Data3D[akey]
#Check output file extension
fpath,fname = checkOCMWExtn(outfname)
save_ocmw_file(fpath,fname,store)
return
[docs]
def saveNodes2ocmw(outfname,glbAttr,nodes,Data2D,Data3D):
"""
Create data dictionary of extracted node data and save to an OCMW file
Parameters
----------
outfname : str
Full path and filename for data output.
glbAttr : str
Global attribute dictionary generated using "generateGlobalAttributes".
nodes : list, float
Indices of nodes where data are extracted.
Data2D : dict
Data dictionary containg the nodes model 2D fields.
Data3D : dict
Data dictionary containg the nodes model 3D fields.
Returns
-------
None.
"""
store = {}
store['globalAttributes'] = glbAttr
store['nodes'] = nodes
store['nodesLoc'] = Data2D['nodesLoc']
store['timeUnits'] = 'days since '+Data2D['TS_epoch']
store['epoch'] = Data2D['TS_epoch']
store['times'] = Data2D['TS_t']
store['depth'] = Data2D['H']
store['height'] = Data2D['SELV']
store['depAvg_U'] = Data2D['U']
store['depAvg_V'] = Data2D['V']
for akey in Data2D.keys():
if akey not in ['U','V','H','SELV','TS_t','TS_epoch','nodesLoc']:
store[akey] = Data2D[akey]
if Data3D is not None:
store['z'] = Data3D['Z']
store['vel_east'] = Data3D['U']
store['vel_north'] = Data3D['V']
store['vel_up'] = Data3D['W']
for akey in Data3D.keys():
if akey not in ['Z','U','V','W']:
store[akey] = Data3D[akey]
#Check output file extension
fpath,fname = checkOCMWExtn(outfname)
save_ocmw_file(fpath,fname,store)
return
# ------------------------------------------------
# Functions for extracting model mesh data
# ------------------------------------------------
[docs]
def enhanceMesh(mesh):
"""
Add element cntroid locations and barycentric weighting to centroid
Parameters
----------
mesh : dict
Mesh data strcuture generated by either extractMesh or extractGeometry
Returns
-------
mesh : dict
Dictionary strcuture contining the mesh data.
"""
X = mesh['meshx'].copy()
Y = mesh['meshy'].copy()
Elems = mesh['elems'].copy()
nElems = np.shape(Elems)[0]
cWghts = np.zeros([nElems,3])
centroid = np.zeros([nElems,2])
area = np.zeros([nElems,])
for ielem in np.arange(nElems):
elem = mesh['elems'][ielem,:].copy()
p1 = [X[elem[0]],Y[elem[0]]]
p2 = [X[elem[1]],Y[elem[1]]]
p3 = [X[elem[2]],Y[elem[2]]]
centroid[ielem,:] = getTriangleCentroid(p1,p2,p3)
area[ielem] = getTriangleArea(p1,p2,p3)
cWghts[ielem,:] = getCentroidWeights(p1,p2,p3)
mesh['area'] = np.float32(area)
mesh['cLoc'] = np.float32(centroid)
mesh['cWghts'] = np.float32(cWghts)
return mesh
# ------------------------------------------------
# Functions for multiple file data extraction
# ------------------------------------------------
[docs]
def get2DElevation(vals,wgts,varNames):
"""
Get time series of surface elevation at a user defined location
Parameters
----------
vals : numpy.array, float
Time sereis of node data extracted from a Selafin 2D model output file.
wgts : list, float
Barycetnric node weightings for interpolation to location.
varNames : list, str
Available model 2D variables.
Returns
-------
elev2D : dict
Dictionary structure contiaining elevation data (nodes and location).
"""
h = getVariableData(vals,varNames,'FREE SURFACE ')
# Interpolate data onto location using Barycentric Weights
TS_h = np.sum(h*wgts,0)
# Store variables in dictionary for return
elev2D = {}
elev2D['TS_NODES_h'] = h
elev2D['TS_h'] = TS_h[:,np.newaxis]
return elev2D
[docs]
def getMultiFileLocElevation(datapath,files,Loc,
glbAttrTitle='Model tidal harmonics at XTRACK locations'):
"""
Gect a multi-file time series of model surface elevation at a user defined
location.
Used for validation of model against X-Track surface elevation data
Parameters
----------
datapath : str
Path to set of Selafin 2D model output files.
files : list, str
List of model 2D Selafin files to process.
Loc : list, [float,float]
Coordinates [x,y] of location elevation calculation.
glbAttrTitle : str, optional
Title string describing use of time series.
The default is 'Model tidal harmonics at XTRACK locations'.
Returns
-------
tsData : dict
Dictionary structure contain time series of elevation at nodes and location.
"""
nFiles = len(files)
TS_t = []
TS_NODES_h = []
TS_h = []
for ifile, f in enumerate(files):
print('Processing: '+f)
slf = InitialiseFile(os.path.join(datapath,f))
idate = slf.datetime
#idate[3] = idate[3]+idate[6]
times = slf.tags['times']
nstps = np.size(times)
varNames = slf.varnames
if ifile == 0:
TS_t = np.zeros((nFiles*nstps,))
TS_NODES_h = np.zeros((3,nFiles*nstps))
TS_h = np.zeros((nFiles*nstps,))
indices = np.arange(nstps)+ifile*nstps
TS_t[indices] = OTMDate2Num(idate) + (times / (24.*60.0*60.))
Elem,NNodes,NWghts,NLocs = getSlfNodesWeights(slf,Loc[0],Loc[1])
wgts = weightsMatrix(NWghts,nstps)
vals = slf.get_series(NNodes+1,showbar=False) # Have to offset Node indices as code assumes numbering from 1 not 0.
elev = get2DElevation(vals,wgts,varNames)
TS_NODES_h[:,indices] = elev['TS_NODES_h']
TS_h[indices] = np.squeeze(elev['TS_h'])
slf.__del__()
v,indx = np.unique(TS_t,return_index=True)
glbattr = generateGlobalAttributes(glbAttrTitle)
tsData = {}
tsData['globalAttributes'] = glbattr
tsData['location'] = Loc
tsData['meshElem'] = Elem
tsData['nearestNodes'] = NNodes
tsData['nodeWeights'] = NWghts
tsData['nodeLocs'] = NLocs
tsData['t'] = TS_t[indx,np.newaxis]
tsData['h_nodes'] = TS_NODES_h[:,indx]
tsData['h'] = TS_h[indx,np.newaxis]
return tsData
# ------------------------------------------------
# Functions for processing Location data
# ------------------------------------------------
[docs]
def getVariableData(slfData,slfVarNames,varName):
"""
Get variable time series extracted from Selafin file
Parameters
----------
slfData : numpy.array, float
Extracted time series of model fields.
slfVarNames : list, str
List of available variables.
varName : str
Name of variable to extract from Selafin file data.
Returns
-------
varData : numpy.array, float
Extracted time series for requested variable.
"""
if type(slfVarNames[0]) == 'bytes':
varIndex = slfVarNames.index(bytes(varName,'utf-8'))
else:
varIndex = slfVarNames.index(varName)
varData = np.squeeze(slfData[varIndex,:,:])
return varData
[docs]
def getLoc2DVars(vals,wgts,varNames):
"""
Extract time series for default set of 2D model fields
Parameters
----------
vals : numpy.array, float
Extracted time series of 2D field data.
wgts : numpy.array, float
Barycentric weighting factors for interpolation from nodes.
varNames : list, str
List of available model 2D fields
Returns
-------
vars2D : dict
Data dictionary containg the interpolated default model 2D fields
"""
U = getVariableData(vals,varNames,'VELOCITY U ')
V = getVariableData(vals,varNames,'VELOCITY V ')
d = getVariableData(vals,varNames,'WATER DEPTH ')
h = getVariableData(vals,varNames,'FREE SURFACE ')
# Interpolate data onto location using Barycentric Weights
TS_U = np.sum(U*wgts,0)
TS_V = np.sum(V*wgts,0)
TS_d = np.sum(d*wgts,0)
TS_h = np.sum(h*wgts,0)
# Store variables in dictionary for return
vars2D = {}
vars2D['TS_NODES_U'] = U
vars2D['TS_NODES_V'] = V
vars2D['TS_NODES_d'] = d
vars2D['TS_NODES_h'] = h
vars2D['TS_U'] = TS_U[:,np.newaxis]
vars2D['TS_V'] = TS_V[:,np.newaxis]
vars2D['TS_d'] = TS_d[:,np.newaxis]
vars2D['TS_h'] = TS_h[:,np.newaxis]
return vars2D
[docs]
def getLoc3DVars(vals,wgts,varNames):
"""
Extract time series for default set of 3D model fields
Parameters
----------
vals : numpy.array, float
Extracted time series of 3D field data.
wgts : numpy.array, float
Barycentric weighting factors for interpolation from nodes.
varNames : list, str
List of available model 3D fields
Returns
-------
vars3D : dict
Data dictionary containg the interpolated default model 3D fields
"""
z = getVariableData(vals,varNames,'ELEVATION Z ')
U = getVariableData(vals,varNames,'VELOCITY U ')
V = getVariableData(vals,varNames,'VELOCITY V ')
W = getVariableData(vals,varNames,'VELOCITY W ')
# Interpolate data onto location using Barycentric Weights
TS_z = np.sum(z*wgts,0)
TS_U = np.sum(U*wgts,0)
TS_V = np.sum(V*wgts,0)
TS_W = np.sum(W*wgts,0)
# Store variables in dictionary for return
vars3D = {}
vars3D['TS_NODES_z'] = z
vars3D['TS_NODES_U'] = U
vars3D['TS_NODES_V'] = V
vars3D['TS_NODES_W'] = W
vars3D['TS_z'] = TS_z[:,np.newaxis]
vars3D['TS_U'] = TS_U[:,np.newaxis]
vars3D['TS_V'] = TS_V[:,np.newaxis]
vars3D['TS_W'] = TS_W[:,np.newaxis]
return vars3D
[docs]
def getLocTS(slfHnd,Elem,NNodes,NWghts):
"""
Get timeseries of model data for nodes around a location and interpolate onto location
Parameters
----------
slfHnd : Selafin obj
Handle to the Selafin object for the file being processed.
Elem : int64
Element ID.
NNodes : list, int64
Nearest node indices.
NWghts : list, float
Barycentric node weights for interpolation to location.
Returns
-------
tsData : dict
Dictionary structure containing the data at the nodes and requested location.
"""
idate = slfHnd.datetime
times = slfHnd.tags['times']
nstps = np.size(times)
npts = slfHnd.npoin2;
nnds = slfHnd.npoin3
nlays = slfHnd.nplan
varNames = slfHnd.varnames
wgts = weightsMatrix(NWghts,nstps)
vals = slfHnd.get_series(NNodes+1,showbar=False) # Have to offset Node indices as code assumes numbering from 1 not 0.
# Get Variables
tsData = {}
if npts == nnds:
# 2D Variables
TS_t = genOTMTimeStamps(idate, times)
vars = getLoc2DVars(vals,wgts,varNames)
tsData['TS_elem'] = Elem
tsData['TS_t'] = TS_t
tsData['TS_epoch'] = OTMDate2Str(default_epoch)
if 'TS_NODES_U' in vars:
tsData['TS_NODES_U'] = vars['TS_NODES_U']
if 'TS_NODES_V' in vars:
tsData['TS_NODES_V'] = vars['TS_NODES_V']
if 'TS_NODES_d' in vars:
tsData['TS_NODES_d'] = vars['TS_NODES_d']
if 'TS_NODES_h' in vars:
tsData['TS_NODES_h'] = vars['TS_NODES_h']
if 'TS_U' in vars:
tsData['TS_U'] = vars['TS_U']
if 'TS_V' in vars:
tsData['TS_V'] = vars['TS_V']
if 'TS_d' in vars:
tsData['TS_d'] = vars['TS_d']
if 'TS_h' in vars:
tsData['TS_h'] = vars['TS_h']
else:
# 3D Variables
TS_NODES_z = np.zeros((3,nstps,nlays))
TS_NODES_U = np.zeros((3,nstps,nlays))
TS_NODES_V = np.zeros((3,nstps,nlays))
TS_NODES_W = np.zeros((3,nstps,nlays))
TS_z = np.zeros((nstps,nlays))
TS_U = np.zeros((nstps,nlays))
TS_V = np.zeros((nstps,nlays))
TS_W = np.zeros((nstps,nlays))
for ilay in np.arange(nlays):
lpts = np.arange(3)+ilay*3
ldat = np.squeeze(vals[:,lpts,:].copy())
vars = getLoc3DVars(ldat,wgts,varNames)
TS_NODES_z[:,:,ilay] = np.squeeze(vars['TS_NODES_z'])
TS_NODES_U[:,:,ilay] = np.squeeze(vars['TS_NODES_U'])
TS_NODES_V[:,:,ilay] = np.squeeze(vars['TS_NODES_V'])
TS_NODES_W[:,:,ilay] = np.squeeze(vars['TS_NODES_W'])
TS_z[:,ilay] = np.squeeze(vars['TS_z'])
TS_U[:,ilay] = np.squeeze(vars['TS_U'])
TS_V[:,ilay] = np.squeeze(vars['TS_V'])
TS_W[:,ilay] = np.squeeze(vars['TS_W'])
del ldat
tsData['TS_NODES_z'] = TS_NODES_z
tsData['TS_NODES_U'] = TS_NODES_U
tsData['TS_NODES_V'] = TS_NODES_V
tsData['TS_NODES_W'] = TS_NODES_W
tsData['TS_z'] = TS_z
tsData['TS_U'] = TS_U
tsData['TS_V'] = TS_V
tsData['TS_W'] = TS_W
return tsData
'''
def extract2DLocData(filePath,slf2DFile,Loc):
"""
Extract OTM model 2D field data interpolated onto a user defined location
Parameters
----------
filePath : str
Path to the data file to be processed.
slf2DFile : str
Filename of the OTM Selafin file containing the 2D model fields.
Loc : list, float
Location [x,y] where data are to be interpolated to.
Returns
-------
list
Set of output data.
Elem = index of element containing the location
NNodes = indices of the corner nodes of the element
NWghts = Barycentric weights for interpolation from nodes to location
Data2D = dictionary containing the 2D data for the nodes and location
"""
tic = time.time()
# Load 2D data file into SELAFIN structure
file2D = os.path.join(filePath,slf2DFile).replace('\\','/')
print(' Initialize 2D SELAFIN file')
slf2D = InitialiseFile(file2D)
# Get model mesh
mesh = extractMesh(file2D, enhance=False, saveMesh=False)
# Get Triangle Nodes and Barycenbtric Weights
Elem,NNodes,NWghts,NLocs = getMeshNodesWeights(mesh,Loc[0],Loc[1])
# Extract time series
print(' Extracting 2D time series')
Data2D = getLocTS(slf2D,Elem,NNodes,NWghts)
toc = time.time()
print('Total Time elapsed = ' + str(toc-tic))
return [Elem, NNodes, NWghts, NLocs, Data2D]
def extract3DLocData(filePath,fileStr,Loc):
"""
Extract OTM 2D and 3D field data interpolated onto a user defined location
Parameters
----------
filePath : str
Path to the data file to be processed.
fileStr : str
Prefix identifying file to be processes
i.e. filename excluding "_2D.slf" or "_3D.slf"
Loc : list, float
Location [x,y] where data are to be interpolated to.
Returns
-------
list
Set of output data.
Elem = index of element containing the location
NNodes = indices of the corner nodes of the element
NWghts = Barycentric weights for interpolation from nodes to location
Data2D = dictionary containing the 2D data for the nodes and location
Data3D = dictionary containing the 3D data for the nodes and location
"""
tic = time.time()
# Load 2D data file into SELAFIN structure
file2D = os.path.join(filePath,(fileStr+'_2D.slf')).replace('\\','/')
print(' Initialize 2D SELAFIN file')
slf2D = InitialiseFile(file2D)
# Get Triangle Nodes and Barycenbtric Weights
Elem,NNodes,NWghts,NLocs = getSlfNodesWeights(slf2D,Loc[0],Loc[1])
# Extract time series
print(' Extracting 2D time series')
Data2D = getLocTS(slf2D,Elem,NNodes,NWghts)
toc = time.time()
print('Time elapsed (2D) = ' + str(toc-tic))
tic1 = time.time()
# Load 3D data file into SELAFIN structure
file3D = os.path.join(filePath,(fileStr+'_3D.slf')).replace('\\','/')
fileExists = os.path.isfile(file3D)
if fileExists:
print(' Initialize 3D SELAFIN file')
slf3D = InitialiseFile(file3D)
# Clear local variables
del Elem,NNodes,NWghts,NLocs
# Get Layer Nodes and Barycenbtric Weights
Elem, NNodes, NWghts, NLocs = getSlfNodesWeights(slf3D,Loc[0],Loc[1])
# Extract time series
print(' Extracting 3D time series')
Data3D = getLocTS(slf3D,Elem,NNodes,NWghts)
else:
Data3D = None
toc = time.time()
print('Time elapsed (3D) = ' + str(toc-tic1))
print('Total time elapsed = ' + str(toc-tic))
return [NNodes, NWghts, NLocs, Data2D, Data3D]
'''
[docs]
def processLocsDataFile(filePath,runStr,locIDs,procStr,locs):
"""
Extract OTM model data for a set of user defined locations
Parameters
----------
filePath : str
Path to Selfin files
runStr : str
Prefix identifying file to be processes
i.e. filename excluding "_2D.slf" or "_3D.slf"
locIDs : list, str
Identifiers used to name extract file for each location.
modelVer : TYPE
DESCRIPTION.
procStr : TYPE
DESCRIPTION.
locs : TYPE
DESCRIPTION.
Returns
-------
None.
"""
crs = 'EPSG:32630'
file2D = filePath+'/'+runStr+'_2D.slf'
fileExists = os.path.isfile(file2D)
if fileExists:
#slf2D = InitialiseFile(file2D)
slf2D = slfFile(file2D)
else:
slf2D = None
file3D = filePath+'/'+runStr+'_3D.slf'
fileExists = os.path.isfile(file3D)
if fileExists:
#slf3D = InitialiseFile(file3D)
slf3D = slfFile(file3D)
else:
slf3D = None
for iloc,loc in enumerate(locs):
# Extract 2D data from run file
if slf2D is not None:
#[Elem, NNodes, NWghts, NLocs, Data2D] = extract2DLocData(slf2D, loc)
Data2D = extract2DLocData(slf2D, loc)
# Extract 3D data from run file
if slf3D is not None:
#[NNodes, NWghts, NLocs, Data3D] = extract3DLocData(slf3D,loc)
Data3D = extract3DLocData(slf3D,loc)
else:
Data3D = None
# Save results
runIDStr = runStr+'_'+locIDs[iloc]
outName = generateOutputFileStr(filePath,runIDStr)
outPath, outFile = os.path.split(outName)
try:
os.mkdir(outPath)
except:
pass
glbAttr = generateGlobalAttributes(procStr+' - '+os.path.split(os.path.splitext(outName)[0])[1])
saveLoc2ocmw(outName,glbAttr,loc,crs,Data2D,Data3D)
if slf2D is not None:
del slf2D
if slf3D is not None:
del slf3D
print(' Processing of '+runStr+' complete.')
print()
return
[docs]
def singleProcess(filePath,runStr,idStr,procStr,Loc):
# Extract data from run files
file2D = filePath+'/'+runStr+'_2D.slf'
slf2D = InitialiseFile(file2D)
[Elem, NNodes, NWghts, NLocs, Data2D] = extract2DLocData(slf2D,Loc)
file3D = filePath+'/'+runStr+'_3D.slf'
slf3D = InitialiseFile(file3D)
[NNodes, NWghts, NLocs, Data3D] = extract3DLocData(slf3D,Loc)
nodes = [NNodes, NWghts, NLocs]
del slf2D
del slf3D
# Save results
idStr = idStr+'_VEL'
outName = generateOutputFileStr(filePath,runStr+idStr)
outPath, outFile = os.path.split(outName)
try:
os.mkdir(outPath)
except:
pass
glbAttr = generateGlobalAttributes(procStr+' - '+os.path.split(os.path.splitext(outName)[0])[1])
save2ocmw(outName,glbAttr,Loc,nodes,Data2D,Data3D)
return
[docs]
def batchProcess(filePath,searchStr,idStr,procStr,Loc):
files = getListOfFiles(filePath,'*'+searchStr+'*_3D.slf')
for fname in files:
print('... Processing '+fname)
runStr = fname[0:-7]
singleProcess(filePath,runStr,idStr,procStr,Loc)
del runStr
return
# ------------------------------------------------
# New functions for extracting all variables
# ------------------------------------------------
[docs]
def dataDict(slfHnd,nlocs):
nstps = slfHnd.nsteps
nlays = slfHnd.nplan
data = {}
for vnam in slfHnd.varnames:
if nlays == 1:
if vnam in OTM2DVarNames:
ivar = OTM2DVarNames.index(vnam)
data[OCMW2DVarNames[ivar]] = np.squeeze(np.zeros((nlocs,nlays,nstps)))
else:
if vnam in OTM3DVarNames:
ivar = OTM3DVarNames.index(vnam)
data[OCMW3DVarNames[ivar]] = np.squeeze(np.zeros((nlocs,nlays,nstps)))
return data
[docs]
def getLocsTS(slf,locs):
locs = np.asarray(locs)
nlocs = locs.shape[0]
locsElem = np.zeros((nlocs,),'int')
data = dataDict(slf,nlocs)
nlays = slf.nplan
nstps = slf.nsteps
if nlays == 1:
idate = slf.datetime
times = slf.tags['times']
TS_t = genOTMTimeStamps(idate, times)
data['TS_t'] = TS_t
data['TS_epoch'] = OTMDate2Str(default_epoch)
for iloc in np.arange(nlocs):
locx = locs[iloc,0]
locy = locs[iloc,1]
Elem, NNodes, NWghts, NLocs = getSlfNodesWeights(slf,locx,locy)
locsElem[iloc] = Elem
wgts = weightsMatrix(NWghts,nstps)
vals = slf.get_series(NNodes+1,showbar=False)
nvars = vals.shape[0]
for ivar in np.arange(nvars):
if nlays == 1:
varindx = OTM2DVarNames.index(slf.varnames[ivar])
varstr = OCMW2DVarNames[varindx]
else:
varindx = OTM3DVarNames.index(slf.varnames[ivar])
varstr = OCMW3DVarNames[varindx]
var = vals[ivar,:]
for ilay in np.arange(nlays):
indx = np.asarray([0,1,2])+ilay*3
v = np.sum(var[indx,:]*wgts,0)
if nlays == 1:
if nlocs == 1:
data[varstr][:] = v
else:
data[varstr][iloc,:] = v
else:
if nlocs == 1:
data[varstr][ilay,:] = v
else:
data[varstr][iloc,ilay,:] = v
if nlays == 1:
data['locsElem'] = locsElem
return data
# ------------------------------------------------
# Functions for processing a linear transect
# ------------------------------------------------
[docs]
def getTransectPoints(loc01,loc02,dL):
"""
Define transect data locations based on start and end points, and spacing
Parameters
----------
loc01 : list[float,float]
Transect start point coordinates [x,y]
loc02 : list[float,float]
Transect end point coordinates [x,y]
dL : float
Transect sample spacing
Returns
-------
Locs : numpy.array, [float,float]
Coordinates of transect data points.
Theta : float
Direction of transect - mathematical convention.
"""
Rx = loc02[0]-loc01[0]
Ry = loc02[1]-loc01[1]
R = np.sqrt(Rx**2. + Ry**2.)
Theta = np.rad2deg(np.arctan2(Rx,Ry)) # opposite to normal as assume x-axis perpendicular to transect
npts = int(np.ceil(R/dL))
R2 = npts*dL
Rx = (R2/R)*Rx
Ry = (R2/R)*Ry
dX = Rx/npts
dY = Ry/npts
X0 = loc01[0]
Y0 = loc01[1]
Locs = np.zeros((npts,2))
for ipt in np.arange(npts):
Locs[ipt,0] = X0 + ipt*dX
Locs[ipt,1] = Y0 + ipt*dY
return Locs,Theta
'''
def extract2DTransect(filePath,runStr,tranStr,loc01,loc02,dL,crs):
"""
Extract a transect of 2-D data from Selafin file
Parameters
----------
filePath : str
Path to Selfin files
runStr : str
Prefix identifying file to be processes
i.e. filename excluding "_2D.slf"
tranStr : str
String containing transect description for setting the Global Attributes title.
loc01 : list[float]
Transect start point location [x,y]
loc02 : list[float]
Transect end point location [x,y]
dL : float
Transect sample separation in m
Returns
-------
Trans_Data : dict
Data dictionary containing transect data
"""
# Get transect points and transect orientation
Locs,Theta = getTransectPoints(loc01,loc02,dL)
npts = Locs.shape[0]
# Load 3D data file into SELAFIN structure
slfFile = os.path.join(filePath,(runStr+'_2D.slf'))
slfHnd = InitialiseFile(slfFile)
# Get array dimensions
nstps = np.size(slfHnd.tags['times'])
# Initialize storage space
Trans_d = np.zeros((npts,nstps))
Trans_h = np.zeros((npts,nstps))
Trans_U = np.zeros((npts,nstps))
Trans_V = np.zeros((npts,nstps))
Trans_Vpp = np.zeros((npts,nstps))
Trans_Vtr = np.zeros((npts,nstps))
# Extract time series
for ipt in range(npts):
Loc = np.squeeze(Locs[ipt,:])
#[Elems, NNodes, NWghts, NLocs, Data2D] = extract2DLocData(filePath,slfFile,Loc)
[Elems, NNodes, NWghts, NLocs, Data2D] = extract2DLocData(slfHnd,Loc)
U = Data2D['TS_U']
V = Data2D['TS_V']
W = np.zeros(U.shape)
[Vpp,Vtr,Up] = rotateVectorField(U,V,W,Theta)
Trans_d[ipt,:] = np.squeeze(Data2D['TS_d'])
Trans_h[ipt,:] = np.squeeze(Data2D['TS_h'])
Trans_U[ipt,:] = np.squeeze(Data2D['TS_U'])
Trans_V[ipt,:] = np.squeeze(Data2D['TS_V'])
Trans_Vtr[ipt,:] = np.squeeze(Vtr)
Trans_Vpp[ipt,:] = np.squeeze(Vpp)
# Generate global attribute record
glbAttr = generateGlobalAttributes('Transect Data - '+tranStr+' from '+slfFile)
# Create data dictionary
Trans_Data = {}
Trans_Data['globalAttributes'] = glbAttr
Trans_Data['transStart'] = loc01
Trans_Data['transEnd'] = loc02
Trans_Data['transSpcing'] = dL
Trans_Data['crs'] = crs
Trans_Data['only2D'] = True
Trans_Data['timeUnits'] = 'days since '+Data2D['TS_epoch']
Trans_Data['epoch'] = Data2D['TS_epoch']
Trans_Data['times'] = Data2D['TS_t']
Trans_Data['depth'] = Trans_d
Trans_Data['height'] = Trans_h
Trans_Data['dataLoc'] = Locs
Trans_Data['depAvg_U'] = Trans_U
Trans_Data['depAvg_V'] = Trans_V
Trans_Data['depAvg_Vpp'] = Trans_Vpp
Trans_Data['depAvg_Vtr'] = Trans_Vtr
# Close slfFile
slfHnd.__del__()
# Return results
return Trans_Data
def extract3DTransect(filePath,runStr,tranStr,loc01,loc02,dL,crs):
"""
Extract a transect of 3-D data from Selafin file
Parameters
----------
filePath : str
Path to Selfin files
runStr : str
Prefix identifying file to be processes
i.e. filename excluding "_3D.slf"
tranStr : str
String containing transect description for setting the Global Attributes title.
loc01 : list[float]
Transect start point location [x,y]
loc02 : list[float]
Transect end point location [x,y]
dL : float
Transect sample separation in m
Returns
-------
Trans_Data : dict
Data dictionary containing transect data
"""
# Get time stamps from 2D file
slfFile = os.path.join(filePath,(runStr+'_2D.slf'))
slfHnd = InitialiseFile(slfFile)
idate = slfHnd.datetime
times = slfHnd.tags['times']
Data2D = {}
Data2D['TS_t'] = genOTMTimeStamps(idate, times)
Data2D['TS_epoch'] = OTMDate2Str(default_epoch)
slfHnd.__del__()
# Get transect points and transect orientation
Locs,Theta = getTransectPoints(loc01,loc02,dL)
npts = Locs.shape[0]
# Load 3D data file into SELAFIN structure
slfFile = os.path.join(filePath,(runStr+'_3D.slf'))
slfHnd = InitialiseFile(slfFile)
# Get array dimensions
nstps = np.size(slfHnd.tags['times'])
nlays = slfHnd.nplan
# Initialize storage space
Trans_z = np.zeros((npts,nstps,nlays))
Trans_U = np.zeros((npts,nstps,nlays))
Trans_V = np.zeros((npts,nstps,nlays))
Trans_W = np.zeros((npts,nstps,nlays))
Trans_Vpp = np.zeros((npts,nstps,nlays))
Trans_Vtr = np.zeros((npts,nstps,nlays))
# Extract time series
for ipt in range(npts):
Loc = np.squeeze(Locs[ipt,:])
[NNodes, NWghts, NLocs, Data3D] = extract3DLocData(slfHnd,Loc)
U = Data3D['TS_U']
V = Data3D['TS_V']
W = Data3D['TS_W']
[Vpp,Vtr,Up] = rotateVectorField(U,V,W,Theta)
Trans_z[ipt,:,:] = Data3D['TS_z']
Trans_U[ipt,:,:] = Data3D['TS_U']
Trans_V[ipt,:,:] = Data3D['TS_V']
Trans_W[ipt,:,:] = Data3D['TS_W']
Trans_Vtr[ipt,:,:] = Vtr
Trans_Vpp[ipt,:,:] = Vpp
# Generate global attribute record
glbAttr = generateGlobalAttributes('Transect Data - '+tranStr+' from '+runStr)
# Create data dictionary
Trans_Data = {}
Trans_Data['globalAttributes'] = glbAttr
Trans_Data['transStart'] = loc01
Trans_Data['transEnd'] = loc02
Trans_Data['transSpcing'] = dL
Trans_Data['crs'] = crs
Trans_Data['only2D'] = False
Trans_Data['timeUnits'] = 'days since '+Data2D['TS_epoch']
Trans_Data['epoch'] = Data2D['TS_epoch']
Trans_Data['times'] = Data2D['TS_t']
Trans_Data['z'] = Trans_z
Trans_Data['dataLoc'] = Locs
Trans_Data['vel_east'] = Trans_U
Trans_Data['vel_north'] = Trans_V
Trans_Data['vel_up'] = Trans_W
Trans_Data['vel_pp'] = Trans_Vpp
Trans_Data['vel_tr'] = Trans_Vtr
# Close slfFile
slfHnd.__del__()
# Return results
return Trans_Data
def extractTransect(filePath,runStr,tranStr,loc01,loc02,dL,crs,only2D=False):
"""
Extract a transect of 2-D or 3-D data from Selafin file
Parameters
----------
filePath : str
Path to Selfin files
runStr : str
Prefix identifying file to be processes
i.e. filename excluding "_2D.slf" or "_3D.slf"
tranStr : str
String containing transect description for setting the Global Attributes title.
loc01 : list[float]
Transect start point location [x,y]
loc02 : list[float]
Transect end point location [x,y]
dL : float
Transect sample separation in m
only2D : bool, optional
Toggle to extract only 2D data. The default is False.
Returns
-------
Trans_Data : dict
Data dictionary containing transect data
"""
if only2D:
slfFile = os.path.join(filePath,(runStr+'_2D.slf'))
fileExists = os.path.isfile(slfFile)
if fileExists:
Trans_Data = extract2DTransect(filePath, runStr, tranStr, loc01, loc02, dL, crs)
else:
print('WARNING: Data for '+runStr+' not available...')
print(' Check data path '+filePath+' exists.')
print(' Check data file(s) '+runStr+'_*D.slf exist.')
Trans_Data = None
else:
# Check if 3D data are available
slfFile = os.path.join(filePath,(runStr+'_3D.slf'))
fileExists = os.path.isfile(slfFile)
if fileExists:
Trans_Data = extract3DTransect(filePath, runStr, tranStr, loc01, loc02, dL, crs)
else:
slfFile = os.path.join(filePath,(runStr+'_2D.slf'))
fileExists = os.path.isfile(slfFile)
if fileExists:
Trans_Data = extract2DTransect(filePath, runStr, tranStr, loc01, loc02, dL, crs)
else:
print('WARNING: Data for '+runStr+' not available...')
print(' Check data path '+filePath+' exists.')
print(' Check data file(s) '+runStr+'_*D.slf exist.')
Trans_Data = None
# Return results
return Trans_Data
'''
# ------------------------------------------------
# Functions for processing Node data
# ------------------------------------------------
[docs]
def getVarTimeSeries(slfHnd,nodes,varStr):
"""
Extract timeseries of data for a particular model variable a set of model nodes.
"""
nnodes = len(nodes)
nstps = slfHnd.tags['times'].size
npts = slfHnd.npoin2
nlays = slfHnd.nplan
varNames = slfHnd.varnames
varStr = varStr.ljust(16,' ')
if varStr in varNames:
varTS = np.zeros((nnodes,nstps,nlays))
for ilay in np.arange(nlays):
nodeset = nodes+npts*ilay
vals = slfHnd.get_series(nodeset+1,showbar=False)
var = getVariableData(vals,varNames,varStr)
varTS[:,:,ilay] = var
'''
if nlays == 1:
vals = slfHnd.get_series(nodes+1,showbar=False)
var = getVariableData(vals,varNames,varStr)
varTS[:,:,0] = var
else:
indx = varNames.index(varStr)
for ts in np.arange(nstps):
vals = slfHnd.get_variables_at(ts,[indx]).reshape((npts,nlays))
varTS[:,ts,:] = vals[nodes,:]
'''
else:
varTS = None
return np.squeeze(varTS)
'''
def getNodesTS(slfHnd,nodes):
"""
Extract timeseries of data for a specific model node.
"""
idate = slfHnd.datetime
times = slfHnd.tags['times']
npts = slfHnd.npoin2;
nnds = slfHnd.npoin3
nodesLoc = np.asarray([slfHnd.meshx[nodes],slfHnd.meshy[nodes]]).transpose()
# Get Variables
tsData = {}
if npts == nnds:
# 2D Variables
TS_t = genOTMTimeStamps(idate, times)
print(' Extracting 2D time series: VELOCITY U')
TS_U = getVarTimeSeries(slfHnd,nodes,'VELOCITY U')
print(' Extracting 2D time series: VELOCITY V')
TS_V = getVarTimeSeries(slfHnd,nodes,'VELOCITY V')
print(' Extracting 2D time series: WATER DEPTH')
TS_d = getVarTimeSeries(slfHnd,nodes,'WATER DEPTH')
print(' Extracting 2D time series: FREE SURFACE')
TS_h = getVarTimeSeries(slfHnd,nodes,'FREE SURFACE')
tsData['TS_t'] = TS_t
tsData['TS_epoch'] = OTMDate2Str(default_epoch)
tsData['nodesLoc'] = nodesLoc
if TS_U is not None:
tsData['TS_U'] = TS_U
if TS_V is not None:
tsData['TS_V'] = TS_V
if TS_d is not None:
tsData['TS_d'] = TS_d
if TS_h is not None:
tsData['TS_h'] = TS_h
else:
# 3D Variables
print(' Extracting 3D time series: ELEVATION Z')
TS_z = getVarTimeSeries(slfHnd,nodes,'ELEVATION Z')
print(' Extracting 3D time series: VELOCITY U')
TS_U = getVarTimeSeries(slfHnd,nodes,'VELOCITY U')
print(' Extracting 3D time series: VELOCITY V')
TS_V = getVarTimeSeries(slfHnd,nodes,'VELOCITY V')
print(' Extracting 3D time series: VELOCITY W')
TS_W = getVarTimeSeries(slfHnd,nodes,'VELOCITY W')
tsData['TS_z'] = TS_z
tsData['TS_U'] = TS_U
tsData['TS_V'] = TS_V
tsData['TS_W'] = TS_W
return tsData
def extractNodeTSData(filePath,fileStr,nodes):
"""
Extract timeseries of data for a set of model nodes.
"""
tic = time.time()
# Load 2D data file into SELAFIN structure
file2D = os.path.join(filePath,(fileStr+'_2D.slf'))
print(' Initialize 2D SELAFIN file')
slf2D = InitialiseFile(file2D)
# Extract time series
Data2D = getNodesTS(slf2D,nodes)
# Close slfFile
slf2D.__del__()
toc = time.time()
print('Time elapsed (2D) = ' + str(toc-tic))
tic1 = time.time()
# Load 3D data file into SELAFIN structure
file3D = os.path.join(filePath,(fileStr+'_3D.slf'))
fileExists = os.path.isfile(file3D)
if fileExists:
print(' Initialize 3D SELAFIN file')
slf3D = InitialiseFile(file3D)
print(' Extracting 3D time series')
Data3D = getNodesTS(slf3D,nodes)
# Close slfFile
slf3D.__del__()
else:
Data3D = None
toc = time.time()
print('Time elapsed (3D) = ' + str(toc-tic1))
print('Total time elapsed = ' + str(toc-tic))
return [Data2D, Data3D]
'''
[docs]
def getNodesTS(slfHnd,nodes):
"""
Extract timeseries of data for a specific model node.
"""
nodesLoc = np.asarray([slfHnd.meshx[nodes],slfHnd.meshy[nodes]]).transpose()
nlocs = np.asarray(nodes).shape[0]
tsData = dataDict(slfHnd,nlocs)
nnodes = slfHnd.npoin2
nlays = slfHnd.nplan
nstps = slfHnd.nsteps
nvars = len(slfHnd.varnames)
if nlays == 1:
idate = slfHnd.datetime
times = slfHnd.tags['times']
TS_t = genOTMTimeStamps(idate, times)
tsData['TS_t'] = TS_t
tsData['TS_epoch'] = OTMDate2Str(default_epoch)
tsData['nodesLoc'] = nodesLoc
for ilay in np.arange(nlays):
dnodes = nodes + ilay*nnodes
vals = slfHnd.get_series(dnodes+1,showbar=False)
for ivar in np.arange(nvars):
if nlays == 1:
varindx = OTM2DVarNames.index(slfHnd.varnames[ivar])
varstr = OCMW2DVarNames[varindx]
else:
varindx = OTM3DVarNames.index(slfHnd.varnames[ivar])
varstr = OCMW3DVarNames[varindx]
var = vals[ivar,:]
if nlays == 1:
tsData[varstr][:,:] = var
else:
tsData[varstr][:,ilay,:] = var
return tsData
[docs]
def processNodeDataFile(filePath,runStr,idStr,procStr,nodes):
"""
Extract data from a set of nodes based on user supplied run file.
"""
# Extract data from run files
[Data2D, Data3D] = extractNodeTSData(filePath,runStr,nodes)
# Save results
outName = generateOutputFileStr(filePath,runStr+idStr)
outPath, outFile = os.path.split(outName)
try:
os.mkdir(outPath)
except:
pass
glbAttr = generateGlobalAttributes(procStr+' - '+os.path.split(os.path.splitext(outName)[0])[1])
saveNodes2ocmw(outName,glbAttr,nodes,Data2D,Data3D)
print(' Processing of '+runStr+' complete.')
print()
return
# ------------------------------------------------
# Functions for processing Polygon data
# ------------------------------------------------
[docs]
def getMeshNodeLocs(filePath,runStr,setStr=None):
"""
Extract mesh node indices and node coordinates from a 2D Selafine file
Parameters
----------
filePath : TYPE
DESCRIPTION.
runStr : TYPE
DESCRIPTION.
setStr : TYPE, optional
DESCRIPTION. The default is None.
Returns
-------
nnodes : TYPE
DESCRIPTION.
meshx : TYPE
DESCRIPTION.
meshy : TYPE
DESCRIPTION.
"""
if setStr is None:
slfFile = os.path.join(filePath,(runStr+'_2D.slf'))
else:
slfFile = os.path.join(filePath,(runStr+'_'+setStr+'_2D.slf'))
slf = InitialiseFile(slfFile)
nnodes = slf.npoin2
meshx = slf.meshx
meshy = slf.meshy
# Close slfFile
slf.__del__()
return nnodes, meshx, meshy
[docs]
def processPolygonDataFile(filePath,runStr,idStr,procStr,polygon):
"""
Extract data from nodes within a polygon based on user supplied run file.
"""
nnodes, meshx, meshy = getMeshNodeLocs(filePath, runStr)
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]
[Data2D, Data3D] = extractNodeTSData(filePath,runStr,nodes)
# Save results
outName = generateOutputFileStr(filePath,runStr+idStr)
outPath, outFile = os.path.split(outName)
try:
os.mkdir(outPath)
except:
pass
glbAttr = generateGlobalAttributes(procStr+' - '+os.path.split(os.path.splitext(outName)[0])[1])
saveNodes2ocmw(outName,glbAttr,nodes,Data2D,Data3D)
print(' Processing of '+runStr+' complete.')
print()
return
# ------------------------------------------------
# Functions for parsing ASCII run files
# ------------------------------------------------
[docs]
def parseLocsRunFile(fileName):
"""
Parser to read run parameters from a configuration file for an extract at
user defined locations
Parameters
----------
fileName : str
Full path and file name for run configuration file.
Returns
-------
params : dict
Dicttionary containing the run paramters.
"""
params = {}
nData = 0
with open(fileName) as fin:
for line in fin:
astr = line.strip()
if '=' in astr:
if 'Data Path' in astr:
val = astr.split(' = ')[1]
params['dataPath'] = val
if 'Search String' in astr:
val = astr.split(' = ')[1]
params['srchStr'] = val
if 'Start Date' in astr:
val = astr.split(' = ')[1]
params['startDate'] = val
if 'Number of Days' in astr:
val = astr.split(' = ')[1]
params['nDays'] = int(val)
if 'Locations CRS' in astr:
val = astr.split(' = ')[1]
params['CRS'] = val
if 'Number of Locations' in astr:
val = astr.split(' = ')[1]
nLocs = int(val)
params['nLocs'] = nLocs
xloc = np.zeros([nLocs,])
yloc = np.zeros([nLocs,])
locid = []
else:
vals = astr.split()
xloc[nData] = float(vals[0])
yloc[nData] = float(vals[1])
locid.append(vals[2])
nData += 1
params['xLoc'] = xloc
params['yLoc'] = yloc
params['locID'] = locid
return params
[docs]
def parseNodesRunFile(fileName):
"""
Parser to read run parameters from a configuration file for an extract at
model nodes
Parameters
----------
fileName : str
Full path and file name for run configuration file.
Returns
-------
params : dict
Dicttionary containing the run paramters.
"""
params = {}
nData = 0
with open(fileName) as fin:
for line in fin:
astr = line.strip()
if '=' in astr:
if 'Data Path' in astr:
val = astr.split(' = ')[1]
params['dataPath'] = val
if 'Search String' in astr:
val = astr.split(' = ')[1]
params['srchStr'] = val
if 'Start Date' in astr:
val = astr.split(' = ')[1]
params['startDate'] = val
if 'Number of Days' in astr:
val = astr.split(' = ')[1]
params['nDays'] = int(val)
if 'Number of Nodes' in astr:
val = astr.split(' = ')[1]
nNodes = int(val)
params['nNodes'] = nNodes
nodes = np.zeros([nNodes,],dtype=int)
else:
vals = astr.split()
nodes[nData] = int(vals[0])
nData += 1
params['nodes'] = nodes
return params
[docs]
def parseTransectRunFile(fileName):
"""
Parser to read run parameters from a configuration file for extract along
a trnasect
Parameters
----------
fileName : str
Full path and file name for run configuration file.
Returns
-------
params : dict
Dicttionary containing the run paramters.
"""
params = {}
with open(fileName) as fin:
for line in fin:
astr = line.strip()
if '=' in astr:
if 'Data Path' in astr:
val = astr.split(' = ')[1]
params['dataPath'] = val
if 'Output Path' in astr:
val = astr.split(' = ')[1]
params['outPath'] = val
if 'Run String' in astr:
val = astr.split(' = ')[1]
params['runStr'] = val
if 'Output File' in astr:
val = astr.split(' = ')[1]
params['outFile'] = val
if 'Title String' in astr:
val = astr.split(' = ')[1]
params['titleStr'] = val
if 'Position CRS' in astr:
val = astr.split(' = ')[1]
params['CRS'] = val
if 'Start Point' in astr:
val = astr.split(' = ')[1]
vals = np.asarray(val.strip('[').strip(']').split(','))
params['loc01'] = [float(vals[0]),float(vals[1])]
if 'End Point' in astr:
val = astr.split(' = ')[1]
vals = np.asarray(val.strip('[').strip(']').split(','))
params['loc02'] = [float(vals[0]),float(vals[1])]
if 'Transect Sampling' in astr:
val = astr.split(' = ')[1]
params['dL'] = float(val)
return params
[docs]
def parsePolygonRunFile(fileName):
"""
Parser to read run parameters from a configuration file for an extract of
nodes within a polygon
Parameters
----------
fileName : str
Full path and file name for run configuration file.
Returns
-------
params : dict
Dicttionary containing the run paramters.
"""
params = {}
nData = 0
with open(fileName) as fin:
for line in fin:
astr = line.strip()
if '=' in astr:
if 'Data Path' in astr:
val = astr.split(' = ')[1]
params['dataPath'] = val
if 'Search String' in astr:
val = astr.split(' = ')[1]
params['srchStr'] = val
if 'Start Date' in astr:
val = astr.split(' = ')[1]
params['startDate'] = val
if 'Number of Days' in astr:
val = astr.split(' = ')[1]
params['nDays'] = int(val)
if 'Locations CRS' in astr:
val = astr.split(' = ')[1]
params['CRS'] = val
if 'Number of Vertices' in astr:
val = astr.split(' = ')[1]
nVert = int(val)
params['nVert'] = nVert
vertX = np.zeros([nVert,])
vertY = np.zeros([nVert,])
else:
vals = astr.split()
vertX[nData] = float(vals[0])
vertY[nData] = float(vals[1])
nData += 1
params['vertX'] = vertX
params['vertY'] = vertY
return params