Merge pull request #5 from tbenthompson/projection_transform

EFCS <--> TDCS and projection transformations.

.github/workflows/test.yml

@@ -8,7 +8,7 @@ jobs:
       fail-fast: false
       matrix:
         os: ["ubuntu-latest", "macos-latest"]
-        python-version: ["3.6", "3.7", "3.8", "3.9"]
+        python-version: ["3.8", "3.9"]
     name: Test (${{ matrix.python-version }}, ${{ matrix.os }})
     runs-on:  ${{ matrix.os }}
     defaults:

.gitignore

@@ -1,6 +1,5 @@
 *.swp
 *.swo
-*.png
 tags
 __pycache__
 .pytest_cache

cutde/__init__.py

@@ -6,3 +6,8 @@
     strain_all_pairs,
     strain_to_stress,
 )
+from .geometry import (  # noqa: F401
+    compute_efcs_to_tdcs_rotations,
+    compute_normal_vectors,
+    compute_projection_transforms,
+)

cutde/geometry.py

@@ -0,0 +1,123 @@
+import numpy as np
+import pyproj
+
+
+def compute_normal_vectors(tri_pts) -> np.ndarray:
+    """
+    Compute normal vectors for each triangle.
+
+    Parameters
+    ----------
+    tri_pts : {array-like}, shape (n_triangles, 3, 3)
+        The vertices of the triangles.
+
+    Returns
+    -------
+    normals : np.ndarray, shape (n_triangles, 3)
+        The normal vectors for each triangle.
+    """
+    leg1 = tri_pts[:, 1] - tri_pts[:, 0]
+    leg2 = tri_pts[:, 2] - tri_pts[:, 0]
+    # The normal vector is one axis of the TDCS and can be
+    # computed as the cross product of the two corner-corner tangent vectors.
+    Vnormal = np.cross(leg1, leg2, axis=1)
+    # And it should be normalized to have unit length of course!
+    Vnormal /= np.linalg.norm(Vnormal, axis=1)[:, None]
+    return Vnormal
+
+
+def compute_projection_transforms(
+    origins, transformer: pyproj.Transformer
+) -> np.ndarray:
+    """
+    Convert vectors from one coordinate system to another. Unlike positions,
+    this cannot be done with a simple pyproj call. We first need to set up a
+    vector start and end point, convert those into the new coordinate system
+    and then recompute the direction/distance between the start and end point.
+
+    The output matrices are not pure rotation matrices because there is also
+    a scaling of vector lengths. For example, converting from latitude to
+    meters will result in a large scale factor.
+
+    You can obtain the inverse transformation either by computing the inverse
+    of the matrix or by passing an inverse pyproj.Transformer.
+
+    Parameters
+    ----------
+    origins : {array-like}, shape (N, 3)
+        The points at which we will compute rotation matrices
+    transformer : pyproj.Transformer
+        A pyproj.Transformer that will perform the necessary projection step.
+
+    Returns
+    -------
+    transform_mats : np.ndarray, shape (n_triangles, 3, 3)
+        The 3x3 rotation and scaling matrices that transform vectors from the
+        EFCS to TDCS.
+    """
+
+    out = np.empty((origins.shape[0], 3, 3), dtype=origins.dtype)
+    for d in range(3):
+        eps = 1.0
+        targets = origins.copy()
+        targets[:, d] += eps
+        proj_origins = np.array(
+            transformer.transform(origins[:, 0], origins[:, 1], origins[:, 2])
+        ).T.copy()
+        proj_targets = np.array(
+            transformer.transform(targets[:, 0], targets[:, 1], targets[:, 2])
+        ).T.copy()
+        out[:, :, d] = proj_targets - proj_origins
+    return out
+
+
+def compute_efcs_to_tdcs_rotations(tri_pts) -> np.ndarray:
+    """
+    Build rotation matrices that convert from an Earth-fixed coordinate system
+    (EFCS) to a triangular dislocation coordinate system (TDCS).
+
+    In the EFCS, the vectors will be directions/length in a map projection or
+    an elliptical coordinate system.
+
+    In the TDCS, the coordinates/vectors will be separated into:
+    `(along-strike-distance, along-dip-distance, tensile-distance)`
+
+    Note that in the Nikhoo and Walter 2015 and the Okada convention, the dip
+    vector points upwards. This is different from the standard geologic
+    convention where the dip vector points downwards.
+
+    It may be useful to extract normal, dip or strike vectors from the rotation
+    matrices that are returned by this function. The strike vectors are:
+    `rot_mats[:, 0, :]`, the dip vectors are `rot_mats[:, 1, :]` and the normal
+    vectors are `rot_mats[:, 2, :]`.
+
+    To transform from TDCS back to EFCS, we simply need the transpose of the
+    rotation matrices because the inverse of an orthogonal matrix is its
+    transpose. To get this you can run `np.transpose(rot_mats, (0, 2, 1))`.
+
+    Parameters
+    ----------
+    tri_pts : {array-like}, shape (n_triangles, 3, 3)
+        The vertices of the triangles.
+
+    Returns
+    -------
+    rot_mats : np.ndarray, shape (n_triangles, 3, 3)
+        The 3x3 rotation matrices that transform vectors from the EFCS to TDCS.
+    """
+    Vnormal = compute_normal_vectors(tri_pts)
+    eY = np.array([0, 1, 0])
+    eZ = np.array([0, 0, 1])
+    # The strike vector is defined as orthogonal to both the (0,0,1) vector and
+    # the normal.
+    Vstrike_raw = np.cross(eZ[None, :], Vnormal, axis=1)
+    Vstrike_length = np.linalg.norm(Vstrike_raw, axis=1)
+
+    # If eZ == Vnormal, we will get Vstrike = (0,0,0). In this case, just set
+    # Vstrike equal to (0,±1,0).
+    Vstrike = np.where(
+        Vstrike_length[:, None] == 0, eY[None, :] * Vnormal[:, 2, None], Vstrike_raw
+    )
+    Vstrike /= np.linalg.norm(Vstrike, axis=1)[:, None]
+    Vdip = np.cross(Vnormal, Vstrike, axis=1)
+    return np.transpose(np.array([Vstrike, Vdip, Vnormal]), (1, 0, 2))

