Source code for scarplet.WindowedTemplate

# -*- coding: utf-8
""" Class for windowed template matching over a spatial grid """

import numexpr
import numpy as np

from scipy.special import erf, erfinv

np.seterr(divide='ignore', invalid='ignore')


[docs]class WindowedTemplate(object): """Base class for windowed template function Attributes ---------- d : float Scale of windowed template function in data projection units alpha : float Orientation of windowed template function in radians c : float Curvature limit of template nx : int Number of columns in template array ny : int Number of rows in template array de : float Spacing of template grid cells in dat projection units Methods ------- get_coordinates(): Get arrays of coordinates for template grid points get_mask(): Get mask array giving curvature extent of template window get_window_limits(): Get mask array giving window extent """ def __init__(self): self.d = None self.alpha = None self.nx = None self.ny = None self.c = self.nx / 2. self.de = None
[docs] def get_coordinates(self): x = self.de * np.linspace(1, self.nx, num=self.nx) y = self.de * np.linspace(1, self.ny, num=self.ny) x = x - np.mean(x) y = y - np.mean(y) x, y = np.meshgrid(x, y) xr = x * np.cos(self.alpha) + y * np.sin(self.alpha) yr = -x * np.sin(self.alpha) + y * np.cos(self.alpha) return xr, yr
[docs] def get_mask(self): xr, yr = self.get_coordinates() mask = (abs(xr) < self.c) & (abs(yr) < self.d) return mask
[docs] def get_window_limits(self): x4 = self.d*np.cos(self.alpha - np.pi/2) y4 = self.d*np.sin(self.alpha - np.pi/2) x1 = self.d*np.cos(self.alpha) y1 = self.d*np.sin(self.alpha) an_y = abs((x4 - x1) + 2 * self.c * np.cos(self.alpha - np.pi/2)) an_x = abs((y1 - y4) + 2 * self.c * np.sin(self.alpha - np.pi/2)) x = self.de*np.linspace(1, self.nx, num=self.nx) y = self.de*np.linspace(1, self.ny, num=self.ny) x = x - np.mean(x) y = y - np.mean(y) X, Y = np.meshgrid(x, y) mask = ((X < (min(x) + an_x)) | (X > (max(x) - an_x)) | (Y < (min(y) + an_y)) | (Y > (max(y) - an_y))) return mask
[docs]class Scarp(WindowedTemplate): """Curvature template for vertical scarp Attributes ---------- d : float Scale of windowed template function in data projection units alpha : float Orientation of windowed template function in radians kt : float Morphologic age of template in m2 nx : int Number of columns in template array ny : int Number of rows in template array de : float Spacing of template grid cells in dat projection units Methods ------- template(): Returns array of windowed template function template_numexpr(): Returns array of windowed template function optimized using numexpr References ---------- Adapted from template derived in Hilley, G.E., DeLong, S., Prentice, C., Blisniuk, K. and Arrowsmith, J.R., 2010. Morphologic dating of fault scarps using airborne laser swath mapping (ALSM) data. Geophysical Research Letters, 37(4). https://dx.doi.org/10.1029/2009GL042044 Based on solutions to the diffusion equation published in Hanks, T.C., 2000. The age of scarplike landforms from diffusion‐equation analysis. Quaternary geochronology, 4, pp.313-338. and many references therein. """ def __init__(self, d, kt, alpha, nx, ny, de): """Constructor method for scarp template Attributes ---------- d : float Scale of windowed template function in data projection units kt : float Morphologic age of template in m2 alpha : float Orientation of windowed template function in radians nx : int Number of columns in template array ny : int Number of rows in template array de : float Spacing of template grid cells in dat projection units """ self.d = d self.kt = kt self.alpha = -alpha self.nx = nx self.ny = ny self.de = de frac = 0.9 self.c = abs(2 * np.sqrt(self.kt) * erfinv(frac))
[docs] def template(self): """Return template function Returns ------- W : numpy array Windowed template function """ x = self.de * np.linspace(1, self.nx, num=self.nx) y = self.de * np.linspace(1, self.ny, num=self.ny) x = x - np.mean(x) y = y - np.mean(y) x, y = np.meshgrid(x, y) xr = x * np.cos(self.alpha) + y * np.sin(self.alpha) yr = -x * np.sin(self.alpha) + y * np.cos(self.alpha) W = (-xr / (2. * self.kt ** (3 / 2.) * np.sqrt(np.pi))) \ * np.exp(-xr ** 2. / (4. * self.kt)) mask = self.get_mask() W = W * mask return W
[docs] def template_numexpr(self): """Return template function (uses numexpr where possible) Returns ------- W : numpy array Windowed template function """ alpha = self.alpha kt = self.kt c = self.c d = self.d x = self.de * np.linspace(1, self.nx, num=self.nx) y = self.de * np.linspace(1, self.ny, num=self.ny) x = x - np.mean(x) y = y - np.mean(y) x, y = np.meshgrid(x, y) xr = numexpr.evaluate("x * cos(alpha) + y * sin(alpha)") yr = numexpr.evaluate("-x * sin(alpha) + y * cos(alpha)") pi = np.pi W = numexpr.evaluate("(-xr / (2 * kt ** (3/2) * sqrt(pi))) * \ exp(-xr ** 2 / (4 * kt))") mask = numexpr.evaluate("(abs(xr) < c) & (abs(yr) < d)") W = numexpr.evaluate("W * mask") return W
[docs]class RightFacingUpperBreakScarp(Scarp): """Template for upper slope break of vertical scarp (right-facting) Overrides template function to correct facign direction Attributes ---------- d : float Scale of windowed template function in data projection units alpha : float Orientation of windowed template function in radians kt : float Morphologic age of template in m2 nx : int Number of columns in template array ny : int Number of rows in template array de : float Spacing of template grid cells in dat projection units Methods ------- get_error_mask(): Return mask array that masks the lower slope break of scarp template(): Returns array of windowed template function """
[docs] def template(self): """Return template function (uses numexpr where possible) Returns ------- W : numpy array Windowed template function """ W = super().template_numexpr() return -W
[docs] def get_err_mask(self): """Return mask array masking the lower half of scarp Returns ------- mask : numpy array Mask array for lower half of scarp """ xr, _ = self.get_coordinates() mask = xr <= 0 return mask
[docs]class LeftFacingUpperBreakScarp(Scarp): """Template for upper slope break of vertical scarp (left-facting) Attributes ---------- d : float Scale of windowed template function in data projection units alpha : float Orientation of windowed template function in radians kt : float Morphologic age of template in m2 nx : int Number of columns in template array ny : int Number of rows in template array de : float Spacing of template grid cells in dat projection units Methods ------- get_error_mask(): Return mask array that masks the lower slope break of scarp """
[docs] def get_err_mask(self): """Return mask array masking the lower half of scarp Returns ------- mask : numpy array Mask array for lower hald of scarp """ xr, _ = self.get_coordinates() mask = xr >= 0 return mask
[docs]class ShiftedTemplateMixin(WindowedTemplate): """Mix-in for template that is offset from the window center Overrides template function to shift template Attributes ---------- d : float Scale of windowed template function in data projection units alpha : float Orientation of windowed template function in radians kt : float Morphologic age of template in m2 nx : int Number of columns in template array ny : int Number of rows in template array de : float Spacing of template grid cells in dat projection units dx : float X Offset of template center in data projection units dy : float Y Offset of template center data projection units Methods ------- set_offset(dx, dy): Set offset attrivutes t odx and dy shift_template(W, dx, dy): Shift template array W by dx and dy template(): Returns array of windowed template function """ def __init__(self, *args, **kwargs): """Constructor for shifted template Parameters ---------- dx : float X Offset of template center in data projection units dy : float Y Offset of template center data projection units """ super().__init__(*args) self.set_offset(kwargs['dx'], kwargs['dy'])
[docs] def set_offset(self, dx, dy): """Set offset values Parameters ---------- dx : float X Offset of template center in data projection units dy : float Y Offset of template center data projection units """ self.dx = dx self.dy = dy
[docs] def shift_template(self, W, dx, dy): """Shift template Parameters ---------- W : numpy array Windowed template function dx : float X Offset of template center in data projection units dy : float Y Offset of template center data projection units Returns ------- W : numpy array Shifted windowed template function """ ny, nx = W.shape if dx > 0: left = np.zeros((ny, dx)) W = W[:, 0:-dx] W = np.hstack([left, W]) else: dx = abs(dx) right = np.zeros((ny, dx)) W = W[:, dx:] W = np.hstack([W, right]) if dy > 0: bottom = np.zeros((dy, nx)) W = W[0:-dy, :] W = np.vstack([W, bottom]) else: dy = abs(dy) top = np.zeros((dy, nx)) W = W[dy:, :] W = np.vstack([top, W]) return W
[docs] def template(self): """Template function for shifted template Returns ------- W : numpy array Shifted windowed template function """ W = super().template() W = self.shift_template(W, self.dx, self.dy) return W
[docs]class ShiftedLeftFacingUpperBreakScarp(ShiftedTemplateMixin, LeftFacingUpperBreakScarp): pass
[docs]class ShiftedRightFacingUpperBreakScarp(ShiftedTemplateMixin, RightFacingUpperBreakScarp): pass
[docs]class Ricker(WindowedTemplate): """Template using 2D Ricker wavelet Attributes ---------- d : float Scale of windowed template function in data projection units alpha : float Orientation of windowed template function in radians kt : float Morphologic age of template in m2 nx : int Number of columns in template array ny : int Number of rows in template array de : float Spacing of template grid cells in dat projection units Methods ------- template(): Returns array of windowed template function References ---------- This implements a Ricker wavelet similar to thet used in the following work Lashermes, B., Foufoula‐Georgiou, E., and Dietrich, W. E., 2007, Channel network extraction from high resolution topography using wavelets. Geophysical Research Letters, 34(23). https://doi.org/10.1029/2007GL031140 """ def __init__(self, d, f, alpha, nx, ny, de): """Constructor method for Ricker template Paramters --------- d : float Scale of windowed template function in data projection units kt : float Morphologic age of template in m2 alpha : float Orientation of windowed template function in radians nx : int Number of columns in template array ny : int Number of rows in template array de : float Spacing of template grid cells in dat projection units """ self.d = d self.f = f self.alpha = -alpha self.nx = nx self.ny = ny self.c = nx self.de = de
[docs] def get_window_limits(self): return np.zeros((self.ny, self.nx), dtype=bool)
[docs] def template(self): """Template function for windowed Ricker wavelet Returns ------- W : numpy array Windowed template function """ alpha = self.alpha d = self.d f = self.f xr, yr = self.get_coordinates() pi = np.pi W = numexpr.evaluate("(1. - 2. * (pi * f * xr) ** 2.) * \ exp(-(pi * f * xr) ** 2.)") mask = self.get_mask() W = W * mask return W
[docs]class Channel(Ricker): """Duplicate class for Ricker wavelet used for fluvial channels""" pass
[docs]class Crater(WindowedTemplate): """Template for radially symmetric crater Attributes ---------- r : float Radius of crater in pixels kt : float Morphologic age of template crater rim in m2 nx : int Number of columns in template array ny : int Number of rows in template array de : float Spacing of template grid cells in dat projection units """ def __init__(self, r, kt, nx, ny, de): """Constructor methodfor radially symmetric crater Parameters ---------- r : float Radius of crater in pixels kt : float Morphologic age of template crater rim in m2 nx : int Number of columns in template array ny : int Number of rows in template array de : float Spacing of template grid cells in dat projection units """ self.r = r / de self.kt = kt self.nx = nx self.ny = ny self.de = de
[docs] def template(self): """Template function for radially symmetric crater Returns ------- W : numpy array Windowed template function """ x = self.de * np.linspace(1, self.nx, num=self.nx) y = self.de * np.linspace(1, self.ny, num=self.ny) x = x - np.mean(x) y = y - np.mean(y) x, y = np.meshgrid(x, y) W = np.zeros_like(x) thetas = np.linspace(0, 2 * np.pi, num=359, endpoint=False) for theta in thetas: alpha = -theta dx = self.r * np.cos(theta) dy = self.r * np.sin(theta) xr = (x - dx) * np.cos(alpha) + (y + dy) * np.sin(alpha) yr = -(x - dx) * np.sin(alpha) + (y + dy) * np.cos(alpha) this_W = (-xr / (2. * self.kt ** (3 / 2.) * np.sqrt(np.pi))) \ * np.exp(-xr ** 2. / (4. * self.kt)) mask = (abs(xr) < 1) & (abs(yr) < 5 / self.de) this_W *= mask if theta > np.pi / 2 and theta < 3 * np.pi / 2: this_W *= -1 W += this_W return W