source_modelling.sources

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.

  1"""Module for representing the geometry of seismic sources: point sources, fault planes and faults.
  2
  3This module provides classes and functions for representing fault planes and
  4faults, along with methods for calculating various properties such as
  5dimensions, orientation, and coordinate transformations.
  6
  7Classes
  8-------
  9Point:
 10    A representation of a point source.
 11
 12Plane:
 13    A representation of a single plane of a Fault.
 14
 15Fault:
 16    A representation of a fault, consisting of one or more Planes.
 17"""
 18
 19import dataclasses
 20from typing import Optional, Protocol
 21
 22import numpy as np
 23import scipy as sp
 24import shapely
 25from qcore import coordinates, geo, grid
 26
 27_KM_TO_M = 1000
 28
 29
 30@dataclasses.dataclass
 31class Point:
 32    """A representation of point source.
 33
 34    Attributes
 35    ----------
 36    bounds : np.ndarray
 37        The coordinates (NZTM) of the point source.
 38    length_m : float
 39        Length used to approximate the point source as a small planar patch (metres).
 40    strike : float
 41        The strike angle of the point source in degrees.
 42    dip : float
 43        The dip angle of the point source in degrees.
 44    dip_dir : float
 45        The dip direction of the point source in degrees.
 46    """
 47
 48    # The bounds of a point source are just the coordinates of the point
 49    bounds: np.ndarray
 50    # used to approximate point source as a small planar patch (metres).
 51    length_m: float
 52    # The usual strike, dip, dip direction, etc cannot be calculated
 53    # from a point source and so must be provided by the user.
 54    strike: float
 55    dip: float
 56    dip_dir: float
 57
 58    @staticmethod
 59    def from_lat_lon_depth(point_coordinates: np.ndarray, **kwargs) -> "Point":
 60        """Construct a point source from a lat, lon, depth format.
 61
 62        Parameters
 63        ----------
 64        point_coordinates : np.ndarray
 65            The coordinates of the point in lat, lon, depth format.
 66        kwargs
 67            The remaining point source arguments (see the class-level docstring)
 68
 69        Returns
 70        -------
 71        Point
 72            The Point source representing this geometry.
 73
 74        """
 75        return Point(bounds=coordinates.wgs_depth_to_nztm(point_coordinates), **kwargs)
 76
 77    @property
 78    def coordinates(self) -> np.ndarray:
 79        """np.ndarray: The coordinates of the point in (lat, lon, depth) format. Depth is in metres."""
 80        return coordinates.nztm_to_wgs_depth(self.bounds)
 81
 82    @property
 83    def length(self) -> float:
 84        """float: The length of the approximating planar patch (in kilometres)."""
 85        return self.length_m / _KM_TO_M
 86
 87    @property
 88    def width_m(self) -> float:
 89        """float: The width of the approximating planar patch (in metres)."""
 90        return self.length_m
 91
 92    @property
 93    def width(self) -> float:
 94        """float: The width of the approximating planar patch (in kilometres)."""
 95        return self.width_m / _KM_TO_M
 96
 97    @property
 98    def centroid(self) -> np.ndarray:
 99        """np.ndarray: The centroid of the point source (which is just the point's coordinates)."""
