From 8cb956adf3492a766d30a378a34d4a2d389e150b Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Tue, 25 Jun 2024 15:50:17 +1200 Subject: [PATCH 01/14] add sources module --- requirements.txt | 3 + source_modelling/sources/sources.py | 714 ++++++++++++++++++++++++++++ wiki/Sources.md | 213 +++++++++ wiki/images/source_coordinates.svg | 449 +++++++++++++++++ wiki/images/sources.ipe | 378 +++++++++++++++ 5 files changed, 1757 insertions(+) create mode 100644 source_modelling/sources/sources.py create mode 100644 wiki/Sources.md create mode 100644 wiki/images/source_coordinates.svg create mode 100644 wiki/images/sources.ipe diff --git a/requirements.txt b/requirements.txt index e69de29..398ca04 100644 --- a/requirements.txt +++ b/requirements.txt @@ -0,0 +1,3 @@ +git+https://github.com/ucgmsim/qcore.git +numpy +scipy diff --git a/source_modelling/sources/sources.py b/source_modelling/sources/sources.py new file mode 100644 index 0000000..0a78b25 --- /dev/null +++ b/source_modelling/sources/sources.py @@ -0,0 +1,714 @@ +"""Module for representing the geometry seismic sources: point sources, fault planes and faults. + +This module provides classes and functions for representing fault planes and +faults, along with methods for calculating various properties such as +dimensions, orientation, and coordinate transformations. + +Classes +------- +PointSource: + A representation of a point source. + +Plane: + A representation of a single plane of a Fault. + +Fault: + A representation of a fault, consisting of one or more Planes. +""" + +import dataclasses +from typing import Optional, Protocol + +import numpy as np +import scipy as sp + +from qcore import coordinates, geo, grid + +_KM_TO_M = 1000 + + +@dataclasses.dataclass +class Point: + """A representation of point source.""" + + point_coordinates: np.ndarray + # used to approximate point source as a small planar patch (metres). + length_m: float + # The usual strike, dip, dip direction, etc cannot be calculated + # from a point source and so must be provided by the user. + strike: float + dip: float + dip_dir: float + + @property + def length(self) -> float: + """ + Returns + ------- + float + The length of the approximating planar patch (in kilometres). + """ + return self.length_m / _KM_TO_M + + @property + def width_m(self) -> float: + """ + Returns + ------- + float + The width of the approximating planar patch (in metres). + """ + return self.length_m + + @property + def width(self) -> float: + """ + Returns + ------- + float + The width of the approximating planar patch (in kilometres). + """ + return self.width_m / _KM_TO_M + + def fault_coordinates_to_wgs_depth_coordinates( + self, fault_coordinates: np.ndarray + ) -> np.ndarray: + """Convert fault-local coordinates to global (lat, lon, depth) coordinates. + + Parameters + ---------- + fault_coordinates : np.ndarray + The local fault coordinates + + Returns + ------- + np.ndarray + The global coordinates for these fault-local + coordinates. Because this is a point-source, the global + coordinates are just the location of the point source. + """ + + return self.point_coordinates + + def wgs_depth_coordinates_to_fault_coordinates( + self, wgs_depth_coordinates: np.ndarray + ) -> np.ndarray: + """Convert global coordinates into fault-local coordinates. + + Parameters + ---------- + wgs_depth_coordinates : np.ndarray + The global coordinates to convert. + + Returns + ------- + np.ndarray + The fault-local coordinates. Because this is a + point-source, the local coordinates are simply (1/2, 1/2) + near the source point and undefined elsewhere. + + Raises + ------ + ValueError + If the point is not near the source point. + """ + nztm_coordinates = coordinates.wgs_depth_to_nztm(wgs_depth_coordinates) + if np.all( + np.abs(nztm_coordinates - self.point_coordinates)[:2] / _KM_TO_M + < self.length + ): + return np.array([1 / 2, 1 / 2]) # Point is in the centre of the small patch + raise ValueError("Given global coordinates out of bounds for point source.") + + +@dataclasses.dataclass +class Plane: + """A representation of a single plane of a Fault. + + This class represents a single plane of a fault, providing various + properties and methods for calculating its dimensions, orientation, and + converting coordinates between different reference frames. + + Attributes + ---------- + 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 + + @property + def corners(self) -> np.ndarray: + """ + Returns + ------- + np.ndarray + The corners of the fault plane in (lat, lon, depth) format. The + corners are the same as in corners_nztm. + """ + return coordinates.nztm_to_wgs_depth(self.corners_nztm) + + @property + def length_m(self) -> float: + """ + Returns + ------- + float + The length of the fault plane (in metres). + """ + return np.linalg.norm(self.corners_nztm[1] - self.corners_nztm[0]) + + @property + def width_m(self) -> float: + """ + Returns + ------- + float + The width of the fault plane (in metres). + """ + return np.linalg.norm(self.corners_nztm[-1] - self.corners_nztm[0]) + + @property + def bottom_m(self) -> float: + """ + Returns + ------- + float + The bottom depth (in metres). + """ + return self.corners_nztm[-1, -1] + + @property + def width(self) -> float: + """ + Returns + ------- + float + The width of the fault plane (in kilometres). + """ + return self.width_m / _KM_TO_M + + @property + def length(self) -> float: + """ + Returns + ------- + float + The length of the fault plane (in kilometres). + """ + return self.length_m / _KM_TO_M + + @property + def projected_width_m(self) -> float: + """ + Returns + ------- + float + The projected width of the fault plane (in metres). + """ + return self.length_m * np.cos(np.radians(self.dip)) + + @property + def projected_width(self) -> float: + """ + Returns + ------- + float + The projected width of the fault plane (in kilometres). + """ + return self.projected_width_m / _KM_TO_M + + @property + def strike(self) -> float: + """ + Returns + ------- + float + The bearing of the strike direction of the fault + (from north; in degrees) + """ + + north_direction = np.array([1, 0, 0]) + up_direction = np.array([0, 0, 1]) + strike_direction = self.corners_nztm[1] - self.corners_nztm[0] + return geo.oriented_bearing_wrt_normal( + north_direction, strike_direction, up_direction + ) + + @property + def dip_dir(self) -> float: + """ + Returns + ------- + float + The bearing of the dip direction (from north; in degrees). + """ + if np.isclose(self.dip, 90): + return 0 # TODO: Is this right for this case? + north_direction = np.array([1, 0, 0]) + up_direction = np.array([0, 0, 1]) + dip_direction = self.corners_nztm[-1] - self.corners_nztm[0] + dip_direction[-1] = 0 + return geo.oriented_bearing_wrt_normal( + north_direction, dip_direction, up_direction + ) + + @property + def dip(self) -> float: + """ + Returns + ------- + float + The dip angle of the fault. + """ + return np.degrees(np.arcsin(np.abs(self.bottom_m) / self.width_m)) + + @staticmethod + def from_centroid_strike_dip( + centroid: np.ndarray, + strike: float, + dip_dir: Optional[float], + top: float, + bottom: float, + length: float, + width: float, + ) -> "Plane": + """Create a fault plane from the centroid, strike, dip_dir, top, bottom, length, and width + + This is used for older descriptions of sources. Internally + converts everything to corners so self.strike ~ strike (but + not exactly due to rounding errors). + + Parameters + ---------- + centroid : np.ndarray + The centre of the fault plane in lat, lon coordinate.s + strike : float + The strike of the fault (in degrees). + dip_dir : Optional[float] + The dip direction of the fault (in degrees). If None this is assumed to be strike + 90 degrees. + top : float + The top depth of the plane (in km). + bottom : float + The bottom depth of the plane (in km). + length : float + The length of the fault plane (in km). + width : float + The width of the fault plane (in km). + + Returns + ------- + Plane + The fault plane with centre at `centroid`, and where the + parameters strike, dip_dir, top, bottom, length and width + match what is passed to this function. + """ + corners = grid.grid_corners( + centroid, + strike, + dip_dir if dip_dir is not None else (strike + 90), + top, + bottom, + length, + width, + ) + return Plane(corners) + + def fault_coordinates_to_wgs_depth_coordinates( + self, plane_coordinates: np.ndarray + ) -> np.ndarray: + """Convert plane coordinates to nztm global coordinates. + + Parameters + ---------- + plane_coordinates : np.ndarray + Plane coordinates to convert. Plane coordinates are + 2D coordinates (x, y) given for a fault plane (a plane), where x + represents displacement along the strike, and y + displacement along the dip (see diagram below). The + origin for plane coordinates is the centre of the fault. + + +x + 0 0 ─────────────────> + ┌─────────────────────┐ │ + │ < strike > │ │ + │ ^ │ │ + │ dip │ │ +y + │ v │ │ + │ │ │ + └─────────────────────┘ ∨ + 1,1 + + Returns + ------- + np.ndarray + An 3d-vector of (lat, lon, depth) transformed coordinates. + """ + origin = self.corners_nztm[0] + top_right = self.corners_nztm[1] + bottom_left = self.corners_nztm[-1] + frame = np.vstack((top_right - origin, bottom_left - origin)) + + return coordinates.nztm_to_wgs_depth(origin + plane_coordinates @ frame) + + def wgs_depth_coordinates_to_fault_coordinates( + self, + global_coordinates: np.ndarray, + ) -> np.ndarray: + """Convert coordinates (lat, lon, depth) to plane coordinates (x, y). + + See plane_coordinates_to_global_coordinates for a description of + plane coordinates. + + Parameters + ---------- + global_coordinates : np.ndarray + Global coordinates to convert. + + Returns + ------- + np.ndarray + The plane coordinates (x, y) representing the position of + global_coordinates on the fault plane. + + Raises + ------ + ValueError + If the given coordinates do not lie in the fault plane. + """ + origin = self.corners_nztm[0] + top_right = self.corners_nztm[1] + bottom_left = self.corners_nztm[-1] + frame = np.vstack((top_right - origin, bottom_left - origin)) + offset = coordinates.wgs_depth_to_nztm(global_coordinates) - origin + plane_coordinates, residual, _, _ = np.linalg.lstsq(frame.T, offset, rcond=None) + if not np.isclose(residual[0], 0, atol=1e-02): + raise ValueError("Coordinates do not lie in fault plane.") + return np.clip(plane_coordinates, 0, 1) + + def wgs_depth_coordinates_in_plane(self, global_coordinates: np.ndarray) -> bool: + """Test if some global coordinates lie in the bounds of a plane. + + Parameters + ---------- + global_coordinates : np.ndarray + The global coordinates to check + + Returns + ------- + bool + True if the given global coordinates (lat, lon, depth) lie on the + fault plane. + """ + + try: + plane_coordinates = self.wgs_depth_coordinates_to_fault_coordinates( + global_coordinates + ) + return np.all( + np.logical_or( + np.abs(plane_coordinates) < 1 / 2, + np.isclose(np.abs(plane_coordinates), 1 / 2, atol=1e-3), + ) + ) + except ValueError: + return False + + def centroid(self) -> np.ndarray: + """Returns the centre of the fault plane. + + Returns + ------- + np.ndarray + A 1 x 3 dimensional vector representing the centroid of the fault + plane in (lat, lon, depth) format. + + """ + + return coordinates.nztm_to_wgs_depth( + np.mean(self.corners_nztm, axis=0).reshape((1, -1)) + ).ravel() + + +@dataclasses.dataclass +class Fault: + """A representation of a fault, consisting of one or more Planes. + + This class represents a fault, which is composed of one or more Planes. + It provides methods for computing the area of the fault, getting the widths and + lengths of all fault planes, retrieving all corners of the fault, converting + global coordinates to fault coordinates, converting fault coordinates to global + coordinates, generating a random hypocentre location within the fault, and + computing the expected fault coordinates. + + Attributes + ---------- + planes : list[Plane] + A list containing all the Planes that constitute the fault. + + Methods + ------- + area: + Compute the area of a fault. + widths: + Get the widths of all fault planes. + lengths: + Get the lengths of all fault planes. + corners: + Get all corners of a fault. + global_coordinates_to_fault_coordinates: + Convert global coordinates to fault coordinates. + fault_coordinates_to_wgsdepth_coordinates: + Convert fault coordinates to global coordinates. + """ + + planes: list[Plane] + + def area(self) -> float: + """Compute the area of a fault. + + Returns + ------- + float + The area of the fault. + """ + return self.width * np.sum(self.lengths) + + @property + def lengths(self) -> np.ndarray: + """The lengths of each plane in the fault. + + Returns + ------- + np.ndarray + A numpy array of each plane length (in km). + """ + return np.array([fault.length for fault in self.planes]) + + @property + def length(self) -> float: + """The length of the fault. + + Returns + ------- + float + The total length of each fault plane. + """ + + return self.lengths.sum() + + @property + def width(self) -> float: + """The width of the fault. + + Returns + ------- + float + The width of the first fault plane (A fault is assumed to + have planes of constant width). + """ + return self.planes[0].width + + @property + def dip_dir(self) -> float: + """The dip direction of the fault. + + Returns + ------- + float + The dip direction of the first fault plane (A fault is + assumed to have planes of constant dip direction). + """ + return self.planes[0].dip_dir + + def corners(self) -> np.ndarray: + """Get all corners of a fault. + + Returns + ------- + np.ndarray of shape (4n x 3) + The corners in (lat, lon, depth) format of each fault plane in the + fault, stacked vertically. + """ + return np.vstack([plane.corners for plane in self.planes]) + + def corners_nztm(self) -> np.ndarray: + """Get all corners of a fault. + + Returns + ------- + np.ndarray of shape (4n x 3) + The corners in NZTM format of each fault plane in the fault, stacked vertically. + """ + return np.vstack([plane.corners_nztm for plane in self.planes]) + + def wgs_depth_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 + from the top left, and d the distance from the top of the + fault (refer to the diagram). The coordinates are normalised + such that (0, 0) is the top left and (1, 1) the bottom right. + + (0, 0) + ┌──────────────────────┬──────┐ + │ | │ │ + │ | │ │ + │ | d │ │ + │ | │ │ + ├----------* │ │ + │ s ^ │ │ + │ | │ │ + │ | │ │ + │ | │ │ + └──────────|───────────┴──────┘ + + (1, 1) + 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. + + """ + # the right edges as a cumulative proportion of the fault length (e.g. [0.1, ..., 0.8]) + right_edges = self.lengths.cumsum() / self.length + right_edges = np.append(right_edges, 1) + for i, plane in enumerate(self.planes): + if plane.wgs_depth_coordinates_in_plane(global_coordinates): + plane_coordinates = plane.wgs_depth_coordinates_to_fault_coordinates( + global_coordinates + ) + return np.array([right_edges[i], 0]) + plane_coordinates * np.array( + [right_edges[i + 1] - right_edges[i], 1] + ) + raise ValueError("Given coordinates are not on fault.") + + def fault_coordinates_to_wgs_depth_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. + """ + + # the right edges as a cumulative proportion of the fault length (e.g. [0.1, ..., 0.8]) + right_edges = self.lengths.cumsum() / self.length + fault_segment_index = np.searchsorted(right_edges, fault_coordinates[0]) + left_proportion = ( + right_edges[fault_segment_index - 1] if fault_segment_index > 0 else 0 + ) + right_proportion = ( + right_edges[fault_segment_index + 1] + if fault_segment_index < len(right_edges) - 1 + else 1 + ) + segment_proportion = (fault_coordinates[0] - left_proportion) / ( + right_proportion - left_proportion + ) + return self.planes[ + fault_segment_index + ].fault_coordinates_to_wgs_depth_coordinates( + np.array([segment_proportion, fault_coordinates[1]]) + ) + + +class HasCoordinates(Protocol): + """Type definition for a source with local coordinates.""" + + def fault_coordinates_to_wgs_depth_coordinates( + self, + fault_coordinates: np.ndarray, + ) -> np.ndarray: ... + + def wgs_depth_coordinates_to_fault_coordinates( + self, + fault_coordinates: np.ndarray, + ) -> np.ndarray: ... + + +def closest_point_between_sources( + source_a: HasCoordinates, source_b: HasCoordinates +) -> tuple[np.ndarray, np.ndarray]: + """Find the closest point between two sources that have local coordinates. + + Parameters + ---------- + source_a : HasCoordinates + The first source. Must have a two-dimensional fault coordinate system. + source_b : HasCoordinates + The first source. Must have a two-dimensional fault coordinate system. + + Raises + ------ + ValueError + Raised when we are unable to converge on the closest points between sources. + + Returns + ------- + source_a_coordinates : np.ndarray + The source-local coordinates of the closest point on source a. + source_b_coordinates : np.ndarray + The source-local coordinates of the closest point on source b. + """ + + def fault_coordinate_distance(fault_coordinates: np.ndarray) -> float: + source_a_global_coordinates = ( + source_a.fault_coordinates_to_wgs_depth_coordinates(fault_coordinates[:2]) + ) + source_b_global_coordinates = ( + source_b.fault_coordinates_to_wgs_depth_coordinates(fault_coordinates[2:]) + ) + return coordinates.distance_between_wgs_depth_coordinates( + source_a_global_coordinates, source_b_global_coordinates + ) + + res = sp.optimize.minimize( + fault_coordinate_distance, + np.array([1 / 2, 1 / 2, 1 / 2, 1 / 2]), + bounds=[(0, 1)] * 4, + ) + + if not res.success: + raise ValueError( + f"Optimisation failed to converge for provided sources: {res.message}" + ) + + return res.x[:2], res.x[2:] diff --git a/wiki/Sources.md b/wiki/Sources.md new file mode 100644 index 0000000..8897d7b --- /dev/null +++ b/wiki/Sources.md @@ -0,0 +1,213 @@ +# Sources in QCore + +A number of packages need to deal with fault geometry including: + +- **Pre-processing**: Where fault geometry is used extensively to + compute fault discretisation and plane rupture propagation (in + type-5). +- **Velocity Model**: Where fault geometry is used to help approximate + the bounds of a suitable velocity model domain. +- **Visualisation Tools**: Where fault geometry is plotted to help + visualise ruptures or GMM model outputs. + +Over and over again we keep defining common routines for working with the geometry, answering questions like: + +1. How do I find the corners of a fault? +2. What are the closest points between two faults? +3. How do I compute an evenly spaced set of grid points along a fault? +4. What is the _insert property of fault plane here (dip\_dir, strike, dip, ...)_ of this fault? + +The constant redefinition and re-implementation leads to needless +repetition, bugs, and outdated assumptions. Examples of such +assumptions and mistakes include: + +1. Always assuming dip direction is 90 degrees from strike (it isn't in NSHM2022!), and +2. Always assuming the top depth of the fault is zero. + +The constant repetitive definitions, and the differing locations we +implement fault geometries begs for an implementation of source +geometry that we can all agree on in a central location. This is the +aim of the `qcore.sources` module. + +## Design Aims of the Sources Module + +The goal of the `qcore.sources` module is to provide a Pythonic +interface to source geometry. It must be general to accommodate new +kinds of sources, and it should not depend on naming conventions that +can change (type-I, type-II, etc). The definition should minimise +parameter redundancy. That is, instead of providing strike, dip, +length, width, bottom, top, we should let as many parameters as +possible be derived. The bottom and top depth values for example, can +be derived from strike, dip, length, and width values. In fact, all +the essential parameters can be found ideally from the supplied +corners of a fault plane in three dimensions. The fault corners, rather than the standard centroid-dip-strike-length-width-... specification we have use in the past, will now be the privileged information defining a fault. Everything else will be measured from the definition of the corners. This has a number of advantages over the old approach: + +1. It completely minimises parameter redundancy, and ensures that all + the paramaters are geometrically consistent. It will now be + impossible, for example, to specify a plane with inconsistent + length, width and depth parameters since these will be derived from + the corners of the fault. +2. Using corners allows us to frame problems of fault geometry as + problems in linear algebra. The advantage of this is that we can + take advantage of the wealth of tools available in numpy, scipy, + etc. In the past, we would write functions like `geo.ll2gp` and do + everything without any vectorisation. In the future we can use + matrix transformations to manipulate faults in an efficient and + concise manner. + +We also refrain from using inheritance hierarchies and factories to +define sources, instead using simple dataclasses and duck-typing. This +approach more closely matches Python development standards than, for +example, factories and other gang-of-four style patterns common to +Java. Accordingly, there is no Source superclass, and instead a +Protocol (like an interface) that defines the functions that should +exist for any object to be considered a source geometry. + +## What Is a Source Geometry? + +A *source geometry* is an object with two properties: + +1. A geometric definition of its bounds. For a fault plane, this is its corners, for a point-source the bounds are the point itself. +2. A _local coordinate system_. This local coordinate system is + essential for finding points inside the fault. EMOD3D has its own + definitions of fault-local coordinates, for example, which we pass + as `shyp` and `dhyp` parameters to `genslip` and friends. + +Note that this definition does not require the geometry to be flat +like a plane, or connected, or anything. It is simply a closed and +bounded region with coordinates. The choice of a general definition is +to allow for the flexible addition of sources to this framework, such +a rough surfaces. + +## Sources Used in Ground Motion Simulation + +While we have five types of sources (per [Source Modelling for GMSim](https://wiki.canterbury.ac.nz/display/QuakeCore/Source+Modelling+for+GMSim)), there are essentially only three source geometries we work with: + +1. **Point Source Geometry**: This is a 0-dimensional geometry consisting of a single point. The `qcore.sources` module uses the `Point` class to model the source geometry for a point. +2. **Plane Geometry**: This a 2-dimensional source geometry consisting of a single plane. The extents of the geometry are its corners. The `qcore.sources` module uses the `Plane` class to model single plane geometry. +3. **Multi-Planar Geometry**: This is a 2-dimensional source geometry consisting of a number of connected planes. The extents of the geometry are the corners of the end fault planes. The `qcore.sources` module uses the `Fault` class to model multi-planar geometry. + +Type-I fault are an instance of the the first geometry, type-2 faults +are plane geometries, and type-4 and type-5 are multi-planar +geometries. + +Note that the term *2-dimensional* here refers to the dimensions of +the local coordinate system, rather than their appearance as +"flat". Both planar and multi-planar geometry are 2-dimensional +because they can be given a local coordinate system with only two +parameters $(s, d)$, where $s$ is length along the strike direction +and $d$ length along the dip direction. Points are 0-dimensional +because their local coordinate system is just a single point. + +![A demonstration of the three different coordinate systems for sources.](images/source_coordinates.svg) + +To make the source module easy to use, we have elected to normalise +all the coordinate systems so that the coordinate systems are always +points $(s, d)$ where $0 \leq s \leq 1$ and $0 \leq d \leq 1$. Note +that this means the boundary of any geometry always corresponds to the +same set of points $\{(0, t)\,|\, 0 \leq t \leq 1\} \cup \{(s, 0)\,|\, +0 \leq s \leq 1\} \cup \{(1, t) \,|\, 0 \leq t \leq 1\} \cup \{(s, +1)\,|\, 0 \leq s \leq 1\}$. + + +Sources defined in `qcore.sources` will have two methods for converting back and forth between fault local and global coordinate systems: + +1. `fault_coordinates_to_wgs_depth_coordinates`: Going from fault-local coordinates to global coordinates. +2. `wgs_depth_coordinates_to_fault_coordinates`: Going from global + coordinates to fault-local coordinates **if the global coordinates + lie in the source geometry**. Sources will raise a `ValueError` if + the supplied coordinates are not in the domain. + +For functions that expect these coordinates to exist, the +`HasCoordinates` Protocol class allows you to require the existence of +these methods without specifying a superclass. + +```python +def minimum_depth(source: HasCoordinates) -> float: + return sp.optimize.minimize( + lambda x: source.fault_coordinates_to_wgs_depth_coordinates(x)[2], + np.array([1/2, 1/2]), + bounds=(0, 1) + ) +``` + +## Answering Geometry Questions with the Sources Module + +Below is a kind of cookbook demonstrating how to use the new sources module to answer source geometry questions. + +Q: How do I find the corners of a geometry? + +A: Assuming your geometry has corners (which it may not, if the geometry is rough), then the corners can be found simply as + +```python +source = Plane(...) # Or FaultPlane, or even Point! +corners = source.fault_coordinates_to_wgs_depth_coordinates(np.array([ + [0, 0], # top left + [1, 0], # top right + [1, 1], # bottom right + [0, 1] # bottom left +]])) +``` + +Q: How can I find the basic parameters of the geometry (strike, dip, rake, etc.)? + +A: A `Plane` source has these defined as properties computed from the corners you supply to construct the source +```python +plane = Plane(...) +strike = plane.strike +dip = plane.dip +dip_dir = plane.dip_dir +length = plane.length +length_metres = plane.length_m +``` +For some geometries, some of these values are not going to be +well-defined (what is the strike of a multi-planar geometry when it +changes?). + +Q: How can I discretise my geometry into a number of evenly spaced points? + +A: Here fault-local coordinates really shine because they make discretising sources extremely trivial: +```python +complicated_fault_geometry = FaultPlane(predefined_corners) +# define a meshgrid of points, 50 along the strike and 100 along the dip. +xv, yv = np.meshgrid(np.linspace(0, 1, num=50), np.linspace(0, 1, num=100)) +fault_local_meshgrid = np.vstack([xv.ravel(), yv.ravel()]) +# Convert the fault local coordinates into global coordinates +global_point_meshgrid = ( + complicated_fault_geometry.fault_coordinates_to_wgs_depth_coordinates( + fault_local_meshgrid + ) +) +``` + +Q: How can I find the closest points between two fault geometries? + +A: Again, having a local coordinate system turns this into a simple problem. Here is the straight-forward implementation of closest points in the `qcore.sources` module: + +```python +source_a = ... # some source geometry like FaultPlane +source_b = ... # some source geometry like Fault + +def fault_coordinate_distance(fault_coordinates: np.ndarray) -> float: + source_a_global_coordinates = ( + source_a.fault_coordinates_to_wgs_depth_coordinates(fault_coordinates[:2]) + ) + source_b_global_coordinates = ( + source_b.fault_coordinates_to_wgs_depth_coordinates(fault_coordinates[2:]) + ) + return coordinates.distance_between_wgs_depth_coordinates( + source_a_global_coordinates, source_b_global_coordinates + ) + +res = sp.optimize.minimize( + fault_coordinate_distance, + np.array([1 / 2, 1 / 2, 1 / 2, 1 / 2]), + bounds=[(0, 1)] * 4, +) + + +# Now the res.x value contains four elements +# [ s_a, d_a, s_b, d_b ] +# where (s_a, d_a) and (s_b, d_b) are the fault-local coordinates of +# the two closest points between source_a and source_b. +``` diff --git a/wiki/images/source_coordinates.svg b/wiki/images/source_coordinates.svg new file mode 100644 index 0000000..5bc2ceb --- /dev/null +++ b/wiki/images/source_coordinates.svg @@ -0,0 +1,449 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/wiki/images/sources.ipe b/wiki/images/sources.ipe new file mode 100644 index 0000000..0f1c982 --- /dev/null +++ b/wiki/images/sources.ipe @@ -0,0 +1,378 @@ + + + + + + + +0 0 m +-1 0.333 l +-1 -0.333 l +h + + + + +0 0 m +-1 0.333 l +-1 -0.333 l +h + + + + +0 0 m +-1 0.333 l +-0.8 0 l +-1 -0.333 l +h + + + + +0 0 m +-1 0.333 l +-0.8 0 l +-1 -0.333 l +h + + + + +0.6 0 0 0.6 0 0 e +0.4 0 0 0.4 0 0 e + + + + +0.6 0 0 0.6 0 0 e + + + + + +0.5 0 0 0.5 0 0 e + + +0.6 0 0 0.6 0 0 e +0.4 0 0 0.4 0 0 e + + + + + +-0.6 -0.6 m +0.6 -0.6 l +0.6 0.6 l +-0.6 0.6 l +h +-0.4 -0.4 m +0.4 -0.4 l +0.4 0.4 l +-0.4 0.4 l +h + + + + +-0.6 -0.6 m +0.6 -0.6 l +0.6 0.6 l +-0.6 0.6 l +h + + + + + +-0.5 -0.5 m +0.5 -0.5 l +0.5 0.5 l +-0.5 0.5 l +h + + +-0.6 -0.6 m +0.6 -0.6 l +0.6 0.6 l +-0.6 0.6 l +h +-0.4 -0.4 m +0.4 -0.4 l +0.4 0.4 l +-0.4 0.4 l +h + + + + + + +-0.43 -0.57 m +0.57 0.43 l +0.43 0.57 l +-0.57 -0.43 l +h + + +-0.43 0.57 m +0.57 -0.43 l +0.43 -0.57 l +-0.57 0.43 l +h + + + + + +0 0 m +-1 0.333 l +-1 -0.333 l +h + + + + +0 0 m +-1 0.333 l +-0.8 0 l +-1 -0.333 l +h + + + + +0 0 m +-1 0.333 l +-0.8 0 l +-1 -0.333 l +h + + + + +-1 0.333 m +0 0 l +-1 -0.333 l + + + + +0 0 m +-1 0.333 l +-1 -0.333 l +h +-1 0 m +-2 0.333 l +-2 -0.333 l +h + + + + +0 0 m +-1 0.333 l +-1 -0.333 l +h +-1 0 m +-2 0.333 l +-2 -0.333 l +h + + + + +0.5 0 m +-0.5 0.333 l +-0.5 -0.333 l +h + + + + +0.5 0 m +-0.5 0.333 l +-0.5 -0.333 l +h + + + + +0.5 0 m +-0.5 0.333 l +-0.3 0 l +-0.5 -0.333 l +h + + + + +0.5 0 m +-0.5 0.333 l +-0.3 0 l +-0.5 -0.333 l +h + + + + +1 0 m +0 0.333 l +0 -0.333 l +h +0 0 m +-1 0.333 l +-1 -0.333 l +h + + + + +1 0 m +0 0.333 l +0 -0.333 l +h +0 0 m +-1 0.333 l +-1 -0.333 l +h + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +(0,0) = (1/2, 1/2) = (1, 1) + +224 736 m +224 640 l +320 640 l +320 736 l +h + + +224 736 m +320 736 l + + +208 720 m +208 624 l + +(1, 1) +(0,0) ++s ++t + +320 672 m +384 688 l +448 752 l +512 768 l + + +320 688 m +384 704 l +448 768 l +512 784 l + ++s + +320 656 m +368 608 l + ++t +(1, 0) +(0, 0) +(0, 1) +(0.25, 1) +(0.75, 1) +(1, 1) + +320 672 m +384 608 l +448 624 l +512 688 l +576 704 l +512 768 l + + +384 688 m +448 624 l + + +448 752 m +512 688 l + + + From 4da862cf5d3bbaf3ebbc4c90f499f311cfd972dd Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Tue, 25 Jun 2024 16:26:27 +1200 Subject: [PATCH 02/14] Add pyproject metadata and qcore to requirements --- pyproject.toml | 7 +++++-- requirements.txt | 2 +- source_modelling/{sources => }/sources.py | 0 3 files changed, 6 insertions(+), 3 deletions(-) rename source_modelling/{sources => }/sources.py (100%) diff --git a/pyproject.toml b/pyproject.toml index ecc984f..27e3198 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -3,14 +3,17 @@ requires = ["setuptools", "setuptools-scm"] build-backend = "setuptools.build_meta" [project] -name = "PACKAGE NAME HERE" +name = "source_modelling" authors = [ {name = "QuakeCoRE" }, ] -description = "PACKAGE DESCRIPTION HERE" +description = "Source modelling library" readme = "README.md" requires-python = ">=3.9" dynamic = ["version", "dependencies"] +[tool.setuptools.package-dir] +source_modelling = "source_modelling" + [tool.setuptools.dynamic] dependencies = {file = ["requirements.txt"]} diff --git a/requirements.txt b/requirements.txt index 398ca04..aba7c42 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,3 +1,3 @@ -git+https://github.com/ucgmsim/qcore.git +qcore @ git+https://github.com/ucgmsim/qcore.git numpy scipy diff --git a/source_modelling/sources/sources.py b/source_modelling/sources.py similarity index 100% rename from source_modelling/sources/sources.py rename to source_modelling/sources.py From 03af787c1e2ffb2cb0c72096e5de4ed74e9945a6 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Tue, 25 Jun 2024 16:27:46 +1200 Subject: [PATCH 03/14] fix title --- wiki/Sources.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/wiki/Sources.md b/wiki/Sources.md index 8897d7b..c69cd6b 100644 --- a/wiki/Sources.md +++ b/wiki/Sources.md @@ -1,4 +1,4 @@ -# Sources in QCore +# Source Geometry A number of packages need to deal with fault geometry including: From ec759b8339cc4ea23fc00756388b1ca2eb1b56b9 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Thu, 27 Jun 2024 15:49:57 +1200 Subject: [PATCH 04/14] trailing whitespace --- wiki/Sources.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/wiki/Sources.md b/wiki/Sources.md index c69cd6b..61dfc89 100644 --- a/wiki/Sources.md +++ b/wiki/Sources.md @@ -1,4 +1,4 @@ -# Source Geometry +# Source Geometry A number of packages need to deal with fault geometry including: From 69cba25e601e22051feda1ea7d03b91f7f303a7c Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Thu, 27 Jun 2024 16:29:59 +1200 Subject: [PATCH 05/14] s/Type-I/Type-1 --- wiki/Sources.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/wiki/Sources.md b/wiki/Sources.md index 61dfc89..5f3d578 100644 --- a/wiki/Sources.md +++ b/wiki/Sources.md @@ -87,7 +87,7 @@ While we have five types of sources (per [Source Modelling for GMSim](https://wi 2. **Plane Geometry**: This a 2-dimensional source geometry consisting of a single plane. The extents of the geometry are its corners. The `qcore.sources` module uses the `Plane` class to model single plane geometry. 3. **Multi-Planar Geometry**: This is a 2-dimensional source geometry consisting of a number of connected planes. The extents of the geometry are the corners of the end fault planes. The `qcore.sources` module uses the `Fault` class to model multi-planar geometry. -Type-I fault are an instance of the the first geometry, type-2 faults +Type-1 fault are an instance of the the first geometry, type-2 faults are plane geometries, and type-4 and type-5 are multi-planar geometries. From 3a5687f3863166ed8f89258a517e86ef280c7d1d Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Thu, 27 Jun 2024 16:31:21 +1200 Subject: [PATCH 06/14] specify all sources have bounds --- source_modelling/sources.py | 62 +++++++++++++++++++++++-------------- 1 file changed, 38 insertions(+), 24 deletions(-) diff --git a/source_modelling/sources.py b/source_modelling/sources.py index 0a78b25..f4115ea 100644 --- a/source_modelling/sources.py +++ b/source_modelling/sources.py @@ -31,7 +31,8 @@ class Point: """A representation of point source.""" - point_coordinates: np.ndarray + # The bounds of a point source are just the coordinates of the point + bounds: np.ndarray # used to approximate point source as a small planar patch (metres). length_m: float # The usual strike, dip, dip direction, etc cannot be calculated @@ -40,6 +41,18 @@ class Point: dip: float dip_dir: float + @property + def coordinates(self) -> np.ndarray: + """Return the coordinates of the point in (lat, lon, depth) format. + + Returns + ------- + np.ndarray + The coordinates of the point in (lat, lon, depth) format. + Depth is in metres. + """ + return coordinates.wgs_depth_to_nztm(self.bounds) + @property def length(self) -> float: """ @@ -88,7 +101,7 @@ def fault_coordinates_to_wgs_depth_coordinates( coordinates are just the location of the point source. """ - return self.point_coordinates + return self.coordinates def wgs_depth_coordinates_to_fault_coordinates( self, wgs_depth_coordinates: np.ndarray @@ -113,10 +126,7 @@ def wgs_depth_coordinates_to_fault_coordinates( If the point is not near the source point. """ nztm_coordinates = coordinates.wgs_depth_to_nztm(wgs_depth_coordinates) - if np.all( - np.abs(nztm_coordinates - self.point_coordinates)[:2] / _KM_TO_M - < self.length - ): + if np.all(np.abs(nztm_coordinates - self.bounds)[:2] / _KM_TO_M < self.length): return np.array([1 / 2, 1 / 2]) # Point is in the centre of the small patch raise ValueError("Given global coordinates out of bounds for point source.") @@ -149,7 +159,8 @@ class Plane: 3 2 """ - corners_nztm: np.ndarray + # Bounds for plane are just the corners + bounds: np.ndarray @property def corners(self) -> np.ndarray: @@ -160,7 +171,7 @@ def corners(self) -> np.ndarray: The corners of the fault plane in (lat, lon, depth) format. The corners are the same as in corners_nztm. """ - return coordinates.nztm_to_wgs_depth(self.corners_nztm) + return coordinates.nztm_to_wgs_depth(self.bounds) @property def length_m(self) -> float: @@ -170,7 +181,7 @@ def length_m(self) -> float: float The length of the fault plane (in metres). """ - return np.linalg.norm(self.corners_nztm[1] - self.corners_nztm[0]) + return np.linalg.norm(self.bounds[1] - self.bounds[0]) @property def width_m(self) -> float: @@ -180,7 +191,7 @@ def width_m(self) -> float: float The width of the fault plane (in metres). """ - return np.linalg.norm(self.corners_nztm[-1] - self.corners_nztm[0]) + return np.linalg.norm(self.bounds[-1] - self.bounds[0]) @property def bottom_m(self) -> float: @@ -190,7 +201,7 @@ def bottom_m(self) -> float: float The bottom depth (in metres). """ - return self.corners_nztm[-1, -1] + return self.bounds[-1, -1] @property def width(self) -> float: @@ -244,7 +255,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] + strike_direction = self.bounds[1] - self.bounds[0] return geo.oriented_bearing_wrt_normal( north_direction, strike_direction, up_direction ) @@ -261,7 +272,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_nztm[-1] - self.corners_nztm[0] + dip_direction = self.bounds[-1] - self.bounds[0] dip_direction[-1] = 0 return geo.oriented_bearing_wrt_normal( north_direction, dip_direction, up_direction @@ -358,9 +369,9 @@ def fault_coordinates_to_wgs_depth_coordinates( np.ndarray An 3d-vector of (lat, lon, depth) transformed coordinates. """ - origin = self.corners_nztm[0] - top_right = self.corners_nztm[1] - bottom_left = self.corners_nztm[-1] + origin = self.bounds[0] + top_right = self.bounds[1] + bottom_left = self.bounds[-1] frame = np.vstack((top_right - origin, bottom_left - origin)) return coordinates.nztm_to_wgs_depth(origin + plane_coordinates @ frame) @@ -390,9 +401,9 @@ def wgs_depth_coordinates_to_fault_coordinates( ValueError If the given coordinates do not lie in the fault plane. """ - origin = self.corners_nztm[0] - top_right = self.corners_nztm[1] - bottom_left = self.corners_nztm[-1] + origin = self.bounds[0] + top_right = self.bounds[1] + bottom_left = self.bounds[-1] frame = np.vstack((top_right - origin, bottom_left - origin)) offset = coordinates.wgs_depth_to_nztm(global_coordinates) - origin plane_coordinates, residual, _, _ = np.linalg.lstsq(frame.T, offset, rcond=None) @@ -440,7 +451,7 @@ def centroid(self) -> np.ndarray: """ return coordinates.nztm_to_wgs_depth( - np.mean(self.corners_nztm, axis=0).reshape((1, -1)) + np.mean(self.bounds, axis=0).reshape((1, -1)) ).ravel() @@ -546,7 +557,8 @@ def corners(self) -> np.ndarray: """ return np.vstack([plane.corners for plane in self.planes]) - def corners_nztm(self) -> np.ndarray: + @property + def bounds(self) -> np.ndarray: """Get all corners of a fault. Returns @@ -554,7 +566,7 @@ def corners_nztm(self) -> np.ndarray: 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]) + return np.vstack([plane.bounds for plane in self.planes]) def wgs_depth_coordinates_to_fault_coordinates( self, global_coordinates: np.ndarray @@ -650,9 +662,11 @@ def fault_coordinates_to_wgs_depth_coordinates( ) -class HasCoordinates(Protocol): +class IsSource(Protocol): """Type definition for a source with local coordinates.""" + bounds: np.ndarray + def fault_coordinates_to_wgs_depth_coordinates( self, fault_coordinates: np.ndarray, @@ -665,7 +679,7 @@ def wgs_depth_coordinates_to_fault_coordinates( def closest_point_between_sources( - source_a: HasCoordinates, source_b: HasCoordinates + source_a: IsSource, source_b: IsSource ) -> tuple[np.ndarray, np.ndarray]: """Find the closest point between two sources that have local coordinates. From 66cce70ccb2ef4dcfc92f25bea79daaee8a9e6b4 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Mon, 1 Jul 2024 15:25:16 +1200 Subject: [PATCH 07/14] add useful staticmethod constructors --- source_modelling/sources.py | 53 ++++++++++++++++++++++++++++++++++++- 1 file changed, 52 insertions(+), 1 deletion(-) diff --git a/source_modelling/sources.py b/source_modelling/sources.py index f4115ea..0433f8b 100644 --- a/source_modelling/sources.py +++ b/source_modelling/sources.py @@ -41,6 +41,25 @@ class Point: dip: float dip_dir: float + @staticmethod + def from_lat_lon_depth(point_coordinates: np.ndarray, **kwargs) -> "Point": + """Construct a point source from a lat, lon, depth format. + + Parameters + ---------- + point_coordinates : np.ndarray + The coordinates of the point in lat, lon, depth format. + kwargs + The remaining point source arguments (see the class-level docstring) + + Returns + ------- + Point + The Point source representing this geometry. + + """ + return Point(bounds=coordinates.wgs_depth_to_nztm(point_coordinates), **kwargs) + @property def coordinates(self) -> np.ndarray: """Return the coordinates of the point in (lat, lon, depth) format. @@ -141,7 +160,7 @@ class Plane: Attributes ---------- - corners_nztm : np.ndarray + bounds : 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. @@ -162,6 +181,22 @@ class Plane: # Bounds for plane are just the corners bounds: np.ndarray + @staticmethod + def from_lat_lon_depth(corners: np.ndarray) -> "Plane": + """Construct a plane point source from its corners. + + Parameters + ---------- + corners : np.ndarray + The corners in lat, lon, depth format. Has shape (4 x 3). + + Returns + ------- + Plane + The plane source representing this geometry. + """ + return Plane(coordinates.wgs_depth_to_nztm(corners)) + @property def corners(self) -> np.ndarray: """ @@ -489,6 +524,22 @@ class Fault: planes: list[Plane] + @staticmethod + def from_lat_lon_depth(fault_corners: np.ndarray) -> "Fault": + """Construct a plane source geometry from the corners of the plane. + + Parameters + ---------- + corners : np.ndarray + The corners of the plane in lat, lon, depth format. Has shape (n x 4 x 3). + + Returns + ------- + Fault + The fault object representing this geometry. + """ + return Fault([Plane.from_lat_lon_depth(corners) for corners in fault_corners]) + def area(self) -> float: """Compute the area of a fault. From e6bea12d347278245b301b0231714c3e18c3c031 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Mon, 1 Jul 2024 15:36:36 +1200 Subject: [PATCH 08/14] used incorrect coordinate change function --- source_modelling/sources.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source_modelling/sources.py b/source_modelling/sources.py index 0433f8b..707497d 100644 --- a/source_modelling/sources.py +++ b/source_modelling/sources.py @@ -70,7 +70,7 @@ def coordinates(self) -> np.ndarray: The coordinates of the point in (lat, lon, depth) format. Depth is in metres. """ - return coordinates.wgs_depth_to_nztm(self.bounds) + return coordinates.nztm_to_wgs_depth(self.bounds) @property def length(self) -> float: From 8981550d469b3e85a98eb7a37da45fb1c1b40a8e Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Mon, 1 Jul 2024 16:40:55 +1200 Subject: [PATCH 09/14] fix names --- source_modelling/sources.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/source_modelling/sources.py b/source_modelling/sources.py index 707497d..54c362a 100644 --- a/source_modelling/sources.py +++ b/source_modelling/sources.py @@ -182,7 +182,7 @@ class Plane: bounds: np.ndarray @staticmethod - def from_lat_lon_depth(corners: np.ndarray) -> "Plane": + def from_corners(corners: np.ndarray) -> "Plane": """Construct a plane point source from its corners. Parameters @@ -525,7 +525,7 @@ class Fault: planes: list[Plane] @staticmethod - def from_lat_lon_depth(fault_corners: np.ndarray) -> "Fault": + def from_corners(fault_corners: np.ndarray) -> "Fault": """Construct a plane source geometry from the corners of the plane. Parameters @@ -538,7 +538,7 @@ def from_lat_lon_depth(fault_corners: np.ndarray) -> "Fault": Fault The fault object representing this geometry. """ - return Fault([Plane.from_lat_lon_depth(corners) for corners in fault_corners]) + return Fault([Plane.from_corners(corners) for corners in fault_corners]) def area(self) -> float: """Compute the area of a fault. From 4e1c6472467be95e5893fce7a039d3bb1f46f056 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Wed, 3 Jul 2024 10:43:34 +1200 Subject: [PATCH 10/14] fix dip bug, make Fault.corners property for consistency --- source_modelling/sources.py | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/source_modelling/sources.py b/source_modelling/sources.py index 54c362a..4de89ec 100644 --- a/source_modelling/sources.py +++ b/source_modelling/sources.py @@ -238,6 +238,16 @@ def bottom_m(self) -> float: """ return self.bounds[-1, -1] + @property + def top_m(self) -> float: + """ + Returns + ------- + float + The top depth of the fault. + """ + return self.bounds[0, -1] + @property def width(self) -> float: """ @@ -321,7 +331,7 @@ def dip(self) -> float: float The dip angle of the fault. """ - return np.degrees(np.arcsin(np.abs(self.bottom_m) / self.width_m)) + return np.degrees(np.arcsin(np.abs(self.bottom_m - self.top_m) / self.width_m)) @staticmethod def from_centroid_strike_dip( @@ -597,6 +607,7 @@ def dip_dir(self) -> float: """ return self.planes[0].dip_dir + @property def corners(self) -> np.ndarray: """Get all corners of a fault. From 6f81c927fc9d375837af4db3d01e6058f69266d0 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Wed, 3 Jul 2024 12:29:27 +1200 Subject: [PATCH 11/14] use a more reliable method to convert to fault-local coordinates --- source_modelling/sources.py | 46 +++++++++++++++++++++---------------- 1 file changed, 26 insertions(+), 20 deletions(-) diff --git a/source_modelling/sources.py b/source_modelling/sources.py index 4de89ec..c843bd3 100644 --- a/source_modelling/sources.py +++ b/source_modelling/sources.py @@ -446,15 +446,24 @@ def wgs_depth_coordinates_to_fault_coordinates( ValueError If the given coordinates do not lie in the fault plane. """ - origin = self.bounds[0] - top_right = self.bounds[1] - bottom_left = self.bounds[-1] - frame = np.vstack((top_right - origin, bottom_left - origin)) - offset = coordinates.wgs_depth_to_nztm(global_coordinates) - origin - plane_coordinates, residual, _, _ = np.linalg.lstsq(frame.T, offset, rcond=None) - if not np.isclose(residual[0], 0, atol=1e-02): - raise ValueError("Coordinates do not lie in fault plane.") - return np.clip(plane_coordinates, 0, 1) + strike_direction = self.bounds[1, :2] - self.bounds[0, :2] + dip_direction = self.bounds[-1, :2] - self.bounds[0, :2] + offset = ( + coordinates.wgs_depth_to_nztm(global_coordinates[:2]) - self.bounds[0, :2] + ) + projection_matrix = np.array( + [ + strike_direction / (strike_direction**2).sum(), + dip_direction / (dip_direction**2).sum(), + ] + ).T + fault_local_coordinates = offset @ projection_matrix + if not np.all( + (fault_local_coordinates > 0 | np.isclose(fault_local_coordinates, 0)) + & (fault_local_coordinates < 1 | np.isclose(fault_local_coordinates, 1)) + ): + raise ValueError("Specified coordinates do not lie in plane") + return np.clip(fault_local_coordinates, 0, 1) def wgs_depth_coordinates_in_plane(self, global_coordinates: np.ndarray) -> bool: """Test if some global coordinates lie in the bounds of a plane. @@ -475,12 +484,7 @@ def wgs_depth_coordinates_in_plane(self, global_coordinates: np.ndarray) -> bool plane_coordinates = self.wgs_depth_coordinates_to_fault_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), - ) - ) + return True except ValueError: return False @@ -672,16 +676,18 @@ def wgs_depth_coordinates_to_fault_coordinates( """ # the right edges as a cumulative proportion of the fault length (e.g. [0.1, ..., 0.8]) - right_edges = self.lengths.cumsum() / self.length - right_edges = np.append(right_edges, 1) + left_edges = self.lengths.cumsum() / self.length + left_edges = np.insert(left_edges, 0, 0) for i, plane in enumerate(self.planes): - if plane.wgs_depth_coordinates_in_plane(global_coordinates): + try: plane_coordinates = plane.wgs_depth_coordinates_to_fault_coordinates( global_coordinates ) - return np.array([right_edges[i], 0]) + plane_coordinates * np.array( - [right_edges[i + 1] - right_edges[i], 1] + return np.array([left_edges[i], 0]) + plane_coordinates * np.array( + [left_edges[i + 1] - left_edges[i], 1] ) + except ValueError: + continue raise ValueError("Given coordinates are not on fault.") def fault_coordinates_to_wgs_depth_coordinates( From adc2eeee349656bf30958b6ac70f1d273a9caf00 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Wed, 3 Jul 2024 12:40:00 +1200 Subject: [PATCH 12/14] improve docs --- source_modelling/sources.py | 35 +++++++++++++++++------------------ 1 file changed, 17 insertions(+), 18 deletions(-) diff --git a/source_modelling/sources.py b/source_modelling/sources.py index c843bd3..6a84abc 100644 --- a/source_modelling/sources.py +++ b/source_modelling/sources.py @@ -1,4 +1,4 @@ -"""Module for representing the geometry seismic sources: point sources, fault planes and faults. +"""Module for representing the geometry of seismic sources: point sources, fault planes and faults. This module provides classes and functions for representing fault planes and faults, along with methods for calculating various properties such as @@ -6,7 +6,7 @@ Classes ------- -PointSource: +Point: A representation of a point source. Plane: @@ -29,7 +29,21 @@ @dataclasses.dataclass class Point: - """A representation of point source.""" + """A representation of point source. + + Attributes + ---------- + bounds : np.ndarray + The coordinates (NZTM) of the point source. + length_m : float + Length used to approximate the point source as a small planar patch (metres). + strike : float + The strike angle of the point source in degrees. + dip : float + The dip angle of the point source in degrees. + dip_dir : float + The dip direction of the point source in degrees. + """ # The bounds of a point source are just the coordinates of the point bounds: np.ndarray @@ -519,21 +533,6 @@ class Fault: ---------- planes : list[Plane] A list containing all the Planes that constitute the fault. - - Methods - ------- - area: - Compute the area of a fault. - widths: - Get the widths of all fault planes. - lengths: - Get the lengths of all fault planes. - corners: - Get all corners of a fault. - global_coordinates_to_fault_coordinates: - Convert global coordinates to fault coordinates. - fault_coordinates_to_wgsdepth_coordinates: - Convert fault coordinates to global coordinates. """ planes: list[Plane] From e2f8dbef5fa72e75dcb180fdf0e17177c5020261 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Thu, 4 Jul 2024 15:28:11 +1200 Subject: [PATCH 13/14] add centroid property, improve numerical stability --- source_modelling/sources.py | 46 +++++++++++++++++++++++++++++++++---- 1 file changed, 42 insertions(+), 4 deletions(-) diff --git a/source_modelling/sources.py b/source_modelling/sources.py index 6a84abc..879685c 100644 --- a/source_modelling/sources.py +++ b/source_modelling/sources.py @@ -21,7 +21,6 @@ import numpy as np import scipy as sp - from qcore import coordinates, geo, grid _KM_TO_M = 1000 @@ -116,6 +115,18 @@ def width(self) -> float: """ return self.width_m / _KM_TO_M + @property + def centroid(self) -> np.ndarray: + """ + + Returns + ------- + np.ndarray + The centroid of the point source (which is just the + point's coordinates). + """ + return self.coordinates + def fault_coordinates_to_wgs_depth_coordinates( self, fault_coordinates: np.ndarray ) -> np.ndarray: @@ -398,6 +409,17 @@ def from_centroid_strike_dip( ) return Plane(corners) + @property + def centroid(self) -> np.ndarray: + """ + + Returns + ------- + np.ndarray + The centre of the fault plane. + """ + return self.fault_coordinates_to_wgs_depth_coordinates(np.array([1 / 2, 1 / 2])) + def fault_coordinates_to_wgs_depth_coordinates( self, plane_coordinates: np.ndarray ) -> np.ndarray: @@ -633,6 +655,16 @@ def bounds(self) -> np.ndarray: """ return np.vstack([plane.bounds for plane in self.planes]) + @property + def centroid(self) -> np.ndarray: + """ + Returns + ------- + np.ndarray + The centre of the fault. + """ + return self.fault_coordinates_to_wgs_depth_coordinates(np.array([1 / 2, 1 / 2])) + def wgs_depth_coordinates_to_fault_coordinates( self, global_coordinates: np.ndarray ) -> np.ndarray: @@ -674,7 +706,7 @@ def wgs_depth_coordinates_to_fault_coordinates( If the given point does not lie on the fault. """ - # the right edges as a cumulative proportion of the fault length (e.g. [0.1, ..., 0.8]) + # the left edges as a cumulative proportion of the fault length (e.g. [0.1, ..., 0.8]) left_edges = self.lengths.cumsum() / self.length left_edges = np.insert(left_edges, 0, 0) for i, plane in enumerate(self.planes): @@ -710,7 +742,12 @@ def fault_coordinates_to_wgs_depth_coordinates( # the right edges as a cumulative proportion of the fault length (e.g. [0.1, ..., 0.8]) right_edges = self.lengths.cumsum() / self.length - fault_segment_index = np.searchsorted(right_edges, fault_coordinates[0]) + for i in range(len(self.planes)): + if fault_coordinates[0] < right_edges[i] or np.isclose( + fault_coordinates[0], right_edges[i] + ): + break + fault_segment_index = i left_proportion = ( right_edges[fault_segment_index - 1] if fault_segment_index > 0 else 0 ) @@ -722,6 +759,7 @@ def fault_coordinates_to_wgs_depth_coordinates( segment_proportion = (fault_coordinates[0] - left_proportion) / ( right_proportion - left_proportion ) + return self.planes[ fault_segment_index ].fault_coordinates_to_wgs_depth_coordinates( @@ -785,11 +823,11 @@ def fault_coordinate_distance(fault_coordinates: np.ndarray) -> float: fault_coordinate_distance, np.array([1 / 2, 1 / 2, 1 / 2, 1 / 2]), bounds=[(0, 1)] * 4, + options={"ftol": 1e-2}, ) if not res.success: raise ValueError( f"Optimisation failed to converge for provided sources: {res.message}" ) - return res.x[:2], res.x[2:] From 462b4d3d3016712483de4311ee91ef5b10d6a8a8 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Tue, 9 Jul 2024 15:12:15 +1200 Subject: [PATCH 14/14] centroid + numerical stability --- source_modelling/sources.py | 27 ++++++++++++++++++++++++--- 1 file changed, 24 insertions(+), 3 deletions(-) diff --git a/source_modelling/sources.py b/source_modelling/sources.py index 6a84abc..5cc6f79 100644 --- a/source_modelling/sources.py +++ b/source_modelling/sources.py @@ -116,6 +116,10 @@ def width(self) -> float: """ return self.width_m / _KM_TO_M + @property + def centroid(self) -> np.ndarray: + return self.coordinates + def fault_coordinates_to_wgs_depth_coordinates( self, fault_coordinates: np.ndarray ) -> np.ndarray: @@ -396,7 +400,11 @@ def from_centroid_strike_dip( length, width, ) - return Plane(corners) + return Plane(coordinates.wgs_depth_to_nztm(corners)) + + @property + def centroid(self) -> np.ndarray: + return self.fault_coordinates_to_wgs_depth_coordinates(np.array([1 / 2, 1 / 2])) def fault_coordinates_to_wgs_depth_coordinates( self, plane_coordinates: np.ndarray @@ -633,6 +641,10 @@ def bounds(self) -> np.ndarray: """ return np.vstack([plane.bounds for plane in self.planes]) + @property + def centroid(self) -> np.ndarray: + return self.fault_coordinates_to_wgs_depth_coordinates(np.array([1 / 2, 1 / 2])) + def wgs_depth_coordinates_to_fault_coordinates( self, global_coordinates: np.ndarray ) -> np.ndarray: @@ -674,7 +686,7 @@ def wgs_depth_coordinates_to_fault_coordinates( If the given point does not lie on the fault. """ - # the right edges as a cumulative proportion of the fault length (e.g. [0.1, ..., 0.8]) + # the left edges as a cumulative proportion of the fault length (e.g. [0.1, ..., 0.8]) left_edges = self.lengths.cumsum() / self.length left_edges = np.insert(left_edges, 0, 0) for i, plane in enumerate(self.planes): @@ -710,7 +722,12 @@ def fault_coordinates_to_wgs_depth_coordinates( # the right edges as a cumulative proportion of the fault length (e.g. [0.1, ..., 0.8]) right_edges = self.lengths.cumsum() / self.length - fault_segment_index = np.searchsorted(right_edges, fault_coordinates[0]) + for i in range(len(self.planes)): + if fault_coordinates[0] < right_edges[i] or np.isclose( + fault_coordinates[0], right_edges[i] + ): + break + fault_segment_index = i left_proportion = ( right_edges[fault_segment_index - 1] if fault_segment_index > 0 else 0 ) @@ -722,6 +739,8 @@ def fault_coordinates_to_wgs_depth_coordinates( segment_proportion = (fault_coordinates[0] - left_proportion) / ( right_proportion - left_proportion ) + if fault_segment_index >= len(self.planes): + breakpoint() return self.planes[ fault_segment_index ].fault_coordinates_to_wgs_depth_coordinates( @@ -785,9 +804,11 @@ def fault_coordinate_distance(fault_coordinates: np.ndarray) -> float: fault_coordinate_distance, np.array([1 / 2, 1 / 2, 1 / 2, 1 / 2]), bounds=[(0, 1)] * 4, + options={"ftol": 1e-2}, ) if not res.success: + breakpoint() raise ValueError( f"Optimisation failed to converge for provided sources: {res.message}" )