diff --git a/nshmdb/fault.py b/nshmdb/fault.py deleted file mode 100644 index f30db2a..0000000 --- a/nshmdb/fault.py +++ /dev/null @@ -1,507 +0,0 @@ -"""Module for representing fault planes and faults. - -This module provides classes and functions for representing fault planes and -faults, along with methods for calculating various properties such as -dimensions, orientation, and coordinate transformations. - -Classes -------- -TectType: - An enumeration of all the different kinds of fault types. - -FaultPlane: - A representation of a single plane of a Fault. - -Fault: - A representation of a fault, consisting of one or more FaultPlanes. -""" - -import dataclasses -from enum import Enum -from typing import Optional - -import numpy as np -from qcore import coordinates, geo - - -class TectType(Enum): - """An enumeration of all the different kinds of fault types.""" - - ACTIVE_SHALLOW = 1 - VOLCANIC = 2 - SUBDUCTION_INTERFACE = 3 - SUBDUCTION_SLAB = 4 - - -_KM_TO_M = 1000 - - -@dataclasses.dataclass -class FaultPlane: - """A representation of a single plane of a Fault. - - This class represents a single plane of a fault, providing various - properties and methods for calculating its dimensions, orientation, and - converting coordinates between different reference frames. - - Attributes - ---------- - rake : float - The rake angle of the fault plane. - corners_nztm : np.ndarray - The corners of the fault plane, in NZTM format. The order of the - corners is given clockwise from the top left (according to strike - and dip). See the diagram below. - - 0 1 - ┌──────────┐ - │ │ - │ │ - │ │ - │ │ - │ │ - │ │ - │ │ - └──────────┘ - 3 2 - """ - - corners_nztm: np.ndarray - rake: float - - @property - def corners(self) -> np.ndarray: - """ - Returns - ------- - np.ndarray - The corners of the fault plane in (lat, lon, depth) format. The - corners are the same as in corners_nztm. - """ - return coordinates.nztm_to_wgs_depth(self.corners_nztm) - - @property - def length_m(self) -> float: - """ - Returns - ------- - float - The length of the fault plane (in metres). - """ - return np.linalg.norm(self.corners_nztm[1] - self.corners_nztm[0]) - - @property - def width_m(self) -> float: - """ - Returns - ------- - float - The width of the fault plane (in metres). - """ - return np.linalg.norm(self.corners_nztm[-1] - self.corners_nztm[0]) - - @property - def bottom_m(self) -> float: - """ - Returns - ------- - float - The bottom depth (in metres). - """ - return self.corners_nztm[-1, -1] - - @property - def top_m(self) -> float: - return self.corners_nztm[0, -1] - - @property - def width(self) -> float: - """ - Returns - ------- - float - The width of the fault plane (in kilometres). - """ - return self.width_m / _KM_TO_M - - @property - def length(self) -> float: - """ - Returns - ------- - float - The length of the fault plane (in kilometres). - """ - return self.length_m / _KM_TO_M - - @property - def projected_width_m(self) -> float: - """ - Returns - ------- - float - The projected width of the fault plane (in metres). - """ - return self.length_m * np.cos(np.radians(self.dip)) - - @property - def projected_width(self) -> float: - """ - Returns - ------- - float - The projected width of the fault plane (in kilometres). - """ - return self.projected_width_m / _KM_TO_M - - @property - def strike(self) -> float: - """ - Returns - ------- - float - The bearing of the strike direction of the fault - (from north; in degrees) - """ - - north_direction = np.array([1, 0, 0]) - up_direction = np.array([0, 0, 1]) - strike_direction = self.corners_nztm[1] - self.corners_nztm[0] - return geo.oriented_bearing_wrt_normal( - north_direction, strike_direction, up_direction - ) - - @property - def dip_dir(self) -> float: - """ - Returns - ------- - float - The bearing of the dip direction (from north; in degrees). - """ - if np.isclose(self.dip, 90): - 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_nztm[-1] - self.corners_nztm[0] - dip_direction[-1] = 0 - return geo.oriented_bearing_wrt_normal( - north_direction, dip_direction, up_direction - ) - - @property - def dip(self) -> float: - """ - Returns - ------- - float - The dip angle of the fault. - """ - return np.degrees(np.arcsin(np.abs(self.bottom_m - self.top_m) / self.width_m)) - - def plane_coordinates_to_global_coordinates( - self, plane_coordinates: np.ndarray - ) -> np.ndarray: - """Convert plane coordinates to nztm global coordinates. - - Parameters - ---------- - 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 strike, and y - displacement along the dip (see diagram below). The - origin for plane coordinates is the centre of the fault. - - +x - -1/2,-1/2 ─────────────────> - ┌─────────────────────┐ │ - │ < strike > │ │ - │ ^ │ │ - │ dip │ │ +y - │ v │ │ - │ │ │ - └─────────────────────┘ ∨ - 1/2,1/2 - - Returns - ------- - np.ndarray - An 3d-vector of (lat, lon, depth) transformed coordinates. - """ - 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]) - - return coordinates.nztm_to_wgs_depth( - origin + (plane_coordinates + offset) @ frame - ) - - def global_coordinates_to_plane_coordinates( - self, - global_coordinates: np.ndarray, - ) -> np.ndarray: - """Convert coordinates (lat, lon, depth) to plane coordinates (x, y). - - See plane_coordinates_to_global_coordinates for a description of - plane coordinates. - - Parameters - ---------- - global_coordinates : np.ndarray - Global coordinates to convert. - - Returns - ------- - np.ndarray - The plane coordinates (x, y) representing the position of - global_coordinates on the fault plane. - - Raises - ------ - ValueError - If the given coordinates do not lie in the fault plane. - """ - 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 = coordinates.wgs_depth_to_nztm(global_coordinates) - origin - plane_coordinates, residual, _, _ = np.linalg.lstsq(frame.T, offset, rcond=None) - if not np.isclose(residual[0], 0, atol=1e-02): - raise ValueError("Coordinates do not lie in fault plane.") - return np.clip(plane_coordinates - np.array([1 / 2, 1 / 2]), -1 / 2, 1 / 2) - - def global_coordinates_in_plane(self, global_coordinates: np.ndarray) -> bool: - """Test if some global coordinates lie in the bounds of a plane. - - Parameters - ---------- - global_coordinates : np.ndarray - The global coordinates to check - - Returns - ------- - bool - True if the given global coordinates (lat, lon, depth) lie on the - fault plane. - """ - - try: - plane_coordinates = self.global_coordinates_to_plane_coordinates( - global_coordinates - ) - return np.all( - np.logical_or( - np.abs(plane_coordinates) < 1 / 2, - np.isclose(np.abs(plane_coordinates), 1 / 2, atol=1e-3), - ) - ) - except ValueError: - return False - - def centroid(self) -> np.ndarray: - """Returns the centre of the fault plane. - - Returns - ------- - np.ndarray - A 1 x 3 dimensional vector representing the centroid of the fault - plane in (lat, lon, depth) format. - - """ - - return coordinates.nztm_to_wgs_depth( - np.mean(self.corners_nztm, axis=0).reshape((1, -1)) - ).ravel() - - -@dataclasses.dataclass -class Fault: - """A representation of a fault, consisting of one or more FaultPlanes. - - This class represents a fault, which is composed of one or more FaultPlanes. - It provides methods for computing the area of the fault, getting the widths and - lengths of all fault planes, retrieving all corners of the fault, converting - global coordinates to fault coordinates, converting fault coordinates to global - coordinates, generating a random hypocentre location within the fault, and - computing the expected fault coordinates. - - Attributes - ---------- - name : str - The name of the fault. - tect_type : TectType | None - The type of fault this is (e.g. crustal, volcanic, subduction). - planes : list[FaultPlane] - A list containing all the FaultPlanes that constitute the fault. - - Methods - ------- - area: - Compute the area of a fault. - widths: - Get the widths of all fault planes. - lengths: - Get the lengths of all fault planes. - corners: - Get all corners of a fault. - global_coordinates_to_fault_coordinates: - Convert global coordinates to fault coordinates. - fault_coordinates_to_wgsdepth_coordinates: - Convert fault coordinates to global coordinates. - """ - - name: str - tect_type: Optional[TectType] - planes: list[FaultPlane] - - def area(self) -> float: - """Compute the area of a fault. - - Returns - ------- - float - The area of the fault. - """ - return sum(plane.width * plane.length for plane in self.planes) - - def widths(self) -> np.ndarray: - """Get the widths of all fault planes. - - Returns - ------- - np.ndarray of shape (1 x n) - The widths of all fault planes contained in this fault. - """ - return np.array([seg.width for seg in self.planes]) - - def lengths(self) -> np.ndarray: - """Get the lengths of all fault planes. - - Returns - ------- - np.ndarray of shape (1 x n) - The lengths of all fault planes contained in this fault. - """ - return np.array([seg.length for seg in self.planes]) - - def corners(self) -> np.ndarray: - """Get all corners of a fault. - - Returns - ------- - np.ndarray of shape (4n x 3) - 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: - """Convert global coordinates in (lat, lon, depth) format to fault coordinates. - - Fault coordinates are a tuple (s, d) where s is the distance (in - kilometres) from the top centre, and d the distance from the top of the - fault (refer to the diagram). - - ┌─────────┬──────────────┬────┐ - │ │ ╎ │ │ - │ │ ╎ │ │ - │ │ d ╎ │ │ - │ │ ╎ │ │ - │ │ └╶╶╶╶╶╶╶╶╶╶+ │ - │ │ s │ ∧ │ - │ │ │ │ │ - │ │ │ │ │ - └─────────┴──────────────┴──┼─┘ - │ - point: (s, d) - - Parameters - ---------- - global_coordinates : np.ndarray of shape (1 x 3) - The global coordinates to convert. - - Returns - ------- - np.ndarray - The fault coordinates. - - Raises - ------ - ValueError - If the given point does not lie on the fault. - """ - - running_length = 0.0 - midpoint = np.sum(self.lengths()) / 2 - for plane in self.planes: - if plane.global_coordinates_in_plane(global_coordinates): - plane_coordinates = plane.global_coordinates_to_plane_coordinates( - global_coordinates - ) - strike_length = plane_coordinates[0] + 1 / 2 - dip_length = plane_coordinates[1] + 1 / 2 - return np.array( - [ - running_length + strike_length * plane.length - midpoint, - max(dip_length * plane.width, 0), - ] - ) - running_length += plane.length - raise ValueError("Specified coordinates not contained on fault.") - - def fault_coordinates_to_wgsdepth_coordinates( - self, fault_coordinates: np.ndarray - ) -> np.ndarray: - """Convert fault coordinates to global coordinates. - - See global_coordinates_to_fault_coordinates for a description of fault - coordinates. - - Parameters - ---------- - fault_coordinates : np.ndarray - The fault coordinates of the point. - - Returns - ------- - np.ndarray - The global coordinates (lat, lon, depth) for this point. - - Raises - ------ - ValueError - If the fault coordinates are out of bounds. - """ - midpoint = np.sum(self.lengths()) / 2 - remaining_length = fault_coordinates[0] + midpoint - for plane in self.planes: - plane_length = plane.length - if remaining_length < plane_length or np.isclose( - remaining_length, plane_length - ): - return plane.plane_coordinates_to_global_coordinates( - np.array( - [ - remaining_length / plane_length - 1 / 2, - fault_coordinates[1] / plane.width - 1 / 2, - ] - ), - ) - remaining_length -= plane_length - raise ValueError("Specified fault coordinates out of bounds.") diff --git a/nshmdb/nshmdb.py b/nshmdb/nshmdb.py index bcd46fd..f53c2d2 100644 --- a/nshmdb/nshmdb.py +++ b/nshmdb/nshmdb.py @@ -21,20 +21,42 @@ >>> db.get_rupture_faults(0) # Should return two faults in this rupture. """ +import collections import dataclasses import importlib.resources import sqlite3 from pathlib import Path from sqlite3 import Connection +from typing import Optional import numpy as np -import qcore.coordinates - -from nshmdb import fault -from nshmdb.fault import Fault +from qcore import coordinates +from source_modelling.sources import Fault, Plane @dataclasses.dataclass +class FaultInfo: + """Fault metadata stored in the database. + + fault_id : int + The id of the fault. + name : str + The name of the fault. + parent_id : int + The id of the parent fault for this fault. + rake : float + The rake of the fault. + tect_type : int + The tectonic type of the fault. + """ + + fault_id: int + name: str + parent_id: int + rake: float + tect_type: Optional[int] + + class NSHMDB: """Class for interacting with the NSHMDB database. @@ -46,6 +68,9 @@ class NSHMDB: db_filepath: Path + def __init__(self, db_filepath: Path): + self.db_filepath = db_filepath + def create(self): """Create the tables for the NSHMDB database.""" schema_traversable = importlib.resources.files("nshmdb.schema") / "schema.sql" @@ -87,7 +112,13 @@ def insert_parent(self, conn: Connection, parent_id: int, parent_name: str): ) def insert_fault( - self, conn: Connection, fault_id: int, parent_id: int, fault: Fault + self, + conn: Connection, + fault_id: int, + parent_id: int, + fault_name: str, + fault_rake: float, + fault: Fault, ): """Insert fault data into the database. @@ -99,12 +130,16 @@ def insert_fault( ID of the fault. parent_id : int ID of the parent fault. + fault_name : str + The name of the fault. + fault_rake : float + The rake of the fault. fault : Fault Fault object containing fault geometry. """ conn.execute( - """INSERT OR REPLACE INTO fault (fault_id, name, parent_id) VALUES (?, ?, ?)""", - (fault_id, fault.name, parent_id), + """INSERT OR REPLACE INTO fault (fault_id, name, rake, parent_id) VALUES (?, ?, ?, ?)""", + (fault_id, fault_name, fault_rake, parent_id), ) for plane in fault.planes: conn.execute( @@ -119,16 +154,14 @@ def insert_fault( bottom_left_lon, top_depth, bottom_depth, - rake, fault_id ) VALUES ( - ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ? + ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ? )""", ( *plane.corners[:, :2].ravel(), plane.corners[0, 2], plane.corners[-1, 2], - plane.rake, fault_id, ), ) @@ -183,7 +216,6 @@ def get_fault(self, fault_id: int) -> Fault: bottom_left_lon, top, bottom, - rake, _, ) in cursor.fetchall(): corners = np.array( @@ -194,14 +226,17 @@ def get_fault(self, fault_id: int) -> Fault: [bottom_left_lat, bottom_left_lon, bottom], ] ) - planes.append( - fault.FaultPlane(qcore.coordinates.wgs_depth_to_nztm(corners), rake) - ) + planes.append(Plane(coordinates.wgs_depth_to_nztm(corners))) + cursor.execute("SELECT * from fault where fault_id = ?", (fault_id,)) + return Fault(planes) + + def get_fault_info(self, fault_id: int) -> FaultInfo: + with self.connection() as conn: + cursor = conn.cursor() cursor.execute("SELECT * from fault where fault_id = ?", (fault_id,)) - fault_id, name, _, _ = cursor.fetchone() - return Fault(name, None, planes) + return FaultInfo(*cursor.fetchone()) - def get_rupture_faults(self, rupture_id: int) -> list[Fault]: + def get_rupture_faults(self, rupture_id: int) -> dict[str, Fault]: """Retrieve faults involved in a rupture from the database. Parameters @@ -210,7 +245,9 @@ def get_rupture_faults(self, rupture_id: int) -> list[Fault]: Returns ------- - list[Fault] + dict[str, Fault] + A dictionary with fault names as keys, and fault geometry + as values. """ with self.connection() as conn: cursor = conn.cursor() @@ -225,8 +262,7 @@ def get_rupture_faults(self, rupture_id: int) -> list[Fault]: (rupture_id,), ) fault_planes = cursor.fetchall() - cur_parent_id = None - faults = [] + faults = collections.defaultdict(lambda: Fault([])) for ( _, top_left_lat, @@ -239,20 +275,10 @@ def get_rupture_faults(self, rupture_id: int) -> list[Fault]: bottom_left_lon, top, bottom, - rake, _, parent_id, parent_name, ) in fault_planes: - if parent_id != cur_parent_id: - faults.append( - Fault( - name=parent_name, - tect_type=None, - planes=[], - ) - ) - cur_parent_id = parent_id corners = np.array( [ [top_left_lat, top_left_lon, top], @@ -261,7 +287,7 @@ def get_rupture_faults(self, rupture_id: int) -> list[Fault]: [bottom_left_lat, bottom_left_lon, bottom], ] ) - faults[-1].planes.append( - fault.FaultPlane(qcore.coordinates.wgs_depth_to_nztm(corners), rake) + faults[parent_name].planes.append( + Plane(coordinates.wgs_depth_to_nztm(corners)) ) return faults diff --git a/nshmdb/plotting/rupture.py b/nshmdb/plotting/rupture.py index e703ebd..5e6e2e9 100644 --- a/nshmdb/plotting/rupture.py +++ b/nshmdb/plotting/rupture.py @@ -10,8 +10,8 @@ from pathlib import Path import numpy as np -from nshmdb.fault import Fault from pygmt_helper import plotting +from source_modelling.sources import Fault def plot_rupture(title: str, faults: list[Fault], output_filepath: Path): @@ -26,7 +26,7 @@ def plot_rupture(title: str, faults: list[Fault], output_filepath: Path): output_filepath : Path The filepath to output the figure to. """ - corners = np.vstack([fault.corners() for fault in faults]) + corners = np.vstack([fault.corners for fault in faults]) region = ( corners[:, 1].min() - 0.5, corners[:, 1].max() + 0.5, diff --git a/nshmdb/schema/schema.sql b/nshmdb/schema/schema.sql index a566c19..16bba0e 100644 --- a/nshmdb/schema/schema.sql +++ b/nshmdb/schema/schema.sql @@ -2,6 +2,7 @@ CREATE TABLE IF NOT EXISTS fault ( fault_id INTEGER PRIMARY KEY NOT NULL, name TEXT NOT NULL, parent_id INTEGER NOT NULL, + rake REAL NOT NULL, tect_type INT, FOREIGN KEY(parent_id) REFERENCES parent_fault(parent_id) ); @@ -23,7 +24,6 @@ CREATE TABLE IF NOT EXISTS fault_plane ( bottom_left_lon REAL NOT NULL, top_depth REAL NOT NULL, bottom_depth REAL NOT NULL, - rake REAL NOT NULL, fault_id INTEGER NOT NULL, FOREIGN KEY(fault_id) REFERENCES fault(fault_id) ); diff --git a/nshmdb/scripts/nshm_db_generator.py b/nshmdb/scripts/nshm_db_generator.py index f845e22..bb2d1cd 100644 --- a/nshmdb/scripts/nshm_db_generator.py +++ b/nshmdb/scripts/nshm_db_generator.py @@ -22,20 +22,20 @@ Example: python generate_nshm2022_data.py data/cru_solutions.zip output/nshm2022.sqlite """ + import zipfile from pathlib import Path from typing import Annotated import geojson -import nshmdb.fault import numpy as np import pandas as pd import qcore.coordinates import tqdm import typer from geojson import FeatureCollection -from nshmdb.fault import Fault from nshmdb.nshmdb import NSHMDB +from source_modelling.sources import Fault, Plane app = typer.Typer() @@ -46,7 +46,7 @@ def extract_faults_from_info( fault_info_list: FeatureCollection, -) -> list[Fault]: +) -> dict[str, Fault]: """Extract the fault geometry from the fault information description. Parameters @@ -56,10 +56,11 @@ def extract_faults_from_info( Returns ------- - list[Fault] - The list of extracted faults. + dict[str, Fault] + A dictionary of extracted faults. The key is the name of the + fault. """ - faults = [] + faults = {} for i in range(len(fault_info_list.features)): fault_feature = fault_info_list[i] fault_trace = list(geojson.utils.coords(fault_feature)) @@ -71,7 +72,6 @@ def extract_faults_from_info( projected_width = 0 else: projected_width = bottom / np.tan(np.radians(dip)) - rake = fault_feature.properties["Rake"] planes = [] for i in range(len(fault_trace) - 1): top_left = qcore.coordinates.wgs_depth_to_nztm( @@ -93,8 +93,8 @@ def extract_faults_from_info( bottom_left = top_left + dip_dir_direction bottom_right = top_right + dip_dir_direction corners = np.array([top_left, top_right, bottom_right, bottom_left]) - planes.append(nshmdb.fault.FaultPlane(corners, rake)) - faults.append(Fault(name, None, planes)) + planes.append(Plane(corners)) + faults[name] = Fault(planes) return faults @@ -125,7 +125,6 @@ def main( with zipfile.ZipFile( cru_solutions_zip_path, "r" ) as cru_solutions_zip_file, db.connection() as conn: - with cru_solutions_zip_file.open( str(FAULT_INFORMATION_PATH) ) as fault_info_handle: @@ -134,12 +133,15 @@ def main( faults = extract_faults_from_info(faults_info) if not skip_faults_creation: - for i, fault in enumerate(faults): + for i, fault in enumerate(faults.values()): fault_info = faults_info[i] + fault_id = fault_info.properties["FaultID"] parent_id = fault_info.properties["ParentID"] + fault_name = fault_info.properties["FaultName"] + fault_rake = fault_info.properties["Rake"] db.insert_parent(conn, parent_id, fault_info.properties["ParentName"]) db.insert_fault( - conn, fault_info.properties["FaultID"], parent_id, fault + conn, fault_id, parent_id, fault_name, fault_rake, fault ) if not skip_rupture_creation: with cru_solutions_zip_file.open(