# -*- coding: utf-8 -*-
"""
Bathymetry manipulation tools.
Chris Old
IES, School of Engineering, University of Edinburgh
Nov 2023
"""
# ----------------------------------------------------------------------------
# IMPORTS
# ----------------------------------------------------------------------------
# Standard Python Dependencies
import os
import numpy as np
# Non-Standard Python Dependencies
import metpy.interpolate as interp
# Local Module Dependencies
from ocmw.dataman.dataReaders import ncread
# Other Dependencies
#--------------------------------------------------------------------------
# GLOBAL CONSTANTS
#--------------------------------------------------------------------------
# ----------------------------------------------------------------------------
# CLASS DEFINITIONS
# ----------------------------------------------------------------------------
#--------------------------------------------------------------------------
# FUNCTION DEFINITIONS
#--------------------------------------------------------------------------
[docs]
def load_bathy(bathyfile,boffsets=[0.0,0.0,0.0]):
"""
Load bathymetry data from ASCII .xyz or.csv files, or netCDF files
Parameters
----------
bathyfile : str
Full path and name to the bathymery data file.
boffsets : list, float, optional
Translation offsets if required. The default is [0.0,0.0,0.0].
Returns
-------
bathy : TYPE
DESCRIPTION.
"""
fname, fextn = os.path.splitext(bathyfile)
if fextn.lower() in ['.xyz','.csv']:
bathy = np.genfromtxt(bathyfile, delimiter=',')
elif fextn.lower() == '.nc':
x = ncread(bathyfile,'x')[:].data
dx = np.nanmean(np.diff(x))
x = x + boffsets[0]*dx
print('Bathy x range: ',np.min(x),np.max(x))
y = ncread(bathyfile,'y')[:].data
dy = np.nanmean(np.diff(y))
y = y + boffsets[1]*dy
print('Bathy y range: ',np.min(y),np.max(y))
#z = ncread(bathyfile,'Band1')[:].data.transpose()
z = ncread(bathyfile,'Band1')[:].data
z = z + boffsets[2]
print('Bathy data shape: ',z.shape)
print('Bathy data range: ', np.min(z),np.max(z))
bx, by = np.meshgrid(x,y)
px = bx.reshape((np.prod(bx.shape),))
py = by.reshape((np.prod(by.shape),))
pz = z.reshape((np.prod(z.shape),))
bathy = np.asarray((px,py,pz)).transpose()
print('New bathy data shape: ',bathy.shape)
print('New bathy data range: ', np.min(bathy[:,2]),np.max(bathy[:,2]))
return bathy
[docs]
def bathy_msl_offset(bathyfile,offset,outfile,delimiter=',',headerSkip=0):
"""
Function for adding MSL correction to ASCII bathymetr data.
Currently only works for *.xyz and *.csv files.
Parameters
----------
bathyfile : str
Full path and filename of ASCII bathymetry data file.
offset : float
Offset in m to be applied to correct bathymetry data to MSL.
outfile : str
Full path and filename of output data file for corrected bathymetry.
delimiter : str, optional
Delimiter used in ASCII file. The default is ','.
headerSkip : int, optional
Number of header lines to copy and skip. The default is 0.
Returns
-------
None.
"""
fname, fextn = os.path.splitext(bathyfile)
if fextn.lower() in ['.xyz','.csv']:
print('get header lines...')
if headerSkip > 0:
fhnd = open(bathyfile,'r')
headerlines = ''
for line in np.arange(headerSkip):
headerlines = headerlines+fhnd.readline()
headerlines = headerlines.strip()
fhnd.close()
print('load bathymetry data...')
bathy = np.genfromtxt(bathyfile, delimiter=delimiter,skip_header=headerSkip)
print('added MSL correction...')
bathy[:,2] = bathy[:,2] + offset
print('save corrected bathymetry to file...')
np.savetxt(outfile,bathy,fmt='%.6f',delimiter=delimiter,header=headerlines,comments='')
print('processing complete.')
else:
print('*'+fextn+' not supported')
return
[docs]
def interp_bathy(nodes,bathy,radius):
"""
Interpolation of bathymetry data onto nodel nodes.
Parameters
----------
nodes : numpy.ndarray, float
Mesh node coordinates [x,y].
bathy : numpy.ndarray, float
Bathyemtry data [x,y,z].
radius : float
Radius of influence for interpolation.
Returns
-------
z : numpy.array, float
Bathymetry values at node locations.
"""
z = interp.inverse_distance_to_points(bathy[:,0:2],
bathy[:,2],
nodes,
radius,
min_neighbors=2)
return z
#--------------------------------------------------------------------------
# MAIN PROCESS
#--------------------------------------------------------------------------
#--------------------------------------------------------------------------
# INTERFACE
#--------------------------------------------------------------------------