# -*- coding: utf-8 -*-
"""
@authors: Rochlitz.R & Günther.T
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import RegularGridInterpolator
from custEM.meshgen import meshgen_utils as mu
import matplotlib.tri as mtri
[docs]
class DEM:
"""
Interpolation class for digital elevation models (DEM). Automatically
creates interpolation objects during the class initialization, based on
DEM information stored as list of coordinates in *txt* format
or as matrix in *asc* format.
"""
def __init__(self, demfile, x=None, y=None, centering=False,
easting_shift=None, northing_shift=None,
revert_z=False):
"""
Initialize DEM (regular grid) interpolation object.
Required arguments
------------------
- demfile, 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 or irregular points
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
- revert_z = False, type bool
Revert z-values of the DEM if the original height information was
provided in a donward-oriented coorindate system
"""
self.x = None
self.y = None
self.rotate = None
self.revert_z = revert_z
if isinstance(demfile, str):
if demfile[-4:].lower() == '.asc':
self.loadASC(demfile, centering, easting_shift, northing_shift)
else:
self.loadTXT(demfile, centering, easting_shift, northing_shift)
elif x is not None and y is not None:
self.z = demfile
else:
raise Exception("Either DEM file or z with x and y must be given!")
def __call__(self, x, y=None, rotation=None):
"""
Interpolation function for either 2D or 3D coordinates.
Required arguments
------------------
- x, type iterable (list, ndarray)
list of target x-coordinates for 2D interpolation
Keyword arguments
-----------------
- y = None, type iterable (list, ndarray)
list of target y-coordinates for 3D interpolation
- rotation (RAD) = None, type float ()
rotation of the DEM around the center coordinate (set via the
*centering* flag or *easting/northing_shift*) relative to the
computational domain, given in rad, not degree
"""
if y is None:
if rotation is not None:
X = np.hstack((x.reshape(-1, 1), np.zeros((len(x), 2))))
x = mu.rotate(X, np.deg2rad(rotation))[:, 0]
return self.dem(x)
else:
if rotation is not None:
X = np.hstack((np.hstack((x.reshape(-1, 1), y.reshape(-1, 1))),
np.zeros((len(y), 1))))
X_rot = mu.rotate(X, np.deg2rad(rotation))
x, y = X_rot[:, 0], X_rot[:, 1]
if hasattr(self, 'tri'):
return self.dem(x, y)
else:
return self.dem((x, y))
[docs]
def loadTXT(self, demfile, centering, easting_shift, northing_shift):
"""
Load *txt* file with column-based list of DEM coordinates.
Required arguments
------------------
- see **DEM** class description
"""
xp, yp, zp = np.loadtxt(demfile, unpack=True)
if len(np.unique(xp)) * len(np.unique(yp)) > len(xp):
# Perform irregular grid interpolation based on triangulation
self.x = xp
self.y = yp
self.z = zp
self.translate(centering, easting_shift, northing_shift)
self.tri = mtri.Triangulation(self.x, self.y)
self.dem = mtri.LinearTriInterpolator(self.tri, self.z)
return
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)).T
if y[1] < y[0]:
y.sort()
zp = np.flipud(zp)
if x[1] < x[0]:
x.sort()
zp = np.fliplr(zp)
self.z = zp
if self.revert_z:
self.z *= -1.
self.x = x
self.y = y
print('DEM original extent (x min, x max, y min, y max):')
print(np.min(self.x), np.max(self.x), np.min(self.y), np.max(self.y))
self.translate(centering, easting_shift, northing_shift)
self.dem = RegularGridInterpolator((self.x, self.y), self.z.T)
[docs]
def loadASC(self, ascfile, centering, easting_shift, northing_shift):
"""
Load ASC file (DEM matrix with location header).
Required arguments
------------------
- see **DEM** class description
"""
with open(ascfile) 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
self.z = np.flipud(np.genfromtxt(ascfile, skip_header=nheader))
if 'NODATA_value' in header:
self.z[self.z == header['NODATA_value']] = np.nan
if self.revert_z:
self.z *= -1.
dx = header.pop('cellsize', 1.0) # just a guess
self.x = np.arange(header['ncols']) * dx + header['xllcorner']
self.y = np.arange(header['nrows']) * dx + header['yllcorner']
print('DEM original extent (x min, x max, y min, y max):')
print(np.min(self.x), np.max(self.x), np.min(self.y), np.max(self.y))
self.translate(centering, easting_shift, northing_shift)
self.dem = RegularGridInterpolator((self.x, self.y), self.z.T)
[docs]
def translate(self, centering=False, easting_shift=None,
northing_shift=None, back=False):
"""
Function to translate coordinates in horizontal direction.
Keyword arguments
-----------------
- see **DEM** class description
- back = False, type bool
set True to translate back in the opposite direction; useful for
rotating coordinates around a specific origin
"""
if centering:
if easting_shift is None and northing_shift is None:
print("shifted by:")
print("x-dir: ", np.mean([np.min(self.x), np.max(self.x)]))
print("y-dir: ", np.mean([np.min(self.y), np.max(self.y)]))
if not back:
self.x -= np.mean([np.min(self.x), np.max(self.x)])
self.y -= np.mean([np.min(self.y), np.max(self.y)])
else:
self.x += np.mean([np.min(self.x), np.max(self.x)])
self.y += np.mean([np.min(self.y), np.max(self.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:
self.x -= (np.mean(self.x) - easting_shift)
if northing_shift is not None:
self.y -= (np.mean(self.y) - northing_shift)
elif not centering:
if easting_shift is not None:
self.x += easting_shift
if northing_shift is not None:
self.y += northing_shift
[docs]
def show(self, cmap="terrain", cbar=True, ax=None, **kwargs):
"""
Show digital elevation model. Needs an update.
Keyword arguments
-----------------
- cmap = "terrain", type str ()
matplotlib colormap definiton
- cbar = True, type bool
add colorbar to the plot or not
- ax = None, type matplotlib figure axes object
add the plot to a given axes object or create a new one
- **kwargs, type keyword arguments
add additional keyword arguments for the plot style (e.g., *lw*)
"""
if ax is None:
fig, ax = plt.subplots(figsize=kwargs.pop('figsize', (15, 15)))
# extract some kwargs for axis setting and colorbar
orientation = kwargs.pop('orientation', 'vertical')
xlim = kwargs.pop('xlim', (-9e99, 9e99))
ylim = kwargs.pop('ylim', (-9e99, 9e99))
x = self.dem.grid[0]
y = self.dem.grid[1]
ix0 = np.argmin(x < xlim[0])
ix1 = np.argmax(x > xlim[1]) - 1
if ix1 < 0:
ix1 = len(x)
iy0 = np.argmin(y < ylim[0])
iy1 = np.argmax(y > ylim[1]) - 1
if iy1 < 0:
iy1 = len(y)
im = ax.pcolormesh(x[ix0:ix1], y[iy0:iy1],
self.dem.values[ix0:ix1, iy0:iy1].T,
cmap=cmap, **kwargs)
if cbar:
plt.colorbar(im, ax=ax, orientation=orientation)
ax.set_aspect(1.0)
return ax