From 1fff878f1da7de5b6a182d6edca42595009b5ef8 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Thu, 16 May 2024 14:09:16 +1200 Subject: [PATCH 01/39] package and scriptify nshm db generator --- nshmdb/__init__.py | 0 .../scripts/nshm_db_generator.py | 0 pyproject.toml | 9 +++++---- 3 files changed, 5 insertions(+), 4 deletions(-) create mode 100644 nshmdb/__init__.py rename nshm_geojson_fault_parser.py => nshmdb/scripts/nshm_db_generator.py (100%) diff --git a/nshmdb/__init__.py b/nshmdb/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/nshm_geojson_fault_parser.py b/nshmdb/scripts/nshm_db_generator.py similarity index 100% rename from nshm_geojson_fault_parser.py rename to nshmdb/scripts/nshm_db_generator.py diff --git a/pyproject.toml b/pyproject.toml index f78d61f..f377b2a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,5 +1,5 @@ [tool.poetry] -name = "nshm2022db" +name = "nshmdb" version = "0.1.0" description = "Fault database for the GNS National Seismic Hazard Model." authors = ["QuakeCoRE"] @@ -7,9 +7,8 @@ readme = "README.md" [tool.poetry.dependencies] python = "^3.12" -numpy = "*" qcore = { git = "https://github.com/ucgmsim/qcore" } -typer = "*" +numpy = "^1.26.4" [tool.poetry.group.dev.dependencies] python-lsp-server = "*" @@ -22,7 +21,9 @@ mccabe = "*" pycodestyle = "*" doq = "*" - [build-system] requires = ["poetry-core"] build-backend = "poetry.core.masonry.api" + +[scripts] +nshm-db-generator = "nshmdb.scripts.nshm_db_generator:app" From 9daf496c1f106708b6d397c8382c358cb69b1af5 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Thu, 16 May 2024 14:20:32 +1200 Subject: [PATCH 02/39] add mypy type-checking --- pyproject.toml | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/pyproject.toml b/pyproject.toml index f377b2a..30cb685 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -20,6 +20,8 @@ pyflakes = "*" mccabe = "*" pycodestyle = "*" doq = "*" +mypy = "^1.10.0" +pylsp-mypy = "^0.6.8" [build-system] requires = ["poetry-core"] @@ -27,3 +29,9 @@ build-backend = "poetry.core.masonry.api" [scripts] nshm-db-generator = "nshmdb.scripts.nshm_db_generator:app" + +[tool.pylsp-mypy] +enabled = true +live_mode = true +strict = false +exclude = ["tests/*"] From 48a669f60f37072c4ccd25ffc3af79578e9dc95a Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Thu, 16 May 2024 14:20:54 +1200 Subject: [PATCH 03/39] Add fault module --- nshmdb/fault.py | 483 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 483 insertions(+) create mode 100644 nshmdb/fault.py diff --git a/nshmdb/fault.py b/nshmdb/fault.py new file mode 100644 index 0000000..507e146 --- /dev/null +++ b/nshmdb/fault.py @@ -0,0 +1,483 @@ +"""Module for representing fault segments and faults. + +This module provides classes and functions for representing fault segments 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. + +FaultSegment: + A representation of a single segment of a Fault. + +Fault: + A representation of a fault, consisting of one or more FaultSegments. +""" + +import dataclasses +from enum import Enum + +import numpy as np +import qcore.coordinates +import qcore.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 FaultSegment: + """A representation of a single segment of a Fault. + + This class represents a single segment 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 segment. + """ + + _corners: np.ndarray + rake: float + + def __init__(self, corners: np.ndarray, rake: float): + self._corners = qcore.coordinates.wgs_depth_to_nztm(corners) + self.rake = rake + + @property + def corners(self) -> np.ndarray: + """ + + Returns + ------- + np.ndarray + The corners of the fault segment in (lat, lon, depth) format. + """ + return qcore.coordinates.nztm_to_wgs_depth(self._corners) + + @property + def length_m(self) -> float: + """ + Returns + ------- + float + The length of the fault segment (in metres). + """ + return np.linalg.norm(self._corners[1] - self._corners[0]) + + @property + def width_m(self) -> float: + """ + Returns + ------- + float + The width of the fault segment (in metres). + """ + return np.linalg.norm(self._corners[-1] - self._corners[0]) + + @property + def bottom_m(self) -> float: + """ + Returns + ------- + float + The bottom depth (in metres). + """ + return self._corners[-1, -1] + + @property + def width(self) -> float: + """ + Returns + ------- + float + The width of the fault segment (in kilometres). + """ + return self.width_m / _KM_TO_M + + @property + def length(self) -> float: + """ + Returns + ------- + float + The length of the fault segment (in kilometres). + """ + return self.length_m / _KM_TO_M + + @property + def projected_width_m(self) -> float: + """ + Returns + ------- + float + The projected width of the fault segment (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 segment (in kilometres). + """ + return self.projected_width / _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[1] - self._corners[0] + return ( + np.degrees( + qcore.geo.oriented_angle_wrt_normal( + north_direction, strike_direction, up_direction + ) + ) + % 360 + ) + + @property + def dip_dir(self) -> float: + """ + Returns + ------- + float + The bearing of the dip direction (from north; in degrees). + """ + north_direction = np.array([1, 0, 0]) + up_direction = np.array([0, 0, 1]) + dip_direction = self._corners[-1] - self._corners[0] + dip_direction[-1] = 0 + return ( + np.degrees( + qcore.geo.oriented_angle_wrt_normal( + north_direction, dip_direction, up_direction + ) + ) + % 360 + ) + + @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 segment_coordinates_to_global_coordinates( + self, segment_coordinates: np.ndarray + ) -> np.ndarray: + """Convert segment coordinates to nztm global coordinates. + + Parameters + ---------- + segment_coordinates : np.ndarray + Segment coordinates to convert. Segment coordinates are + 2D coordinates (x, y) given for a fault segment (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 + origin for segment coordinates is the centre of the fault. + + +x + -1/2,-1/2 ─────────────────> + ┌─────────────────────┐ │ + │ < width > │ │ + │ ^ │ │ + │ length│ │ +y + │ v │ │ + │ │ │ + └─────────────────────┘ ∨ + 1/2,1/2 + + Returns + ------- + 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] + frame = np.vstack((top_right - origin, bottom_left - origin)) + offset = np.array([1 / 2, 1 / 2]) + + return qcore.coordinates.nztm_to_wgs_depth( + origin + (segment_coordinates + offset) @ frame + ) + + def global_coordinates_to_segment_coordinates( + self, + global_coordinates: np.ndarray, + ) -> np.ndarray: + """Convert coordinates (lat, lon, depth) to segment coordinates (x, y). + + See segment_coordinates_to_global_coordinates for a description of + segment coordinates. + + Parameters + ---------- + global_coordinates : np.ndarray + Global coordinates to convert. + + Returns + ------- + np.ndarray + The segment coordinates (x, y) representing the position of + global_coordinates on the fault segment. + + Raises + ------ + 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] + frame = np.vstack((top_right - origin, bottom_left - origin)) + offset = qcore.coordinates.wgs_depth_to_nztm(global_coordinates) - origin + segment_coordinates, residual, _, _ = np.linalg.lstsq(frame.T, offset) + if not np.isclose(residual[0], 0): + raise ValueError("Coordinates do not lie in fault plane.") + return segment_coordinates - np.array([1 / 2, 1 / 2]) + + def global_coordinates_in_segment(self, global_coordinates: np.ndarray) -> bool: + """Test if some global coordinates lie in the bounds of a segment. + + 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 segment. + """ + + segment_coordinates = self.global_coordinates_to_segment_coordinates( + global_coordinates + ) + return np.all( + np.logical_or( + np.abs(segment_coordinates) < 1 / 2, + np.isclose(np.abs(segment_coordinates), 1 / 2, atol=1e-4), + ) + ) + + def centroid(self) -> np.ndarray: + """Returns the centre of the fault segment. + + Returns + ------- + np.ndarray + A 1 x 3 dimensional vector representing the centroid of the fault + plane in (lat, lon, depth) format. + + """ + + return qcore.coordinates.nztm_to_wgs_depth( + np.mean(self._corners, axis=0).reshape((1, -1)) + ).ravel() + + +@dataclasses.dataclass +class Fault: + """A representation of a fault, consisting of one or more FaultSegments. + + This class represents a fault, which is composed of one or more FaultSegments. + It provides methods for computing the area of the fault, getting the widths and + lengths of all fault segments, 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 : str + The type of fault this is (e.g. crustal, volcanic, subduction). + segments : list[FaultSegment] + A list containing all the FaultSegments that constitute the fault. + + Methods + ------- + area: + Compute the area of a fault. + widths: + Get the widths of all fault segments. + lengths: + Get the lengths of all fault segments. + 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: str + segments: list[FaultSegment] + + def area(self) -> float: + """Compute the area of a fault. + + Returns + ------- + float + The area of the fault. + """ + return sum(segment.width * segment.length for segment in self.segments) + + def widths(self) -> np.ndarray: + """Get the widths of all fault segments. + + Returns + ------- + np.ndarray of shape (1 x n) + The widths of all fault segments contained in this fault. + """ + return np.array([seg.width for seg in self.segments]) + + def lengths(self) -> np.ndarray: + """Get the lengths of all fault segments. + + Returns + ------- + np.ndarray of shape (1 x n) + The lengths of all fault segments contained in this fault. + """ + return np.array([seg.length for seg in self.segments]) + + def corners(self) -> np.ndarray: + """Get all corners of a fault. + + Returns + ------- + np.ndarray of shape (4n x 3) + The corners of each fault segment in the fault, stacked vertically. + """ + return np.vstack([segment.corners() for segment in self.segments]) + + 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 segment in self.segments: + if segment.global_coordinates_in_segment(global_coordinates): + segment_coordinates = segment.global_coordinates_to_segment_coordinates( + global_coordinates + ) + strike_length = segment_coordinates[0] + 1 / 2 + dip_length = segment_coordinates[1] + 1 / 2 + return np.array( + [ + running_length + strike_length * segment.length - midpoint, + max(dip_length * segment.width, 0), + ] + ) + running_length += segment.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 segment in self.segments: + segment_length = segment.length + if remaining_length < segment_length: + return segment.segment_coordinates_to_global_coordinates( + np.array( + [ + remaining_length / segment_length - 1 / 2, + fault_coordinates[1] / segment.width - 1 / 2, + ] + ), + ) + remaining_length -= segment_length + raise ValueError("Specified fault coordinates out of bounds.") From de2663f800106b823206daa762aabc3fc5ee2e12 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Thu, 16 May 2024 15:01:49 +1200 Subject: [PATCH 04/39] Add nshmdb interface --- nshmdb/nshmdb.py | 170 +++++++++++++++++++++++++++++++++++++++++++++++ schema.sql | 19 +++--- 2 files changed, 180 insertions(+), 9 deletions(-) create mode 100644 nshmdb/nshmdb.py diff --git a/nshmdb/nshmdb.py b/nshmdb/nshmdb.py new file mode 100644 index 0000000..7013d04 --- /dev/null +++ b/nshmdb/nshmdb.py @@ -0,0 +1,170 @@ +""" +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 sqlite3 +from pathlib import Path +from sqlite3 import Connection + +import numpy as np + +import nshmdb.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 connection(self) -> Connection: + """Establish a connection to the SQLite database. + + Returns + ------- + Connection + """ + return sqlite3.connect(self.db_filepath) + + def insert_fault(self, fault_id: int, parent_id: int, fault: Fault): + """Insert fault data into the database. + + Parameters + ---------- + fault_id : int + ID of the fault. + parent_id : int + ID of the parent fault. + fault : Fault + Fault object containing fault geometry. + """ + with self.connection() as conn: + conn.execute( + """INSERT INTO fault (fault_id, name, parent_id) VALUES (?, ?, ?)""", + (fault_id, fault.name, parent_id), + ) + for segment in fault.segments: + conn.execute( + """INSERT INTO fault_segment ( + top_left_lat, + top_left_lon, + top_right_lat, + top_right_lon, + bottom_right_lat, + bottom_right_lon, + top_depth, + bottom_depth, + rake, + fault_id + ) VALUES ( + ?, ?, ?, ?, ?, ?, ?, ?, ?, ? + )""", + *segment.corners[:, :2].ravel(), + segment.corners[0, 2], + segment.corners[-1, 2] + ) + + def insert_rupture(self, rupture_id: int, fault_ids: list[int]): + """Insert rupture data into the database. + + Parameters + ---------- + rupture_id : int + ID of the rupture. + fault_ids : list[int] + List of faults involved in the rupture. + """ + with self.connection() as conn: + conn.execute("INSERT INTO rupture (rupture_id) VALUES (?)", (rupture_id,)) + for fault_id in fault_ids: + conn.execute( + "INSERT INTO rupture_faults VALUES (?, ?)", (rupture_id, fault_id) + ) + + 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_segment 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_segments = 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_segments: + if parent_id != cur_parent_id: + faults.append( + Fault( + name=parent_name, + tect_type=None, + segments=[], + ) + ) + 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_right_lon, bottom], + ] + ) + faults[-1].segments.append(nshmdb.fault.FaultSegment(corners, rake)) + return faults diff --git a/schema.sql b/schema.sql index 9837b3e..4d431fc 100644 --- a/schema.sql +++ b/schema.sql @@ -13,16 +13,17 @@ CREATE TABLE IF NOT EXISTS parent_fault ( CREATE TABLE IF NOT EXISTS fault_segment ( segment_id INTEGER PRIMARY KEY, - strike REAL NOT NULL, + 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) ); From 60a2a10e8c0054fb16fb085b9f306d24f4ef18bb Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Thu, 16 May 2024 16:57:08 +1200 Subject: [PATCH 05/39] add insert parent method --- nshmdb/nshmdb.py | 16 ++++++++++++++++ schema.sql => nshmdb/schema/schema.sql | 0 2 files changed, 16 insertions(+) rename schema.sql => nshmdb/schema/schema.sql (100%) diff --git a/nshmdb/nshmdb.py b/nshmdb/nshmdb.py index 7013d04..e73558c 100644 --- a/nshmdb/nshmdb.py +++ b/nshmdb/nshmdb.py @@ -53,6 +53,22 @@ def connection(self) -> Connection: """ return sqlite3.connect(self.db_filepath) + def insert_parent(self, parent_id: int, parent_name: str): + """Insert parent fault data into the database. + + Parameters + ---------- + parent_id : int + ID of the parent fault. + name : str + Name of the parent fault. + """ + with self.connection() as conn: + conn.execute( + """INSERT OR REPLACE INTO parent_fault (parent_id, name) VALUES (?, ?)""", + (parent_id, parent_name), + ) + def insert_fault(self, fault_id: int, parent_id: int, fault: Fault): """Insert fault data into the database. diff --git a/schema.sql b/nshmdb/schema/schema.sql similarity index 100% rename from schema.sql rename to nshmdb/schema/schema.sql From 229f2994f80e48faeb768b88d2b5b03caedb472d Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Thu, 16 May 2024 16:57:18 +1200 Subject: [PATCH 06/39] add schema creation method --- nshmdb/nshmdb.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/nshmdb/nshmdb.py b/nshmdb/nshmdb.py index e73558c..e28cc0c 100644 --- a/nshmdb/nshmdb.py +++ b/nshmdb/nshmdb.py @@ -22,13 +22,14 @@ """ import dataclasses +import importlib.resources import sqlite3 from pathlib import Path from sqlite3 import Connection import numpy as np -import nshmdb.fault +from nshmdb import fault from nshmdb.fault import Fault @@ -44,6 +45,14 @@ class NSHMDB: db_filepath: Path + def create(self): + 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. From a374553a21b68fe26016ee1bf3587b82bd9b7d9f Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Thu, 16 May 2024 16:57:36 +1200 Subject: [PATCH 07/39] fix segment reference --- nshmdb/nshmdb.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nshmdb/nshmdb.py b/nshmdb/nshmdb.py index e28cc0c..89221e6 100644 --- a/nshmdb/nshmdb.py +++ b/nshmdb/nshmdb.py @@ -191,5 +191,5 @@ def get_rupture_faults(self, rupture_id: int) -> list[Fault]: [bottom_left_lat, bottom_right_lon, bottom], ] ) - faults[-1].segments.append(nshmdb.fault.FaultSegment(corners, rake)) + faults[-1].segments.append(fault.FaultSegment(corners, rake)) return faults From 12902248bad54f1a728f0ad1ea040e9a5e49ea81 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Thu, 16 May 2024 16:57:47 +1200 Subject: [PATCH 08/39] begin script refactor --- nshmdb/scripts/nshm_db_generator.py | 295 ++++++++-------------------- pyproject.toml | 2 + 2 files changed, 84 insertions(+), 213 deletions(-) diff --git a/nshmdb/scripts/nshm_db_generator.py b/nshmdb/scripts/nshm_db_generator.py index ced1aa2..4518125 100644 --- a/nshmdb/scripts/nshm_db_generator.py +++ b/nshmdb/scripts/nshm_db_generator.py @@ -1,236 +1,105 @@ #!/usr/bin/env python3 import csv -import json +import functools import sqlite3 +import zipfile from pathlib import Path -from sqlite3 import Connection -from typing import Annotated, Any +from typing import Annotated +import geojson +import nshmdb.fault import numpy as np +import pandas as pd import qcore.geo import typer +from geojson import FeatureCollection +from nshmdb.fault import Fault +from nshmdb.nshmdb import NSHMDB 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), - ) +POLYGONS_PATH = Path("ruptures") / "sect_polygons.geojson" +FAULT_INFORMATION_PATH = Path("ruptures") / "fault_sections.geojson" +RUPTURE_FAULT_JOIN_PATH = Path("ruptures") / "fast_indices.csv" + + +def extract_faults_from_polygons( + fault_info_list: FeatureCollection, +) -> list[Fault]: + faults = [] + for fault_feature in fault_info_list: + fault_trace = geojson.utils.coords(fault_feature) + name = fault_feature.properties["FaultName"] + dip_dir = fault_feature.properties["DipDir"] + dip = fault_feature.properties["DipDeg"] + width = fault_feature.properties["Width"] + projected_width = width * np.cos(np.radians(dip)) + bottom = fault_feature.properties["LowDepth"] + rake = fault_feature.properties["Rake"] + + segments = [] + for i in range(len(fault_trace) - 1): + top_left = fault_trace[i][::-1] + top_right = fault_trace[i + 1][::-1] + bottom_left = qcore.geo.ll_shift(*top_left, dip_dir, projected_width) + bottom_right = qcore.geo.ll_shift(*top_right, dip_dir, projected_width) + corners = np.array([top_left, top_right, bottom_right, bottom_left]) + corners = np.append(corners, np.array([0, 0, bottom, bottom]), axis=1) + segments.append(nshmdb.fault.FaultSegment(corners, rake)) + faults.append(Fault(name, None, segments)) + return faults + + +def insert_rupture(db: NSHMDB, rupture_df: pd.DataFrame): + rupture_id = rupture_df["rupture"].iloc[0] + rupture_faults = rupture_df["section"].to_list() + db.insert_rupture(rupture_id, rupture_faults) @app.command() def main( - fault_sections_geojson_filepath: Annotated[ + cru_solutions_zip_path: 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", + typer.Argument( + help="CRU solutions zip file", readable=True, dir_ok=False, exists=True + ), + ], sqlite_db_path: Annotated[ - Path, typer.Option(help="Output SQLite DB path", writable=True, exists=True) - ] = "nshm2022.db", + Path, typer.Argument(help="Output SQLite DB path", writable=True, exists=True) + ], ): """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) + + db = NSHMDB(sqlite_db_path) + db.create() + + with zipfile.ZipFile(cru_solutions_zip_path, "r") as cru_solutions_zip_file: + + 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) + + 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 i, fault in enumerate(faults): + fault_info = faults_info[i] + parent_id = fault_info.proprties["ParentID"] + db.insert_parent(parent_id, fault_info.properties["ParentName"]) + db.insert_fault(fault_info.properties["FaultID"], parent_id, fault) + + rupture_fault_join_df.groupby("rupture").apply( + functools.partial(insert_rupture, db=db) + ) if __name__ == "__main__": diff --git a/pyproject.toml b/pyproject.toml index 30cb685..74e95a6 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -9,6 +9,8 @@ readme = "README.md" python = "^3.12" qcore = { git = "https://github.com/ucgmsim/qcore" } numpy = "^1.26.4" +geojson = "^3.1.0" +pandas = "^2.2.2" [tool.poetry.group.dev.dependencies] python-lsp-server = "*" From da52082cb50e2803bd1472c4df4c5c34da025a51 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Fri, 17 May 2024 10:19:07 +1200 Subject: [PATCH 09/39] s/segment/plane --- nshmdb/fault.py | 142 ++++++++++++++-------------- nshmdb/scripts/nshm_db_generator.py | 6 +- 2 files changed, 74 insertions(+), 74 deletions(-) diff --git a/nshmdb/fault.py b/nshmdb/fault.py index 507e146..20f7141 100644 --- a/nshmdb/fault.py +++ b/nshmdb/fault.py @@ -1,6 +1,6 @@ -"""Module for representing fault segments and faults. +"""Module for representing fault planes and faults. -This module provides classes and functions for representing fault segments and +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. @@ -9,11 +9,11 @@ TectType: An enumeration of all the different kinds of fault types. -FaultSegment: - A representation of a single segment of a Fault. +FaultPlane: + A representation of a single plane of a Fault. Fault: - A representation of a fault, consisting of one or more FaultSegments. + A representation of a fault, consisting of one or more FaultPlanes. """ import dataclasses @@ -37,17 +37,17 @@ class TectType(Enum): @dataclasses.dataclass -class FaultSegment: - """A representation of a single segment of a Fault. +class FaultPlane: + """A representation of a single plane of a Fault. - This class represents a single segment of a fault, providing various + 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 segment. + The rake angle of the fault plane. """ _corners: np.ndarray @@ -64,7 +64,7 @@ def corners(self) -> np.ndarray: Returns ------- np.ndarray - The corners of the fault segment in (lat, lon, depth) format. + The corners of the fault plane in (lat, lon, depth) format. """ return qcore.coordinates.nztm_to_wgs_depth(self._corners) @@ -74,7 +74,7 @@ def length_m(self) -> float: Returns ------- float - The length of the fault segment (in metres). + The length of the fault plane (in metres). """ return np.linalg.norm(self._corners[1] - self._corners[0]) @@ -84,7 +84,7 @@ def width_m(self) -> float: Returns ------- float - The width of the fault segment (in metres). + The width of the fault plane (in metres). """ return np.linalg.norm(self._corners[-1] - self._corners[0]) @@ -104,7 +104,7 @@ def width(self) -> float: Returns ------- float - The width of the fault segment (in kilometres). + The width of the fault plane (in kilometres). """ return self.width_m / _KM_TO_M @@ -114,7 +114,7 @@ def length(self) -> float: Returns ------- float - The length of the fault segment (in kilometres). + The length of the fault plane (in kilometres). """ return self.length_m / _KM_TO_M @@ -124,7 +124,7 @@ def projected_width_m(self) -> float: Returns ------- float - The projected width of the fault segment (in metres). + The projected width of the fault plane (in metres). """ return self.length_m * np.cos(np.radians(self.dip)) @@ -134,7 +134,7 @@ def projected_width(self) -> float: Returns ------- float - The projected width of the fault segment (in kilometres). + The projected width of the fault plane (in kilometres). """ return self.projected_width / _KM_TO_M @@ -191,19 +191,19 @@ def dip(self) -> float: """ return np.degrees(np.arcsin(np.abs(self.bottom_m) / self.width_m)) - def segment_coordinates_to_global_coordinates( - self, segment_coordinates: np.ndarray + def plane_coordinates_to_global_coordinates( + self, plane_coordinates: np.ndarray ) -> np.ndarray: - """Convert segment coordinates to nztm global coordinates. + """Convert plane coordinates to nztm global coordinates. Parameters ---------- - segment_coordinates : np.ndarray - Segment coordinates to convert. Segment coordinates are - 2D coordinates (x, y) given for a fault segment (a plane), where x + 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 - origin for segment coordinates is the centre of the fault. + origin for plane coordinates is the centre of the fault. +x -1/2,-1/2 ─────────────────> @@ -228,17 +228,17 @@ def segment_coordinates_to_global_coordinates( offset = np.array([1 / 2, 1 / 2]) return qcore.coordinates.nztm_to_wgs_depth( - origin + (segment_coordinates + offset) @ frame + origin + (plane_coordinates + offset) @ frame ) - def global_coordinates_to_segment_coordinates( + def global_coordinates_to_plane_coordinates( self, global_coordinates: np.ndarray, ) -> np.ndarray: - """Convert coordinates (lat, lon, depth) to segment coordinates (x, y). + """Convert coordinates (lat, lon, depth) to plane coordinates (x, y). - See segment_coordinates_to_global_coordinates for a description of - segment coordinates. + See plane_coordinates_to_global_coordinates for a description of + plane coordinates. Parameters ---------- @@ -248,8 +248,8 @@ def global_coordinates_to_segment_coordinates( Returns ------- np.ndarray - The segment coordinates (x, y) representing the position of - global_coordinates on the fault segment. + The plane coordinates (x, y) representing the position of + global_coordinates on the fault plane. Raises ------ @@ -261,13 +261,13 @@ def global_coordinates_to_segment_coordinates( bottom_left = self._corners[-1] frame = np.vstack((top_right - origin, bottom_left - origin)) offset = qcore.coordinates.wgs_depth_to_nztm(global_coordinates) - origin - segment_coordinates, residual, _, _ = np.linalg.lstsq(frame.T, offset) + plane_coordinates, residual, _, _ = np.linalg.lstsq(frame.T, offset) if not np.isclose(residual[0], 0): raise ValueError("Coordinates do not lie in fault plane.") - return segment_coordinates - np.array([1 / 2, 1 / 2]) + return plane_coordinates - np.array([1 / 2, 1 / 2]) - def global_coordinates_in_segment(self, global_coordinates: np.ndarray) -> bool: - """Test if some global coordinates lie in the bounds of a segment. + def global_coordinates_in_plane(self, global_coordinates: np.ndarray) -> bool: + """Test if some global coordinates lie in the bounds of a plane. Parameters ---------- @@ -278,21 +278,21 @@ def global_coordinates_in_segment(self, global_coordinates: np.ndarray) -> bool: ------- bool True if the given global coordinates (lat, lon, depth) lie on the - fault segment. + fault plane. """ - segment_coordinates = self.global_coordinates_to_segment_coordinates( + plane_coordinates = self.global_coordinates_to_plane_coordinates( global_coordinates ) return np.all( np.logical_or( - np.abs(segment_coordinates) < 1 / 2, - np.isclose(np.abs(segment_coordinates), 1 / 2, atol=1e-4), + np.abs(plane_coordinates) < 1 / 2, + np.isclose(np.abs(plane_coordinates), 1 / 2, atol=1e-4), ) ) def centroid(self) -> np.ndarray: - """Returns the centre of the fault segment. + """Returns the centre of the fault plane. Returns ------- @@ -309,11 +309,11 @@ def centroid(self) -> np.ndarray: @dataclasses.dataclass class Fault: - """A representation of a fault, consisting of one or more FaultSegments. + """A representation of a fault, consisting of one or more FaultPlanes. - This class represents a fault, which is composed of one or more FaultSegments. + 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 segments, retrieving all corners of the fault, converting + 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. @@ -324,17 +324,17 @@ class Fault: The name of the fault. tect_type : str The type of fault this is (e.g. crustal, volcanic, subduction). - segments : list[FaultSegment] - A list containing all the FaultSegments that constitute the fault. + 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 segments. + Get the widths of all fault planes. lengths: - Get the lengths of all fault segments. + Get the lengths of all fault planes. corners: Get all corners of a fault. global_coordinates_to_fault_coordinates: @@ -345,7 +345,7 @@ class Fault: name: str tect_type: str - segments: list[FaultSegment] + planes: list[FaultPlane] def area(self) -> float: """Compute the area of a fault. @@ -355,27 +355,27 @@ def area(self) -> float: float The area of the fault. """ - return sum(segment.width * segment.length for segment in self.segments) + return sum(plane.width * plane.length for plane in self.planes) def widths(self) -> np.ndarray: - """Get the widths of all fault segments. + """Get the widths of all fault planes. Returns ------- np.ndarray of shape (1 x n) - The widths of all fault segments contained in this fault. + The widths of all fault planes contained in this fault. """ - return np.array([seg.width for seg in self.segments]) + return np.array([seg.width for seg in self.planes]) def lengths(self) -> np.ndarray: - """Get the lengths of all fault segments. + """Get the lengths of all fault planes. Returns ------- np.ndarray of shape (1 x n) - The lengths of all fault segments contained in this fault. + The lengths of all fault planes contained in this fault. """ - return np.array([seg.length for seg in self.segments]) + return np.array([seg.length for seg in self.planes]) def corners(self) -> np.ndarray: """Get all corners of a fault. @@ -383,9 +383,9 @@ def corners(self) -> np.ndarray: Returns ------- np.ndarray of shape (4n x 3) - The corners of each fault segment in the fault, stacked vertically. + The corners of each fault plane in the fault, stacked vertically. """ - return np.vstack([segment.corners() for segment in self.segments]) + return np.vstack([plane.corners for plane in self.planes]) def global_coordinates_to_fault_coordinates( self, global_coordinates: np.ndarray @@ -427,20 +427,20 @@ def global_coordinates_to_fault_coordinates( running_length = 0.0 midpoint = np.sum(self.lengths()) / 2 - for segment in self.segments: - if segment.global_coordinates_in_segment(global_coordinates): - segment_coordinates = segment.global_coordinates_to_segment_coordinates( + 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 = segment_coordinates[0] + 1 / 2 - dip_length = segment_coordinates[1] + 1 / 2 + strike_length = plane_coordinates[0] + 1 / 2 + dip_length = plane_coordinates[1] + 1 / 2 return np.array( [ - running_length + strike_length * segment.length - midpoint, - max(dip_length * segment.width, 0), + running_length + strike_length * plane.length - midpoint, + max(dip_length * plane.width, 0), ] ) - running_length += segment.length + running_length += plane.length raise ValueError("Specified coordinates not contained on fault.") def fault_coordinates_to_wgsdepth_coordinates( @@ -468,16 +468,16 @@ def fault_coordinates_to_wgsdepth_coordinates( """ midpoint = np.sum(self.lengths()) / 2 remaining_length = fault_coordinates[0] + midpoint - for segment in self.segments: - segment_length = segment.length - if remaining_length < segment_length: - return segment.segment_coordinates_to_global_coordinates( + for plane in self.planes: + plane_length = plane.length + if remaining_length < plane_length: + return plane.plane_coordinates_to_global_coordinates( np.array( [ - remaining_length / segment_length - 1 / 2, - fault_coordinates[1] / segment.width - 1 / 2, + remaining_length / plane_length - 1 / 2, + fault_coordinates[1] / plane.width - 1 / 2, ] ), ) - remaining_length -= segment_length + remaining_length -= plane_length raise ValueError("Specified fault coordinates out of bounds.") diff --git a/nshmdb/scripts/nshm_db_generator.py b/nshmdb/scripts/nshm_db_generator.py index 4518125..558c950 100644 --- a/nshmdb/scripts/nshm_db_generator.py +++ b/nshmdb/scripts/nshm_db_generator.py @@ -38,7 +38,7 @@ def extract_faults_from_polygons( bottom = fault_feature.properties["LowDepth"] rake = fault_feature.properties["Rake"] - segments = [] + planes = [] for i in range(len(fault_trace) - 1): top_left = fault_trace[i][::-1] top_right = fault_trace[i + 1][::-1] @@ -46,8 +46,8 @@ def extract_faults_from_polygons( bottom_right = qcore.geo.ll_shift(*top_right, dip_dir, projected_width) corners = np.array([top_left, top_right, bottom_right, bottom_left]) corners = np.append(corners, np.array([0, 0, bottom, bottom]), axis=1) - segments.append(nshmdb.fault.FaultSegment(corners, rake)) - faults.append(Fault(name, None, segments)) + planes.append(nshmdb.fault.FaultPlane(corners, rake)) + faults.append(Fault(name, None, planes)) return faults From 215fa67f636050a1e9f4b7bdb6f16d9c65dd6a8d Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Fri, 17 May 2024 10:20:02 +1200 Subject: [PATCH 10/39] fix type of tect_type --- nshmdb/fault.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/nshmdb/fault.py b/nshmdb/fault.py index 20f7141..c5cb82a 100644 --- a/nshmdb/fault.py +++ b/nshmdb/fault.py @@ -322,7 +322,7 @@ class Fault: ---------- name : str The name of the fault. - tect_type : str + 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. @@ -344,7 +344,7 @@ class Fault: """ name: str - tect_type: str + tect_type: TectType | None planes: list[FaultPlane] def area(self) -> float: From f7ccabe53324e8e0d381eb5359e6c19fbaf28027 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Fri, 17 May 2024 10:22:20 +1200 Subject: [PATCH 11/39] remove unused imports --- nshmdb/scripts/nshm_db_generator.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/nshmdb/scripts/nshm_db_generator.py b/nshmdb/scripts/nshm_db_generator.py index 558c950..b143b55 100644 --- a/nshmdb/scripts/nshm_db_generator.py +++ b/nshmdb/scripts/nshm_db_generator.py @@ -1,7 +1,5 @@ #!/usr/bin/env python3 -import csv import functools -import sqlite3 import zipfile from pathlib import Path from typing import Annotated From 76ac658b8a83cd188454b6d17c49c03064fe7cef Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Fri, 17 May 2024 10:22:29 +1200 Subject: [PATCH 12/39] fix function name --- nshmdb/scripts/nshm_db_generator.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nshmdb/scripts/nshm_db_generator.py b/nshmdb/scripts/nshm_db_generator.py index b143b55..45b30d8 100644 --- a/nshmdb/scripts/nshm_db_generator.py +++ b/nshmdb/scripts/nshm_db_generator.py @@ -22,7 +22,7 @@ RUPTURE_FAULT_JOIN_PATH = Path("ruptures") / "fast_indices.csv" -def extract_faults_from_polygons( +def extract_faults_from_info( fault_info_list: FeatureCollection, ) -> list[Fault]: faults = [] From 7d5242763febd66d585e1f47c22090365292041e Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Fri, 17 May 2024 10:22:35 +1200 Subject: [PATCH 13/39] add comment to explain section casting --- nshmdb/scripts/nshm_db_generator.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/nshmdb/scripts/nshm_db_generator.py b/nshmdb/scripts/nshm_db_generator.py index 45b30d8..5a1bc8c 100644 --- a/nshmdb/scripts/nshm_db_generator.py +++ b/nshmdb/scripts/nshm_db_generator.py @@ -85,6 +85,9 @@ def main( str(RUPTURE_FAULT_JOIN_PATH) ) as rupture_fault_join_handle: rupture_fault_join_df = pd.read_csv(rupture_fault_join_handle) + # The fast_indices.csv file has two columns, rupture and section. + # Both are ids but the section id is a float for some reason, so we + # need to cast to an integer. rupture_fault_join_df["section"] = rupture_fault_join_df["section"].astype( "Int64" ) From f33d38543220545ff77550190041b786dc7cf743 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Fri, 17 May 2024 17:11:24 +1200 Subject: [PATCH 14/39] s/segment/plane & replace reuse sqlite connection --- nshmdb/nshmdb.py | 137 +++++++++++++++++++++++++++------------ nshmdb/schema/schema.sql | 4 +- 2 files changed, 97 insertions(+), 44 deletions(-) diff --git a/nshmdb/nshmdb.py b/nshmdb/nshmdb.py index 89221e6..aace944 100644 --- a/nshmdb/nshmdb.py +++ b/nshmdb/nshmdb.py @@ -62,27 +62,37 @@ def connection(self) -> Connection: """ return sqlite3.connect(self.db_filepath) - def insert_parent(self, parent_id: int, parent_name: str): + # 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. """ - with self.connection() as conn: - conn.execute( - """INSERT OR REPLACE INTO parent_fault (parent_id, name) VALUES (?, ?)""", - (parent_id, parent_name), - ) - - def insert_fault(self, fault_id: int, parent_id: int, fault: 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 @@ -90,48 +100,89 @@ def insert_fault(self, fault_id: int, parent_id: int, fault: Fault): fault : Fault Fault object containing fault geometry. """ - with self.connection() as conn: + 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 (fault_id, name, parent_id) VALUES (?, ?, ?)""", - (fault_id, fault.name, parent_id), + """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, + ), ) - for segment in fault.segments: - conn.execute( - """INSERT INTO fault_segment ( - top_left_lat, - top_left_lon, - top_right_lat, - top_right_lon, - bottom_right_lat, - bottom_right_lon, - top_depth, - bottom_depth, - rake, - fault_id - ) VALUES ( - ?, ?, ?, ?, ?, ?, ?, ?, ?, ? - )""", - *segment.corners[:, :2].ravel(), - segment.corners[0, 2], - segment.corners[-1, 2] - ) - def insert_rupture(self, rupture_id: int, fault_ids: list[int]): + 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) -> Fault: with self.connection() as conn: - conn.execute("INSERT INTO rupture (rupture_id) VALUES (?)", (rupture_id,)) - for fault_id in fault_ids: - conn.execute( - "INSERT INTO rupture_faults VALUES (?, ?)", (rupture_id, fault_id) + 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(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. @@ -148,7 +199,7 @@ def get_rupture_faults(self, rupture_id: int) -> list[Fault]: cursor = conn.cursor() cursor.execute( """SELECT fs.*, p.parent_id, p.name - FROM fault_segment fs + 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 @@ -156,10 +207,11 @@ def get_rupture_faults(self, rupture_id: int) -> list[Fault]: ORDER BY f.parent_id""", (rupture_id,), ) - fault_segments = cursor.fetchall() + fault_planes = cursor.fetchall() cur_parent_id = None faults = [] for ( + _, top_left_lat, top_left_lon, top_right_lat, @@ -171,15 +223,16 @@ def get_rupture_faults(self, rupture_id: int) -> list[Fault]: top, bottom, rake, + _, parent_id, parent_name, - ) in fault_segments: + ) in fault_planes: if parent_id != cur_parent_id: faults.append( Fault( name=parent_name, tect_type=None, - segments=[], + planes=[], ) ) cur_parent_id = parent_id @@ -188,8 +241,8 @@ def get_rupture_faults(self, rupture_id: int) -> list[Fault]: [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_right_lon, bottom], + [bottom_left_lat, bottom_left_lon, bottom], ] ) - faults[-1].segments.append(fault.FaultSegment(corners, rake)) + faults[-1].planes.append(fault.FaultPlane(corners, rake)) return faults diff --git a/nshmdb/schema/schema.sql b/nshmdb/schema/schema.sql index 4d431fc..a566c19 100644 --- a/nshmdb/schema/schema.sql +++ b/nshmdb/schema/schema.sql @@ -11,8 +11,8 @@ CREATE TABLE IF NOT EXISTS parent_fault ( name TEXT NOT NULL ); -CREATE TABLE IF NOT EXISTS fault_segment ( - segment_id INTEGER PRIMARY KEY, +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, From 38eaebb55174afad6eb29979f18bce09e832b9af Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Fri, 17 May 2024 17:12:30 +1200 Subject: [PATCH 15/39] speed up db gen --- nshmdb/scripts/nshm_db_generator.py | 42 +++++++++++++++-------------- 1 file changed, 22 insertions(+), 20 deletions(-) diff --git a/nshmdb/scripts/nshm_db_generator.py b/nshmdb/scripts/nshm_db_generator.py index 5a1bc8c..914dc3e 100644 --- a/nshmdb/scripts/nshm_db_generator.py +++ b/nshmdb/scripts/nshm_db_generator.py @@ -81,26 +81,28 @@ def main( faults = extract_faults_from_info(faults_info) - 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) - # The fast_indices.csv file has two columns, rupture and section. - # Both are ids but the section id is a float for some reason, so we - # need to cast to an integer. - rupture_fault_join_df["section"] = rupture_fault_join_df["section"].astype( - "Int64" - ) - - for i, fault in enumerate(faults): - fault_info = faults_info[i] - parent_id = fault_info.proprties["ParentID"] - db.insert_parent(parent_id, fault_info.properties["ParentName"]) - db.insert_fault(fault_info.properties["FaultID"], parent_id, fault) - - rupture_fault_join_df.groupby("rupture").apply( - functools.partial(insert_rupture, db=db) - ) + if not skip_faults_insertion: + 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__": From c772631a996b2c917e2b3fccf845ddc760b65da8 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Tue, 21 May 2024 11:14:33 +1200 Subject: [PATCH 16/39] add rupture plotting functionality --- nshmdb/plotting/rupture.py | 44 ++++++++++++++++++++++++++++++++++++++ pyproject.toml | 6 ++++++ 2 files changed, 50 insertions(+) create mode 100644 nshmdb/plotting/rupture.py diff --git a/nshmdb/plotting/rupture.py b/nshmdb/plotting/rupture.py new file mode 100644 index 0000000..85623dc --- /dev/null +++ b/nshmdb/plotting/rupture.py @@ -0,0 +1,44 @@ +""" +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. +""" + +import numpy as np +from nshmdb.fault import Fault +from pygmt_helper import plotting + + +def plot_rupture(title: str, faults: list[Fault]): + """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. + """ + 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.show() diff --git a/pyproject.toml b/pyproject.toml index 74e95a6..e8ef87d 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -11,6 +11,12 @@ qcore = { git = "https://github.com/ucgmsim/qcore" } numpy = "^1.26.4" geojson = "^3.1.0" pandas = "^2.2.2" +typer = "^0.12.3" +pyproj = "^3.6.1" +tqdm = "^4.66.4" +pygmt-helper = {git = "https://github.com/ucgmsim/pygmt_helper.git"} +pygmt = "^0.12.0" +geopandas = "^0.14.4" [tool.poetry.group.dev.dependencies] python-lsp-server = "*" From 11f83d6f05da0dbb6b80cb97319a4d56719bae3d Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Tue, 21 May 2024 11:23:59 +1200 Subject: [PATCH 17/39] Fix CLI scripts --- nshmdb/scripts/nshm_db_generator.py | 87 ++++++++++++++++++++++------- 1 file changed, 66 insertions(+), 21 deletions(-) diff --git a/nshmdb/scripts/nshm_db_generator.py b/nshmdb/scripts/nshm_db_generator.py index 914dc3e..40170e9 100644 --- a/nshmdb/scripts/nshm_db_generator.py +++ b/nshmdb/scripts/nshm_db_generator.py @@ -1,18 +1,42 @@ #!/usr/bin/env python3 -import functools +""" +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 +from nshmdb.fault import Fault +from nshmdb.nshmdb import NSHMDB + import geojson import nshmdb.fault import numpy as np import pandas as pd import qcore.geo +import tqdm import typer from geojson import FeatureCollection -from nshmdb.fault import Fault -from nshmdb.nshmdb import NSHMDB app = typer.Typer() @@ -25,54 +49,75 @@ 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 fault_feature in fault_info_list: - fault_trace = geojson.utils.coords(fault_feature) + 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"] - width = fault_feature.properties["Width"] - projected_width = width * np.cos(np.radians(dip)) 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 = fault_trace[i][::-1] top_right = fault_trace[i + 1][::-1] - bottom_left = qcore.geo.ll_shift(*top_left, dip_dir, projected_width) - bottom_right = qcore.geo.ll_shift(*top_right, dip_dir, projected_width) + bottom_left = qcore.geo.ll_shift(*top_left, projected_width, dip_dir) + bottom_right = qcore.geo.ll_shift(*top_right, projected_width, dip_dir) corners = np.array([top_left, top_right, bottom_right, bottom_left]) - corners = np.append(corners, np.array([0, 0, bottom, bottom]), axis=1) + corners = np.append( + corners, + np.array([0, 0, bottom * 1000, bottom * 1000]).reshape((-1, 1)), + axis=1, + ) planes.append(nshmdb.fault.FaultPlane(corners, rake)) faults.append(Fault(name, None, planes)) return faults -def insert_rupture(db: NSHMDB, rupture_df: pd.DataFrame): - rupture_id = rupture_df["rupture"].iloc[0] - rupture_faults = rupture_df["section"].to_list() - db.insert_rupture(rupture_id, rupture_faults) - - @app.command() def main( cru_solutions_zip_path: Annotated[ Path, typer.Argument( - help="CRU solutions zip file", readable=True, dir_ok=False, exists=True + 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, exists=True) + 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: + 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) @@ -81,7 +126,7 @@ def main( faults = extract_faults_from_info(faults_info) - if not skip_faults_insertion: + if not skip_faults_creation: for i, fault in enumerate(faults): fault_info = faults_info[i] parent_id = fault_info.properties["ParentID"] From d00556e27a0ccee73bd0053dff0d5ff7200282a0 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Tue, 21 May 2024 11:58:51 +1200 Subject: [PATCH 18/39] update strike and dip_dir calculation for new api --- nshmdb/fault.py | 18 ++++-------------- 1 file changed, 4 insertions(+), 14 deletions(-) diff --git a/nshmdb/fault.py b/nshmdb/fault.py index c5cb82a..9d1b6b3 100644 --- a/nshmdb/fault.py +++ b/nshmdb/fault.py @@ -151,13 +151,8 @@ 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] - return ( - np.degrees( - qcore.geo.oriented_angle_wrt_normal( - north_direction, strike_direction, up_direction - ) - ) - % 360 + return qcore.geo.oriented_bearing_wrt_normal( + north_direction, strike_direction, up_direction ) @property @@ -172,13 +167,8 @@ def dip_dir(self) -> float: up_direction = np.array([0, 0, 1]) dip_direction = self._corners[-1] - self._corners[0] dip_direction[-1] = 0 - return ( - np.degrees( - qcore.geo.oriented_angle_wrt_normal( - north_direction, dip_direction, up_direction - ) - ) - % 360 + return qcore.geo.oriented_bearing_wrt_normal( + north_direction, dip_direction, up_direction ) @property From 0f75bacd24c5dfc36be8ee7b288dbaef0c7e71dc Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Tue, 21 May 2024 12:01:38 +1200 Subject: [PATCH 19/39] add docstring for NSHMDB.get_fault --- nshmdb/nshmdb.py | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/nshmdb/nshmdb.py b/nshmdb/nshmdb.py index aace944..7795b29 100644 --- a/nshmdb/nshmdb.py +++ b/nshmdb/nshmdb.py @@ -151,7 +151,20 @@ def add_fault_to_rupture(self, conn: Connection, rupture_id: int, fault_id: int) (rupture_id, fault_id), ) - def get_fault(self, fault_id) -> Fault: + 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,)) From 9b9e6924b87b9760c5c1894b1af128769fe57ca0 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Tue, 21 May 2024 12:02:46 +1200 Subject: [PATCH 20/39] add docstring for NSHMDB.create --- nshmdb/nshmdb.py | 1 + 1 file changed, 1 insertion(+) diff --git a/nshmdb/nshmdb.py b/nshmdb/nshmdb.py index 7795b29..04358a4 100644 --- a/nshmdb/nshmdb.py +++ b/nshmdb/nshmdb.py @@ -46,6 +46,7 @@ class NSHMDB: 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: From c12f98de5352ddaa643abd62e820f80d8853d11f Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Tue, 21 May 2024 12:05:39 +1200 Subject: [PATCH 21/39] raise error if db does not exist --- nshmdb/nshmdb.py | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/nshmdb/nshmdb.py b/nshmdb/nshmdb.py index 04358a4..16091f2 100644 --- a/nshmdb/nshmdb.py +++ b/nshmdb/nshmdb.py @@ -45,6 +45,23 @@ class NSHMDB: db_filepath: Path + def __init__(self, db_filepath: Path): + """Create a NSHMDB object + + Parameters + ---------- + db_filepath : Path + The path to the nshm database. + + Raises + ------ + ValueError + If the database at db_filepath does not exist. + """ + if not db_filepath.exists(): + raise ValueError(f"DB filepath {str(db_filepath)} does not exist.") + self.db_filepath = db_filepath + def create(self): """Create the tables for the NSHMDB database.""" schema_traversable = importlib.resources.files("nshmdb.schema") / "schema.sql" From 87f7b715016615f81eb4a025e0378c388fbd8513 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Tue, 21 May 2024 16:38:23 +1200 Subject: [PATCH 22/39] sensitivity tweaks for fault coordinates --- nshmdb/fault.py | 25 ++++++++++++++----------- 1 file changed, 14 insertions(+), 11 deletions(-) diff --git a/nshmdb/fault.py b/nshmdb/fault.py index 9d1b6b3..b698ec2 100644 --- a/nshmdb/fault.py +++ b/nshmdb/fault.py @@ -251,10 +251,10 @@ def global_coordinates_to_plane_coordinates( bottom_left = self._corners[-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) - if not np.isclose(residual[0], 0): + 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 plane_coordinates - np.array([1 / 2, 1 / 2]) + 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. @@ -271,15 +271,18 @@ def global_coordinates_in_plane(self, global_coordinates: np.ndarray) -> bool: fault plane. """ - 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-4), + 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. From 7be42c7f903f64de727ec4306ab1692297321a0d Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Wed, 22 May 2024 11:24:39 +1200 Subject: [PATCH 23/39] convert to nztm --- nshmdb/fault.py | 6 ++--- nshmdb/nshmdb.py | 26 ++++++---------------- nshmdb/scripts/nshm_db_generator.py | 34 ++++++++++++++++++----------- 3 files changed, 30 insertions(+), 36 deletions(-) diff --git a/nshmdb/fault.py b/nshmdb/fault.py index b698ec2..7b1b329 100644 --- a/nshmdb/fault.py +++ b/nshmdb/fault.py @@ -53,10 +53,6 @@ class FaultPlane: _corners: np.ndarray rake: float - def __init__(self, corners: np.ndarray, rake: float): - self._corners = qcore.coordinates.wgs_depth_to_nztm(corners) - self.rake = rake - @property def corners(self) -> np.ndarray: """ @@ -163,6 +159,8 @@ def dip_dir(self) -> float: 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[-1] - self._corners[0] diff --git a/nshmdb/nshmdb.py b/nshmdb/nshmdb.py index 16091f2..bcd46fd 100644 --- a/nshmdb/nshmdb.py +++ b/nshmdb/nshmdb.py @@ -28,6 +28,7 @@ from sqlite3 import Connection import numpy as np +import qcore.coordinates from nshmdb import fault from nshmdb.fault import Fault @@ -45,23 +46,6 @@ class NSHMDB: db_filepath: Path - def __init__(self, db_filepath: Path): - """Create a NSHMDB object - - Parameters - ---------- - db_filepath : Path - The path to the nshm database. - - Raises - ------ - ValueError - If the database at db_filepath does not exist. - """ - if not db_filepath.exists(): - raise ValueError(f"DB filepath {str(db_filepath)} does not exist.") - self.db_filepath = db_filepath - def create(self): """Create the tables for the NSHMDB database.""" schema_traversable = importlib.resources.files("nshmdb.schema") / "schema.sql" @@ -210,7 +194,9 @@ def get_fault(self, fault_id: int) -> Fault: [bottom_left_lat, bottom_left_lon, bottom], ] ) - planes.append(fault.FaultPlane(corners, rake)) + 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) @@ -275,5 +261,7 @@ def get_rupture_faults(self, rupture_id: int) -> list[Fault]: [bottom_left_lat, bottom_left_lon, bottom], ] ) - faults[-1].planes.append(fault.FaultPlane(corners, rake)) + faults[-1].planes.append( + fault.FaultPlane(qcore.coordinates.wgs_depth_to_nztm(corners), rake) + ) return faults diff --git a/nshmdb/scripts/nshm_db_generator.py b/nshmdb/scripts/nshm_db_generator.py index 40170e9..30394d7 100644 --- a/nshmdb/scripts/nshm_db_generator.py +++ b/nshmdb/scripts/nshm_db_generator.py @@ -26,17 +26,16 @@ from pathlib import Path from typing import Annotated -from nshmdb.fault import Fault -from nshmdb.nshmdb import NSHMDB - import geojson import nshmdb.fault import numpy as np import pandas as pd -import qcore.geo +import qcore.coordinates import tqdm import typer from geojson import FeatureCollection +from nshmdb.fault import Fault +from nshmdb.nshmdb import NSHMDB app = typer.Typer() @@ -76,16 +75,25 @@ def extract_faults_from_info( rake = fault_feature.properties["Rake"] planes = [] for i in range(len(fault_trace) - 1): - top_left = fault_trace[i][::-1] - top_right = fault_trace[i + 1][::-1] - bottom_left = qcore.geo.ll_shift(*top_left, projected_width, dip_dir) - bottom_right = qcore.geo.ll_shift(*top_right, projected_width, dip_dir) - corners = np.array([top_left, top_right, bottom_right, bottom_left]) - corners = np.append( - corners, - np.array([0, 0, bottom * 1000, bottom * 1000]).reshape((-1, 1)), - axis=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 From 9c010a00bbd40badb62e2b9009a41df337672222 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Wed, 22 May 2024 16:48:29 +1200 Subject: [PATCH 24/39] fix boundary case fault coord -> global coord bug --- nshmdb/fault.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/nshmdb/fault.py b/nshmdb/fault.py index 7b1b329..d3d1ebb 100644 --- a/nshmdb/fault.py +++ b/nshmdb/fault.py @@ -461,7 +461,9 @@ def fault_coordinates_to_wgsdepth_coordinates( remaining_length = fault_coordinates[0] + midpoint for plane in self.planes: plane_length = plane.length - if remaining_length < plane_length: + if remaining_length < plane_length or np.isclose( + remaining_length, plane_length + ): return plane.plane_coordinates_to_global_coordinates( np.array( [ From 90cf1e488395c4d1197cce8e83045e0672d26250 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Wed, 22 May 2024 16:48:48 +1200 Subject: [PATCH 25/39] downgrade used version of python --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index e8ef87d..6d54e6b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -6,7 +6,7 @@ authors = ["QuakeCoRE"] readme = "README.md" [tool.poetry.dependencies] -python = "^3.12" +python = "^3.11" qcore = { git = "https://github.com/ucgmsim/qcore" } numpy = "^1.26.4" geojson = "^3.1.0" From 139318c20ab784e22b03299b1b67cedc15a52915 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Wed, 29 May 2024 11:33:44 +1200 Subject: [PATCH 26/39] make _corners not private --- nshmdb/fault.py | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/nshmdb/fault.py b/nshmdb/fault.py index d3d1ebb..184902f 100644 --- a/nshmdb/fault.py +++ b/nshmdb/fault.py @@ -50,7 +50,7 @@ class FaultPlane: The rake angle of the fault plane. """ - _corners: np.ndarray + corners_nztm: np.ndarray rake: float @property @@ -62,7 +62,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 +72,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 +82,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 +92,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: @@ -146,7 +146,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 +163,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 @@ -209,9 +209,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 +244,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 +294,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() From 55deac60c3424af6a4c6f80e1d668c481f5a28ba Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Wed, 29 May 2024 11:34:05 +1200 Subject: [PATCH 27/39] improve docs --- nshmdb/fault.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/nshmdb/fault.py b/nshmdb/fault.py index 184902f..37416fd 100644 --- a/nshmdb/fault.py +++ b/nshmdb/fault.py @@ -189,8 +189,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 @@ -374,7 +374,8 @@ 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]) From 8be0ace7e3a32c779a57c92f4427396a81496881 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Wed, 29 May 2024 11:34:16 +1200 Subject: [PATCH 28/39] add fault-level nztm corner calculations --- nshmdb/fault.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/nshmdb/fault.py b/nshmdb/fault.py index 37416fd..8a28159 100644 --- a/nshmdb/fault.py +++ b/nshmdb/fault.py @@ -379,6 +379,16 @@ def corners(self) -> np.ndarray: """ 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: From f1a56e499882f5b698aaeeed2c2ead9213182d8a Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Wed, 29 May 2024 11:35:34 +1200 Subject: [PATCH 29/39] write plot to file instead of fig.show --- nshmdb/plotting/rupture.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/nshmdb/plotting/rupture.py b/nshmdb/plotting/rupture.py index 85623dc..e703ebd 100644 --- a/nshmdb/plotting/rupture.py +++ b/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) From 142afb32de9ace884799687e0b1750777abf9fc1 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Wed, 29 May 2024 11:39:06 +1200 Subject: [PATCH 30/39] use requirements.txt --- pyproject.toml | 45 --------------------------------------------- requirements.txt | 10 ++++++++++ 2 files changed, 10 insertions(+), 45 deletions(-) delete mode 100644 pyproject.toml create mode 100644 requirements.txt diff --git a/pyproject.toml b/pyproject.toml deleted file mode 100644 index 6d54e6b..0000000 --- a/pyproject.toml +++ /dev/null @@ -1,45 +0,0 @@ -[tool.poetry] -name = "nshmdb" -version = "0.1.0" -description = "Fault database for the GNS National Seismic Hazard Model." -authors = ["QuakeCoRE"] -readme = "README.md" - -[tool.poetry.dependencies] -python = "^3.11" -qcore = { git = "https://github.com/ucgmsim/qcore" } -numpy = "^1.26.4" -geojson = "^3.1.0" -pandas = "^2.2.2" -typer = "^0.12.3" -pyproj = "^3.6.1" -tqdm = "^4.66.4" -pygmt-helper = {git = "https://github.com/ucgmsim/pygmt_helper.git"} -pygmt = "^0.12.0" -geopandas = "^0.14.4" - -[tool.poetry.group.dev.dependencies] -python-lsp-server = "*" -python-lsp-black = "*" -python-lsp-isort = "*" -black = "*" -rope = "*" -pyflakes = "*" -mccabe = "*" -pycodestyle = "*" -doq = "*" -mypy = "^1.10.0" -pylsp-mypy = "^0.6.8" - -[build-system] -requires = ["poetry-core"] -build-backend = "poetry.core.masonry.api" - -[scripts] -nshm-db-generator = "nshmdb.scripts.nshm_db_generator:app" - -[tool.pylsp-mypy] -enabled = true -live_mode = true -strict = false -exclude = ["tests/*"] 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 From e7190a4c6566aa4b7e9ebd54068ecc094eeba135 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Wed, 29 May 2024 11:40:12 +1200 Subject: [PATCH 31/39] remove unused variable --- nshmdb/scripts/nshm_db_generator.py | 1 - 1 file changed, 1 deletion(-) diff --git a/nshmdb/scripts/nshm_db_generator.py b/nshmdb/scripts/nshm_db_generator.py index 30394d7..f845e22 100644 --- a/nshmdb/scripts/nshm_db_generator.py +++ b/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" From 23c560fc2a55810e29727c36ee515ba9574a4a94 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Wed, 29 May 2024 11:41:39 +1200 Subject: [PATCH 32/39] fix typo --- nshmdb/fault.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nshmdb/fault.py b/nshmdb/fault.py index 8a28159..e88eda8 100644 --- a/nshmdb/fault.py +++ b/nshmdb/fault.py @@ -132,7 +132,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: From a6adccad00dab56e207c745fe3d0becfd8b91373 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Wed, 29 May 2024 11:42:01 +1200 Subject: [PATCH 33/39] optimise imports --- nshmdb/fault.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/nshmdb/fault.py b/nshmdb/fault.py index e88eda8..a092313 100644 --- a/nshmdb/fault.py +++ b/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): From aa9f1ff921b8a2140610ea23fdda970a0dad2f84 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Thu, 30 May 2024 15:26:33 +1200 Subject: [PATCH 34/39] update names and docstrings --- nshmdb/fault.py | 38 +++++++++++++++++++++++++++----------- 1 file changed, 27 insertions(+), 11 deletions(-) diff --git a/nshmdb/fault.py b/nshmdb/fault.py index a092313..013f9fe 100644 --- a/nshmdb/fault.py +++ b/nshmdb/fault.py @@ -47,7 +47,23 @@ class FaultPlane: ---------- 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 @@ -55,13 +71,13 @@ class FaultPlane: @property def corners(self) -> np.ndarray: """ - Returns ------- np.ndarray - The corners of the fault plane in (lat, lon, depth) format. + The corners of the fault plane in (lat, lon, depth) format. The + corners are the same as in corners_nztm. """ - return qcore.coordinates.nztm_to_wgs_depth(self.corners_nztm) + return coordinates.nztm_to_wgs_depth(self.corners_nztm) @property def length_m(self) -> float: @@ -146,7 +162,7 @@ def strike(self) -> float: 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 qcore.geo.oriented_bearing_wrt_normal( + return geo.oriented_bearing_wrt_normal( north_direction, strike_direction, up_direction ) @@ -164,7 +180,7 @@ def dip_dir(self) -> float: up_direction = np.array([0, 0, 1]) dip_direction = self.corners_nztm[-1] - self.corners_nztm[0] dip_direction[-1] = 0 - return qcore.geo.oriented_bearing_wrt_normal( + return geo.oriented_bearing_wrt_normal( north_direction, dip_direction, up_direction ) @@ -195,9 +211,9 @@ def plane_coordinates_to_global_coordinates( +x -1/2,-1/2 ─────────────────> ┌─────────────────────┐ │ - │ < width > │ │ + │ < strike > │ │ │ ^ │ │ - │ length│ │ +y + │ dip │ │ +y │ v │ │ │ │ │ └─────────────────────┘ ∨ @@ -214,7 +230,7 @@ def plane_coordinates_to_global_coordinates( frame = np.vstack((top_right - origin, bottom_left - origin)) offset = np.array([1 / 2, 1 / 2]) - return qcore.coordinates.nztm_to_wgs_depth( + return coordinates.nztm_to_wgs_depth( origin + (plane_coordinates + offset) @ frame ) @@ -247,7 +263,7 @@ def global_coordinates_to_plane_coordinates( 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 + 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.") @@ -292,7 +308,7 @@ def centroid(self) -> np.ndarray: """ - return qcore.coordinates.nztm_to_wgs_depth( + return coordinates.nztm_to_wgs_depth( np.mean(self.corners_nztm, axis=0).reshape((1, -1)) ).ravel() From 03d84c0cc1dd62b03b648b1271babab3e13f5020 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Thu, 30 May 2024 15:27:52 +1200 Subject: [PATCH 35/39] black --- nshmdb/fault.py | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/nshmdb/fault.py b/nshmdb/fault.py index 013f9fe..bd29317 100644 --- a/nshmdb/fault.py +++ b/nshmdb/fault.py @@ -51,19 +51,19 @@ class FaultPlane: 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 - """ + + 0 1 + ┌──────────┐ + │ │ + │ │ + │ │ + │ │ + │ │ + │ │ + │ │ + └──────────┘ + 3 2 + """ corners_nztm: np.ndarray rake: float From 48ba6a8f48726c9d17009a796cbad86bde16d1d4 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Fri, 31 May 2024 14:07:46 +1200 Subject: [PATCH 36/39] use optional --- nshmdb/fault.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/nshmdb/fault.py b/nshmdb/fault.py index bd29317..4057259 100644 --- a/nshmdb/fault.py +++ b/nshmdb/fault.py @@ -20,6 +20,7 @@ from enum import Enum import numpy as np +from typing import Optional from qcore import coordinates, geo @@ -350,7 +351,7 @@ class Fault: """ name: str - tect_type: TectType | None + tect_type: Optional[TectType] planes: list[FaultPlane] def area(self) -> float: From debf646a94cd29c32cd18554c88a57eb390a8695 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Fri, 31 May 2024 15:58:02 +1200 Subject: [PATCH 37/39] update readme --- README.md | 69 ++++--------------------------------------------------- 1 file changed, 4 insertions(+), 65 deletions(-) diff --git a/README.md b/README.md index f999533..fc35cc2 100644 --- a/README.md +++ b/README.md @@ -63,70 +63,9 @@ 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. +*YOU DO NOT NEED TO EXTRACT THIS FILE ANYWHERE* +2. After cloning this repository and installing the depedencies, run the following ```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 - } - - 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" +python nshmdb/scripts/nshm_db_generator.py nshmdb.db ``` +This will take some time. From 75ab5ffa2bc72c6ba7b4478cb665b6011de33075 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Fri, 31 May 2024 15:58:43 +1200 Subject: [PATCH 38/39] spacing --- README.md | 1 + 1 file changed, 1 insertion(+) diff --git a/README.md b/README.md index fc35cc2..7850640 100644 --- a/README.md +++ b/README.md @@ -64,6 +64,7 @@ To integrate with PyCharm (if you use it) see [PyCharm's Documentation](https:// └── WARNING.md ``` *YOU DO NOT NEED TO EXTRACT THIS FILE ANYWHERE* + 2. After cloning this repository and installing the depedencies, run the following ```bash python nshmdb/scripts/nshm_db_generator.py nshmdb.db From f2f4fbc56602a92b01201ea417b51046ce96147e Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Fri, 31 May 2024 15:58:57 +1200 Subject: [PATCH 39/39] bold text --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 7850640..3ad7bd2 100644 --- a/README.md +++ b/README.md @@ -63,7 +63,7 @@ To integrate with PyCharm (if you use it) see [PyCharm's Documentation](https:// │   └── sub_seismo_on_fault_mfds.csv └── WARNING.md ``` -*YOU DO NOT NEED TO EXTRACT THIS FILE ANYWHERE* +**YOU DO NOT NEED TO EXTRACT THIS FILE ANYWHERE** 2. After cloning this repository and installing the depedencies, run the following ```bash