100        return self.coordinates
101
102    @property
103    def geometry(self) -> shapely.Point:
104        """shapely.Point: A shapely geometry for the point (projected onto the surface)."""
105        return shapely.Point(self.bounds)
106
107    def fault_coordinates_to_wgs_depth_coordinates(
108        self, fault_coordinates: np.ndarray
109    ) -> np.ndarray:
110        """Convert fault-local coordinates to global (lat, lon, depth) coordinates.
111
112        Parameters
113        ----------
114        fault_coordinates : np.ndarray
115            The local fault coordinates
116
117        Returns
118        -------
119        np.ndarray
120            The global coordinates for these fault-local
121            coordinates. Because this is a point-source, the global
122            coordinates are just the location of the point source.
123        """
124        return self.coordinates
125
126    def wgs_depth_coordinates_to_fault_coordinates(
127        self, wgs_depth_coordinates: np.ndarray
128    ) -> np.ndarray:
129        """Convert global coordinates into fault-local coordinates.
130
131        Parameters
132        ----------
133        wgs_depth_coordinates : np.ndarray
134            The global coordinates to convert.
135
136        Returns
137        -------
138        np.ndarray
139            The fault-local coordinates. Because this is a
140            point-source, the local coordinates are simply (1/2, 1/2)
141            near the source point and undefined elsewhere.
142
143        Raises
144        ------
145        ValueError
146            If the point is not near the source point.
147        """
148        nztm_coordinates = coordinates.wgs_depth_to_nztm(wgs_depth_coordinates)
149        distance = np.abs(nztm_coordinates - self.bounds).max() / _KM_TO_M
150        if distance < self.length / 2 or np.isclose(distance, self.length / 2):
151            return np.array([1 / 2, 1 / 2])  # Point is in the centre of the small patch
152        raise ValueError("Given global coordinates out of bounds for point source.")
153
154
155@dataclasses.dataclass
156class Plane:
157    """A representation of a single plane of a Fault.
158
159    This class represents a single plane of a fault, providing various
160    properties and methods for calculating its dimensions, orientation, and
161    converting coordinates between different reference frames.
162
163    Attributes
164    ----------
165    bounds : np.ndarray
166        The corners of the fault plane, in NZTM format. The order of the
167        corners is given clockwise from the top left (according to strike
168        and dip). See the diagram below.
169
170         0            1
171          ┌──────────┐
172          │          │
173          │          │
174          │          │
175          │          │
176          │          │
177          │          │
178          │          │
179          └──────────┘
180         3            2
181    """
182
183    # Bounds for plane are just the corners
184    bounds: np.ndarray
185
186    @staticmethod
187    def from_corners(corners: np.ndarray) -> "Plane":
188        """Construct a plane point source from its corners.
189
190        Parameters
191        ----------
192        corners : np.ndarray
193            The corners in lat, lon, depth format. Has shape (4 x 3).
194
195        Returns
196        -------
197        Plane
198            The plane source representing this geometry.
199        """
200        return Plane(coordinates.wgs_depth_to_nztm(corners))
201
202    @property
203    def corners(self) -> np.ndarray:
204        """np.ndarray: The corners of the fault plane in (lat, lon, depth) format. The corners are the same as in corners_nztm."""
205        return coordinates.nztm_to_wgs_depth(self.bounds)
206
207    @property
208    def length_m(self) -> float:
209        """float: The length of the fault plane (in metres)."""
210        return np.linalg.norm(self.bounds[1] - self.bounds[0])
211
212    @property
213    def width_m(self) -> float:
214        """float: The width of the fault plane (in metres)."""
215        return np.linalg.norm(self.bounds[-1] - self.bounds[0])
216
217    @property
218    def bottom_m(self) -> float:
219        """float: The bottom depth (in metres)."""
220        return self.bounds[-1, -1]
221
222    @property
223    def top_m(self) -> float:
224        """float: The top depth of the fault."""
225        return self.bounds[0, -1]
226
227    @property
228    def width(self) -> float:
229        """float: The width of the fault plane (in kilometres)."""
230        return self.width_m / _KM_TO_M
231
232    @property
233    def length(self) -> float:
234        """float: The length of the fault plane (in kilometres)."""
235        return self.length_m / _KM_TO_M
236
237    @property
238    def area(self) -> float:
239        """float: The area of the plane (in km^2)."""
240        return self.length * self.width
241
242    @property
243    def projected_width_m(self) -> float:
244        """float: The projected width of the fault plane (in metres)."""
245        return self.width_m * np.cos(np.radians(self.dip))
246
247    @property
248    def projected_width(self) -> float:
249        """float: The projected width of the fault plane (in kilometres)."""
250        return self.projected_width_m / _KM_TO_M
251
252    @property
253    def strike(self) -> float:
254        """float: The bearing of the strike direction of the fault (from north; in degrees)."""
255        north_direction = np.array([1, 0, 0])
256        up_direction = np.array([0, 0, 1])
257        strike_direction = self.bounds[1] - self.bounds[0]
258        return geo.oriented_bearing_wrt_normal(
259            north_direction, strike_direction, up_direction
260        )
261
262    @property
263    def dip_dir(self) -> float:
264        """float: The bearing of the dip direction (from north; in degrees)."""
265        if np.isclose(self.dip, 90):
266            return 0  # TODO: Is this right for this case?
267        north_direction = np.array([1, 0, 0])
268        up_direction = np.array([0, 0, 1])
269        dip_direction = self.bounds[-1] - self.bounds[0]
270        dip_direction[-1] = 0
271        return geo.oriented_bearing_wrt_normal(
272            north_direction, dip_direction, up_direction
273        )
274
275    @property
276    def dip(self) -> float:
277        """float: The dip angle of the fault."""
278        return np.degrees(np.arcsin(np.abs(self.bottom_m - self.top_m) / self.width_m))
279
280    @property
281    def geometry(self) -> shapely.Polygon:
282        """shapely.Polygon: A shapely geometry for the plane (projected onto the surface)."""
283        return shapely.Polygon(self.bounds)
284
285    @staticmethod
286    def from_centroid_strike_dip(
287        centroid: np.ndarray,
288        strike: float,
289        dip_dir: Optional[float],
290        dtop: float,
291        dbottom: float,
292        length: float,
293        projected_width: float,
294    ) -> "Plane":
295        """Create a fault plane from the centroid, strike, dip_dir, top, bottom, length, and width.
296
297        This is used for older descriptions of sources. Internally
298        converts everything to corners so self.strike ~ strike (but
299        not exactly due to rounding errors).
300
301        Parameters
302        ----------
303        centroid : np.ndarray
304            The centre of the fault plane in lat, lon coordinates.
305        strike : float
306            The strike of the fault (in degrees).
307        dip_dir : Optional[float]
308            The dip direction of the fault (in degrees). If None this is assumed to be strike + 90 degrees.
309        top : float
310            The top depth of the plane (in km).
311        bottom : float
312            The bottom depth of the plane (in km).
313        length : float
314            The length of the fault plane (in km).
315        projected_width : float
316            The projected width of the fault plane (in km).
317
318        Returns
319        -------
320        Plane
321            The fault plane with centre at `centroid`, and where the
322            parameters strike, dip_dir, top, bottom, length and width
323            match what is passed to this function.
324        """
325        corners = grid.grid_corners(
326            centroid,
327            strike,
328            dip_dir if dip_dir is not None else (strike + 90),
329            dtop,
330            dbottom,
331            length,
332            projected_width,
333        )
334        corners[[2, 3]] = corners[[3, 2]]
335        return Plane(coordinates.wgs_depth_to_nztm(np.array(corners)))
336
337    @property
338    def centroid(self) -> np.ndarray:
339        """np.ndarray: The centre of the fault plane."""
340        return self.fault_coordinates_to_wgs_depth_coordinates(np.array([1 / 2, 1 / 2]))
341
342    def fault_coordinates_to_wgs_depth_coordinates(
343        self, plane_coordinates: np.ndarray
344    ) -> np.ndarray:
345        """Convert plane coordinates to nztm global coordinates.
346
347        Parameters
348        ----------
349        plane_coordinates : np.ndarray
350            Plane coordinates to convert. Plane coordinates are
351            2D coordinates (x, y) given for a fault plane (a plane), where x
352            represents displacement along the strike, and y
353            displacement along the dip (see diagram below). The
354            origin for plane coordinates is the centre of the fault.
355
356                          +x
357             0 0   ─────────────────>
358                ┌─────────────────────┐ │
359                │      < strike >     │ │
360                │                 ^   │ │
361                │                dip  │ │ +y
362                │                 v   │ │
363                │                     │ │
364                └─────────────────────┘ ∨
365                                       1,1
366
367        Returns
368        -------
369        np.ndarray
370            An 3d-vector of (lat, lon, depth) transformed coordinates.
371        """
372        origin = self.bounds[0]
373        top_right = self.bounds[1]
374        bottom_left = self.bounds[-1]
375        frame = np.vstack((top_right - origin, bottom_left - origin))
376
377        return coordinates.nztm_to_wgs_depth(origin + plane_coordinates @ frame)
378
379    def wgs_depth_coordinates_to_fault_coordinates(
380        self,
381        global_coordinates: np.ndarray,
382    ) -> np.ndarray:
383        """Convert coordinates (lat, lon, depth) to plane coordinates (x, y).
384
385        See plane_coordinates_to_global_coordinates for a description of
386        plane coordinates.
387
388        Parameters
389        ----------
390        global_coordinates : np.ndarray
391            Global coordinates to convert.
392
393        Returns
394        -------
395        np.ndarray
396            The plane coordinates (x, y) representing the position of
397            global_coordinates on the fault plane.
398
399        Raises
400        ------
401        ValueError
402            If the given coordinates do not lie in the fault plane.
403        """
404        strike_direction = self.bounds[1, :2] - self.bounds[0, :2]
405        dip_direction = self.bounds[-1, :2] - self.bounds[0, :2]
406        offset = (
407            coordinates.wgs_depth_to_nztm(global_coordinates[:2]) - self.bounds[0, :2]
408        )
409        fault_local_coordinates, _, _, _ = np.linalg.lstsq(
410            np.array([strike_direction, dip_direction]).T, offset, rcond=None
411        )
412        if not np.all(
413            (
414                (fault_local_coordinates > 0)
415                | np.isclose(fault_local_coordinates, 0, atol=1e-6)
416            )
417            & (
418                (fault_local_coordinates < 1)
419                | np.isclose(fault_local_coordinates, 1, atol=1e-6)
420            )
421        ):
422            raise ValueError("Specified coordinates do not lie in plane")
423        return np.clip(fault_local_coordinates, 0, 1)
424
425
426@dataclasses.dataclass
427class Fault:
428    """A representation of a fault, consisting of one or more Planes.
429
430    This class represents a fault, which is composed of one or more Planes.
431    It provides methods for computing the area of the fault, getting the widths and
432    lengths of all fault planes, retrieving all corners of the fault, converting
433    global coordinates to fault coordinates, converting fault coordinates to global
434    coordinates, generating a random hypocentre location within the fault, and
435    computing the expected fault coordinates.
436
437    Attributes
438    ----------
439    planes : list[Plane]
440        A list containing all the Planes that constitute the fault.
441    """
442
443    planes: list[Plane]
444
445    @staticmethod
446    def from_corners(fault_corners: np.ndarray) -> "Fault":
447        """Construct a plane source geometry from the corners of the plane.
448
449        Parameters
450        ----------
451        corners : np.ndarray
452            The corners of the plane in lat, lon, depth format. Has shape (n x 4 x 3).
453
454        Returns
455        -------
456        Fault
457            The fault object representing this geometry.
458        """
459        return Fault([Plane.from_corners(corners) for corners in fault_corners])
460
461    def area(self) -> float:
462        """Compute the area of a fault.
463
464        Returns
465        -------
466        float
467            The area of the fault.
468        """
469        return self.width * np.sum(self.lengths)
470
471    @property
472    def lengths(self) -> np.ndarray:
473        """np.ndarray: A numpy array of each plane length (in km)."""
474        return np.array([fault.length for fault in self.planes])
475
476    @property
477    def length(self) -> float:
478        """float: The total length of each fault plane."""
479        return self.lengths.sum()
480
481    @property
482    def width(self) -> float:
483        """The width of the fault.
484
485        Returns
486        -------
487        float
488            The width of the first fault plane (A fault is assumed to
489            have planes of constant width).
490        """
491        return self.planes[0].width
492
493    @property
494    def dip_dir(self) -> float:
495        """float: The dip direction of the first fault plane (A fault is assumed to have planes of constant dip direction)."""
496        return self.planes[0].dip_dir
497
498    @property
499    def corners(self) -> np.ndarray:
500        """np.ndarray of shape (4n x 3): The corners in (lat, lon, depth) format of each fault plane in the fault, stacked vertically."""
501        return np.vstack([plane.corners for plane in self.planes])
502
503    @property
504    def bounds(self) -> np.ndarray:
505        """np.ndarray of shape (4n x 3): The corners in NZTM format of each fault plane in the fault, stacked vertically."""
506        return np.vstack([plane.bounds for plane in self.planes])
507
508    @property
509    def centroid(self) -> np.ndarray:
510        """np.ndarray: The centre of the fault."""
511        return self.fault_coordinates_to_wgs_depth_coordinates(np.array([1 / 2, 1 / 2]))
512
513    @property
514    def geometry(self) -> shapely.Polygon:
515        """shapely.Geometry: A shapely geometry for the fault (projected onto the surface)."""
516        return shapely.normalize(
517            shapely.union_all([plane.geometry for plane in self.planes])
518        )
519
520    def wgs_depth_coordinates_to_fault_coordinates(
521        self, global_coordinates: np.ndarray
522    ) -> np.ndarray:
523        """Convert global coordinates in (lat, lon, depth) format to fault coordinates.
524
525        Fault coordinates are a tuple (s, d) where s is the distance
526        from the top left, and d the distance from the top of the
527        fault (refer to the diagram). The coordinates are normalised
528        such that (0, 0) is the top left and (1, 1) the bottom right.
529
530        (0, 0)
531          ┌──────────────────────┬──────┐
532          │          |           │      │
533          │          |           │      │
534          │          | d         │      │
535          │          |           │      │
536          ├----------*           │      │
537          │    s     ^           │      │
538          │          |           │      │
539          │          |           │      │
540          │          |           │      │
541          └──────────|───────────┴──────┘
542                     +                    (1, 1)
543                  point: (s, d)
544
545        Parameters
546        ----------
547        global_coordinates : np.ndarray of shape (1 x 3)
548            The global coordinates to convert.
549
550        Returns
551        -------
552        np.ndarray
553            The fault coordinates.
554
555        Raises
556        ------
557        ValueError
558            If the given point does not lie on the fault.
559
560        """
561        # the left edges as a cumulative proportion of the fault length (e.g. [0.1, ..., 0.8])
562        left_edges = self.lengths.cumsum() / self.length
563        left_edges = np.insert(left_edges, 0, 0)
564        for i, plane in enumerate(self.planes):
565            try:
566                plane_coordinates = plane.wgs_depth_coordinates_to_fault_coordinates(
567                    global_coordinates
568                )
569                return np.array([left_edges[i], 0]) + plane_coordinates * np.array(
570                    [left_edges[i + 1] - left_edges[i], 1]
571                )
572            except ValueError:
573                continue
574        raise ValueError("Given coordinates are not on fault.")
575
576    def fault_coordinates_to_wgs_depth_coordinates(
577        self, fault_coordinates: np.ndarray
578    ) -> np.ndarray:
579        """Convert fault coordinates to global coordinates.
580
581        See global_coordinates_to_fault_coordinates for a description of fault
582        coordinates.
583
584        Parameters
585        ----------
586        fault_coordinates : np.ndarray
587            The fault coordinates of the point.
588
589        Returns
590        -------
591        np.ndarray
592            The global coordinates (lat, lon, depth) for this point.
593        """
594        # the right edges as a cumulative proportion of the fault length (e.g. [0.1, ..., 0.8])
595        right_edges = self.lengths.cumsum() / self.length
596        for i in range(len(self.planes)):
597            if fault_coordinates[0] < right_edges[i] or np.isclose(
598                fault_coordinates[0], right_edges[i]
599            ):
600                break
601        fault_segment_index = i
602        left_proportion = (
603            right_edges[fault_segment_index - 1] if fault_segment_index > 0 else 0
604        )
605        right_proportion = right_edges[fault_segment_index]
606        segment_proportion = (fault_coordinates[0] - left_proportion) / (
607            right_proportion - left_proportion
608        )
609
610        return self.planes[
611            fault_segment_index
612        ].fault_coordinates_to_wgs_depth_coordinates(
613            np.array([segment_proportion, fault_coordinates[1]])
614        )
615
616
617class IsSource(Protocol):
618    """Type definition for a source with local coordinates."""
619
620    bounds: np.ndarray
621
622    def fault_coordinates_to_wgs_depth_coordinates(
623        self,
624        fault_coordinates: np.ndarray,
625    ) -> np.ndarray: ...
626
627    def wgs_depth_coordinates_to_fault_coordinates(
628        self,
629        fault_coordinates: np.ndarray,
630    ) -> np.ndarray: ...
631
632
633def closest_point_between_sources(
634    source_a: IsSource, source_b: IsSource
635) -> tuple[np.ndarray, np.ndarray]:
636    """Find the closest point between two sources that have local coordinates.
637
638    Parameters
639    ----------
640    source_a : HasCoordinates
641        The first source. Must have a two-dimensional fault coordinate system.
642    source_b : HasCoordinates
643        The first source. Must have a two-dimensional fault coordinate system.
644
645    Raises
646    ------
647    ValueError
648        Raised when we are unable to converge on the closest points between sources.
649
650    Returns
651    -------
652    source_a_coordinates : np.ndarray
653        The source-local coordinates of the closest point on source a.
654    source_b_coordinates : np.ndarray
655        The source-local coordinates of the closest point on source b.
656    """
657
658    def fault_coordinate_distance(fault_coordinates: np.ndarray) -> float:
659        source_a_global_coordinates = (
660            source_a.fault_coordinates_to_wgs_depth_coordinates(fault_coordinates[:2])
661        )
662        source_b_global_coordinates = (
663            source_b.fault_coordinates_to_wgs_depth_coordinates(fault_coordinates[2:])
664        )
665        return coordinates.wgs_depth_to_nztm(
666            source_a_global_coordinates
667        ) - coordinates.wgs_depth_to_nztm(source_b_global_coordinates)
668
669    res = sp.optimize.least_squares(
670        fault_coordinate_distance,
671        np.array([1 / 2, 1 / 2, 1 / 2, 1 / 2]),
672        bounds=([0] * 4, [1] * 4),
673        gtol=1e-5,
674        ftol=1e-5,
675        xtol=1e-5,
676    )
677
678    if not res.success and res.status != 0:
679        raise ValueError(
680            f"Optimisation failed to converge for provided sources: {res.message} with x = {res.x}"
681        )
682    return res.x[:2], res.x[2:]
@dataclasses.dataclass
class Point: 
 31@dataclasses.dataclass
 32class Point:
 33    """A representation of point source.
 34
 35    Attributes
 36    ----------
 37    bounds : np.ndarray
 38        The coordinates (NZTM) of the point source.
 39    length_m : float
 40        Length used to approximate the point source as a small planar patch (metres).
 41    strike : float
 42        The strike angle of the point source in degrees.
 43    dip : float
 44        The dip angle of the point source in degrees.
 45    dip_dir : float
 46        The dip direction of the point source in degrees.
 47    """
 48
 49    # The bounds of a point source are just the coordinates of the point
 50    bounds: np.ndarray
 51    # used to approximate point source as a small planar patch (metres).
 52    length_m: float
 53    # The usual strike, dip, dip direction, etc cannot be calculated
 54    # from a point source and so must be provided by the user.
 55    strike: float
 56    dip: float
 57    dip_dir: float
 58
 59    @staticmethod
 60    def from_lat_lon_depth(point_coordinates: np.ndarray, **kwargs) -> "Point":
 61        """Construct a point source from a lat, lon, depth format.
 62
 63        Parameters
 64        ----------
 65        point_coordinates : np.ndarray
 66            The coordinates of the point in lat, lon, depth format.
 67        kwargs
 68            The remaining point source arguments (see the class-level docstring)
 69
 70        Returns
 71        -------
 72        Point
 73            The Point source representing this geometry.
 74
 75        """
 76        return Point(bounds=coordinates.wgs_depth_to_nztm(point_coordinates), **kwargs)
 77
 78    @property
 79    def coordinates(self) -> np.ndarray:
 80        """np.ndarray: The coordinates of the point in (lat, lon, depth) format. Depth is in metres."""
 81        return coordinates.nztm_to_wgs_depth(self.bounds)
 82
 83    @property
 84    def length(self) -> float:
 85        """float: The length of the approximating planar patch (in kilometres)."""
 86        return self.length_m / _KM_TO_M
 87
 88    @property
 89    def width_m(self) -> float:
 90        """float: The width of the approximating planar patch (in metres)."""
 91        return self.length_m
 92
 93    @property
 94    def width(self) -> float:
 95        """float: The width of the approximating planar patch (in kilometres)."""
 96        return self.width_m / _KM_TO_M
 97
 98    @property
 99    def centroid(self) -> np.ndarray:
