diff --git a/README.md b/README.md index f999533..3ad7bd2 100644 --- a/README.md +++ b/README.md @@ -63,70 +63,10 @@ To integrate with PyCharm (if you use it) see [PyCharm's Documentation](https:// │   └── sub_seismo_on_fault_mfds.csv └── WARNING.md ``` -2. Copy the `fault_sections.geojson` and `fast_indices.csv` file into the directory contaning the script. -3. For some reason (pandas shenanigans?) the `fast_indices.csv` which maps between rupture ids and fault segment ids uses ints for the ruptures and *floats* for the fault segment ids. Run the following in the command line to turn the floats into integers. -```bash -$ sed -i 's/\.0//' fast_indices.csv -``` -4. Creating an empty `nshm2022.db` database file. -```bash -$ sqlite3 nshm2022.db < schema.db -``` -5. Populate the database using the `nshm_geojson_fault_parser.py` script. -```bash -$ python nshm_geojson_fault_parser.py -``` -This step will take some time. - -# Database Schema -```mermaid -erDiagram - fault { - fault_id INT PK - name TEXT - parent_id INT FK - tect_type INT - } +**YOU DO NOT NEED TO EXTRACT THIS FILE ANYWHERE** - parent_fault { - parent_id INT PK - name TEXT - } - - fault_segment { - segment_id INT PK - strike REAL - rake REAL - dip REAL - dtop REAL - dbottom REAL - length REAL - width REAL - dip_dir REAL - clon REAL - clat REAL - } - - rupture { - rupture_id INT PK - } - - rupture_faults { - rupture_fault_id INT PK - rupture_id INT FK - fault_id INT FK - } - - magnitude_frequency_distribution { - entry_id INT PK - fault_id INT FK - magnitude REAL - rate REAL - } - - rupture }o--o{ rupture_faults : "has many" - fault }o--o{ rupture_faults : "belonging to many" - fault ||--o{ fault_segment : "containing many" - parent_fault ||--o{ fault : "are groups of" - fault ||--o{ magnitude_frequency_distribution : "erupts with annual rate at given magnitude" +2. After cloning this repository and installing the depedencies, run the following +```bash +python nshmdb/scripts/nshm_db_generator.py nshmdb.db ``` +This will take some time. diff --git a/nshm_geojson_fault_parser.py b/nshm_geojson_fault_parser.py deleted file mode 100644 index ced1aa2..0000000 --- a/nshm_geojson_fault_parser.py +++ /dev/null @@ -1,237 +0,0 @@ -#!/usr/bin/env python3 -import csv -import json -import sqlite3 -from pathlib import Path -from sqlite3 import Connection -from typing import Annotated, Any - -import numpy as np -import qcore.geo -import typer - -app = typer.Typer() - - -def centre_point_of_fault( - point_a_coords: (float, float), - point_b_coords: (float, float), - dip: float, - dip_dir: float, - width: float, -) -> (float, float): - """Find the centre point of a fault defined by two points and a dip direction. - - This function calculates the centre point of a fault given the coordinates - of two points on the fault trace, the dip angle, dip direction, and the - width of the fault. The centre point is defined as the midpoint of the - line connecting points a and a, shifted along the dip direction by half - the *projected* fault width. See the diagram. - - point a point b - +───────────────+ - │ │ │ │ - │ │ │ │ - │ │dip dir│ │ - │ │ │ │ - │ ∨ │ │ width - │ *centre │ │ - │ │ │ - │ │ │ - │ │ │ - └───────────────┘ - - The dip angle should be provided in degrees, with 0 degrees indicating - an (impossible) horizontal fault and 90 degrees indicating a vertical - fault. The dip direction is a bearing in degrees from north. - - Parameters - ---------- - a : (float, float) - The (lat, lon) coordinates of the point a. - b : (float, float) - The (lat, lon) coordinates of the point b. - dip : float - The dip angle of the fault. - dip_dir : float - The dip direction of the fault. - width : float - The width of the fault. Note this is the actual width, not the - projected width. - - Returns - ------- - (float, float) - The (lat, lon) coordinates of the (projected) centre point of - the fault. - """ - a_lat, a_lon = point_a_coords - b_lat, b_lon = point_b_coords - c_lon, c_lat = qcore.geo.ll_mid(a_lon, a_lat, b_lon, b_lat) - projected_width = width * np.cos(np.radians(dip)) - return qcore.geo.ll_shift(c_lat, c_lon, projected_width / 2, dip_dir) - - -def insert_magnitude_frequency_distribution( - conn: Connection, magnitude_frequency_distribution: list[dict[str, float | str]] -): - """Inserts magnitude-frequency distribution data into a database. - - Parameters - ---------- - conn : sqlite3.Connection - A connection object to the SQLite database. - magnitude_frequency_distribution : list[dict[str, Union[float, str]]] - A list of dictionaries containing magnitude-frequency distribution - data. Each dictionary should have keys 'Section Index', representing - the fault section index, and keys representing magnitudes associated - with their respective annual occurrence rate. - """ - for section_distribution in magnitude_frequency_distribution: - segment_id = int(section_distribution["Section Index"]) - for magnitude_key, rate_raw in section_distribution.items(): - if magnitude_key == "Section Index": - continue - magnitude = float(magnitude_key) - rate = float(rate_raw) - conn.execute( - """ - INSERT INTO magnitude_frequency_distribution ( - fault_id, magnitude, rate - ) VALUES (?, ?, ?) - """, - (segment_id, magnitude, rate), - ) - - -def insert_faults(conn: Connection, fault_map: dict[str, Any]): - """Inserts fault data into the database. - - Parameters - ---------- - conn : sqlite3.Connection - A connection object to the SQLite database. - fault_map : list[dict[str, Any]] - A list of dictionaries containing fault data. Each dictionary should - represent a fault feature and contain necessary properties like dip, - rake, depth, etc. - """ - for feature in fault_map: - properties = feature["properties"] - dip = properties["DipDeg"] - rake = properties["Rake"] - dbottom = properties["LowDepth"] - dtop = properties["UpDepth"] - dip_dir = properties["DipDir"] - leading_edge = feature["geometry"]["coordinates"] - width = float(dbottom / np.sin(np.radians(dip))) - fault_name = properties["FaultName"] - fault_id = properties["FaultID"] - parent_id = properties["ParentID"] - parent_name = properties["ParentName"] - conn.execute( - """INSERT OR REPLACE INTO parent_fault (parent_id, name) - VALUES (?, ?)""", - (parent_id, parent_name), - ) - conn.execute( - "INSERT INTO fault (fault_id, name, parent_id) VALUES (?, ?, ?)", - (fault_id, fault_name, parent_id), - ) - for i in range(len(leading_edge) - 1): - left = tuple(reversed(leading_edge[i])) - right = tuple(reversed(leading_edge[i + 1])) - c_lat, c_lon = centre_point_of_fault(left, right, dip, dip_dir, dbottom) - strike = qcore.geo.ll_bearing(*left[::-1], *right[::-1]) - length = qcore.geo.ll_dist(*left[::-1], *right[::-1]) - conn.execute( - """ - INSERT INTO fault_segment ( - strike, - rake, - dip, - dtop, - dbottom, - length, - width, - dip_dir, - clon, - clat, - fault_id - ) VALUES (?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?) - """, - ( - strike, - rake, - dip, - dtop, - dbottom, - length, - width, - dip_dir, - c_lon, - c_lat, - fault_id, - ), - ) - - -def insert_ruptures(conn: Connection, indices: dict[int, int]): - """Inserts rupture data into the database. - - Parameters - ---------- - conn : sqlite3.Connection - A connection object to the SQLite database. - indices : dict[int, int] - A dictionary containing rupture indices mapped to fault indices. If - indices[0] = 1, then the fault 1 is involved in rupture 0. - """ - for row in indices: - rupture_idx, fault_idx = [int(value) for value in row.values()] - conn.execute( - "INSERT OR REPLACE INTO rupture (rupture_id) VALUES (?)", (rupture_idx,) - ) - - conn.execute( - "INSERT INTO rupture_faults (rupture_id, fault_id) VALUES (?, ?)", - (rupture_idx, fault_idx), - ) - - -@app.command() -def main( - fault_sections_geojson_filepath: Annotated[ - Path, - typer.Option(help="Fault sections geojson file", readable=True, exists=True), - ] = "fault_sections.geojson", - fast_indices_filepath: Annotated[ - Path, typer.Option(help="Fast indices csv file", readable=True, exists=True) - ] = "fast_indices.csv", - mfds_filepath: Annotated[ - Path, typer.Option(help="MFDS filepath", readable=True, exists=True) - ] = "sub_seismo_on_fault_mfds.csv", - sqlite_db_path: Annotated[ - Path, typer.Option(help="Output SQLite DB path", writable=True, exists=True) - ] = "nshm2022.db", -): - """Generate the NSHM2022 rupture data from a CRU system solution package.""" - with open(fault_sections_geojson_filepath, "r", encoding="utf-8") as fault_file: - geojson_object = json.load(fault_file) - with open(fast_indices_filepath, "r", encoding="utf-8") as csv_file_handle: - csv_reader = csv.DictReader(csv_file_handle) - indices = list(csv_reader) - with open(mfds_filepath, "r", encoding="utf-8") as mfds_file_handle: - csv_reader = csv.DictReader(mfds_file_handle) - magnitude_frequency_distribution = list(csv_reader) - with sqlite3.connect(sqlite_db_path) as conn: - # To enforce foreign key constraints. - conn.execute("PRAGMA foreign_keys = 1") - - insert_faults(conn, geojson_object["features"]) - insert_ruptures(conn, indices) - insert_magnitude_frequency_distribution(conn, magnitude_frequency_distribution) - - -if __name__ == "__main__": - app() diff --git a/nshmdb/__init__.py b/nshmdb/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/nshmdb/fault.py b/nshmdb/fault.py new file mode 100644 index 0000000..4057259 --- /dev/null +++ b/nshmdb/fault.py @@ -0,0 +1,503 @@ +"""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 + +import numpy as np +from typing import Optional +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 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.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 new file mode 100644 index 0000000..bcd46fd --- /dev/null +++ b/nshmdb/nshmdb.py @@ -0,0 +1,267 @@ +""" +Module to interact with the NSHMDB (National Seismic Hazard Model Database). + +This module provides classes and functions to interact with an SQLite database +containing national seismic hazard model data. It includes functionalities to +insert fault and rupture data into the database, as well as retrieve fault +information associated with ruptures. + +Classes +------- +NSHMDB + Class for interacting with the NSHMDB database. + +Usage +----- +Initialize an instance of NSHMDB with the path to the SQLite database file. +Use the methods of the NSHMDB class to interact with fault and rupture data +in the database. + +>>> db = NSHMDB('path/to/nshm.db') +>>> db.get_rupture_faults(0) # Should return two faults in this rupture. +""" + +import dataclasses +import importlib.resources +import sqlite3 +from pathlib import Path +from sqlite3 import Connection + +import numpy as np +import qcore.coordinates + +from nshmdb import fault +from nshmdb.fault import Fault + + +@dataclasses.dataclass +class NSHMDB: + """Class for interacting with the NSHMDB database. + + Parameters + ---------- + db_filepath : Path + Path to the SQLite database file. + """ + + db_filepath: Path + + def create(self): + """Create the tables for the NSHMDB database.""" + schema_traversable = importlib.resources.files("nshmdb.schema") / "schema.sql" + with importlib.resources.as_file(schema_traversable) as schema_path: + with open(schema_path, "r", encoding="utf-8") as schema_file_handle: + schema = schema_file_handle.read() + with self.connection() as conn: + conn.executescript(schema) + + def connection(self) -> Connection: + """Establish a connection to the SQLite database. + + Returns + ------- + Connection + """ + return sqlite3.connect(self.db_filepath) + + # The functions `insert_parent`, `insert_fault`, and `add_fault_to_rupture` + # reuse a connection for efficiency (rather than use db.connection()). There + # are thousands of faults and tens of millions of rupture, fault binding + # pairs. Without reusing a connection it takes hours to setup the database. + + def insert_parent(self, conn: Connection, parent_id: int, parent_name: str): + """Insert parent fault data into the database. + + Parameters + ---------- + conn : Connection + The db connection object. + parent_id : int + ID of the parent fault. + name : str + Name of the parent fault. + """ + conn.execute( + """INSERT OR REPLACE INTO parent_fault (parent_id, name) VALUES (?, ?)""", + (parent_id, parent_name), + ) + + def insert_fault( + self, conn: Connection, fault_id: int, parent_id: int, fault: Fault + ): + """Insert fault data into the database. + + Parameters + ---------- + conn : Connection + The db connection object. + fault_id : int + ID of the fault. + parent_id : int + ID of the parent 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), + ) + for plane in fault.planes: + conn.execute( + """INSERT INTO fault_plane ( + top_left_lat, + top_left_lon, + top_right_lat, + top_right_lon, + bottom_right_lat, + bottom_right_lon, + bottom_left_lat, + 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, + ), + ) + + def add_fault_to_rupture(self, conn: Connection, rupture_id: int, fault_id: int): + """Insert rupture data into the database. + + Parameters + ---------- + conn : Connection + The db connection object. + rupture_id : int + ID of the rupture. + fault_ids : list[int] + List of faults involved in the rupture. + """ + conn.execute( + "INSERT OR REPLACE INTO rupture (rupture_id) VALUES (?)", (rupture_id,) + ) + conn.execute( + "INSERT INTO rupture_faults (rupture_id, fault_id) VALUES (?, ?)", + (rupture_id, fault_id), + ) + + def get_fault(self, fault_id: int) -> Fault: + """Get a specific fault definition from a database. + + Parameters + ---------- + fault_id : int + The id of the fault to retreive. + + Returns + ------- + Fault + The fault geometry. + """ + + with self.connection() as conn: + cursor = conn.cursor() + cursor.execute("SELECT * from fault_plane where fault_id = ?", (fault_id,)) + planes = [] + for ( + _, + top_left_lat, + top_left_lon, + top_right_lat, + top_right_lon, + bottom_right_lat, + bottom_right_lon, + bottom_left_lat, + bottom_left_lon, + top, + bottom, + rake, + _, + ) in cursor.fetchall(): + corners = np.array( + [ + [top_left_lat, top_left_lon, top], + [top_right_lat, top_right_lon, top], + [bottom_right_lat, bottom_right_lon, bottom], + [bottom_left_lat, bottom_left_lon, bottom], + ] + ) + planes.append( + fault.FaultPlane(qcore.coordinates.wgs_depth_to_nztm(corners), rake) + ) + cursor.execute("SELECT * from fault where fault_id = ?", (fault_id,)) + fault_id, name, _, _ = cursor.fetchone() + return Fault(name, None, planes) + + def get_rupture_faults(self, rupture_id: int) -> list[Fault]: + """Retrieve faults involved in a rupture from the database. + + Parameters + ---------- + rupture_id : int + + Returns + ------- + list[Fault] + """ + with self.connection() as conn: + cursor = conn.cursor() + cursor.execute( + """SELECT fs.*, p.parent_id, p.name + FROM fault_plane fs + JOIN rupture_faults rf ON fs.fault_id = rf.fault_id + JOIN fault f ON fs.fault_id = f.fault_id + JOIN parent_fault p ON f.parent_id = p.parent_id + WHERE rf.rupture_id = ? + ORDER BY f.parent_id""", + (rupture_id,), + ) + fault_planes = cursor.fetchall() + cur_parent_id = None + faults = [] + for ( + _, + top_left_lat, + top_left_lon, + top_right_lat, + top_right_lon, + bottom_right_lat, + bottom_right_lon, + bottom_left_lat, + 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], + [top_right_lat, top_right_lon, top], + [bottom_right_lat, bottom_right_lon, bottom], + [bottom_left_lat, bottom_left_lon, bottom], + ] + ) + faults[-1].planes.append( + fault.FaultPlane(qcore.coordinates.wgs_depth_to_nztm(corners), rake) + ) + return faults diff --git a/nshmdb/plotting/rupture.py b/nshmdb/plotting/rupture.py new file mode 100644 index 0000000..e703ebd --- /dev/null +++ b/nshmdb/plotting/rupture.py @@ -0,0 +1,48 @@ +""" +Module to plot ruptures using PyGMT. + +This module provides functionality to visualize fault ruptures on a map using PyGMT. + +Functions: + 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], output_filepath: Path): + """Plot faults involved in a rupture scenario using pygmt. + + Parameters + ---------- + title : str + 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 = ( + corners[:, 1].min() - 0.5, + corners[:, 1].max() + 0.5, + corners[:, 0].min() - 0.25, + corners[:, 0].max() + 0.25, + ) + fig = plotting.gen_region_fig(title, region=region) + + for rupture_fault in faults: + for plane in rupture_fault.planes: + corners = plane.corners + fig.plot( + x=corners[:, 1].tolist() + [corners[0, 1]], + y=corners[:, 0].tolist() + [corners[0, 0]], + pen="1p", + fill="red", + ) + + fig.savefig(output_filepath) diff --git a/schema.sql b/nshmdb/schema/schema.sql similarity index 73% rename from schema.sql rename to nshmdb/schema/schema.sql index 9837b3e..a566c19 100644 --- a/schema.sql +++ b/nshmdb/schema/schema.sql @@ -11,18 +11,19 @@ CREATE TABLE IF NOT EXISTS parent_fault ( name TEXT NOT NULL ); -CREATE TABLE IF NOT EXISTS fault_segment ( - segment_id INTEGER PRIMARY KEY, - strike REAL NOT NULL, +CREATE TABLE IF NOT EXISTS fault_plane ( + plane_id INTEGER PRIMARY KEY, + top_left_lat REAL NOT NULL, + top_left_lon REAL NOT NULL, + top_right_lat REAL NOT NULL, + top_right_lon REAL NOT NULL, + bottom_right_lat REAL NOT NULL, + bottom_right_lon REAL NOT NULL, + bottom_left_lat REAL NOT NULL, + bottom_left_lon REAL NOT NULL, + top_depth REAL NOT NULL, + bottom_depth REAL NOT NULL, rake REAL NOT NULL, - dip REAL NOT NULL, - dtop REAL NOT NULL, - dbottom REAL NOT NULL, - length REAL NOT NULL, - width REAL NOT NULL, - dip_dir REAL NOT NULL, - clon REAL NOT NULL, - clat 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 new file mode 100644 index 0000000..f845e22 --- /dev/null +++ b/nshmdb/scripts/nshm_db_generator.py @@ -0,0 +1,161 @@ +#!/usr/bin/env python3 +""" +NSHM2022 Rupture Data Generation Script + +This script generates NSHM2022 rupture data from a CRU system solution package. + +Usage: + python script_name.py [OPTIONS] CRU_SOLUTIONS_ZIP_PATH SQLITE_DB_PATH + +Arguments: + CRU_SOLUTIONS_ZIP_PATH : str + Path to the CRU solutions zip file. + SQLITE_DB_PATH : str + Output SQLite DB path. + +Options: + --skip-faults-creation : bool, optional + If flag is set, skip fault creation. + --skip-rupture-creation : bool, optional + If flag is set, skip rupture creation. + +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 + +app = typer.Typer() + + +FAULT_INFORMATION_PATH = Path("ruptures") / "fault_sections.geojson" +RUPTURE_FAULT_JOIN_PATH = Path("ruptures") / "fast_indices.csv" + + +def extract_faults_from_info( + fault_info_list: FeatureCollection, +) -> list[Fault]: + """Extract the fault geometry from the fault information description. + + Parameters + ---------- + fault_info_list : FeatureCollection + The GeoJson object containing the fault definitions. + + Returns + ------- + list[Fault] + The list of extracted 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)) + name = fault_feature.properties["FaultName"] + dip_dir = fault_feature.properties["DipDir"] + dip = fault_feature.properties["DipDeg"] + bottom = fault_feature.properties["LowDepth"] + if dip == 90: + 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( + np.append(fault_trace[i][::-1], 0) + ) + top_right = qcore.coordinates.wgs_depth_to_nztm( + np.append(fault_trace[i + 1][::-1], 0) + ) + dip_dir_direction = ( + np.array( + [ + projected_width * np.cos(np.radians(dip_dir)), + projected_width * np.sin(np.radians(dip_dir)), + bottom, + ] + ) + * 1000 + ) + 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)) + return faults + + +@app.command() +def main( + cru_solutions_zip_path: Annotated[ + Path, + typer.Argument( + help="CRU solutions zip file", readable=True, dir_okay=False, exists=True + ), + ], + sqlite_db_path: Annotated[ + Path, + typer.Argument(help="Output SQLite DB path", writable=True, dir_okay=False), + ], + skip_faults_creation: Annotated[ + bool, typer.Option(help="If flag is set, skip fault creation.") + ] = False, + skip_rupture_creation: Annotated[ + bool, typer.Option(help="If flag is set, skip rupture creation.") + ] = False, +): + """Generate the NSHM2022 rupture data from a CRU system solution package.""" + + db = NSHMDB(sqlite_db_path) + db.create() + + 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: + faults_info = geojson.load(fault_info_handle) + + faults = extract_faults_from_info(faults_info) + + if not skip_faults_creation: + for i, fault in enumerate(faults): + fault_info = faults_info[i] + parent_id = fault_info.properties["ParentID"] + db.insert_parent(conn, parent_id, fault_info.properties["ParentName"]) + db.insert_fault( + conn, fault_info.properties["FaultID"], parent_id, fault + ) + if not skip_rupture_creation: + with cru_solutions_zip_file.open( + str(RUPTURE_FAULT_JOIN_PATH) + ) as rupture_fault_join_handle: + rupture_fault_join_df = pd.read_csv(rupture_fault_join_handle) + rupture_fault_join_df["section"] = rupture_fault_join_df[ + "section" + ].astype("Int64") + for _, row in tqdm.tqdm( + rupture_fault_join_df.iterrows(), + desc="Binding ruptures to faults", + total=len(rupture_fault_join_df), + ): + db.add_fault_to_rupture(conn, row["rupture"], row["section"]) + + +if __name__ == "__main__": + app() diff --git a/pyproject.toml b/pyproject.toml deleted file mode 100644 index f78d61f..0000000 --- a/pyproject.toml +++ /dev/null @@ -1,28 +0,0 @@ -[tool.poetry] -name = "nshm2022db" -version = "0.1.0" -description = "Fault database for the GNS National Seismic Hazard Model." -authors = ["QuakeCoRE"] -readme = "README.md" - -[tool.poetry.dependencies] -python = "^3.12" -numpy = "*" -qcore = { git = "https://github.com/ucgmsim/qcore" } -typer = "*" - -[tool.poetry.group.dev.dependencies] -python-lsp-server = "*" -python-lsp-black = "*" -python-lsp-isort = "*" -black = "*" -rope = "*" -pyflakes = "*" -mccabe = "*" -pycodestyle = "*" -doq = "*" - - -[build-system] -requires = ["poetry-core"] -build-backend = "poetry.core.masonry.api" diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..9a9e03f --- /dev/null +++ b/requirements.txt @@ -0,0 +1,10 @@ +dataclasses +geojson +numpy +pandas +pygments +pygmt +qcore +scipy +tqdm +typer