From 560c778ca50ad13cb6dca802ebb66ae7503e80e5 Mon Sep 17 00:00:00 2001 From: Ben Thompson Date: Mon, 12 Apr 2021 14:21:02 -0400 Subject: [PATCH 1/3] EFCS <--> TDCS and projection transformations. --- .gitignore | 1 - cutde/__init__.py | 5 ++ cutde/geometry.py | 123 ++++++++++++++++++++++++++++++++++++++++++++++ environment.yml | 1 + 4 files changed, 129 insertions(+), 1 deletion(-) create mode 100644 cutde/geometry.py diff --git a/.gitignore b/.gitignore index bc7200e..12de820 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,5 @@ *.swp *.swo -*.png tags __pycache__ .pytest_cache diff --git a/cutde/__init__.py b/cutde/__init__.py index 6e372c7..61d3323 100644 --- a/cutde/__init__.py +++ b/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_rotations, +) diff --git a/cutde/geometry.py b/cutde/geometry.py new file mode 100644 index 0000000..69def51 --- /dev/null +++ b/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 = 0.1 + 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)) diff --git a/environment.yml b/environment.yml index e731632..83d09ee 100644 --- a/environment.yml +++ b/environment.yml @@ -14,3 +14,4 @@ dependencies: - pre-commit - ipdb - pip + - sphinx==3.5.4 From 6012edea66cce526e3ef6f4fb5df1c66d55e2d5e Mon Sep 17 00:00:00 2001 From: Ben Thompson Date: Mon, 12 Apr 2021 14:55:10 -0400 Subject: [PATCH 2/3] Add pyproj dependency, tests for geometry. --- cutde/__init__.py | 2 +- cutde/geometry.py | 2 +- environment.yml | 3 +-- tests/test_geometry.py | 37 +++++++++++++++++++++++++++++++++++++ 4 files changed, 40 insertions(+), 4 deletions(-) create mode 100644 tests/test_geometry.py diff --git a/cutde/__init__.py b/cutde/__init__.py index 61d3323..9ab6785 100644 --- a/cutde/__init__.py +++ b/cutde/__init__.py @@ -9,5 +9,5 @@ from .geometry import ( # noqa: F401 compute_efcs_to_tdcs_rotations, compute_normal_vectors, - compute_projection_rotations, + compute_projection_transforms, ) diff --git a/cutde/geometry.py b/cutde/geometry.py index 69def51..4a46ba2 100644 --- a/cutde/geometry.py +++ b/cutde/geometry.py @@ -58,7 +58,7 @@ def compute_projection_transforms( out = np.empty((origins.shape[0], 3, 3), dtype=origins.dtype) for d in range(3): - eps = 0.1 + eps = 1.0 targets = origins.copy() targets[:, d] += eps proj_origins = np.array( diff --git a/environment.yml b/environment.yml index 83d09ee..a885f36 100644 --- a/environment.yml +++ b/environment.yml @@ -5,7 +5,7 @@ dependencies: - numpy - scipy - matplotlib - - jinja2 + - pyproj - black - flake8 - isort @@ -14,4 +14,3 @@ dependencies: - pre-commit - ipdb - pip - - sphinx==3.5.4 diff --git a/tests/test_geometry.py b/tests/test_geometry.py new file mode 100644 index 0000000..f729038 --- /dev/null +++ b/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) From 0fe466f7315d5394fb0a4834d1176caeefe90722 Mon Sep 17 00:00:00 2001 From: Ben Thompson Date: Tue, 13 Apr 2021 11:47:17 -0400 Subject: [PATCH 3/3] Remove python 3.6/3.7 from github CI workflow. --- .github/workflows/test.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index dc785e8..fe5ceaa 100644 --- a/.github/workflows/test.yml +++ b/.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: