Nshm organisational refactor

nshmdb/fault.py

@@ -20,8 +20,7 @@
 from enum import Enum
 
 import numpy as np
-import qcore.coordinates
-import qcore.geo
+from qcore import coordinates, geo
 
 
 class TectType(Enum):
@@ -50,7 +49,7 @@ class FaultPlane:
         The rake angle of the fault plane.
     """
 
-    _corners: np.ndarray
+    corners_nztm: np.ndarray
     rake: float
 
     @property
@@ -62,7 +61,7 @@ def corners(self) -> np.ndarray:
         np.ndarray
             The corners of the fault plane in (lat, lon, depth) format.
         """
-        return qcore.coordinates.nztm_to_wgs_depth(self._corners)
+        return qcore.coordinates.nztm_to_wgs_depth(self.corners_nztm)
 
     @property
     def length_m(self) -> float:
@@ -72,7 +71,7 @@ def length_m(self) -> float:
         float
             The length of the fault plane (in metres).
         """
-        return np.linalg.norm(self._corners[1] - self._corners[0])
+        return np.linalg.norm(self.corners_nztm[1] - self.corners_nztm[0])
 
     @property
     def width_m(self) -> float:
@@ -82,7 +81,7 @@ def width_m(self) -> float:
         float
             The width of the fault plane (in metres).
         """
-        return np.linalg.norm(self._corners[-1] - self._corners[0])
+        return np.linalg.norm(self.corners_nztm[-1] - self.corners_nztm[0])
 
     @property
     def bottom_m(self) -> float:
@@ -92,7 +91,7 @@ def bottom_m(self) -> float:
         float
             The bottom depth (in metres).
         """
-        return self._corners[-1, -1]
+        return self.corners_nztm[-1, -1]
 
     @property
     def width(self) -> float:
@@ -132,7 +131,7 @@ def projected_width(self) -> float:
         float
             The projected width of the fault plane (in kilometres).
         """
-        return self.projected_width / _KM_TO_M
+        return self.projected_width_m / _KM_TO_M
 
     @property
     def strike(self) -> float:
@@ -146,7 +145,7 @@ def strike(self) -> float:
 
         north_direction = np.array([1, 0, 0])
         up_direction = np.array([0, 0, 1])
-        strike_direction = self._corners[1] - self._corners[0]
+        strike_direction = self.corners_nztm[1] - self.corners_nztm[0]
         return qcore.geo.oriented_bearing_wrt_normal(
             north_direction, strike_direction, up_direction
         )
@@ -163,7 +162,7 @@ def dip_dir(self) -> float:
             return 0  # TODO: Is this right for this case?
         north_direction = np.array([1, 0, 0])
         up_direction = np.array([0, 0, 1])
-        dip_direction = self._corners[-1] - self._corners[0]
+        dip_direction = self.corners_nztm[-1] - self.corners_nztm[0]
         dip_direction[-1] = 0
         return qcore.geo.oriented_bearing_wrt_normal(
             north_direction, dip_direction, up_direction
@@ -189,8 +188,8 @@ def plane_coordinates_to_global_coordinates(
         plane_coordinates : np.ndarray
             Plane coordinates to convert. Plane coordinates are
             2D coordinates (x, y) given for a fault plane (a plane), where x
-            represents displacement along the length of the fault, and y
-            displacement along the width of the fault (see diagram below). The
+            represents displacement along the strike, and y
+            displacement along the dip (see diagram below). The
             origin for plane coordinates is the centre of the fault.
 
                           +x
@@ -209,9 +208,9 @@ def plane_coordinates_to_global_coordinates(
         np.ndarray
             An 3d-vector of (lat, lon, depth) transformed coordinates.
         """
-        origin = self._corners[0]
-        top_right = self._corners[1]
-        bottom_left = self._corners[-1]
+        origin = self.corners_nztm[0]
+        top_right = self.corners_nztm[1]
+        bottom_left = self.corners_nztm[-1]
         frame = np.vstack((top_right - origin, bottom_left - origin))
         offset = np.array([1 / 2, 1 / 2])
 
@@ -244,9 +243,9 @@ def global_coordinates_to_plane_coordinates(
         ValueError
             If the given coordinates do not lie in the fault plane.
         """
-        origin = self._corners[0]
-        top_right = self._corners[1]
-        bottom_left = self._corners[-1]
+        origin = self.corners_nztm[0]
+        top_right = self.corners_nztm[1]
+        bottom_left = self.corners_nztm[-1]
         frame = np.vstack((top_right - origin, bottom_left - origin))
         offset = qcore.coordinates.wgs_depth_to_nztm(global_coordinates) - origin
         plane_coordinates, residual, _, _ = np.linalg.lstsq(frame.T, offset, rcond=None)
@@ -294,7 +293,7 @@ def centroid(self) -> np.ndarray:
         """
 
         return qcore.coordinates.nztm_to_wgs_depth(
-            np.mean(self._corners, axis=0).reshape((1, -1))
+            np.mean(self.corners_nztm, axis=0).reshape((1, -1))
         ).ravel()
 
 
@@ -374,10 +373,21 @@ def corners(self) -> np.ndarray:
         Returns
         -------
         np.ndarray of shape (4n x 3)
-            The corners of each fault plane in the fault, stacked vertically.
+            The corners in (lat, lon, depth) format of each fault plane in the
+            fault, stacked vertically.
         """
         return np.vstack([plane.corners for plane in self.planes])
 
+    def corners_nztm(self) -> np.ndarray:
+        """Get all corners of a fault.
+
+        Returns
+        -------
+        np.ndarray of shape (4n x 3)
+            The corners in NZTM format of each fault plane in the fault, stacked vertically.
+        """
+        return np.vstack([plane.corners_nztm for plane in self.planes])
+
     def global_coordinates_to_fault_coordinates(
         self, global_coordinates: np.ndarray
     ) -> np.ndarray:

nshmdb/plotting/rupture.py

@@ -7,12 +7,14 @@
     plot_rupture: Plot ruptures on faults.
 """
 
+from pathlib import Path
+
 import numpy as np
 from nshmdb.fault import Fault
 from pygmt_helper import plotting
 
 
-def plot_rupture(title: str, faults: list[Fault]):
+def plot_rupture(title: str, faults: list[Fault], output_filepath: Path):
     """Plot faults involved in a rupture scenario using pygmt.
 
     Parameters
@@ -21,6 +23,8 @@ def plot_rupture(title: str, faults: list[Fault]):
         The title of the figure.
     faults : list[Fault]
         The list of faults involved in the rupture.
+    output_filepath : Path
+        The filepath to output the figure to.
     """
     corners = np.vstack([fault.corners() for fault in faults])
     region = (
@@ -41,4 +45,4 @@ def plot_rupture(title: str, faults: list[Fault]):
                 fill="red",
             )
 
-    fig.show()
+    fig.savefig(output_filepath)

nshmdb/scripts/nshm_db_generator.py

@@ -40,7 +40,6 @@
 app = typer.Typer()
 
 
-POLYGONS_PATH = Path("ruptures") / "sect_polygons.geojson"
 FAULT_INFORMATION_PATH = Path("ruptures") / "fault_sections.geojson"
 RUPTURE_FAULT_JOIN_PATH = Path("ruptures") / "fast_indices.csv"
 

requirements.txt

@@ -0,0 +1,10 @@
+dataclasses
+geojson
+numpy
+pandas
+pygments
+pygmt
+qcore
+scipy
+tqdm
+typer