Source code for fastwater.telemac.extraction
# -*- coding: utf-8 -*-
"""
OpenTelemac output data extraction using HR Wallingford parser tools.
"""
# ----------------------------------------------------------------------------
# IMPORTS
# ----------------------------------------------------------------------------
# Standard Python Dependencies
import os
import time
import numpy as np
import scipy.io as sio
# Non-Standard Python Dependencies
# Local Module Dependencies
from fastwater.general.timeFuncs import datetime, date2num, num2date, dateStr, dateNumFromStr
from fastwater.general.timeFuncs import stdTimeUnits, stdTimeFmt, stdCalendar
from fastwater.general.fileTools import getListOfFiles, loadMatFile, generateGlobalAttributes
from fastwater.general.geometry import rotateVectorField
from fastwater.general.meshTools import isInsideTriangle, getTriangleArea
from fastwater.general.meshTools import getTriangleCentroid, getCentroidWeights
from fastwater.general.meshTools import getBarycentricWeights, weightsMatrix
from fastwater.general.meshTools import getElemAroundLoc, getNearestNode, getNodesWeights
from fastwater.otm_gpl.selafin import Selafin as slfObj
# Other Dependencies
#--------------------------------------------------------------------------
# GLOBAL CONSTANTS
#--------------------------------------------------------------------------
default_epoch = np.array([1900, 1, 1, 0, 0, 0])
default_epoch_str = '1900-01-01 00:00:00'
# ----------------------------------------------------------------------------
# CLASS DEFINITIONS
# ----------------------------------------------------------------------------
#--------------------------------------------------------------------------
# FUNCTION DEFINITIONS
#--------------------------------------------------------------------------
[docs]def InitialiseFile(FileName):
SLF = slfObj(FileName)
SLF.set_mpl_tri()
SLF.set_kd_tree()
return SLF
[docs]def Date2Num(timestamp):
dt = datetime(timestamp[0],timestamp[1],timestamp[2],timestamp[3],timestamp[4],timestamp[5])
datenum = date2num(dt,stdTimeUnits,stdCalendar)
return datenum
[docs]def Date2Str(timestamp):
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 genTimeStamps(idate,times,nstps,epoch=default_epoch):
epochOffset = Date2Num(epoch)
startDate = Date2Num(idate)-epochOffset
timeStamp = np.reshape(startDate + (times / (24.*60.0*60.)),[nstps,1])
return timeStamp
[docs]def getTriangle(slfHnd,LocX,LocY):
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, Elem
[docs]def getSlfNodesWeights(slfHnd,LocX,LocY):
NNodes, Elem = getTriangle(slfHnd,LocX,LocY)
p1 = [slfHnd.meshx[NNodes[0]],slfHnd.meshy[NNodes[0]]]
p2 = [slfHnd.meshx[NNodes[1]],slfHnd.meshy[NNodes[1]]]
p3 = [slfHnd.meshx[NNodes[2]],slfHnd.meshy[NNodes[2]]]
p = [LocX,LocY]
NLocs = [p1,p2,p3]
NWghts = getBarycentricWeights(p1,p2,p3,p)
return Elem,NNodes,NWghts,NLocs
[docs]def getVariableData(slfData,slfNames,varName):
if type(slfNames[0]) == 'bytes':
varIndex = slfNames.index(bytes(varName,'utf-8'))
else:
varIndex = slfNames.index(varName)
varData = np.squeeze(slfData[varIndex,:,:])
return varData
[docs]def get2DElevation(vals,wgts,varNames):
h = getVariableData(vals,varNames,'FREE SURFACE ')
# Interpolate data onto location using Barycentric Weights
TS_INTRP_h = np.sum(h*wgts,0)
# Store variables in dictionary for return
elev2D = {}
elev2D['TS_h'] = h
elev2D['TS_INTRP_h'] = TS_INTRP_h[:,np.newaxis]
return elev2D
[docs]def getMultiFileSeries(datapath,files,Loc,showbar=False):
TS_t = []
TS_h = []
TS_INTRP_h = []
for f in files:
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
TS_t = np.concatenate((TS_t , (Date2Num(idate) + (times / (24.*60.0*60.)))),axis=0)
Elem,NNodes,NWghts,NLocs = getSlfNodesWeights(slf,Loc[0],Loc[1])
wgts = weightsMatrix(NWghts,nstps)
vals = slf.get_series(NNodes+1,showbar=showbar) # Have to offset Node indices as code assumes numbering from 1 not 0.
elev = get2DElevation(vals,wgts,varNames)
if TS_h == []:
TS_h = elev['TS_h'].tolist()
else:
TS_h = np.concatenate((TS_h,elev['TS_h'].tolist()),axis=1)
TS_INTRP_h = np.concatenate((TS_INTRP_h,np.squeeze(elev['TS_INTRP_h'].tolist())),axis=0)
v,indx = np.unique(TS_t,return_index=True)
glbattr = generateGlobalAttributes('Model tidal harmonics at XTRACK locations')
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'] = TS_h[:,indx]
tsData['h_loc'] = TS_INTRP_h[indx,np.newaxis]
return tsData
[docs]def get2DVars(vals,wgts,varNames):
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_INTRP_U = np.sum(U*wgts,0)
TS_INTRP_V = np.sum(V*wgts,0)
TS_INTRP_d = np.sum(d*wgts,0)
TS_INTRP_h = np.sum(h*wgts,0)
# Store variables in dictionary for return
vars2D = {}
vars2D['TS_U'] = U
vars2D['TS_V'] = V
vars2D['TS_d'] = d
vars2D['TS_h'] = h
vars2D['TS_INTRP_U'] = TS_INTRP_U[:,np.newaxis]
vars2D['TS_INTRP_V'] = TS_INTRP_V[:,np.newaxis]
vars2D['TS_INTRP_d'] = TS_INTRP_d[:,np.newaxis]
vars2D['TS_INTRP_h'] = TS_INTRP_h[:,np.newaxis]
return vars2D
[docs]def get3DVars(vals,wgts,varNames):
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_INTRP_z = np.sum(z*wgts,0)
TS_INTRP_U = np.sum(U*wgts,0)
TS_INTRP_V = np.sum(V*wgts,0)
TS_INTRP_W = np.sum(W*wgts,0)
# Store variables in dictionary for return
vars3D = {}
vars3D['TS_z'] = z
vars3D['TS_U'] = U
vars3D['TS_V'] = V
vars3D['TS_W'] = W
vars3D['TS_INTRP_z'] = TS_INTRP_z[:,np.newaxis]
vars3D['TS_INTRP_U'] = TS_INTRP_U[:,np.newaxis]
vars3D['TS_INTRP_V'] = TS_INTRP_V[:,np.newaxis]
vars3D['TS_INTRP_W'] = TS_INTRP_W[:,np.newaxis]
return vars3D
[docs]def extractTimeSeries(slfHnd,NNodes,NWghts,loc,showbar=False):
idate = slfHnd.datetime
times = slfHnd.tags['times']
nstps = np.size(times)
npts = slfHnd.npoin2;
nnds = slfHnd.npoin3
nlays = slfHnd.nplan
#nvars = slfHnd.nbv1
varNames = slfHnd.varnames
wgts = weightsMatrix(NWghts,nstps)
vals = slfHnd.get_series(NNodes+1,showbar=showbar) # Have to offset Node indices as code assumes numbering from 1 not 0.
# Get Variables
tsData = {}
if npts == nnds:
# 2D Variables
#TS_t = np.reshape(Date2Num(idate) + (times / (24.*60.0*60.)),[nstps,1])
TS_t = genTimeStamps(idate, times, nstps)
vars = get2DVars(vals,wgts,varNames)
tsData['TS_t'] = TS_t
tsData['TS_epoch'] = Date2Str(default_epoch)
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']
if 'TS_INTRP_U' in vars:
tsData['TS_INTRP_U'] = vars['TS_INTRP_U']
if 'TS_INTRP_V' in vars:
tsData['TS_INTRP_V'] = vars['TS_INTRP_V']
if 'TS_INTRP_d' in vars:
tsData['TS_INTRP_d'] = vars['TS_INTRP_d']
if 'TS_INTRP_h' in vars:
tsData['TS_INTRP_h'] = vars['TS_INTRP_h']
else:
# 3D Variables
TS_z = np.zeros((3,nstps,nlays))
TS_U = np.zeros((3,nstps,nlays))
TS_V = np.zeros((3,nstps,nlays))
TS_W = np.zeros((3,nstps,nlays))
TS_INTRP_z = np.zeros((nstps,nlays))
TS_INTRP_U = np.zeros((nstps,nlays))
TS_INTRP_V = np.zeros((nstps,nlays))
TS_INTRP_W = np.zeros((nstps,nlays))
for ilay in np.arange(nlays):
lpts = np.arange(3)+ilay*3
ldat = np.squeeze(vals[:,lpts,:].copy())
vars = get3DVars(ldat,wgts,varNames)
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'])
TS_INTRP_z[:,ilay] = np.squeeze(vars['TS_INTRP_z'])
TS_INTRP_U[:,ilay] = np.squeeze(vars['TS_INTRP_U'])
TS_INTRP_V[:,ilay] = np.squeeze(vars['TS_INTRP_V'])
TS_INTRP_W[:,ilay] = np.squeeze(vars['TS_INTRP_W'])
del ldat
tsData['TS_z'] = TS_z
tsData['TS_U'] = TS_U
tsData['TS_V'] = TS_V
tsData['TS_W'] = TS_W
tsData['TS_INTRP_z'] = TS_INTRP_z
tsData['TS_INTRP_U'] = TS_INTRP_U
tsData['TS_INTRP_V'] = TS_INTRP_V
tsData['TS_INTRP_W'] = TS_INTRP_W
return tsData
[docs]def extract2DRunData(filePath,fileStr,Loc,showbar=False):
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)
# 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 = extractTimeSeries(slf2D,NNodes,NWghts,Loc,showbar)
toc = time.time()
print('Total Time elapsed = ' + str(toc-tic))
return [Elem, NNodes, NWghts, NLocs, Data2D]
[docs]def extractRunData(filePath,fileStr,Loc,showbar=False):
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)
# 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 = extractTimeSeries(slf2D,NNodes,NWghts,Loc,showbar)
# Clear local variables
del NNodes,NWghts,NLocs
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'))
slf3D = InitialiseFile(file3D)
# Get Triangle Nodes and Barycenbtric Weights
Elem,NNodes,NWghts,NLocs = getSlfNodesWeights(slf3D,Loc[0],Loc[1])
# Extract time series
print(' Extracting 3D time series')
Data3D = extractTimeSeries(slf3D,NNodes,NWghts,Loc,showbar)
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 getTransectPoints(loc01,loc02,dL):
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
[docs]def extractTransect(filePath,runStr,tranStr,loc01,loc02,dL,showbar=False):
# 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, Data2D, Data3D] = extractRunData(filePath,runStr,Loc,showbar)
U = Data3D['TS_INTRP_U']
V = Data3D['TS_INTRP_V']
W = Data3D['TS_INTRP_W']
[Vpp,Vtr,Up] = rotateVectorField(U,V,W,Theta)
Trans_z[ipt,:,:] = Data3D['TS_INTRP_z']
Trans_U[ipt,:,:] = Data3D['TS_INTRP_U']
Trans_V[ipt,:,:] = Data3D['TS_INTRP_V']
Trans_W[ipt,:,:] = Data3D['TS_INTRP_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['timeUnits'] = 'days since '+Data2D['TS_epoch']
Trans_Data['epoch'] = Data2D['TS_epoch']
Trans_Data['t'] = Data2D['TS_t']
Trans_Data['z'] = Trans_z
Trans_Data['loc'] = Locs
Trans_Data['U'] = Trans_U
Trans_Data['V'] = Trans_V
Trans_Data['W'] = Trans_W
Trans_Data['Vpp'] = Trans_Vpp
Trans_Data['Vtr'] = Trans_Vtr
# Return results
return Trans_Data
[docs]def enhanceMesh(mesh):
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
[docs]def extractMesh(filePath,runStr,setStr=None,modelVersion=None,showbar=False):
# Load 2D data file into SELAFIN structure
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)
idate = slf.datetime
times = slf.tags['times']
nsteps = len(times)
nnodes = slf.npoin2
meshx = slf.meshx
meshy = slf.meshy
elems = slf.ikle3
nelems = slf.nelem3
# Construct time variable
t = np.reshape(Date2Num(idate) + (times / (24.*60.0*60.)),[nsteps,1])
# Construct data dictionary
meshData = {}
meshData['times'] = t
meshData['nNodes'] = nnodes
meshData['meshx'] = meshx
meshData['meshy'] = meshy
meshData['nElems'] = nelems
meshData['elems'] = elems
meshData = enhanceMesh(meshData)
# Save mesh data
if modelVersion is None:
meshFile = runStr+'_mesh.mat'
else:
vstr = 'v'+str(modelVersion).zfill(2)
meshFile = runStr+'_mesh_'+vstr+'.mat'
extractPath = filePath.replace('RESULTS','EXTRACT')
try:
os.mkdir(extractPath)
except:
pass
sio.savemat(os.path.join(extractPath,meshFile).replace('\\','/'),meshData)
return
[docs]def extractGeometry(filePath,geomFile,modelVersion=None,saveMesh=False,showbar=False):
# Load Geometry data file into SELAFIN structure
slfFile = os.path.join(filePath,geomFile)
slf = InitialiseFile(slfFile)
nnodes = slf.npoin2
meshx = slf.meshx
meshy = slf.meshy
elems = slf.ikle3
nelems = slf.nelem3
varNames = slf.varnames
# Get variables
vars = slf.get_values(0)
nvars = np.shape(vars)[0]
bathy = vars[0]
if nvars > 1:
if varNames[1] == 'BOTTOM FRICTION ':
ssmap = vars[1]
has_ssmap = True
else:
has_ssmap = False
# Generate global attribute record
glbAttr = generateGlobalAttributes('Model Geometry - '+geomFile)
# Construct data dictionary
meshData = {}
meshData['globalAttributes'] = glbAttr
meshData['nNodes'] = nnodes
meshData['meshx'] = meshx
meshData['meshy'] = meshy
meshData['nElems'] = nelems
meshData['elems'] = elems
# Enhance mesh data - include areas, centroids, etc.
meshData = enhanceMesh(meshData)
# Add bathymetry and bottom friction (if it exists)
meshData['bathy'] = bathy
if has_ssmap:
# NOTE: ssmap = seabed substrate map (converted to friction coefficient).
meshData['fcoeff'] = ssmap
# Save mesh data
if saveMesh:
if modelVersion is None:
meshFile = 'MESH.mat'
else:
vstr = 'v'+str(modelVersion).zfill(2)
meshFile = 'MESH_'+vstr+'.mat'
extractPath = filePath.replace('GEO','EXTRACT')
sio.savemat(os.path.join(extractPath,meshFile),meshData)
return meshData
[docs]def extractElevation(filePath,runStr,modelVersion,showbar=False):
# Load 2D data file into SELAFIN structure
slfFile = os.path.join(filePath,(runStr+'_2D.slf'))
slf = InitialiseFile(slfFile)
times = slf.tags['times']
nsteps = len(times)
varNames = slf.varnames
varIndex = varNames.index('FREE SURFACE ')
extractPath = filePath.replace('RESULTS','EXTRACT')
for istp in np.arange(nsteps):
print(istp)
vars = slf.get_values(istp)
elevData = {}
elevData['elev'] = vars[varIndex,:]
vstr = 'v'+str(modelVersion).zfill(2)
elevFile = runStr+'_elev_frame'+str(istp).zfill(4)+'_'+vstr+'.mat'
sio.savemat(os.path.join(extractPath,elevFile),elevData)
del vars, elevData, elevFile
return
[docs]def extractMultiFileElevation(filePath,modStr,startStr,nFiles,saveFile=True):
files = getListOfFiles(filePath,'*'+modStr+'*_2D.slf')
ifile = 0
irec = 0
epoch_num = Date2Num(default_epoch)
for f in files:
if (startStr in f) & (ifile == 0):
ifile = 1
if (ifile > 0) & (ifile <= nFiles):
# Load 2D data file into SELAFIN structure
print('File: ',f)
slfFile = os.path.join(filePath,f)
slf = InitialiseFile(slfFile)
varNames = slf.varnames
varIndex = varNames.index('FREE SURFACE ')
times = slf.tags['times']
idate = slf.datetime
t_start = Date2Num(idate)
nsteps = len(times)
if ifile == 1:
nNodes = slf.npoin2
nRecs = nsteps*nFiles-(nFiles-1)
print('Number of Nodes = ',nNodes)
print('Number of records = ',nRecs)
tims = np.empty((nRecs,))
elev = np.empty((nRecs,nNodes))
for istp in np.arange(nsteps):
if ((istp == 0) and (ifile == 1)) or (istp > 0):
tims[irec] = t_start + (times[istp] / (24.*60.0*60.))
vals = slf.get_variables_at(istp,[varIndex])
elev[irec,:] = vals
irec += 1
del vals
ifile += 1
elevData = {}
elevData['time_units'] = 'days since '+default_epoch_str
elevData['epoch'] = default_epoch_str
elevData['times'] = tims-epoch_num
elevData['elev'] = elev
if saveFile:
extractPath = filePath.replace('RESULTS','EXTRACT')
elevFile = modStr+'_elev_multi_day.mat'
sio.savemat(os.path.join(extractPath,elevFile),elevData)
return elevData
[docs]def singleProcess(filePath,runStr,idStr,modelVer,Loc,showbar=False):
# Extract data from run files
[NNodes, NWghts, NLocs, Data2D, Data3D] = extractRunData(filePath,runStr,Loc,showbar)
nodes = [NNodes, NWghts, NLocs]
# Save results
idStr = idStr+'_VEL'
outName = generateOutputFileName(filePath,runStr+idStr,modelVer)
#glbAttr = generateGlobalAttributes('Iroise Sea Model Calibration - '+os.path.split(os.path.splitext(outName)[0])[1])
glbAttr = generateGlobalAttributes('Model Scale Impact Analysis - '+os.path.split(os.path.splitext(outName)[0])[1])
save2mat(outName,glbAttr,Loc,nodes,Data2D,Data3D)
return
[docs]def batchProcess(filePath,searchStr,modelVer,idStr,Loc,showbar=False):
files = getListOfFiles(filePath,'*'+searchStr+'*_3D.slf')
for fname in files:
print('... Processing '+fname)
runStr = fname[0:-7]
singleProcess(filePath,runStr,modelVer,idStr,Loc,showbar)
del runStr
return
[docs]def generateOutputFileName(filePath,runStr,modelVer):
outPath = filePath.replace('RESULTS','EXTRACT')
mvStr = 'v'+str(modelVer).zfill(2)
outName = runStr+'_'+mvStr+'.mat'
outFName = os.path.join(outPath,outName)
return outFName
[docs]def save2mat(outfname,glbAttr,Loc,nodes,Data2D,Data3D):
dict = {}
dict['globalAttributes'] = glbAttr
dict['location'] = Loc
dict['nearestNodes'] = nodes[0]
dict['nodeWeights'] = nodes[1]
dict['nodeLocs'] = nodes[2]
dict['timeUnits'] = 'days since '+Data2D['TS_epoch']
dict['epoch'] = Data2D['TS_epoch']
dict['times'] = Data2D['TS_t']
dict['depth'] = Data2D['TS_d']
dict['height'] = Data2D['TS_h']
dict['z'] = Data3D['TS_z']
dict['DepAvg_U'] = Data2D['TS_U']
dict['DepAvg_V'] = Data2D['TS_V']
dict['vel_east'] = Data3D['TS_U']
dict['vel_north'] = Data3D['TS_V']
dict['vel_up'] = Data3D['TS_W']
dict['depth_loc'] = Data2D['TS_INTRP_d']
dict['height_loc'] = Data2D['TS_INTRP_h']
dict['z_loc'] = Data3D['TS_INTRP_z']
dict['DepAvg_U_loc'] = Data2D['TS_INTRP_U']
dict['DepAvg_V_loc'] = Data2D['TS_INTRP_V']
dict['vel_east_loc'] = Data3D['TS_INTRP_U']
dict['vel_north_loc'] = Data3D['TS_INTRP_V']
dict['vel_up_loc'] = Data3D['TS_INTRP_W']
sio.savemat(outfname,dict)
return
[docs]def interpolateModelProfile(zint,z,val):
nstep = len(z)
interpVal = np.zeros((nstep,))
for istep in range(nstep):
interpVal[istep] = np.interp(zint,z[istep,:],val[istep,:])
return interpVal
[docs]def extractModelHubHeight(filePath,searchStr,z_hub):
files = getListOfFiles(filePath,searchStr)
times = []
depth = []
vel_east = []
vel_north = []
vel_up = []
for fname in files:
tsdata = loadMatFile(os.path.join(filePath,fname))
if 'epoch' not in locals():
epoch = tsdata['epoch']
times = np.concatenate((times,np.squeeze(tsdata['times'])),axis=0)
depth = np.concatenate((depth,np.squeeze(tsdata['depth'])),axis=0)
z = tsdata['z']-tsdata['z'][0,0]
vel = interpolateModelProfile(z_hub,z,tsdata['vel_east'])
vel_east = np.concatenate((vel_east,vel),axis=0)
del vel
vel = interpolateModelProfile(z_hub,z,tsdata['vel_north'])
vel_north = np.concatenate((vel_north,vel),axis=0)
del vel
vel = interpolateModelProfile(z_hub,z,tsdata['vel_up'])
vel_up = np.concatenate((vel_up,vel),axis=0)
del tsdata, vel
data = {}
data['epoch'] = epoch
data['times'] = times
data['depth'] = depth
data['hub_hght'] = z_hub
data['vel_east'] = vel_east
data['vel_north'] = vel_north
data['vel_up'] = vel_up
return data
[docs]def extractModelAtHAB(filePath,searchStr,z_hab):
files = getListOfFiles(filePath,searchStr)
times = []
depth = []
vel_east = []
vel_north = []
vel_up = []
for fname in files:
tsdata = loadMatFile(os.path.join(filePath,fname))
if 'epoch' not in locals():
epoch = tsdata['epoch']
times = np.concatenate((times,np.squeeze(tsdata['times'])),axis=0)
depth = np.concatenate((depth,np.squeeze(tsdata['depth'])),axis=0)
z = tsdata['z']-tsdata['z'][0,0]
vel = interpolateModelProfile(z_hab,z,tsdata['vel_east'])
vel_east = np.concatenate((vel_east,vel),axis=0)
del vel
vel = interpolateModelProfile(z_hab,z,tsdata['vel_north'])
vel_north = np.concatenate((vel_north,vel),axis=0)
del vel
vel = interpolateModelProfile(z_hab,z,tsdata['vel_up'])
vel_up = np.concatenate((vel_up,vel),axis=0)
del tsdata, vel
data = {}
data['epoch'] = epoch
data['times'] = times
data['depth'] = depth
data['hab'] = z_hab
data['vel_east'] = vel_east
data['vel_north'] = vel_north
data['vel_up'] = vel_up
return data
# ------------------------------------------------
# Functions for processing Node data
# ------------------------------------------------
[docs]def saveNodes2mat(outfname,glbAttr,nodes,Data2D,Data3D):
dict = {}
dict['globalAttributes'] = glbAttr
dict['nodes'] = nodes
dict['nodesLoc'] = Data2D['nodesLoc']
dict['timeUnits'] = 'days since '+Data2D['TS_epoch']
dict['epoch'] = Data2D['TS_epoch']
dict['times'] = Data2D['TS_t']
dict['depth'] = Data2D['TS_d']
dict['height'] = Data2D['TS_h']
dict['DepAvg_U'] = Data2D['TS_U']
dict['DepAvg_V'] = Data2D['TS_V']
dict['z'] = Data3D['TS_z']
dict['vel_east'] = Data3D['TS_U']
dict['vel_north'] = Data3D['TS_V']
dict['vel_up'] = Data3D['TS_W']
sio.savemat(outfname,dict)
return
[docs]def getNode2DVars(vals,varNames):
U = getVariableData(vals,varNames,'VELOCITY U ')
V = getVariableData(vals,varNames,'VELOCITY V ')
d = getVariableData(vals,varNames,'WATER DEPTH ')
h = getVariableData(vals,varNames,'FREE SURFACE ')
# Store variables in dictionary for return
vars2D = {}
vars2D['TS_U'] = U
vars2D['TS_V'] = V
vars2D['TS_d'] = d
vars2D['TS_h'] = h
return vars2D
[docs]def getNode3DVars(vals,varNames):
z = getVariableData(vals,varNames,'ELEVATION Z ')
U = getVariableData(vals,varNames,'VELOCITY U ')
V = getVariableData(vals,varNames,'VELOCITY V ')
W = getVariableData(vals,varNames,'VELOCITY W ')
# Store variables in dictionary for return
vars3D = {}
vars3D['TS_z'] = z
vars3D['TS_U'] = U
vars3D['TS_V'] = V
vars3D['TS_W'] = W
return vars3D
[docs]def extractNodeTimeSeries(slfHnd,nodes,showbar=False):
nnodes = len(nodes)
idate = slfHnd.datetime
times = slfHnd.tags['times']
nstps = np.size(times)
npts = slfHnd.npoin2;
nnds = slfHnd.npoin3
nlays = slfHnd.nplan
nodesLoc = np.asarray([slfHnd.meshx[nodes],slfHnd.meshy[nodes]]).transpose()
varNames = slfHnd.varnames
# Get Variables
tsData = {}
if npts == nnds:
# 2D Variables
TS_t = genTimeStamps(idate, times, nstps)
# NOTE: The Selafin get_series function assumes node numbering
# starts at 1, hence the nodes+1 being passed.
vals = slfHnd.get_series(nodes+1,showbar=showbar)
vars = getNode2DVars(vals,varNames)
tsData['TS_t'] = TS_t
tsData['TS_epoch'] = Date2Str(default_epoch)
tsData['nodesLoc'] = nodesLoc
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_z = np.zeros((nnodes,nstps,nlays))
TS_U = np.zeros((nnodes,nstps,nlays))
TS_V = np.zeros((nnodes,nstps,nlays))
TS_W = np.zeros((nnodes,nstps,nlays))
for ilay in np.arange(nlays):
print(' Extracting layer '+str(ilay+1)+' of '+str(nlays))
lnodes = nodes + ilay * npts
vals = slfHnd.get_series(lnodes+1,showbar=showbar)
vars = getNode3DVars(vals,varNames)
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 vals, lnodes, vars
tsData['TS_z'] = TS_z
tsData['TS_U'] = TS_U
tsData['TS_V'] = TS_V
tsData['TS_W'] = TS_W
return tsData
[docs]def extractNodeData(filePath,fileStr,nodes,showbar=False):
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 = extractNodeTimeSeries(slf2D,nodes,showbar)
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'))
print(' Initialize 3D SELAFIN file')
slf3D = InitialiseFile(file3D)
print(' Extracting 3D time series')
Data3D = extractNodeTimeSeries(slf3D,nodes,showbar)
toc = time.time()
print('Time elapsed (3D) = ' + str(toc-tic1))
print('Total time elapsed = ' + str(toc-tic))
return [Data2D, Data3D]
[docs]def processNodeDataFile(filePath,runStr,idStr,modelVer,procStr,nodes,showbar=False):
# Extract data from run files
[Data2D, Data3D] = extractNodeData(filePath,runStr,nodes,showbar)
# Save results
outName = generateOutputFileName(filePath,runStr+idStr,modelVer)
glbAttr = generateGlobalAttributes(procStr+' - '+os.path.split(os.path.splitext(outName)[0])[1])
saveNodes2mat(outName,glbAttr,nodes,Data2D,Data3D)
print(' Processing of '+runStr+' complete.')
print()
return
# ------------------------------------------------
# Functions for processing Location data
# ------------------------------------------------
[docs]def saveLoc2mat(outfname,glbAttr,loc,crs,Data2D,Data3D):
store = {}
store['globalAttributes'] = glbAttr
store['dataLoc'] = loc
store['locCRS'] = 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']
store['z'] = Data3D['TS_z']
store['vel_east'] = Data3D['TS_U']
store['vel_north'] = Data3D['TS_V']
store['vel_up'] = Data3D['TS_W']
sio.savemat(outfname,store)
return
[docs]def getLocs2DVars(vals,wgts,varNames):
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_INTRP_U = np.sum(U*wgts,0)
TS_INTRP_V = np.sum(V*wgts,0)
TS_INTRP_d = np.sum(d*wgts,0)
TS_INTRP_h = np.sum(h*wgts,0)
# Store variables in dictionary for return
vars2D = {}
vars2D['TS_U'] = U
vars2D['TS_V'] = V
vars2D['TS_d'] = d
vars2D['TS_h'] = h
vars2D['TS_INTRP_U'] = TS_INTRP_U[:,np.newaxis]
vars2D['TS_INTRP_V'] = TS_INTRP_V[:,np.newaxis]
vars2D['TS_INTRP_d'] = TS_INTRP_d[:,np.newaxis]
vars2D['TS_INTRP_h'] = TS_INTRP_h[:,np.newaxis]
return vars2D
[docs]def getLocs3DVars(vals,wgts,varNames):
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_INTRP_z = np.sum(z*wgts,0)
TS_INTRP_U = np.sum(U*wgts,0)
TS_INTRP_V = np.sum(V*wgts,0)
TS_INTRP_W = np.sum(W*wgts,0)
# Store variables in dictionary for return
vars3D = {}
vars3D['TS_z'] = z
vars3D['TS_U'] = U
vars3D['TS_V'] = V
vars3D['TS_W'] = W
vars3D['TS_INTRP_z'] = TS_INTRP_z[:,np.newaxis]
vars3D['TS_INTRP_U'] = TS_INTRP_U[:,np.newaxis]
vars3D['TS_INTRP_V'] = TS_INTRP_V[:,np.newaxis]
vars3D['TS_INTRP_W'] = TS_INTRP_W[:,np.newaxis]
return vars3D
[docs]def extractLocTimeSeries(slfHnd,Elem,NNodes,NWghts,Loc,showbar=False):
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=showbar) # Have to offset Node indices as code assumes numbering from 1 not 0.
# Get Variables
tsData = {}
if npts == nnds:
# 2D Variables
TS_t = genTimeStamps(idate, times, nstps)
vars = getLocs2DVars(vals,wgts,varNames)
tsData['TS_elem'] = Elem
tsData['TS_t'] = TS_t
tsData['TS_epoch'] = Date2Str(default_epoch)
if 'TS_U' in vars:
tsData['TS_NODES_U'] = vars['TS_U']
if 'TS_V' in vars:
tsData['TS_NODES_V'] = vars['TS_V']
if 'TS_d' in vars:
tsData['TS_NODES_d'] = vars['TS_d']
if 'TS_h' in vars:
tsData['TS_NODES_h'] = vars['TS_h']
if 'TS_INTRP_U' in vars:
tsData['TS_U'] = vars['TS_INTRP_U']
if 'TS_INTRP_V' in vars:
tsData['TS_V'] = vars['TS_INTRP_V']
if 'TS_INTRP_d' in vars:
tsData['TS_d'] = vars['TS_INTRP_d']
if 'TS_INTRP_h' in vars:
tsData['TS_h'] = vars['TS_INTRP_h']
else:
# 3D Variables
TS_z = np.zeros((3,nstps,nlays))
TS_U = np.zeros((3,nstps,nlays))
TS_V = np.zeros((3,nstps,nlays))
TS_W = np.zeros((3,nstps,nlays))
TS_INTRP_z = np.zeros((nstps,nlays))
TS_INTRP_U = np.zeros((nstps,nlays))
TS_INTRP_V = np.zeros((nstps,nlays))
TS_INTRP_W = np.zeros((nstps,nlays))
for ilay in np.arange(nlays):
lpts = np.arange(3)+ilay*3
ldat = np.squeeze(vals[:,lpts,:].copy())
vars = getLocs3DVars(ldat,wgts,varNames)
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'])
TS_INTRP_z[:,ilay] = np.squeeze(vars['TS_INTRP_z'])
TS_INTRP_U[:,ilay] = np.squeeze(vars['TS_INTRP_U'])
TS_INTRP_V[:,ilay] = np.squeeze(vars['TS_INTRP_V'])
TS_INTRP_W[:,ilay] = np.squeeze(vars['TS_INTRP_W'])
del ldat
tsData['TS_NODES_z'] = TS_z
tsData['TS_NODES_U'] = TS_U
tsData['TS_NODES_V'] = TS_V
tsData['TS_NODES_W'] = TS_W
tsData['TS_z'] = TS_INTRP_z
tsData['TS_U'] = TS_INTRP_U
tsData['TS_V'] = TS_INTRP_V
tsData['TS_W'] = TS_INTRP_W
return tsData
[docs]def extract3DLocData(filePath,fileStr,Loc,showbar=False):
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 mesh definition
meshx = np.transpose(slf2D.meshx.copy())
meshy = np.transpose(slf2D.meshy.copy())
elems = slf2D.ikle2
# Get Triangle Nodes and Barycenbtric Weights
Elem, NNodes,NWghts,NLocs = getNodesWeights(meshx,meshy,elems,Loc[0],Loc[1])
# Extract time series
print(' Extracting 2D time series')
Data2D = extractLocTimeSeries(slf2D,Elem,NNodes,NWghts,Loc,showbar)
# Clear local variables
del Elem,NNodes,NWghts,NLocs
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('\\','/')
slf3D = InitialiseFile(file3D)
# 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 = extractLocTimeSeries(slf3D,Elem,NNodes,NWghts,Loc,showbar)
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,modelVer,procStr,locs,showbar=False):
crs = 'EPSG:32630'
for iloc,loc in enumerate(locs):
# Extract data from run files
[NNodes, NWghts, NLocs, Data2D, Data3D] = extract3DLocData(filePath,runStr,loc,showbar)
# Save results
runIDStr = runStr+'_'+locIDs[iloc]
outName = generateOutputFileName(filePath,runIDStr,modelVer)
glbAttr = generateGlobalAttributes(procStr+' - '+os.path.split(os.path.splitext(outName)[0])[1])
saveLoc2mat(outName,glbAttr,loc,crs,Data2D,Data3D)
print(' Processing of '+runStr+' complete.')
print()
return
# ------------------------------------------------
# Functions for parsing ASCII run files
# ------------------------------------------------
[docs]def parseLocsRunFile(fileName):
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):
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):
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