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 e69de29..aba7c42 100644 --- a/requirements.txt +++ b/requirements.txt @@ -0,0 +1,3 @@ +qcore @ git+https://github.com/ucgmsim/qcore.git +numpy +scipy diff --git a/source_modelling/sources.py b/source_modelling/sources.py new file mode 100644 index 0000000..b98d1f9 --- /dev/null +++ b/source_modelling/sources.py @@ -0,0 +1,839 @@ +"""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 +dimensions, orientation, and coordinate transformations. + +Classes +------- +Point: + 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. + + 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 + # 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 + + @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. + + Returns + ------- + np.ndarray + The coordinates of the point in (lat, lon, depth) format. + Depth is in metres. + """ + return coordinates.nztm_to_wgs_depth(self.bounds) + + @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 + + @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: + """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.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.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.") + + +@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 + ---------- + 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. + + 0 1 + ┌──────────┐ + │ │ + │ │ + │ │ + │ │ + │ │ + │ │ + │ │ + └──────────┘ + 3 2 + """ + + # Bounds for plane are just the corners + bounds: np.ndarray + + @staticmethod + def from_corners(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: + """ + 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.bounds) + + @property + def length_m(self) -> float: + """ + Returns + ------- + float + The length of the fault plane (in metres). + """ + return np.linalg.norm(self.bounds[1] - self.bounds[0]) + + @property + def width_m(self) -> float: + """ + Returns + ------- + float + The width of the fault plane (in metres). + """ + return np.linalg.norm(self.bounds[-1] - self.bounds[0]) + + @property + def bottom_m(self) -> float: + """ + Returns + ------- + float + The bottom depth (in metres). + """ + 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: + """ + 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.bounds[1] - self.bounds[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.bounds[-1] - self.bounds[0] + dip_direction[-1] = 0 + return geo.oriented_bearing_wrt_normal( + north_direction, dip_direction, up_direction + ) + + @property + def dip(self) -> float: + """ + Returns + ------- + float + The dip angle of the fault. + """ + return np.degrees(np.arcsin(np.abs(self.bottom_m - self.top_m) / self.width_m)) + + @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(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])) + + @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: + """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.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) + + 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. + """ + 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. + + 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 True + 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.bounds, 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. + """ + + planes: list[Plane] + + @staticmethod + def from_corners(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_corners(corners) for corners in fault_corners]) + + 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 + + @property + 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]) + + @property + def bounds(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.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: + """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 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): + try: + plane_coordinates = plane.wgs_depth_coordinates_to_fault_coordinates( + global_coordinates + ) + 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( + 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 + 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 + ) + 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 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, + ) -> np.ndarray: ... + + def wgs_depth_coordinates_to_fault_coordinates( + self, + fault_coordinates: np.ndarray, + ) -> np.ndarray: ... + + +def closest_point_between_sources( + source_a: IsSource, source_b: IsSource +) -> 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, + options={"ftol": 1e-2}, + ) + + if not res.success: + breakpoint() + 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..5f3d578 --- /dev/null +++ b/wiki/Sources.md @@ -0,0 +1,213 @@ +# Source Geometry + +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-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. + +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 + + +