100        """np.ndarray: The centroid of the point source (which is just the point's coordinates)."""
101        return self.coordinates
102
103    @property
104    def geometry(self) -> shapely.Point:
105        """shapely.Point: A shapely geometry for the point (projected onto the surface)."""
106        return shapely.Point(self.bounds)
107
108    def fault_coordinates_to_wgs_depth_coordinates(
109        self, fault_coordinates: np.ndarray
110    ) -> np.ndarray:
111        """Convert fault-local coordinates to global (lat, lon, depth) coordinates.
112
113        Parameters
114        ----------
115        fault_coordinates : np.ndarray
116            The local fault coordinates
117
118        Returns
119        -------
120        np.ndarray
121            The global coordinates for these fault-local
122            coordinates. Because this is a point-source, the global
123            coordinates are just the location of the point source.
124        """
125        return self.coordinates
126
127    def wgs_depth_coordinates_to_fault_coordinates(
128        self, wgs_depth_coordinates: np.ndarray
129    ) -> np.ndarray:
130        """Convert global coordinates into fault-local coordinates.
131
132        Parameters
133        ----------
134        wgs_depth_coordinates : np.ndarray
135            The global coordinates to convert.
136
137        Returns
138        -------
139        np.ndarray
140            The fault-local coordinates. Because this is a
141            point-source, the local coordinates are simply (1/2, 1/2)
142            near the source point and undefined elsewhere.
143
144        Raises
145        ------
146        ValueError
147            If the point is not near the source point.
148        """
149        nztm_coordinates = coordinates.wgs_depth_to_nztm(wgs_depth_coordinates)
150        distance = np.abs(nztm_coordinates - self.bounds).max() / _KM_TO_M
151        if distance < self.length / 2 or np.isclose(distance, self.length / 2):
152            return np.array([1 / 2, 1 / 2])  # Point is in the centre of the small patch
153        raise ValueError("Given global coordinates out of bounds for 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.
Point( bounds: numpy.ndarray, length_m: float, strike: float, dip: float, dip_dir: float)
bounds: numpy.ndarray
length_m: float
strike: float
dip: float
dip_dir: float
@staticmethod
def from_lat_lon_depth( point_coordinates: numpy.ndarray, **kwargs) -> Point: 
59    @staticmethod
60    def from_lat_lon_depth(point_coordinates: np.ndarray, **kwargs) -> "Point":
61        """Construct a point source from a lat, lon, depth format.
62
63        Parameters
64        ----------
65        point_coordinates : np.ndarray
66            The coordinates of the point in lat, lon, depth format.
67        kwargs
68            The remaining point source arguments (see the class-level docstring)
69
70        Returns
71        -------
72        Point
73            The Point source representing this geometry.
74
75        """
76        return Point(bounds=coordinates.wgs_depth_to_nztm(point_coordinates), **kwargs)

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.
coordinates: numpy.ndarray 
78    @property
79    def coordinates(self) -> np.ndarray:
80        """np.ndarray: The coordinates of the point in (lat, lon, depth) format. Depth is in metres."""
81        return coordinates.nztm_to_wgs_depth(self.bounds)

np.ndarray: The coordinates of the point in (lat, lon, depth) format. Depth is in metres.

length: float 
83    @property
84    def length(self) -> float:
85        """float: The length of the approximating planar patch (in kilometres)."""
86        return self.length_m / _KM_TO_M

float: The length of the approximating planar patch (in kilometres).

width_m: float 
88    @property
89    def width_m(self) -> float:
90        """float: The width of the approximating planar patch (in metres)."""
91        return self.length_m

float: The width of the approximating planar patch (in metres).

width: float 
93    @property
94    def width(self) -> float:
95        """float: The width of the approximating planar patch (in kilometres)."""
96        return self.width_m / _KM_TO_M

float: The width of the approximating planar patch (in kilometres).

centroid: numpy.ndarray 
 98    @property
 99    def centroid(self) -> np.ndarray:
100        """np.ndarray: The centroid of the point source (which is just the point's coordinates)."""
101        return self.coordinates

np.ndarray: The centroid of the point source (which is just the point's coordinates).

geometry: shapely.geometry.point.Point 
103    @property
104    def geometry(self) -> shapely.Point:
105        """shapely.Point: A shapely geometry for the point (projected onto the surface)."""
106        return shapely.Point(self.bounds)

shapely.Point: A shapely geometry for the point (projected onto the surface).

def fault_coordinates_to_wgs_depth_coordinates(self, fault_coordinates: numpy.ndarray) -> numpy.ndarray: 
108    def fault_coordinates_to_wgs_depth_coordinates(
109        self, fault_coordinates: np.ndarray
110    ) -> np.ndarray:
111        """Convert fault-local coordinates to global (lat, lon, depth) coordinates.
112
113        Parameters
114        ----------
115        fault_coordinates : np.ndarray
116            The local fault coordinates
117
118        Returns
119        -------
120        np.ndarray
121            The global coordinates for these fault-local
122            coordinates. Because this is a point-source, the global
123            coordinates are just the location of the point source.
124        """
125        return self.coordinates

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.
def wgs_depth_coordinates_to_fault_coordinates(self, wgs_depth_coordinates: numpy.ndarray) -> numpy.ndarray: 
127    def wgs_depth_coordinates_to_fault_coordinates(
128        self, wgs_depth_coordinates: np.ndarray
129    ) -> np.ndarray:
130        """Convert global coordinates into fault-local coordinates.
131
132        Parameters
133        ----------
134        wgs_depth_coordinates : np.ndarray
135            The global coordinates to convert.
136
137        Returns
138        -------
139        np.ndarray
140            The fault-local coordinates. Because this is a
141            point-source, the local coordinates are simply (1/2, 1/2)
142            near the source point and undefined elsewhere.
143
144        Raises
145        ------
146        ValueError
147            If the point is not near the source point.
148        """
149        nztm_coordinates = coordinates.wgs_depth_to_nztm(wgs_depth_coordinates)
150        distance = np.abs(nztm_coordinates - self.bounds).max() / _KM_TO_M
151        if distance < self.length / 2 or np.isclose(distance, self.length / 2):
152            return np.array([1 / 2, 1 / 2])  # Point is in the centre of the small patch
153        raise ValueError("Given global coordinates out of bounds for point source.")

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.
@dataclasses.dataclass
class Plane: 
156@dataclasses.dataclass
157class Plane:
158    """A representation of a single plane of a Fault.
159
160    This class represents a single plane of a fault, providing various
161    properties and methods for calculating its dimensions, orientation, and
162    converting coordinates between different reference frames.
163
164    Attributes
165    ----------
166    bounds : np.ndarray
167        The corners of the fault plane, in NZTM format. The order of the
168        corners is given clockwise from the top left (according to strike
169        and dip). See the diagram below.
170
171         0            1
172          ┌──────────┐
173          │          │
174          │          │
175          │          │
176          │          │
177          │          │
178          │          │
179          │          │
180          └──────────┘
181         3            2
182    """
183
184    # Bounds for plane are just the corners
185    bounds: np.ndarray
186
187    @staticmethod
188    def from_corners(corners: np.ndarray) -> "Plane":
189        """Construct a plane point source from its corners.
190
191        Parameters
192        ----------
193        corners : np.ndarray
194            The corners in lat, lon, depth format. Has shape (4 x 3).
195
196        Returns
197        -------
198        Plane
199            The plane source representing this geometry.
200        """
201        return Plane(coordinates.wgs_depth_to_nztm(corners))
202
203    @property
204    def corners(self) -> np.ndarray:
205        """np.ndarray: The corners of the fault plane in (lat, lon, depth) format. The corners are the same as in corners_nztm."""
206        return coordinates.nztm_to_wgs_depth(self.bounds)
207
208    @property
209    def length_m(self) -> float:
210        """float: The length of the fault plane (in metres)."""
211        return np.linalg.norm(self.bounds[1] - self.bounds[0])
212
213    @property
214    def width_m(self) -> float:
215        """float: The width of the fault plane (in metres)."""
216        return np.linalg.norm(self.bounds[-1] - self.bounds[0])
217
218    @property
219    def bottom_m(self) -> float:
220        """float: The bottom depth (in metres)."""
221        return self.bounds[-1, -1]
222
223    @property
224    def top_m(self) -> float:
225        """float: The top depth of the fault."""
226        return self.bounds[0, -1]
227
228    @property
229    def width(self) -> float:
230        """float: The width of the fault plane (in kilometres)."""
231        return self.width_m / _KM_TO_M
232
233    @property
234    def length(self) -> float:
235        """float: The length of the fault plane (in kilometres)."""
236        return self.length_m / _KM_TO_M
237
238    @property
239    def area(self) -> float:
240        """float: The area of the plane (in km^2)."""
241        return self.length * self.width
242
243    @property
244    def projected_width_m(self) -> float:
245        """float: The projected width of the fault plane (in metres)."""
246        return self.width_m * np.cos(np.radians(self.dip))
247
248    @property
249    def projected_width(self) -> float:
250        """float: The projected width of the fault plane (in kilometres)."""
251        return self.projected_width_m / _KM_TO_M
252
253    @property
254    def strike(self) -> float:
255        """float: The bearing of the strike direction of the fault (from north; in degrees)."""
256        north_direction = np.array([1, 0, 0])
257        up_direction = np.array([0, 0, 1])
258        strike_direction = self.bounds[1] - self.bounds[0]
259        return geo.oriented_bearing_wrt_normal(
260            north_direction, strike_direction, up_direction
261        )
262
263    @property
264    def dip_dir(self) -> float:
265        """float: The bearing of the dip direction (from north; in degrees)."""
266        if np.isclose(self.dip, 90):
267            return 0  # TODO: Is this right for this case?
268        north_direction = np.array([1, 0, 0])
269        up_direction = np.array([0, 0, 1])
270        dip_direction = self.bounds[-1] - self.bounds[0]
271        dip_direction[-1] = 0
272        return geo.oriented_bearing_wrt_normal(
273            north_direction, dip_direction, up_direction
274        )
275
276    @property
277    def dip(self) -> float:
278        """float: The dip angle of the fault."""
279        return np.degrees(np.arcsin(np.abs(self.bottom_m - self.top_m) / self.width_m))
280
281    @property
282    def geometry(self) -> shapely.Polygon:
283        """shapely.Polygon: A shapely geometry for the plane (projected onto the surface)."""
284        return shapely.Polygon(self.bounds)
285
286    @staticmethod
287    def from_centroid_strike_dip(
288        centroid: np.ndarray,
289        strike: float,
290        dip_dir: Optional[float],
291        dtop: float,
292        dbottom: float,
293        length: float,
294        projected_width: float,
295    ) -> "Plane":
296        """Create a fault plane from the centroid, strike, dip_dir, top, bottom, length, and width.
297
298        This is used for older descriptions of sources. Internally
299        converts everything to corners so self.strike ~ strike (but
300        not exactly due to rounding errors).
301
302        Parameters
303        ----------
304        centroid : np.ndarray
305            The centre of the fault plane in lat, lon coordinates.
306        strike : float
307            The strike of the fault (in degrees).
308        dip_dir : Optional[float]
309            The dip direction of the fault (in degrees). If None this is assumed to be strike + 90 degrees.
310        top : float
311            The top depth of the plane (in km).
312        bottom : float
313            The bottom depth of the plane (in km).
314        length : float
315            The length of the fault plane (in km).
316        projected_width : float
317            The projected width of the fault plane (in km).
318
319        Returns
320        -------
321        Plane
322            The fault plane with centre at `centroid`, and where the
323            parameters strike, dip_dir, top, bottom, length and width
324            match what is passed to this function.
325        """
326        corners = grid.grid_corners(
327            centroid,
328            strike,
329            dip_dir if dip_dir is not None else (strike + 90),
330            dtop,
331            dbottom,
332            length,
333            projected_width,
334        )
335        corners[[2, 3]] = corners[[3, 2]]
336        return Plane(coordinates.wgs_depth_to_nztm(np.array(corners)))
337
338    @property
339    def centroid(self) -> np.ndarray:
340        """np.ndarray: The centre of the fault plane."""
341        return self.fault_coordinates_to_wgs_depth_coordinates(np.array([1 / 2, 1 / 2]))
342
343    def fault_coordinates_to_wgs_depth_coordinates(
344        self, plane_coordinates: np.ndarray
345    ) -> np.ndarray:
346        """Convert plane coordinates to nztm global coordinates.
347
348        Parameters
349        ----------
350        plane_coordinates : np.ndarray
351            Plane coordinates to convert. Plane coordinates are
352            2D coordinates (x, y) given for a fault plane (a plane), where x
353            represents displacement along the strike, and y
354            displacement along the dip (see diagram below). The
355            origin for plane coordinates is the centre of the fault.
356
357                          +x
358             0 0   ─────────────────>
359                ┌─────────────────────┐ │
360                │      < strike >     │ │
361                │                 ^   │ │
362                │                dip  │ │ +y
363                │                 v   │ │
364                │                     │ │
365                └─────────────────────┘ ∨
366                                       1,1
367
368        Returns
369        -------
370        np.ndarray
371            An 3d-vector of (lat, lon, depth) transformed coordinates.
372        """
373        origin = self.bounds[0]
374        top_right = self.bounds[1]
375        bottom_left = self.bounds[-1]
376        frame = np.vstack((top_right - origin, bottom_left - origin))
377
378        return coordinates.nztm_to_wgs_depth(origin + plane_coordinates @ frame)
379
380    def wgs_depth_coordinates_to_fault_coordinates(
381        self,
382        global_coordinates: np.ndarray,
383    ) -> np.ndarray:
384        """Convert coordinates (lat, lon, depth) to plane coordinates (x, y).
385
386        See plane_coordinates_to_global_coordinates for a description of
387        plane coordinates.
388
389        Parameters
390        ----------
391        global_coordinates : np.ndarray
392            Global coordinates to convert.
393
394        Returns
395        -------
396        np.ndarray
397            The plane coordinates (x, y) representing the position of
398            global_coordinates on the fault plane.
399
400        Raises
401        ------
402        ValueError
403            If the given coordinates do not lie in the fault plane.
404        """
405        strike_direction = self.bounds[1, :2] - self.bounds[0, :2]
406        dip_direction = self.bounds[-1, :2] - self.bounds[0, :2]
407        offset = (
408            coordinates.wgs_depth_to_nztm(global_coordinates[:2]) - self.bounds[0, :2]
409        )
410        fault_local_coordinates, _, _, _ = np.linalg.lstsq(
411            np.array([strike_direction, dip_direction]).T, offset, rcond=None
412        )
413        if not np.all(
414            (
415                (fault_local_coordinates > 0)
416                | np.isclose(fault_local_coordinates, 0, atol=1e-6)
417            )
418            & (
419                (fault_local_coordinates < 1)
420                | np.isclose(fault_local_coordinates, 1, atol=1e-6)
421            )
422        ):
423            raise ValueError("Specified coordinates do not lie in plane")
424        return np.clip(fault_local_coordinates, 0, 1)

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

Plane(bounds: numpy.ndarray)
bounds: numpy.ndarray
@staticmethod
def from_corners(corners: numpy.ndarray) -> Plane: 
187    @staticmethod
188    def from_corners(corners: np.ndarray) -> "Plane":
189        """Construct a plane point source from its corners.
190
191        Parameters
192        ----------
193        corners : np.ndarray
194            The corners in lat, lon, depth format. Has shape (4 x 3).
195
196        Returns
197        -------
198        Plane
199            The plane source representing this geometry.
200        """
201        return Plane(coordinates.wgs_depth_to_nztm(corners))

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.
corners: numpy.ndarray 
203    @property
204    def corners(self) -> np.ndarray:
205        """np.ndarray: The corners of the fault plane in (lat, lon, depth) format. The corners are the same as in corners_nztm."""
206        return coordinates.nztm_to_wgs_depth(self.bounds)

np.ndarray: The corners of the fault plane in (lat, lon, depth) format. The corners are the same as in corners_nztm.

length_m: float 
208    @property
209    def length_m(self) -> float:
210        """float: The length of the fault plane (in metres)."""
211        return np.linalg.norm(self.bounds[1] - self.bounds[0])

float: The length of the fault plane (in metres).

width_m: float 
213    @property
214    def width_m(self) -> float:
215        """float: The width of the fault plane (in metres)."""
216        return np.linalg.norm(self.bounds[-1] - self.bounds[0])

float: The width of the fault plane (in metres).

bottom_m: float 
218    @property
219    def bottom_m(self) -> float:
220        """float: The bottom depth (in metres)."""
221        return self.bounds[-1, -1]

float: The bottom depth (in metres).

top_m: float 
223    @property
224    def top_m(self) -> float:
225        """float: The top depth of the fault."""
226        return self.bounds[0, -1]

float: The top depth of the fault.

width: float 
228    @property
229    def width(self) -> float:
230        """float: The width of the fault plane (in kilometres)."""
231        return self.width_m / _KM_TO_M

float: The width of the fault plane (in kilometres).

length: float 
233    @property
234    def length(self) -> float:
235        """float: The length of the fault plane (in kilometres)."""
236        return self.length_m / _KM_TO_M

float: The length of the fault plane (in kilometres).

area: float 
238    @property
239    def area(self) -> float:
240        """float: The area of the plane (in km^2)."""
241        return self.length * self.width

float: The area of the plane (in km^2).

projected_width_m: float 
243    @property
244    def projected_width_m(self) -> float:
245        """float: The projected width of the fault plane (in metres)."""
246        return self.width_m * np.cos(np.radians(self.dip))

float: The projected width of the fault plane (in metres).

projected_width: float 
248    @property
249    def projected_width(self) -> float:
250        """float: The projected width of the fault plane (in kilometres)."""
251        return self.projected_width_m / _KM_TO_M

float: The projected width of the fault plane (in kilometres).

strike: float 
253    @property
254    def strike(self) -> float:
255        """float: The bearing of the strike direction of the fault (from north; in degrees)."""
256        north_direction = np.array([1, 0, 0])
257        up_direction = np.array([0, 0, 1])
258        strike_direction = self.bounds[1] - self.bounds[0]
259        return geo.oriented_bearing_wrt_normal(
260            north_direction, strike_direction, up_direction
261        )

float: The bearing of the strike direction of the fault (from north; in degrees).

dip_dir: float 
263    @property
264    def dip_dir(self) -> float:
265        """float: The bearing of the dip direction (from north; in degrees)."""
266        if np.isclose(self.dip, 90):
267            return 0  # TODO: Is this right for this case?
268        north_direction = np.array([1, 0, 0])
269        up_direction = np.array([0, 0, 1])
270        dip_direction = self.bounds[-1] - self.bounds[0]
271        dip_direction[-1] = 0
272        return geo.oriented_bearing_wrt_normal(
273            north_direction, dip_direction, up_direction
274        )

float: The bearing of the dip direction (from north; in degrees).

dip: float 
276    @property
277    def dip(self) -> float:
278        """float: The dip angle of the fault."""
279        return np.degrees(np.arcsin(np.abs(self.bottom_m - self.top_m) / self.width_m))

float: The dip angle of the fault.

geometry: shapely.geometry.polygon.Polygon 
281    @property
282    def geometry(self) -> shapely.Polygon:
283        """shapely.Polygon: A shapely geometry for the plane (projected onto the surface)."""
284        return shapely.Polygon(self.bounds)

shapely.Polygon: A shapely geometry for the plane (projected onto the surface).

@staticmethod
def from_centroid_strike_dip( centroid: numpy.ndarray, strike: float, dip_dir: Optional[float], dtop: float, dbottom: float, length: float, projected_width: float) -> Plane: 
286    @staticmethod
287    def from_centroid_strike_dip(
288        centroid: np.ndarray,
289        strike: float,
290        dip_dir: Optional[float],
291        dtop: float,
292        dbottom: float,
293        length: float,
294        projected_width: float,
295    ) -> "Plane":
296        """Create a fault plane from the centroid, strike, dip_dir, top, bottom, length, and width.
297
298        This is used for older descriptions of sources. Internally
299        converts everything to corners so self.strike ~ strike (but
300        not exactly due to rounding errors).
301
302        Parameters
303        ----------
304        centroid : np.ndarray
305            The centre of the fault plane in lat, lon coordinates.
306        strike : float
307            The strike of the fault (in degrees).
308        dip_dir : Optional[float]
309            The dip direction of the fault (in degrees). If None this is assumed to be strike + 90 degrees.
310        top : float
311            The top depth of the plane (in km).
312        bottom : float
313            The bottom depth of the plane (in km).
314        length : float
315            The length of the fault plane (in km).
316        projected_width : float
317            The projected width of the fault plane (in km).
318
319        Returns
320        -------
321        Plane
322            The fault plane with centre at `centroid`, and where the
323            parameters strike, dip_dir, top, bottom, length and width
324            match what is passed to this function.
325        """
326        corners = grid.grid_corners(
327            centroid,
328            strike,
329            dip_dir if dip_dir is not None else (strike + 90),
330            dtop,
331            dbottom,
332            length,
333            projected_width,
334        )
335        corners[[2, 3]] = corners[[3, 2]]
336        return Plane(coordinates.wgs_depth_to_nztm(np.array(corners)))

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 coordinates.
  • 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).
  • projected_width (float): The projected 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.
centroid: numpy.ndarray 
338    @property
339    def centroid(self) -> np.ndarray:
340        """np.ndarray: The centre of the fault plane."""
341        return self.fault_coordinates_to_wgs_depth_coordinates(np.array([1 / 2, 1 / 2]))

np.ndarray: The centre of the fault plane.

def fault_coordinates_to_wgs_depth_coordinates(self, plane_coordinates: numpy.ndarray) -> numpy.ndarray: 
343    def fault_coordinates_to_wgs_depth_coordinates(
344        self, plane_coordinates: np.ndarray
345    ) -> np.ndarray:
346        """Convert plane coordinates to nztm global coordinates.
347
348        Parameters
349        ----------
350        plane_coordinates : np.ndarray
351            Plane coordinates to convert. Plane coordinates are
352            2D coordinates (x, y) given for a fault plane (a plane), where x
353            represents displacement along the strike, and y
354            displacement along the dip (see diagram below). The
355            origin for plane coordinates is the centre of the fault.
356
357                          +x
358             0 0   ─────────────────>
359                ┌─────────────────────┐ │
360                │      < strike >     │ │
361                │                 ^   │ │
362                │                dip  │ │ +y
363                │                 v   │ │
364                │                     │ │
365                └─────────────────────┘ ∨
366                                       1,1
367
368        Returns
369        -------
370        np.ndarray
371            An 3d-vector of (lat, lon, depth) transformed coordinates.
372        """
373        origin = self.bounds[0]
374        top_right = self.bounds[1]
375        bottom_left = self.bounds[-1]
376        frame = np.vstack((top_right - origin, bottom_left - origin))
377
378        return coordinates.nztm_to_wgs_depth(origin + plane_coordinates @ frame)

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.
def wgs_depth_coordinates_to_fault_coordinates(self, global_coordinates: numpy.ndarray) -> numpy.ndarray: 
380    def wgs_depth_coordinates_to_fault_coordinates(
381        self,
382        global_coordinates: np.ndarray,
383    ) -> np.ndarray:
384        """Convert coordinates (lat, lon, depth) to plane coordinates (x, y).
385
386        See plane_coordinates_to_global_coordinates for a description of
387        plane coordinates.
388
389        Parameters
390        ----------
391        global_coordinates : np.ndarray
392            Global coordinates to convert.
393
394        Returns
395        -------
396        np.ndarray
397            The plane coordinates (x, y) representing the position of
398            global_coordinates on the fault plane.
399
400        Raises
401        ------
402        ValueError
403            If the given coordinates do not lie in the fault plane.
404        """
405        strike_direction = self.bounds[1, :2] - self.bounds[0, :2]
406        dip_direction = self.bounds[-1, :2] - self.bounds[0, :2]
407        offset = (
408            coordinates.wgs_depth_to_nztm(global_coordinates[:2]) - self.bounds[0, :2]
409        )
410        fault_local_coordinates, _, _, _ = np.linalg.lstsq(
411            np.array([strike_direction, dip_direction]).T, offset, rcond=None
412        )
413        if not np.all(
414            (
415                (fault_local_coordinates > 0)
416                | np.isclose(fault_local_coordinates, 0, atol=1e-6)
417            )
418            & (
419                (fault_local_coordinates < 1)
420                | np.isclose(fault_local_coordinates, 1, atol=1e-6)
421            )
422        ):
423            raise ValueError("Specified coordinates do not lie in plane")
424        return np.clip(fault_local_coordinates, 0, 1)

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.
@dataclasses.dataclass
class Fault: 
427@dataclasses.dataclass
428class Fault:
429    """A representation of a fault, consisting of one or more Planes.
430
431    This class represents a fault, which is composed of one or more Planes.
432    It provides methods for computing the area of the fault, getting the widths and
433    lengths of all fault planes, retrieving all corners of the fault, converting
434    global coordinates to fault coordinates, converting fault coordinates to global
435    coordinates, generating a random hypocentre location within the fault, and
436    computing the expected fault coordinates.
437
438    Attributes
439    ----------
440    planes : list[Plane]
441        A list containing all the Planes that constitute the fault.
442    """
443
444    planes: list[Plane]
445
446    @staticmethod
447    def from_corners(fault_corners: np.ndarray) -> "Fault":
448        """Construct a plane source geometry from the corners of the plane.
449
450        Parameters
451        ----------
452        corners : np.ndarray
453            The corners of the plane in lat, lon, depth format. Has shape (n x 4 x 3).
454
455        Returns
456        -------
457        Fault
458            The fault object representing this geometry.
459        """
460        return Fault([Plane.from_corners(corners) for corners in fault_corners])
461
462    def area(self) -> float:
463        """Compute the area of a fault.
464
465        Returns
466        -------
467        float
468            The area of the fault.
469        """
470        return self.width * np.sum(self.lengths)
471
472    @property
473    def lengths(self) -> np.ndarray:
474        """np.ndarray: A numpy array of each plane length (in km)."""
475        return np.array([fault.length for fault in self.planes])
476
477    @property
478    def length(self) -> float:
479        """float: The total length of each fault plane."""
480        return self.lengths.sum()
481
482    @property
483    def width(self) -> float:
484        """The width of the fault.
485
486        Returns
487        -------
488        float
489            The width of the first fault plane (A fault is assumed to
490            have planes of constant width).
491        """
492        return self.planes[0].width
493
494    @property
495    def dip_dir(self) -> float:
496        """float: The dip direction of the first fault plane (A fault is assumed to have planes of constant dip direction)."""
497        return self.planes[0].dip_dir
498
499    @property
500    def corners(self) -> np.ndarray:
501        """np.ndarray of shape (4n x 3): The corners in (lat, lon, depth) format of each fault plane in the fault, stacked vertically."""
502        return np.vstack([plane.corners for plane in self.planes])
503
504    @property
505    def bounds(self) -> np.ndarray:
506        """np.ndarray of shape (4n x 3): The corners in NZTM format of each fault plane in the fault, stacked vertically."""
507        return np.vstack([plane.bounds for plane in self.planes])
508
509    @property
510    def centroid(self) -> np.ndarray:
511        """np.ndarray: The centre of the fault."""
512        return self.fault_coordinates_to_wgs_depth_coordinates(np.array([1 / 2, 1 / 2]))
513
514    @property
515    def geometry(self) -> shapely.Polygon:
516        """shapely.Geometry: A shapely geometry for the fault (projected onto the surface)."""
517        return shapely.normalize(
518            shapely.union_all([plane.geometry for plane in self.planes])
519        )
520
521    def wgs_depth_coordinates_to_fault_coordinates(
522        self, global_coordinates: np.ndarray
523    ) -> np.ndarray:
524        """Convert global coordinates in (lat, lon, depth) format to fault coordinates.
525
526        Fault coordinates are a tuple (s, d) where s is the distance
527        from the top left, and d the distance from the top of the
528        fault (refer to the diagram). The coordinates are normalised
529        such that (0, 0) is the top left and (1, 1) the bottom right.
530
531        (0, 0)
532          ┌──────────────────────┬──────┐
533          │          |           │      │
534          │          |           │      │
535          │          | d         │      │
536          │          |           │      │
537          ├----------*           │      │
538          │    s     ^           │      │
539          │          |           │      │
540          │          |           │      │
541          │          |           │      │
542          └──────────|───────────┴──────┘
543                     +                    (1, 1)
544                  point: (s, d)
545
546        Parameters
547        ----------
548        global_coordinates : np.ndarray of shape (1 x 3)
549            The global coordinates to convert.
550
551        Returns
552        -------
553        np.ndarray
554            The fault coordinates.
555
556        Raises
557        ------
558        ValueError
559            If the given point does not lie on the fault.
560
561        """
562        # the left edges as a cumulative proportion of the fault length (e.g. [0.1, ..., 0.8])
563        left_edges = self.lengths.cumsum() / self.length
564        left_edges = np.insert(left_edges, 0, 0)
565        for i, plane in enumerate(self.planes):
566            try:
567                plane_coordinates = plane.wgs_depth_coordinates_to_fault_coordinates(
568                    global_coordinates
569                )
570                return np.array([left_edges[i], 0]) + plane_coordinates * np.array(
571                    [left_edges[i + 1] - left_edges[i], 1]
572                )
573            except ValueError:
574                continue
575        raise ValueError("Given coordinates are not on fault.")
576
577    def fault_coordinates_to_wgs_depth_coordinates(
578        self, fault_coordinates: np.ndarray
579    ) -> np.ndarray:
580        """Convert fault coordinates to global coordinates.
581
582        See global_coordinates_to_fault_coordinates for a description of fault
583        coordinates.
584
585        Parameters
586        ----------
587        fault_coordinates : np.ndarray
588            The fault coordinates of the point.
589
590        Returns
591        -------
592        np.ndarray
593            The global coordinates (lat, lon, depth) for this point.
594        """
595        # the right edges as a cumulative proportion of the fault length (e.g. [0.1, ..., 0.8])
596        right_edges = self.lengths.cumsum() / self.length
597        for i in range(len(self.planes)):
598            if fault_coordinates[0] < right_edges[i] or np.isclose(
599                fault_coordinates[0], right_edges[i]
600            ):
601                break
602        fault_segment_index = i
603        left_proportion = (
604            right_edges[fault_segment_index - 1] if fault_segment_index > 0 else 0
605        )
606        right_proportion = right_edges[fault_segment_index]
607        segment_proportion = (fault_coordinates[0] - left_proportion) / (
608            right_proportion - left_proportion
609        )
610
611        return self.planes[
612            fault_segment_index
613        ].fault_coordinates_to_wgs_depth_coordinates(
614            np.array([segment_proportion, fault_coordinates[1]])
615        )

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.
Fault(planes: list[Plane])
planes: list[Plane]
@staticmethod
def from_corners(fault_corners: numpy.ndarray) -> Fault: 
446    @staticmethod
447    def from_corners(fault_corners: np.ndarray) -> "Fault":
448        """Construct a plane source geometry from the corners of the plane.
449
450        Parameters
451        ----------
452        corners : np.ndarray
453            The corners of the plane in lat, lon, depth format. Has shape (n x 4 x 3).
454
455        Returns
456        -------
457        Fault
458            The fault object representing this geometry.
459        """
460        return Fault([Plane.from_corners(corners) for corners in fault_corners])

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.
def area(self) -> float: 
462    def area(self) -> float:
463        """Compute the area of a fault.
464
465        Returns
466        -------
467        float
468            The area of the fault.
469        """
470        return self.width * np.sum(self.lengths)

Compute the area of a fault.

Returns
  • float: The area of the fault.
lengths: numpy.ndarray 
472    @property
473    def lengths(self) -> np.ndarray:
474        """np.ndarray: A numpy array of each plane length (in km)."""
475        return np.array([fault.length for fault in self.planes])

np.ndarray: A numpy array of each plane length (in km).

length: float 
477    @property
478    def length(self) -> float:
479        """float: The total length of each fault plane."""
480        return self.lengths.sum()

float: The total length of each fault plane.

width: float 
482    @property
483    def width(self) -> float:
484        """The width of the fault.
485
486        Returns
487        -------
488        float
489            The width of the first fault plane (A fault is assumed to
490            have planes of constant width).
491        """
492        return self.planes[0].width

The width of the fault.

Returns
  • float: The width of the first fault plane (A fault is assumed to have planes of constant width).
dip_dir: float 
494    @property
495    def dip_dir(self) -> float:
496        """float: The dip direction of the first fault plane (A fault is assumed to have planes of constant dip direction)."""
497        return self.planes[0].dip_dir

float: The dip direction of the first fault plane (A fault is assumed to have planes of constant dip direction).

corners: numpy.ndarray 
499    @property
500    def corners(self) -> np.ndarray:
501        """np.ndarray of shape (4n x 3): The corners in (lat, lon, depth) format of each fault plane in the fault, stacked vertically."""
502        return np.vstack([plane.corners for plane in self.planes])

np.ndarray of shape (4n x 3): The corners in (lat, lon, depth) format of each fault plane in the fault, stacked vertically.

bounds: numpy.ndarray 
504    @property
505    def bounds(self) -> np.ndarray:
506        """np.ndarray of shape (4n x 3): The corners in NZTM format of each fault plane in the fault, stacked vertically."""
507        return np.vstack([plane.bounds for plane in self.planes])

np.ndarray of shape (4n x 3): The corners in NZTM format of each fault plane in the fault, stacked vertically.

centroid: numpy.ndarray 
509    @property
510    def centroid(self) -> np.ndarray:
511        """np.ndarray: The centre of the fault."""
512        return self.fault_coordinates_to_wgs_depth_coordinates(np.array([1 / 2, 1 / 2]))

np.ndarray: The centre of the fault.

geometry: shapely.geometry.polygon.Polygon 
514    @property
515    def geometry(self) -> shapely.Polygon:
516        """shapely.Geometry: A shapely geometry for the fault (projected onto the surface)."""
517        return shapely.normalize(
518            shapely.union_all([plane.geometry for plane in self.planes])
519        )

shapely.Geometry: A shapely geometry for the fault (projected onto the surface).

def wgs_depth_coordinates_to_fault_coordinates(self, global_coordinates: numpy.ndarray) -> numpy.ndarray: 
521    def wgs_depth_coordinates_to_fault_coordinates(
522        self, global_coordinates: np.ndarray
523    ) -> np.ndarray:
524        """Convert global coordinates in (lat, lon, depth) format to fault coordinates.
525
526        Fault coordinates are a tuple (s, d) where s is the distance
527        from the top left, and d the distance from the top of the
528        fault (refer to the diagram). The coordinates are normalised
529        such that (0, 0) is the top left and (1, 1) the bottom right.
530
531        (0, 0)
532          ┌──────────────────────┬──────┐
533          │          |           │      │
534          │          |           │      │
535          │          | d         │      │
536          │          |           │      │
537          ├----------*           │      │
538          │    s     ^           │      │
539          │          |           │      │
540          │          |           │      │
541          │          |           │      │
542          └──────────|───────────┴──────┘
543                     +                    (1, 1)
544                  point: (s, d)
545
546        Parameters
547        ----------
548        global_coordinates : np.ndarray of shape (1 x 3)
549            The global coordinates to convert.
550
551        Returns
552        -------
553        np.ndarray
554            The fault coordinates.
555
556        Raises
557        ------
558        ValueError
559            If the given point does not lie on the fault.
560
561        """
562        # the left edges as a cumulative proportion of the fault length (e.g. [0.1, ..., 0.8])
563        left_edges = self.lengths.cumsum() / self.length
564        left_edges = np.insert(left_edges, 0, 0)
565        for i, plane in enumerate(self.planes):
566            try:
567                plane_coordinates = plane.wgs_depth_coordinates_to_fault_coordinates(
568                    global_coordinates
569                )
570                return np.array([left_edges[i], 0]) + plane_coordinates * np.array(
571                    [left_edges[i + 1] - left_edges[i], 1]
572                )
573            except ValueError:
574                continue
575        raise ValueError("Given coordinates are not on fault.")

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.
def fault_coordinates_to_wgs_depth_coordinates(self, fault_coordinates: numpy.ndarray) -> numpy.ndarray: 
577    def fault_coordinates_to_wgs_depth_coordinates(
578        self, fault_coordinates: np.ndarray
579    ) -> np.ndarray:
580        """Convert fault coordinates to global coordinates.
581
582        See global_coordinates_to_fault_coordinates for a description of fault
583        coordinates.
584
585        Parameters
586        ----------
587        fault_coordinates : np.ndarray
588            The fault coordinates of the point.
589
590        Returns
591        -------
592        np.ndarray
593            The global coordinates (lat, lon, depth) for this point.
594        """
595        # the right edges as a cumulative proportion of the fault length (e.g. [0.1, ..., 0.8])
596        right_edges = self.lengths.cumsum() / self.length
597        for i in range(len(self.planes)):
598            if fault_coordinates[0] < right_edges[i] or np.isclose(
599                fault_coordinates[0], right_edges[i]
600            ):
601                break
602        fault_segment_index = i
603        left_proportion = (
604            right_edges[fault_segment_index - 1] if fault_segment_index > 0 else 0
605        )
606        right_proportion = right_edges[fault_segment_index]
607        segment_proportion = (fault_coordinates[0] - left_proportion) / (
608            right_proportion - left_proportion
609        )
610
611        return self.planes[
612            fault_segment_index
613        ].fault_coordinates_to_wgs_depth_coordinates(
614            np.array([segment_proportion, fault_coordinates[1]])
615        )

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.
class IsSource(typing.Protocol): 
618class IsSource(Protocol):
619    """Type definition for a source with local coordinates."""
620
621    bounds: np.ndarray
622
623    def fault_coordinates_to_wgs_depth_coordinates(
624        self,
625        fault_coordinates: np.ndarray,
626    ) -> np.ndarray: ...
627
628    def wgs_depth_coordinates_to_fault_coordinates(
629        self,
630        fault_coordinates: np.ndarray,
631    ) -> np.ndarray: ...

Type definition for a source with local coordinates.

IsSource(*args, **kwargs) 
1767def _no_init_or_replace_init(self, *args, **kwargs):
1768    cls = type(self)
1769
1770    if cls._is_protocol:
1771        raise TypeError('Protocols cannot be instantiated')
1772
1773    # Already using a custom `__init__`. No need to calculate correct
1774    # `__init__` to call. This can lead to RecursionError. See bpo-45121.
1775    if cls.__init__ is not _no_init_or_replace_init:
1776        return
1777
1778    # Initially, `__init__` of a protocol subclass is set to `_no_init_or_replace_init`.
1779    # The first instantiation of the subclass will call `_no_init_or_replace_init` which
1780    # searches for a proper new `__init__` in the MRO. The new `__init__`
1781    # replaces the subclass' old `__init__` (ie `_no_init_or_replace_init`). Subsequent
1782    # instantiation of the protocol subclass will thus use the new
1783    # `__init__` and no longer call `_no_init_or_replace_init`.
1784    for base in cls.__mro__:
1785        init = base.__dict__.get('__init__', _no_init_or_replace_init)
1786        if init is not _no_init_or_replace_init:
1787            cls.__init__ = init
1788            break
1789    else:
1790        # should not happen
1791        cls.__init__ = object.__init__
1792
1793    cls.__init__(self, *args, **kwargs)
bounds: numpy.ndarray
def fault_coordinates_to_wgs_depth_coordinates(self, fault_coordinates: numpy.ndarray) -> numpy.ndarray: 
623    def fault_coordinates_to_wgs_depth_coordinates(
624        self,
625        fault_coordinates: np.ndarray,
626    ) -> np.ndarray: ...
def wgs_depth_coordinates_to_fault_coordinates(self, fault_coordinates: numpy.ndarray) -> numpy.ndarray: 
628    def wgs_depth_coordinates_to_fault_coordinates(
629        self,
630        fault_coordinates: np.ndarray,
631    ) -> np.ndarray: ...
def closest_point_between_sources( source_a: IsSource, source_b: IsSource) -> tuple[numpy.ndarray, numpy.ndarray]: 
634def closest_point_between_sources(
635    source_a: IsSource, source_b: IsSource
636) -> tuple[np.ndarray, np.ndarray]:
637    """Find the closest point between two sources that have local coordinates.
638
639    Parameters
640    ----------
641    source_a : HasCoordinates
642        The first source. Must have a two-dimensional fault coordinate system.
643    source_b : HasCoordinates
644        The first source. Must have a two-dimensional fault coordinate system.
645
646    Raises
647    ------
648    ValueError
649        Raised when we are unable to converge on the closest points between sources.
650
651    Returns
652    -------
653    source_a_coordinates : np.ndarray
654        The source-local coordinates of the closest point on source a.
655    source_b_coordinates : np.ndarray
656        The source-local coordinates of the closest point on source b.
657    """
658
659    def fault_coordinate_distance(fault_coordinates: np.ndarray) -> float:
660        source_a_global_coordinates = (
661            source_a.fault_coordinates_to_wgs_depth_coordinates(fault_coordinates[:2])
662        )
663        source_b_global_coordinates = (
664            source_b.fault_coordinates_to_wgs_depth_coordinates(fault_coordinates[2:])
665        )
666        return coordinates.wgs_depth_to_nztm(
667            source_a_global_coordinates
668        ) - coordinates.wgs_depth_to_nztm(source_b_global_coordinates)
669
670    res = sp.optimize.least_squares(
671        fault_coordinate_distance,
672        np.array([1 / 2, 1 / 2, 1 / 2, 1 / 2]),
673        bounds=([0] * 4, [1] * 4),
674        gtol=1e-5,
675        ftol=1e-5,
676        xtol=1e-5,
677    )
678
679    if not res.success and res.status != 0:
680        raise ValueError(
681            f"Optimisation failed to converge for provided sources: {res.message} with x = {res.x}"
682        )
683    return res.x[:2], res.x[2:]

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.