# -*- coding: utf-8 -*-
"""
@author: Rochlitz.R
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import RegularGridInterpolator
from custEM.meshgen import meshgen_utils as mu
[docs]
class Bathy:
"""
Coastline identification class for digital elevation models. Please note
that the basic methods are working properly for a couple of specific
models, but a more general application of this class needs some debugging
and further extensions.
"""
def __init__(self, bathy_file, x=None, y=None, centering=False,
easting_shift=None, northing_shift=None, min_dist=50.,
min_points=10, coast_refinement=None, do_nothing=False):
"""
Initialize DEM (regular grid) interpolation object.
Required arguments
------------------
- topo_file, type str or iterable elevations (list, ndarray)
digital elevation file:
* ASC file with lower left corner and spacing in header OR
* x, y, z list of grid points, no irregular grids supported so far
Keyword arguments
-----------------
- x, y = None, type iterable (list, ndarray)
unique x and y positions
- centering = False, type bool
set **True**, if the given DEM should be centered relative to the
computational domain coordinate frame
- easting_shift, northing_shift = None, type float
if centering is **False**, a custom offset of a given DEM relaive
to the computational domain in x- and y- direction can be set
- min_dist = 50., type float
minimum distance between two positions identified on the coastline
- min_points = 10, type int
minimum number of positions that need to be identified on a
coastline to be incorporated in the mesh; used to avoid
comparatively small water areas
- coast_refinement = None, type int
divide each segment on the coastline in *coast_refinement* smaller
segments; needs to be replaced by a *min_dist* formulation
- do_nothing = False, type bool
flag currently for development purposes only
"""
if do_nothing:
return
self.x = None
self.y = None
self.rotate = None
if isinstance(bathy_file, str):
if bathy_file[-4:].lower() == '.asc':
self.x_bathy, self.y_bathy, self.z_bathy = \
self.loadASC(bathy_file, centering,
easting_shift, northing_shift)
else:
self.x_bathy, self.y_bathy, self.z_bathy = \
self.loadTXT(bathy_file, centering,
easting_shift, northing_shift)
self.x_bathy, self.y_bathy = \
self.translate(self.x_bathy, self.y_bathy,
centering, easting_shift, northing_shift)
else:
raise Exception("correct DEM file (tsype str) must be given!")
print('DEM original extent (x min, x max, y min, y max, z):')
print(np.min(self.x_bathy), np.max(self.x_bathy),
np.min(self.y_bathy), np.max(self.y_bathy),
np.min(self.z_bathy), np.max(self.z_bathy))
# f, ax = plt.subplots(1, figsize=(9, 6))
# obj = ax.imshow(self.z_bathy, vmin=-100., vmax=0.,
# cmap='ocean')
# f.colorbar(obj, ax=ax)
# plt.savefig('bla.png')
# raise SystemExit
self.level = 0.
self.atol = 1e-2
self.dx = np.diff(self.x_bathy[:2])
self.min_dist = min_dist
self.min_points = min_points
self.marine_marker = []
self.eval_coast_line()
self.coast_line = np.array(self.coast_line)
if coast_refinement is not None:
tmp_lines = []
# if type(coast_refinement) is not int:
# print('Error! coast_refinement must be of *type int*!')
# raise SystemExit
# ll = coast_refinement + 1
# for seg in self.coast_line:
# new_line = seg[0]
# arr = np.zeros((ll, 3))
# for j in range(len(seg) - 1):
# if (((seg[j, 0] or seg[j+1, 0]) in
# [np.min(self.x_bathy), np.max(self.x_bathy)]) or
# ((seg[j, 1] or seg[j+1, 1]) in
# [np.min(self.y_bathy), np.max(self.y_bathy)])):
# new_line = np.vstack((new_line, seg[j+1]))
# else:
# print('HACKED!')
# print('HACKED! LINE 99 in bathy tools')
# print('HACKED!')
# if seg[j, 0] > -2500. and seg[j, 1] > -500.:
# ll = 5
# arr = np.zeros((ll, 3))
# else:
# ll = coast_refinement + 1
# arr = np.zeros((ll, 3))
# arr[:, 0] = np.linspace(seg[j, 0], seg[j+1, 0], ll)
# arr[:, 1] = np.linspace(seg[j, 1], seg[j+1, 1], ll)
# new_line = np.vstack((new_line, arr[1:]))
# arr[:, 0] = np.linspace(seg[-1, 0], seg[0, 0], ll)
# arr[:, 1] = np.linspace(seg[-1, 1], seg[0, 1], ll)
# new_line = np.vstack((new_line, arr[1:-1]))
# tmp_lines.append(new_line)
# better version with min_length
for seg in self.coast_line:
print('need to rework coast refinement')
bla
new_line = seg[0]
arr = np.zeros((ll, 3))
for j in range(len(seg) - 1):
if (((seg[j, 0] or seg[j+1, 0]) in
[np.min(self.x_bathy), np.max(self.x_bathy)]) or
((seg[j, 1] or seg[j+1, 1]) in
[np.min(self.y_bathy), np.max(self.y_bathy)])):
new_line = np.vstack((new_line, seg[j+1]))
else:
print('HACKED!')
print('HACKED! LINE 99 in bathy tools')
print('HACKED!')
if seg[j, 0] > -2500. and seg[j, 1] > -500.:
ll = 5
arr = np.zeros((ll, 3))
else:
ll = coast_refinement + 1
arr = np.zeros((ll, 3))
arr[:, 0] = np.linspace(seg[j, 0], seg[j+1, 0], ll)
arr[:, 1] = np.linspace(seg[j, 1], seg[j+1, 1], ll)
new_line = np.vstack((new_line, arr[1:]))
arr[:, 0] = np.linspace(seg[-1, 0], seg[0, 0], ll)
arr[:, 1] = np.linspace(seg[-1, 1], seg[0, 1], ll)
new_line = np.vstack((new_line, arr[1:-1]))
tmp_lines.append(new_line)
self.coast_line = tmp_lines
[docs]
def loadTXT(self, topo_file, centering,
easting_shift, northing_shift):
"""
Load *txt* file with column-based list of DEM coordinates. This
function is redudant to an identical function in the **DEM** class
and will be removed in a major update of the bathy-tools.
Required arguments
------------------
- see **Bathy** class description
"""
xp, yp, zp = np.loadtxt(topo_file, unpack=True)
if len(np.unique(xp)) * len(np.unique(yp)) > len(xp):
print('Only regular grids for coast_line search supported')
raise SystemExit
if np.isclose(yp[0], yp[1], rtol=1e-32, atol=1e-2):
if np.isclose(xp[0], xp[1], rtol=1e-32, atol=1e-2):
print('Fatal error! Neither first two x- nor y- coords are '
'increasing in the specified xyz file! Aborting...')
raise SystemExit
elif xp[1] > xp[0]:
nx = np.argwhere(np.diff(xp) < 0)[0][0] + 1
else:
nx = np.argwhere(np.diff(xp) > 0)[0][0] + 1
ny = len(xp) // nx
x = xp[:nx]
y = yp[::nx]
zp = zp.reshape((ny, nx))
else:
if yp[1] > yp[0]:
ny = np.argwhere(np.diff(yp) < 0)[0][0] + 1
else:
ny = np.argwhere(np.diff(yp) > 0)[0][0] + 1
nx = len(xp) // ny
x = xp[::ny]
y = yp[:ny]
zp = zp.reshape((nx, ny))
if y[1] < y[0]:
y.sort()
zp = np.flipud(zp)
if x[1] < x[0]:
x.sort()
zp = np.fliplr(zp)
self.nx = int(nx)
self.ny = int(ny)
return(x, y, zp)
[docs]
def loadASC(self, topo_file, centering,
easting_shift, northing_shift):
"""
Load *asc* file (DEM matrix with location header). This
function is redudant to an identical function in the **DEM** class
and will be removed in a major update of the bathy-tools.
Required arguments
------------------
- see **Bathy** class description
"""
with open(topo_file) as fid:
header = {}
sp = []
nheader = 0
while len(sp) < 3:
sp = fid.readline().split()
if len(sp) < 3:
header[sp[0]] = float(sp[1])
nheader += 1
zp = np.flipud(np.genfromtxt(topo_file, skip_header=nheader)).T
if 'NODATA_value' in header:
zp[zp == header['NODATA_value']] = np.nan
dx = header.pop('cellsize', 1.0) # just a guess
x = np.arange(header['ncols']) * dx + header['xllcorner']
y = np.arange(header['nrows']) * dx + header['yllcorner']
self.nx = int(header['ncols'])
self.ny = int(header['nrows'])
return(x, y, zp)
[docs]
def translate(self, x, y, centering=False, easting_shift=None,
northing_shift=None, back=False):
"""
Function to translate coordinates in horizontal direction.
Required arguments
------------------
- x, y, type list/array of coordinates
coordinates to be translated horizontally
Keyword arguments
-----------------
- see **Bathy** class description
"""
if centering:
if easting_shift is None and northing_shift is None:
if not back:
x -= np.mean(x)
y -= np.mean(y)
else:
x += np.mean(x)
y += np.mean(y)
else:
print('Warning! Never set *centering=True* while '
'*easting_shift* and *northing_shift* are not None, '
'Aborting before something unexpected happens ...')
raise SystemExit
if easting_shift is not None:
x -= (np.mean(x) - easting_shift)
if northing_shift is not None:
y -= (np.mean(y) - northing_shift)
elif not centering:
if easting_shift is not None:
x += easting_shift
if northing_shift is not None:
y += northing_shift
return(x, y)
[docs]
def eval_coast_line(self):
"""
Main method to find intersection of topo and zero-height level.
"""
print(' - evaluating coast line(s) from DEM data - ')
idxs, d = [], 0
m, n = self.z_bathy.shape[0], self.z_bathy.shape[1]
for k in range(m):
for j in range(n):
if not (k in (0, m - 1) or
j in (0, n - 1)):
if self.z_bathy[k, j] < 0. and \
(self.z_bathy[k - 1, j] >= 0. or
self.z_bathy[k, j - 1] >= 0. or
self.z_bathy[k + 1, j] >= 0. or
self.z_bathy[k, j + 1] >= 0. or
self.z_bathy[k - 1, j - 1] >= 0. or
self.z_bathy[k + 1, j - 1] >= 0. or
self.z_bathy[k + 1, j + 1] >= 0. or
self.z_bathy[k - 1, j + 1] >= 0.):
if 'neighbours' not in locals():
neighbours = np.zeros((1, 4), dtype=bool)
else:
neighbours = np.vstack((
neighbours, np.zeros((1, 4), dtype=bool)))
if j + 1 != n - 1:
neighbours[d, 0] = self.check_neighbours(k, j + 1)
if k + 1 != m - 1:
neighbours[d, 1] = self.check_neighbours(k + 1, j)
if j - 1 != 0:
neighbours[d, 2] = self.check_neighbours(k, j - 1)
if k - 1 != 0:
neighbours[d, 3] = self.check_neighbours(k - 1, j)
d += 1
idxs.append([k, j])
self.all_idxs = np.array(idxs)
self.neighbours = neighbours
self.n_neighbours = np.sum(self.neighbours, axis=1)
print('... need to check for ' + str(len(idxs)) + ' points ...')
self.coast_line, self.line_count, self.counter = [], 0, 0
self.every, found = len(idxs), 0
#
# np.save('ids.npy', self.all_idxs)
# np.save('land.npy', self.z_bathy)
# raise SystemExit
blubb = 0
while self.every > self.counter:
idx = self.find_start_point()
if idx is None:
print(' - skipping remaining single points - ')
return
print('START ID', idx)
self.find_direction(idx)
print('LEN', len(self.coast_line[self.line_count]))
# raise SystemExit
if len(self.coast_line[self.line_count]) < self.min_points:
del self.coast_line[-1]
del self.marine_marker[-1]
else:
self.line_count += 1
print('SUMMARY: ', self.every, self.counter)
print(len(self.all_idxs), 'LC: ', self.line_count,
len(self.marine_marker))
# np.save('ids' + str(blubb) + '.npy', self.all_idxs)
# np.save('land' + str(blubb) + '.npy', self.z_bathy)
blubb += 1
# print(self.coast_line)
print(self.line_count)
print('MIN DIST:', self.min_dist)
print('hacked break')
for elem in self.coast_line[0]:
print(elem)
break
# print('HARD CUT AFTER 2 ITERATION(S); NEED TO BE FIXED!')
# if self.line_count == 1:
# return
# self.marine_marker=np.array([[-2000., -2000., -100.]])
[docs]
def find_direction(self, idx, from_boundary=False):
"""
Identifiy clock- or counterclockwise search path direction based on
regular grid DEM data automatically. Does not work in all cases so far.
Required arguments
------------------
- idx, type list of [int, int]
x and y index of the DEM grid to start the direction search from
Keyword arguments
-----------------
- from_boundary = False, type bool
internally used flag to check if a coastline path was following the
horizontal mesh boundaries before entering the grid again
"""
ref = np.array(idx, dtype=int)
self.coast_line.append([])
dummy, found = 0, []
self.last_dir = None
while True:
blubb = False
pre_idx = np.array(idx, dtype=int)
for jjj in range(len(self.all_idxs)):
if idx[0] == self.all_idxs[jjj, 0] and \
idx[1] == self.all_idxs[jjj, 1]:
pos = jjj
blubb = True
break
if dummy == 0:
ref = self.all_idxs[pos]
print('\nJJJ', self.n_neighbours[jjj])
if dummy == 1:
self.all_idxs = np.vstack((self.all_idxs, ref))
if pos not in found:
found.extend([pos])
if not blubb:
idx = ref
if dummy != 0:
if idx[0] == ref[0] and idx[1] == ref[1]:
self.coast_line[self.line_count] = np.array(
self.coast_line[self.line_count])
self.marine_marker.append([self.x_bathy[idx[0]],
self.y_bathy[idx[1]], 0.])
self.all_idxs = np.delete(self.all_idxs, found, 0)
self.neighbours = np.delete(self.neighbours, found, 0)
self.n_neighbours = np.delete(self.n_neighbours, found, 0)
self.all_idxs = np.delete(self.all_idxs, -1, 0)
self.counter += len(found)
return
n_neighbours = self.n_neighbours[pos]
w, d, s, a = self.neighbours[pos]
if dummy == 0:
# print('WARNING! Clockwise or Counterclockwise direction choice'
# ' is still bugged in some cases! Continuing ...')
# if self.z_bathy[idx[1] - 1, idx[1]] < 0.:
# clockwise = False
# else:
# clockwise = True
clockwise = True
print('HACKED DIRECTION', clockwise)
if not from_boundary and (idx[1] == self.z_bathy.shape[1] - 2 or
idx[0] == 1 or
idx[1] == 1 or
idx[0] == self.z_bathy.shape[0] - 2):
self.build_closing_frame(idx, clockwise)
from_boundary = True
elif n_neighbours == 1 or (n_neighbours == 2 and dummy == 0):
if w:
self.add_point_to_coast_line(idx, 'w', clockwise)
idx[1] += 1
self.last_dir = 'w'
elif d:
self.add_point_to_coast_line(idx, 'd', clockwise)
idx[0] += 1
self.last_dir = 'd'
elif s:
self.add_point_to_coast_line(idx, 's', clockwise)
idx[1] -= 1
self.last_dir = 's'
elif a:
self.add_point_to_coast_line(idx, 'a', clockwise)
idx[0] -= 1
self.last_dir = 'a'
elif (n_neighbours >= 2 and dummy != 0) or \
(n_neighbours == 3 and dummy == 0):
if (from_boundary and not
(idx[1] == self.z_bathy.shape[1] - 2 or idx[0] == 1 or
idx[1] == 1 or idx[0] == self.z_bathy.shape[0] - 2)):
from_boundary = False
if self.last_dir == 'w' and not clockwise:
if a:
self.add_point_to_coast_line(idx, 'a', clockwise)
idx[0] -= 1
elif w:
self.add_point_to_coast_line(idx, 'w', clockwise)
idx[1] += 1
elif d:
self.add_point_to_coast_line(idx, 'd', clockwise)
idx[0] += 1
elif self.last_dir == 'w' and clockwise:
if d:
self.add_point_to_coast_line(idx, 'd', clockwise)
idx[0] += 1
elif w:
self.add_point_to_coast_line(idx, 'w', clockwise)
idx[1] += 1
elif a:
self.add_point_to_coast_line(idx, 'a', clockwise)
idx[0] -= 1
elif self.last_dir == 'd' and not clockwise:
if w:
self.add_point_to_coast_line(idx, 'w', clockwise)
idx[1] += 1
elif d:
self.add_point_to_coast_line(idx, 'd', clockwise)
idx[0] += 1
elif s:
self.add_point_to_coast_line(idx, 's', clockwise)
idx[1] -= 1
elif self.last_dir == 'd' and clockwise:
if s:
self.add_point_to_coast_line(idx, 's', clockwise)
idx[1] -= 1
elif d:
self.add_point_to_coast_line(idx, 'd', clockwise)
idx[0] += 1
elif w:
self.add_point_to_coast_line(idx, 'w', clockwise)
idx[1] += 1
elif self.last_dir == 's' and not clockwise:
if d:
self.add_point_to_coast_line(idx, 'd', clockwise)
idx[0] += 1
elif s:
self.add_point_to_coast_line(idx, 's', clockwise)
idx[1] -= 1
elif a:
self.add_point_to_coast_line(idx, 'a', clockwise)
idx[0] -= 1
elif self.last_dir == 's' and clockwise:
if a:
self.add_point_to_coast_line(idx, 'a', clockwise)
idx[0] -= 1
elif s:
self.add_point_to_coast_line(idx, 's', clockwise)
idx[1] -= 1
elif d:
self.add_point_to_coast_line(idx, 'd', clockwise)
idx[0] += 1
elif self.last_dir == 'a' and not clockwise:
if s:
self.add_point_to_coast_line(idx, 's', clockwise)
idx[1] -= 1
elif a:
self.add_point_to_coast_line(idx, 'a', clockwise)
idx[0] -= 1
elif w:
self.add_point_to_coast_line(idx, 'w', clockwise)
idx[1] += 1
elif self.last_dir == 'a' and clockwise:
if w:
self.add_point_to_coast_line(idx, 'w', clockwise)
idx[1] += 1
elif a:
self.add_point_to_coast_line(idx, 'a', clockwise)
idx[0] -= 1
elif s:
self.add_point_to_coast_line(idx, 's', clockwise)
idx[1] -= 1
self.find_last_dir(pre_idx, idx)
dummy += 1
if dummy % 1 == 0:
print('every 100th id:', idx)
if dummy == 10000:
print('Doing this while loop for 10.000 times now, something '
'is probably not correct here! Aborting ...')
raise SystemExit
[docs]
def add_point_to_coast_line(self, idx, direction, clockwise, tol=None):
"""
Add a coordinate to a coastline path if the distance to the previous
position on the coastline is larger than *min_dist*.
Required arguments
------------------
- idx, type list of [int, int]
x and y index of the DEM grid
- direction, type str
direction string, either 'w', 'a', 's', or 'd'
- clockwise, type bool
go in clock- or counterclockwise direction
Keyword arguments
-----------------
- tol = None, type float
custom search tolerance if required, not recommended to be changed
"""
if tol is None:
tol = self.atol
if direction == 'w':
if clockwise:
if self.z_bathy[idx[0] + 1, idx[1]] > self.level - tol:
self.add_point_towards_east(idx)
if self.z_bathy[idx[0], idx[1] + 1] > self.level - tol:
self.add_point_towards_north(idx)
if self.z_bathy[idx[0] - 1, idx[1]] > self.level - tol:
self.add_point_towards_west(idx)
else:
if self.z_bathy[idx[0] - 1, idx[1]] > self.level - tol:
self.add_point_towards_west(idx)
if self.z_bathy[idx[0], idx[1] + 1] > self.level - tol:
self.add_point_towards_north(idx)
if self.z_bathy[idx[0] + 1, idx[1]] > self.level - tol:
self.add_point_towards_east(idx)
elif direction == 'a':
if clockwise:
if self.z_bathy[idx[0], idx[1] + 1] > self.level - tol:
self.add_point_towards_north(idx)
if self.z_bathy[idx[0] - 1, idx[1]] > self.level - tol:
self.add_point_towards_west(idx)
if self.z_bathy[idx[0], idx[1] - 1] > self.level - tol:
self.add_point_towards_south(idx)
else:
if self.z_bathy[idx[0], idx[1] - 1] > self.level - tol:
self.add_point_towards_south(idx)
if self.z_bathy[idx[0] - 1, idx[1]] > self.level - tol:
self.add_point_towards_west(idx)
if self.z_bathy[idx[0], idx[1] + 1] > self.level - tol:
self.add_point_towards_north(idx)
elif direction == 's':
if clockwise:
if self.z_bathy[idx[0] - 1, idx[1]] > self.level - tol:
self.add_point_towards_west(idx)
if self.z_bathy[idx[0], idx[1] - 1] > self.level - tol:
self.add_point_towards_south(idx)
if self.z_bathy[idx[0] + 1, idx[1]] > self.level - tol:
self.add_point_towards_east(idx)
else:
if self.z_bathy[idx[0] + 1, idx[1]] > self.level - tol:
self.add_point_towards_east(idx)
if self.z_bathy[idx[0], idx[1] - 1] > self.level - tol:
self.add_point_towards_south(idx)
if self.z_bathy[idx[0] - 1, idx[1]] > self.level - tol:
self.add_point_towards_west(idx)
elif direction == 'd':
if clockwise:
if self.z_bathy[idx[0], idx[1] - 1] > self.level - tol:
self.add_point_towards_south(idx)
if self.z_bathy[idx[0] + 1, idx[1]] > self.level - tol:
self.add_point_towards_east(idx)
if self.z_bathy[idx[0], idx[1] + 1] > self.level - tol:
self.add_point_towards_north(idx)
else:
if self.z_bathy[idx[0], idx[1] + 1] > self.level - tol:
self.add_point_towards_north(idx)
if self.z_bathy[idx[0] + 1, idx[1]] > self.level - tol:
self.add_point_towards_east(idx)
if self.z_bathy[idx[0], idx[1] - 1] > self.level - tol:
self.add_point_towards_south(idx)
[docs]
def add_point_towards_north(self, idx):
"""
Add a coordinate to a coastline path in northern direction.
Required arguments
------------------
- idx, type list of [int, int]
x and y index of the DEM grid
"""
grad = (self.z_bathy[idx[0], idx[1]] -
self.z_bathy[idx[0], idx[1] + 1]) / self.dx
mx, my = 0., (self.z_bathy[idx[0], idx[1]] / grad)[0]
if self.check_distance(self.x_bathy[idx[0]] + mx,
self.y_bathy[idx[1]] + my):
self.coast_line[self.line_count].append(
[self.x_bathy[idx[0]] + mx,
self.y_bathy[idx[1]] + my, 0.])
[docs]
def add_point_towards_south(self, idx):
"""
Add a coordinate to a coastline path in southern direction.
Required arguments
------------------
- idx, type list of [int, int]
x and y index of the DEM grid
"""
grad = (self.z_bathy[idx[0], idx[1] - 1] -
self.z_bathy[idx[0], idx[1]]) / self.dx
mx, my = 0., (self.z_bathy[idx[0], idx[1]] / grad)[0]
if self.check_distance(self.x_bathy[idx[0]] + mx,
self.y_bathy[idx[1]] + my):
self.coast_line[self.line_count].append(
[self.x_bathy[idx[0]] + mx,
self.y_bathy[idx[1]] + my, 0.])
[docs]
def add_point_towards_east(self, idx):
"""
Add a coordinate to a coastline path in eastern direction.
Required arguments
------------------
- idx, type list of [int, int]
x and y index of the DEM grid
"""
grad = (self.z_bathy[idx[0], idx[1]] -
self.z_bathy[idx[0] + 1, idx[1]]) / self.dx
mx, my = (self.z_bathy[idx[0], idx[1]] / grad)[0], 0.
if self.check_distance(self.x_bathy[idx[0]] + mx,
self.y_bathy[idx[1]] + my):
self.coast_line[self.line_count].append(
[self.x_bathy[idx[0]] + mx,
self.y_bathy[idx[1]] + my, 0.])
[docs]
def add_point_towards_west(self, idx):
"""
Add a coordinate to a coastline path in western direction.
Required arguments
------------------
- idx, type list of [int, int]
x and y index of the DEM grid
"""
grad = (self.z_bathy[idx[0] - 1, idx[1]] -
self.z_bathy[idx[0], idx[1]]) / self.dx
mx, my = (self.z_bathy[idx[0], idx[1]] / grad)[0], 0.
if self.check_distance(self.x_bathy[idx[0]] + mx,
self.y_bathy[idx[1]] + my):
self.coast_line[self.line_count].append(
[self.x_bathy[idx[0]] + mx,
self.y_bathy[idx[1]] + my, 0.])
[docs]
def check_distance(self, x, y):
"""
Calculate distance between a coordinate [x, y] and the previous
position on the coastline.
Required arguments
------------------
- x, y, typefloat
x and y coordinates for the horizontal distance evaluation
"""
if len(self.coast_line[self.line_count]) > 0:
d = np.sqrt((self.coast_line[self.line_count][-1][0] - x)**2 +
(self.coast_line[self.line_count][-1][1] - y)**2)
d_init = np.sqrt((self.coast_line[self.line_count][0][0] - x)**2 +
(self.coast_line[self.line_count][0][1] - y)**2)
if d < self.min_dist or d_init < self.min_dist:
return(False)
else:
return(True)
else:
return(True)
[docs]
def go_clockwise(self, side, idx, from_corner=False):
"""
Continue coastline path in clockwise direction if the horizontal
mesh boundary is reached.
Required arguments
------------------
- side, type str
string specifying one of the four sides
- idx, type list of [int, int]
x and y index of the DEM grid
- from_corner = False, type bool
internally used flag to continue a path after a corner of the
horizontal mesh extent was reached
"""
print('boundary reached going clockwise ', side, ' ', idx)
if side == 'north':
if not from_corner:
self.coast_line[self.line_count].append(
[self.x_bathy[idx[0]], self.y_bathy[idx[1] + 1], 0.])
while idx[0] < self.z_bathy.shape[0]:
idx[0] += 1
if len(np.where((self.all_idxs ==
(idx[0], self.z_bathy.shape[1] - 2)).all(
axis=1))[0]) == 1:
self.coast_line[self.line_count].append(
[self.x_bathy[idx[0]],
self.y_bathy[idx[1] + 1], 0.])
return([idx[0], self.z_bathy.shape[1] - 2])
idx[0] -= 1
self.coast_line[self.line_count].append([self.x_bathy[idx[0] + 1],
self.y_bathy[idx[1] + 1],
0.])
return(self.go_clockwise('east', idx, True))
if side == 'west':
if not from_corner:
self.coast_line[self.line_count].append(
[self.x_bathy[idx[0] - 1], self.y_bathy[idx[1]], 0.])
while idx[1] < self.z_bathy.shape[1]:
idx[1] += 1
if len(np.where((self.all_idxs ==
(1, idx[1])).all(axis=1))[0]) == 1:
self.coast_line[self.line_count].append(
[self.x_bathy[idx[0] - 1],
self.y_bathy[idx[1]], 0.])
return([1, idx[1]])
idx[1] -= 1
self.coast_line[self.line_count].append([self.x_bathy[idx[0] - 1],
self.y_bathy[idx[1] + 1],
0.])
return(self.go_clockwise('north', idx, True))
if side == 'south':
if not from_corner:
self.coast_line[self.line_count].append(
[self.x_bathy[idx[0]], self.y_bathy[idx[1] - 1], 0.])
while idx[0] > 0:
idx[0] -= 1
if len(np.where((self.all_idxs ==
(idx[0], 1)).all(axis=1))[0]) == 1:
self.coast_line[self.line_count].append(
[self.x_bathy[idx[0]],
self.y_bathy[idx[1] - 1], 0.])
return([idx[0], 1])
idx[0] += 1
self.coast_line[self.line_count].append([self.x_bathy[idx[0] - 1],
self.y_bathy[idx[1] - 1],
0.])
return(self.go_clockwise('west', idx, True))
if side == 'east':
if not from_corner:
self.coast_line[self.line_count].append(
[self.x_bathy[idx[0] + 1], self.y_bathy[idx[1]], 0.])
while idx[1] > 0:
idx[1] -= 1
if len(np.where((self.all_idxs ==
(self.z_bathy.shape[0] - 2,
idx[1])).all(axis=1))[0]) == 1:
self.coast_line[self.line_count].append(
[self.x_bathy[idx[0] + 1],
self.y_bathy[idx[1]], 0.])
return([self.z_bathy.shape[0] - 2, idx[1]])
idx[1] += 1
self.coast_line[self.line_count].append([self.x_bathy[idx[0] + 1],
self.y_bathy[idx[1] - 1],
0.])
return(self.go_clockwise('south', idx, True))
[docs]
def go_counter_clockwise(self, side, idx, from_corner=False):
"""
Continue coastline path in counterclockwise direction if the horizontal
mesh boundary is reached.
Required arguments
------------------
- side, type str
string specifying one of the four sides
- idx, type list of [int, int]
x and y index of the DEM grid
- from_corner = False, type bool
internally used flag to continue a path after a corner of the
horizontal mesh extent was reached
"""
print('boundary reached going counterclockwise ', side, ' ', idx)
if side == 'north':
if not from_corner:
self.coast_line[self.line_count].append(
[self.x_bathy[idx[0]], self.y_bathy[idx[1] + 1], 0.])
while idx[0] > 0:
print(idx[0])
idx[0] -= 1
# if len(np.where((self.all_idxs ==
# (idx[0], self.z_bathy.shape[1] - 2)).all(
# axis=1))[0]) == 1:
# print(idx[0])
# vbs
# self.coast_line[self.line_count].append(
# [self.x_bathy[idx[0]],
# self.y_bathy[idx[1] + 1], 0.])
# return([idx[0], self.z_bathy.shape[1] - 2])
# asd
idx[0] += 1
self.coast_line[self.line_count].append([self.x_bathy[idx[0] - 1],
self.y_bathy[idx[1] + 1],
0.])
return(self.go_counter_clockwise('west', idx, True))
if side == 'west':
if not from_corner:
self.coast_line[self.line_count].append(
[self.x_bathy[idx[0] - 1], self.y_bathy[idx[1]], 0.])
while idx[1] > 0:
idx[1] -= 1
# if len(np.where((self.all_idxs ==
# (1, idx[1])).all(axis=1))[0]) == 1:
# self.coast_line[self.line_count].append(
# [self.x_bathy[idx[0] - 1],
# self.y_bathy[idx[1]], 0.])
# return([1, idx[1]])
idx[1] += 1
self.coast_line[self.line_count].append([self.x_bathy[idx[0] - 1],
self.y_bathy[idx[1] - 1],
0.])
return(self.go_counter_clockwise('south', idx, True))
if side == 'south':
if not from_corner:
self.coast_line[self.line_count].append(
[self.x_bathy[idx[0]], self.y_bathy[idx[1] - 1], 0.])
while idx[0] < self.z_bathy.shape[0]:
idx[0] += 1
if len(np.where((self.all_idxs ==
(idx[0], 1)).all(axis=1))[0]) == 1:
self.coast_line[self.line_count].append(
[self.x_bathy[idx[0]],
self.y_bathy[idx[1] - 1], 0.])
return([idx[0], 1])
idx[0] -= 1
self.coast_line[self.line_count].append([self.x_bathy[idx[0] - 1],
self.y_bathy[idx[1] - 1],
0.])
return(self.go_counter_clockwise('east', idx, True))
if side == 'east':
if not from_corner:
self.coast_line[self.line_count].append(
[self.x_bathy[idx[0] + 1], self.y_bathy[idx[1]], 0.])
while idx[1] < self.z_bathy.shape[1]:
idx[1] += 1
if len(np.where((self.all_idxs ==
(self.z_bathy.shape[0] - 2,
idx[1])).all(axis=1))[0]) == 1:
self.coast_line[self.line_count].append(
[self.x_bathy[idx[0] + 1],
self.y_bathy[idx[1]], 0.])
return([self.z_bathy.shape[0] - 2, idx[1]])
idx[1] += 1
self.coast_line[self.line_count].append([self.x_bathy[idx[0] + 1],
self.y_bathy[idx[1] - 1],
0.])
return(self.go_counter_clockwise('north', idx, True))
[docs]
def build_closing_frame(self, idx, clockwise):
"""
Continue coastline path in clockwise or counterclockwise direction
if the horizontal mesh boundary is reached.
Required arguments
------------------
- idx, type list of [int, int]
x and y index of the DEM grid
- clockwise, type bool
go in clock- or counterclockwise direction
"""
if idx[1] == self.z_bathy.shape[1] - 2:
if not clockwise:
idx = self.go_clockwise('north', idx)
else:
idx = self.go_counter_clockwise('north', idx)
elif idx[0] == 1:
if not clockwise:
idx = self.go_clockwise('west', idx)
else:
idx = self.go_counter_clockwise('west', idx)
elif idx[1] == 1:
if not clockwise:
idx = self.go_clockwise('south', idx)
else:
idx = self.go_counter_clockwise('south', idx)
elif idx[0] == self.z_bathy.shape[0] - 2:
if not clockwise:
idx = self.go_clockwise('east', idx)
else:
idx = self.go_counter_clockwise('east', idx)
return(idx)
[docs]
def check_neighbours(self, kk, jj):
"""
Check z-values of the eight neighbors around a position in the DEM.
Utility function that needs an update.
Required arguments
------------------
- kk, jj, type int
x/y indices of the DEM grid
"""
if self.z_bathy[kk, jj] < 0. and\
(self.z_bathy[kk - 1, jj] >= 0. or
self.z_bathy[kk, jj - 1] >= 0. or
self.z_bathy[kk + 1, jj] >= 0. or
self.z_bathy[kk, jj + 1] >= 0. or
self.z_bathy[kk - 1, jj - 1] >= 0. or
self.z_bathy[kk + 1, jj - 1] >= 0. or
self.z_bathy[kk + 1, jj + 1] >= 0. or
self.z_bathy[kk - 1, jj + 1] >= 0.):
return(True)
else:
return(False)
[docs]
def find_last_dir(self, pre_idx, idx):
"""
Evaluate previous direction during the coastline search.
Required arguments
------------------
- pre_idx, type list of [int, int]
previous x and y index of the DEM grid
- idx, type list of [int, int]
x and y index of the DEM grid
"""
if pre_idx[1] == idx[1] - 1:
self.last_dir = 'w'
elif pre_idx[0] == idx[0] - 1:
self.last_dir = 'd'
elif pre_idx[1] == idx[1] + 1:
self.last_dir = 's'
elif pre_idx[0] == idx[0] + 1:
self.last_dir = 'a'
[docs]
def find_start_point(self, jjj=0):
"""
Evaluate start indices for a new coastline search in the DEM.
Keyword arguments
-----------------
- jjj, type int
counter to continute until all indices in the DEM-grid were
considered
"""
while True:
try:
if self.all_idxs[jjj, 0] not in [1, self.z_bathy.shape[0]] and\
self.all_idxs[jjj, 1] not in [1, self.z_bathy.shape[1]]:
return([self.all_idxs[jjj, 0], self.all_idxs[jjj, 1]])
jjj += 1
except IndexError:
return(None)