Add surface-tangent motion model (Brinkerhoff 2017)

Includes both cartesian and cylindrical versions.

src/glimpse/track/__init__.py

@@ -1,8 +1,21 @@
 """Track image features through time using particle filtering."""
 
-from .motion import CartesianMotion, CylindricalMotion
+from .motion import (
+    CartesianMotion,
+    CylindricalMotion,
+    TangentCartesianMotion,
+    TangentCylindricalMotion,
+)
 from .observer import Observer
 from .tracker import Tracker
 from .tracks import Tracks
 
-__all__ = ["CartesianMotion", "CylindricalMotion", "Observer", "Tracker", "Tracks"]
+__all__ = [
+    "CartesianMotion",
+    "CylindricalMotion",
+    "TangentCartesianMotion",
+    "TangentCylindricalMotion",
+    "Observer",
+    "Tracker",
+    "Tracks",
+]

src/glimpse/track/motion.py

@@ -309,3 +309,214 @@ def evolve_particles(self, particles: np.ndarray, dt: datetime.timedelta) -> Non
             time_units * particles[:, 3:6] + 0.5 * axyz * time_units ** 2
         )
         particles[:, 3:6] += time_units * axyz
+
+
+class TangentCartesianMotion(Motion):
+    """
+    Evolves particles tangent to a mean surface.
+
+    Initial particle positions and velocities, and random accelerations, are
+    specified by independent and normally distributed x and y components.
+    Temporal arguments (e.g. :attr:`vxy`, :attr:`axy`) are assumed to be
+    in :attr:`time_unit` time units.
+
+    Particle heights (z) are initialized based on a mean surface (:attr:`dem`)
+    and its uncertainty (:attr:`dem_sigma`), then maintain the same relative
+    height adjusted by a random walk proportional to the particle's horizontal
+    distance and a characteristic slope of small-scale features
+    (:attr:`slope_sigma`).
+
+    This is the model used in Brinkerhoff (2017): Bayesian methods in glaciology
+    (http://hdl.handle.net/11122/8113), chapter 4.
+
+    Attributes:
+        xy (iterable): Mean initial position (x, y).
+        time_unit (timedelta): Length of time unit for temporal arguments.
+        dem: Elevation of the surface on which to track points.
+        dem_sigma: Elevation standard deviations.
+        n (int): Number of particles.
+        xy_sigma (iterable): Standard deviation of initial position (x, y).
+        vxy (iterable): Mean initial velocity (dx/dt, dy/dt).
+        vxy_sigma (iterable): Standard deviation of initial velocity
+            (dx/dt, dy/dt).
+        axy (iterable): Mean acceleration (d^2x/dt^2, d^2y/dt^2).
+        axy_sigma (iterable): Standard deviation of acceleration
+            (d^2x/dt^2, d^2y/dt^2).
+        slope_sigma (float): Standard deviation of the characteristic slope
+            of small-scale features.
+    """
+
+    def __init__(
+        self,
+        xy: Iterable[Number],
+        time_unit: datetime.timedelta,
+        dem: Union[Number, Raster],
+        dem_sigma: Union[Number, Raster] = 0,
+        n: int = 1000,
+        xy_sigma: Iterable[Number] = (0, 0),
+        vxy: Iterable[Number] = (0, 0),
+        vxy_sigma: Iterable[Number] = (0, 0),
+        axy: Iterable[Number] = (0, 0),
+        axy_sigma: Iterable[Number] = (0, 0),
+        slope_sigma: Number = 0,
+    ) -> None:
+        self.xy = xy
+        self.time_unit = time_unit
+        if not isinstance(dem, Raster):
+            dem = Raster(dem)
+        self.dem = dem
+        if not isinstance(dem_sigma, Raster):
+            dem_sigma = Raster(dem_sigma)
+        self.dem_sigma = dem_sigma
+        self.n = n
+        self.xy_sigma = xy_sigma
+        self.vxy = vxy
+        self.vxy_sigma = vxy_sigma
+        self.axy = axy
+        self.axy_sigma = axy_sigma
+        self.slope_sigma = slope_sigma
+
+    def initialize_particles(self) -> np.ndarray:
+        """
+        Initialize particles around an initial mean position.
+
+        Returns:
+            particles: Particle positions and velocities (x, y, z, vx, vy, vz).
+        """
+        particles = np.zeros((self.n, 6), dtype=float)
+        particles[:, 0:2] = self.xy + self.xy_sigma * np.random.randn(self.n, 2)
+        z_offsets = self.dem_sigma.sample(particles[:, 0:2]) * np.random.randn(self.n)
+        particles[:, 2] = self.dem.sample(particles[:, 0:2]) + z_offsets
+        particles[:, 3:5] = self.vxy + self.vxy_sigma * np.random.randn(self.n, 2)
+        return particles
+
+    def evolve_particles(self, particles: np.ndarray, dt: datetime.timedelta) -> None:
+        """
+        Evolve particles through time by stochastic differentiation.
+
+        Arguments:
+            particles: Particle positions and velocities (x, y, z, vx, vy, vz).
+            dt: Time step to evolve particles forward or backward.
+        """
+        n = len(particles)
+        time_units = dt.total_seconds() / self.time_unit.total_seconds()
+        axy = self.axy + self.axy_sigma * np.random.randn(n, 2)
+        dxy = time_units * particles[:, 3:5] + 0.5 * axy * time_units ** 2
+        # HACK: Recover z_offsets (since particles may have been resampled)
+        z_offsets = particles[:, 2] - self.dem.sample(particles[:, 0:2])
+        z_offsets += (
+            self.slope_sigma * np.random.randn(n) * (dxy ** 2).sum(axis=1) ** 0.5
+        )
+        particles[:, 0:2] += dxy
+        particles[:, 2] = self.dem.sample(particles[:, 0:2]) + z_offsets
+        particles[:, 3:5] += time_units * axy
+
+
+class TangentCylindricalMotion(Motion):
+    """
+    Evolves particles tangent to a mean surface using cylindrical motion.
+
+    Identical to :class:`TangentCartesianMotion`, except that particle motion is
+    specified by independent and normally distributed components of magnitude
+    (radius) and direction (theta). Angular arguments are
+    assumed to be in radians counterclockwise from the +x axis.
+
+    Attributes:
+        xy (iterable): Mean initial position (x, y).
+        time_unit (timedelta): Length of time unit for temporal arguments.
+        dem: Elevation of the surface on which to track points.
+        dem_sigma: Elevation standard deviations.
+        n (int): Number of particles.
+        xy_sigma (iterable): Standard deviation of initial position (x, y).
+        vrth (iterable): Mean initial velocity (d radius/dt, theta).
+        vrth_sigma (iterable): Standard deviation of initial velocity
+            (d radius/dt, theta).
+        arth (iterable): Mean acceleration
+            (d^2 radius/dt^2, d theta/dt).
+        arth_sigma (iterable): Standard deviation of acceleration
+            (d^2 radius/dt^2, d theta/dt).
+    """
+
+    def __init__(
+        self,
+        xy: Iterable[Number],
+        time_unit: datetime.timedelta,
+        dem: Union[Number, Raster],
+        dem_sigma: Union[Number, Raster] = None,
+        n: int = 1000,
+        xy_sigma: Iterable[Number] = (0, 0),
+        vrth: Iterable[Number] = (0, 0),
+        vrth_sigma: Iterable[Number] = (0, 0),
+        arth: Iterable[Number] = (0, 0),
+        arth_sigma: Iterable[Number] = (0, 0),
+        slope_sigma: Number = 0,
+    ) -> None:
+        self.xy = xy
+        self.time_unit = time_unit
+        if not isinstance(dem, Raster):
+            dem = Raster(dem)
+        self.dem = dem
+        if not isinstance(dem_sigma, Raster):
+            dem_sigma = Raster(dem_sigma)
+        self.dem_sigma = dem_sigma
+        self.n = n
+        self.xy_sigma = xy_sigma
+        self.vrth = vrth
+        self.vrth_sigma = vrth_sigma
+        self.arth = arth
+        self.arth_sigma = arth_sigma
+        self.slope_sigma = slope_sigma
+
+    def initialize_particles(self) -> np.ndarray:
+        """
+        Initialize particles around an initial mean position.
+
+        Returns:
+            particles: Particle positions and velocities (x, y, z, vx, vy, vz).
+        """
+        particles = np.zeros((self.n, 6), dtype=float)
+        particles[:, 0:2] = self.xy + self.xy_sigma * np.random.randn(self.n, 2)
+        z_offsets = self.dem_sigma.sample(particles[:, 0:2]) * np.random.randn(self.n)
+        particles[:, 2] = self.dem.sample(particles[:, 0:2]) + z_offsets
+        vrth = self.vrth + self.vrth_sigma * np.random.randn(self.n, 2)
+        particles[:, 3:5] = np.column_stack(
+            (
+                # r' * cos(th)
+                vrth[:, 0] * np.cos(vrth[:, 1]),
+                # r' * sin(th)
+                vrth[:, 0] * np.sin(vrth[:, 1]),
+            )
+        )
+        return particles
+
+    def evolve_particles(self, particles: np.ndarray, dt: datetime.timedelta) -> None:
+        """
+        Evolve particles through time by stochastic differentiation.
+
+        Arguments:
+            particles: Particle positions and velocities (x, y, z, vx, vy, vz).
+            dt: Time step to evolve particles forward or backward.
+        """
+        n = len(particles)
+        time_units = dt.total_seconds() / self.time_unit.total_seconds()
+        vx = particles[:, 3]
+        vy = particles[:, 4]
+        vr = np.sqrt(vx ** 2 + vy ** 2)
+        arth = self.arth + self.arth_sigma * np.random.randn(n, 2)
+        axy = np.column_stack(
+            (
+                # r'' * cos(th) - r' * sin(th) * th'
+                arth[:, 0] * (vx / vr) - vy * arth[:, 1],
+                # r'' * sin(th) - r' * cos(th) * th'
+                arth[:, 0] * (vy / vr) + vx * arth[:, 1],
+            )
+        )
+        dxy = time_units * particles[:, 3:5] + 0.5 * axy * time_units ** 2
+        # HACK: Recover z_offsets (since particles may have been resampled)
+        z_offsets = particles[:, 2] - self.dem.sample(particles[:, 0:2])
+        z_offsets += (
+            self.slope_sigma * np.random.randn(n) * (dxy ** 2).sum(axis=1) ** 0.5
+        )
+        particles[:, 0:2] += dxy
+        particles[:, 2] = self.dem.sample(particles[:, 0:2]) + z_offsets
+        particles[:, 3:5] += time_units * axy