environment.yml

@@ -5,7 +5,7 @@ dependencies:
   - numpy
   - scipy
   - matplotlib
-  - jinja2
+  - pyproj
   - black
   - flake8
   - isort

tests/test_geometry.py

@@ -0,0 +1,37 @@
+import numpy as np
+
+from cutde import compute_efcs_to_tdcs_rotations, compute_projection_transforms
+
+
+def test_efcs_to_tdcs_orthogonal():
+    tp = np.random.random_sample((1, 3, 3))
+    R = compute_efcs_to_tdcs_rotations(tp)
+    for d1 in range(3):
+        for d2 in range(d1, 3):
+            c = np.sum(R[:, d1, :] * R[:, d2, :], axis=1)
+            true = 1 if d1 == d2 else 0
+            np.testing.assert_almost_equal(c, true)
+
+
+def test_projection_transforms():
+    from pyproj import Transformer
+
+    transformer = Transformer.from_crs(
+        "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs",
+        "+proj=geocent +datum=WGS84 +units=m +no_defs",
+    )
+    rand_pt = np.random.random_sample(3)
+    pts = np.array([rand_pt, rand_pt + np.array([0, 0, -1000.0])])
+    transforms = compute_projection_transforms(pts, transformer)
+
+    # The vertical component should transform without change in scale.
+    np.testing.assert_almost_equal(np.linalg.norm(transforms[:, :, 2], axis=1), 1.0)
+
+    # The transformation should be identical for vertical motions
+    np.testing.assert_allclose(
+        transforms[0, :, [0, 1]], transforms[1, :, [0, 1]], rtol=1e-2
+    )
+
+    # Check length of degree of latitude and longitude. Approximate.
+    np.testing.assert_allclose(transforms[0, 1, 0], 1.11e5, rtol=5e-3)
+    np.testing.assert_allclose(transforms[0, 2, 1], 1.11e5, rtol=5e-3)