From eed9805cae24156ca3769e271ba3a174045b1687 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Tue, 2 Jul 2024 12:31:52 +1200 Subject: [PATCH 01/11] add rupture propagation module --- source_modelling/rupture_propagation.py | 210 ++++++++++++++++++++++++ 1 file changed, 210 insertions(+) create mode 100644 source_modelling/rupture_propagation.py diff --git a/source_modelling/rupture_propagation.py b/source_modelling/rupture_propagation.py new file mode 100644 index 0000000..b7a71da --- /dev/null +++ b/source_modelling/rupture_propagation.py @@ -0,0 +1,210 @@ +""" +Rupture Propagation Module + +This module provides functions for computing likely rupture paths from +information about the distances between faults. + +Functions: + - shaw_dieterich_distance_model: + Compute fault jump probabilities using the Shaw-Dieterich distance model. + - prune_distance_graph: Prune the distance graph based on a cutoff value. + - probability_graph: + Convert a graph of distances between faults into a graph of jump + probabilities using the Shaw-Dieterich model. + - probabilistic_minimum_spanning_tree: Generate a probabilistic minimum spanning tree. + +Typing Aliases: + - DistanceGraph: A graph representing distances between faults. + - ProbabilityGraph: A graph representing jump probabilities between faults. + - RuptureCausalityTree: A tree representing the causality of ruptures between faults. +""" + +from collections import defaultdict, namedtuple +from typing import Optional, Tuple + +import numpy as np + +from source_modelling import sources + +DistanceGraph = dict[str, dict[str, int]] +ProbabilityGraph = defaultdict[str, dict[str, float]] +RuptureCausalityTree = dict[str, Optional[str]] + +JumpPair = namedtuple("JumpPair", ["from", "to"]) + + +def shaw_dieterich_distance_model(distance: float, d0, delta) -> float: + """ + Compute fault jump probabilities using the Shaw-Dieterich distance model[0]. + + Parameters + ---------- + distance : float + The distance between two faults. + d0 : float + The characteristic distance parameter. + delta : float + The characteristic slip distance parameter. + + Returns + ------- + float + The calculated probability. + + References + ---------- + [0]: Shaw, B. E., & Dieterich, J. H. (2007). Probabilities for jumping fault + segment stepovers. Geophysical Research Letters, 34(1). + """ + return min(1, np.exp(-(distance - delta) / d0)) + + +def prune_distance_graph(distances: DistanceGraph, cutoff: int) -> DistanceGraph: + """ + Prune the distance graph based on a cutoff value. + + Parameters + ---------- + distances : DistanceGraph + The graph of distances between faults. + cutoff : int + The cutoff distance in metres. + + Returns + ------- + DistanceGraph + A copy of the input distance graph, keeping only edges that are less + than the cutoff. + """ + return { + fault_u: { + fault_v: distance + for fault_v, distance in neighbours_fault_u.items() + if distance < cutoff + } + for fault_u, neighbours_fault_u in distances.items() + } + + +def probability_graph( + distances: DistanceGraph, d0: float = 3, delta: float = 1 +) -> ProbabilityGraph: + """ + Convert a graph of distances between faults into a graph of jump + probablities using the Shaw-Dieterich model. + + Parameters + ---------- + distances : DistanceGraph + The distance graph between faults. + d0 : float, optional + The d0 parameter for the Shaw_Dieterich model. See `shaw_dieterich_distance_model`. + delta : float, optional + The delta parameter for the Shaw_Dieterich model. See `shaw_dieterich_distance_model`. + + Returns + ------- + ProbabilityGraph + The graph with faults as vertices. Each edge (fault_u, fault_v) + has a log-probability -p as a weight. The log-probability -p here + is the negative of the log-probability a rupture propogates from + fault_u to fault_v, relative to the probability it propogates to + any of the other neighbours of fault_u. + """ + probabilities_raw = { + fault_u: { + fault_v: shaw_dieterich_distance_model(distance, d0, delta) + for fault_v, distance in neighbours_fault_u.items() + } + for fault_u, neighbours_fault_u in distances.items() + } + + probabilities_log = defaultdict(dict) + for fault_u, neighbours_fault_u in probabilities_raw.items(): + normalising_constant = sum(neighbours_fault_u.values()) + for fault_v, prob in neighbours_fault_u.items(): + probabilities_log[fault_u][fault_v] = normalising_constant - np.log(prob) + return probabilities_log + + +def probabilistic_minimum_spanning_tree( + probability_graph: ProbabilityGraph, initial_fault: str +) -> RuptureCausalityTree: + """ + Generate a probabilistic minimum spanning tree. + + The minimum spanning tree of the probability graph represents rupture + causality tree with the highest likelihood of occuring assuming that + rupture jumps are independent. + + NOTE: While the overall probability of the minimum spanning tree is high, + the paths from the initial fault may not be the most likely. + + Parameters + ---------- + probability_graph : ProbabilityGraph + The probability graph. + initial_fault : str + The initial fault. + + Returns + ------- + RuptureCausalityTree + The probabilistic minimum spanning tree. + """ + rupture_causality_tree = {initial_fault: None} + path_probabilities = defaultdict(lambda: np.inf) + path_probabilities[initial_fault] = 0 + processed_faults = set() + for _ in range(len(probability_graph)): + current_fault = min( + (fault for fault in probability_graph if fault not in processed_faults), + key=lambda fault: path_probabilities[fault], + ) + processed_faults.add(current_fault) + for fault_neighbour, log_probability in probability_graph[ + current_fault + ].items(): + if ( + fault_neighbour not in processed_faults + and log_probability < path_probabilities[fault_neighbour] + ): + path_probabilities[fault_neighbour] = log_probability + rupture_causality_tree[fault_neighbour] = current_fault + return rupture_causality_tree + + +def estimate_most_likely_rupture_propagation( + sources_map: dict[str, sources.IsSource], + initial_source: str, + jump_impossibility_limit_distance: int = 15000, +) -> RuptureCausalityTree: + distance_graph = { + source_a_name: { + source_b_name: sources.closest_point_between_sources(source_a, source_b) + for source_b_name, source_b in sources_map.items() + } + for source_a_name, source_a in sources_map.items() + } + distance_graph = prune_distance_graph( + distance_graph, jump_impossibility_limit_distance + ) + jump_probability_graph = probability_graph(distance_graph) + return probabilistic_minimum_spanning_tree( + jump_probability_graph, initial_fault=initial_source + ) + + +def jump_points_from_rupture_tree( + source_map: dict[str, sources.IsSource], + rupture_causality_tree: RuptureCausalityTree, +) -> dict[str, JumpPair]: + jump_points = {} + for source, parent in source_map.items(): + if parent is None: + continue + source_point, parent_point = sources.closest_point_between_sources( + source_map[source], source_map[parent] + ) + jump_points[source] = JumpPair(parent_point, source_point) + return jump_points From 1b1a7f7ba7a272f03a8b7416ebfb2511a3a7e821 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Wed, 3 Jul 2024 10:41:09 +1200 Subject: [PATCH 02/11] make Fault.corners a property for consistency --- source_modelling/sources.py | 1 + 1 file changed, 1 insertion(+) diff --git a/source_modelling/sources.py b/source_modelling/sources.py index 54c362a..0faf59a 100644 --- a/source_modelling/sources.py +++ b/source_modelling/sources.py @@ -597,6 +597,7 @@ def dip_dir(self) -> float: """ return self.planes[0].dip_dir + @property def corners(self) -> np.ndarray: """Get all corners of a fault. From 25e9598b1c05b7be601197d0da1a01133ac91fab Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Wed, 3 Jul 2024 10:41:32 +1200 Subject: [PATCH 03/11] fix bug in dip calculation --- source_modelling/sources.py | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/source_modelling/sources.py b/source_modelling/sources.py index 0faf59a..4de89ec 100644 --- a/source_modelling/sources.py +++ b/source_modelling/sources.py @@ -238,6 +238,16 @@ def bottom_m(self) -> float: """ return self.bounds[-1, -1] + @property + def top_m(self) -> float: + """ + Returns + ------- + float + The top depth of the fault. + """ + return self.bounds[0, -1] + @property def width(self) -> float: """ @@ -321,7 +331,7 @@ def dip(self) -> float: float The dip angle of the fault. """ - return np.degrees(np.arcsin(np.abs(self.bottom_m) / self.width_m)) + return np.degrees(np.arcsin(np.abs(self.bottom_m - self.top_m) / self.width_m)) @staticmethod def from_centroid_strike_dip( From 8aba38faeb0730ab8c5f9fac6f8e33392b3cb57b Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Wed, 3 Jul 2024 10:42:01 +1200 Subject: [PATCH 04/11] fix namedtuple using keyword arguments --- source_modelling/rupture_propagation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source_modelling/rupture_propagation.py b/source_modelling/rupture_propagation.py index b7a71da..954fe70 100644 --- a/source_modelling/rupture_propagation.py +++ b/source_modelling/rupture_propagation.py @@ -30,7 +30,7 @@ ProbabilityGraph = defaultdict[str, dict[str, float]] RuptureCausalityTree = dict[str, Optional[str]] -JumpPair = namedtuple("JumpPair", ["from", "to"]) +JumpPair = namedtuple("JumpPair", ["from_point", "to_point"]) def shaw_dieterich_distance_model(distance: float, d0, delta) -> float: From 6e1d66ad8f61166abdcb760036dbd779a3592471 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Wed, 3 Jul 2024 16:29:15 +1200 Subject: [PATCH 05/11] correct comment --- source_modelling/sources.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source_modelling/sources.py b/source_modelling/sources.py index 6a84abc..13e796f 100644 --- a/source_modelling/sources.py +++ b/source_modelling/sources.py @@ -674,7 +674,7 @@ def wgs_depth_coordinates_to_fault_coordinates( If the given point does not lie on the fault. """ - # the right edges as a cumulative proportion of the fault length (e.g. [0.1, ..., 0.8]) + # the left edges as a cumulative proportion of the fault length (e.g. [0.1, ..., 0.8]) left_edges = self.lengths.cumsum() / self.length left_edges = np.insert(left_edges, 0, 0) for i, plane in enumerate(self.planes): From 23c69de83f80ce36be183dae1d26523879039cfc Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Thu, 4 Jul 2024 15:22:46 +1200 Subject: [PATCH 06/11] add centroid, make algorithms more numerically stable --- source_modelling/rupture_propagation.py | 26 ++++++++++++++------- source_modelling/sources.py | 30 ++++++++++++++++++++----- 2 files changed, 43 insertions(+), 13 deletions(-) diff --git a/source_modelling/rupture_propagation.py b/source_modelling/rupture_propagation.py index 954fe70..a85cce2 100644 --- a/source_modelling/rupture_propagation.py +++ b/source_modelling/rupture_propagation.py @@ -21,6 +21,8 @@ from collections import defaultdict, namedtuple from typing import Optional, Tuple +from qcore import coordinates + import numpy as np @@ -113,17 +115,19 @@ def probability_graph( """ probabilities_raw = { fault_u: { - fault_v: shaw_dieterich_distance_model(distance, d0, delta) + fault_v: shaw_dieterich_distance_model(distance / 1000, d0, delta) for fault_v, distance in neighbours_fault_u.items() } for fault_u, neighbours_fault_u in distances.items() } - probabilities_log = defaultdict(dict) for fault_u, neighbours_fault_u in probabilities_raw.items(): normalising_constant = sum(neighbours_fault_u.values()) + if normalising_constant == 0: + for fault_v, _ in neighbours_fault_u.items(): + probabilities_log[fault_u][fault_v] = 1/len(neighbours_fault_u) for fault_v, prob in neighbours_fault_u.items(): - probabilities_log[fault_u][fault_v] = normalising_constant - np.log(prob) + probabilities_log[fault_u][fault_v] = prob / normalising_constant return probabilities_log @@ -162,17 +166,22 @@ def probabilistic_minimum_spanning_tree( key=lambda fault: path_probabilities[fault], ) processed_faults.add(current_fault) - for fault_neighbour, log_probability in probability_graph[ + for fault_neighbour, probability in probability_graph[ current_fault ].items(): if ( fault_neighbour not in processed_faults - and log_probability < path_probabilities[fault_neighbour] + and probability < path_probabilities[fault_neighbour] ): - path_probabilities[fault_neighbour] = log_probability + path_probabilities[fault_neighbour] = probability rupture_causality_tree[fault_neighbour] = current_fault return rupture_causality_tree +def distance_between(source_a: sources.IsSource, source_b: sources.IsSource, source_a_point: np.ndarray, source_b_point: np.ndarray) -> float: + global_point_a = source_a.fault_coordinates_to_wgs_depth_coordinates(source_a_point) + global_point_b = source_b.fault_coordinates_to_wgs_depth_coordinates(source_b_point) + return coordinates.distance_between_wgs_depth_coordinates(global_point_a, global_point_b) + def estimate_most_likely_rupture_propagation( sources_map: dict[str, sources.IsSource], @@ -181,8 +190,9 @@ def estimate_most_likely_rupture_propagation( ) -> RuptureCausalityTree: distance_graph = { source_a_name: { - source_b_name: sources.closest_point_between_sources(source_a, source_b) + source_b_name: distance_between(source_a, source_b, *sources.closest_point_between_sources(source_a, source_b)) for source_b_name, source_b in sources_map.items() + if source_a_name != source_b_name } for source_a_name, source_a in sources_map.items() } @@ -200,7 +210,7 @@ def jump_points_from_rupture_tree( rupture_causality_tree: RuptureCausalityTree, ) -> dict[str, JumpPair]: jump_points = {} - for source, parent in source_map.items(): + for source, parent in rupture_causality_tree.items(): if parent is None: continue source_point, parent_point = sources.closest_point_between_sources( diff --git a/source_modelling/sources.py b/source_modelling/sources.py index 13e796f..b5fca53 100644 --- a/source_modelling/sources.py +++ b/source_modelling/sources.py @@ -21,7 +21,6 @@ import numpy as np import scipy as sp - from qcore import coordinates, geo, grid _KM_TO_M = 1000 @@ -116,6 +115,10 @@ def width(self) -> float: """ return self.width_m / _KM_TO_M + @property + def centroid(self) -> np.ndarray: + return self.coordinates + def fault_coordinates_to_wgs_depth_coordinates( self, fault_coordinates: np.ndarray ) -> np.ndarray: @@ -398,6 +401,10 @@ def from_centroid_strike_dip( ) return Plane(corners) + @property + def centroid(self) -> np.ndarray: + return self.fault_coordinates_to_wgs_depth_coordinates(np.array([1 / 2, 1 / 2])) + def fault_coordinates_to_wgs_depth_coordinates( self, plane_coordinates: np.ndarray ) -> np.ndarray: @@ -633,6 +640,10 @@ def bounds(self) -> np.ndarray: """ return np.vstack([plane.bounds for plane in self.planes]) + @property + def centroid(self) -> np.ndarray: + return self.fault_coordinates_to_wgs_depth_coordinates(np.array([1 / 2, 1 / 2])) + def wgs_depth_coordinates_to_fault_coordinates( self, global_coordinates: np.ndarray ) -> np.ndarray: @@ -707,10 +718,13 @@ def fault_coordinates_to_wgs_depth_coordinates( np.ndarray The global coordinates (lat, lon, depth) for this point. """ - + # the right edges as a cumulative proportion of the fault length (e.g. [0.1, ..., 0.8]) right_edges = self.lengths.cumsum() / self.length - fault_segment_index = np.searchsorted(right_edges, fault_coordinates[0]) + 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 ) @@ -720,8 +734,10 @@ def fault_coordinates_to_wgs_depth_coordinates( else 1 ) segment_proportion = (fault_coordinates[0] - left_proportion) / ( - right_proportion - left_proportion + right_proportion - left_proportion ) + if fault_segment_index >= len(self.planes): + breakpoint() return self.planes[ fault_segment_index ].fault_coordinates_to_wgs_depth_coordinates( @@ -769,7 +785,6 @@ def closest_point_between_sources( 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]) @@ -785,11 +800,16 @@ def fault_coordinate_distance(fault_coordinates: np.ndarray) -> float: fault_coordinate_distance, np.array([1 / 2, 1 / 2, 1 / 2, 1 / 2]), bounds=[(0, 1)] * 4, + options={ + 'ftol': 1e-2 + } ) if not res.success: + breakpoint() raise ValueError( f"Optimisation failed to converge for provided sources: {res.message}" ) + return res.x[:2], res.x[2:] From 76bb1d64973b9ce60add72dc41a7e2b47225c67d Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Fri, 5 Jul 2024 15:41:59 +1200 Subject: [PATCH 07/11] small fix for centroid -> plane --- source_modelling/sources.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source_modelling/sources.py b/source_modelling/sources.py index b5fca53..6a91697 100644 --- a/source_modelling/sources.py +++ b/source_modelling/sources.py @@ -399,7 +399,7 @@ def from_centroid_strike_dip( length, width, ) - return Plane(corners) + return Plane(coordinates.wgs_depth_to_nztm(corners)) @property def centroid(self) -> np.ndarray: From b1126a959deea887fdf3d7743c0ced1a6f85250a Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Fri, 5 Jul 2024 15:42:12 +1200 Subject: [PATCH 08/11] format and sort --- source_modelling/sources.py | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/source_modelling/sources.py b/source_modelling/sources.py index 6a91697..5cc6f79 100644 --- a/source_modelling/sources.py +++ b/source_modelling/sources.py @@ -21,6 +21,7 @@ import numpy as np import scipy as sp + from qcore import coordinates, geo, grid _KM_TO_M = 1000 @@ -718,11 +719,13 @@ def fault_coordinates_to_wgs_depth_coordinates( 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]): + if fault_coordinates[0] < right_edges[i] or np.isclose( + fault_coordinates[0], right_edges[i] + ): break fault_segment_index = i left_proportion = ( @@ -734,7 +737,7 @@ def fault_coordinates_to_wgs_depth_coordinates( else 1 ) segment_proportion = (fault_coordinates[0] - left_proportion) / ( - right_proportion - left_proportion + right_proportion - left_proportion ) if fault_segment_index >= len(self.planes): breakpoint() @@ -785,6 +788,7 @@ def closest_point_between_sources( 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]) @@ -800,9 +804,7 @@ def fault_coordinate_distance(fault_coordinates: np.ndarray) -> float: fault_coordinate_distance, np.array([1 / 2, 1 / 2, 1 / 2, 1 / 2]), bounds=[(0, 1)] * 4, - options={ - 'ftol': 1e-2 - } + options={"ftol": 1e-2}, ) if not res.success: @@ -810,6 +812,5 @@ def fault_coordinate_distance(fault_coordinates: np.ndarray) -> float: raise ValueError( f"Optimisation failed to converge for provided sources: {res.message}" ) - return res.x[:2], res.x[2:] From b3793e9d7775453417f305475bbb4f505caa52dc Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Fri, 5 Jul 2024 15:44:15 +1200 Subject: [PATCH 09/11] formatting --- source_modelling/rupture_propagation.py | 27 ++++++++++++++++--------- 1 file changed, 18 insertions(+), 9 deletions(-) diff --git a/source_modelling/rupture_propagation.py b/source_modelling/rupture_propagation.py index a85cce2..1516797 100644 --- a/source_modelling/rupture_propagation.py +++ b/source_modelling/rupture_propagation.py @@ -21,11 +21,10 @@ from collections import defaultdict, namedtuple from typing import Optional, Tuple -from qcore import coordinates - import numpy as np +from qcore import coordinates from source_modelling import sources DistanceGraph = dict[str, dict[str, int]] @@ -125,7 +124,7 @@ def probability_graph( normalising_constant = sum(neighbours_fault_u.values()) if normalising_constant == 0: for fault_v, _ in neighbours_fault_u.items(): - probabilities_log[fault_u][fault_v] = 1/len(neighbours_fault_u) + probabilities_log[fault_u][fault_v] = 1 / len(neighbours_fault_u) for fault_v, prob in neighbours_fault_u.items(): probabilities_log[fault_u][fault_v] = prob / normalising_constant return probabilities_log @@ -166,9 +165,7 @@ def probabilistic_minimum_spanning_tree( key=lambda fault: path_probabilities[fault], ) processed_faults.add(current_fault) - for fault_neighbour, probability in probability_graph[ - current_fault - ].items(): + for fault_neighbour, probability in probability_graph[current_fault].items(): if ( fault_neighbour not in processed_faults and probability < path_probabilities[fault_neighbour] @@ -177,10 +174,18 @@ def probabilistic_minimum_spanning_tree( rupture_causality_tree[fault_neighbour] = current_fault return rupture_causality_tree -def distance_between(source_a: sources.IsSource, source_b: sources.IsSource, source_a_point: np.ndarray, source_b_point: np.ndarray) -> float: + +def distance_between( + source_a: sources.IsSource, + source_b: sources.IsSource, + source_a_point: np.ndarray, + source_b_point: np.ndarray, +) -> float: global_point_a = source_a.fault_coordinates_to_wgs_depth_coordinates(source_a_point) global_point_b = source_b.fault_coordinates_to_wgs_depth_coordinates(source_b_point) - return coordinates.distance_between_wgs_depth_coordinates(global_point_a, global_point_b) + return coordinates.distance_between_wgs_depth_coordinates( + global_point_a, global_point_b + ) def estimate_most_likely_rupture_propagation( @@ -190,7 +195,11 @@ def estimate_most_likely_rupture_propagation( ) -> RuptureCausalityTree: distance_graph = { source_a_name: { - source_b_name: distance_between(source_a, source_b, *sources.closest_point_between_sources(source_a, source_b)) + source_b_name: distance_between( + source_a, + source_b, + *sources.closest_point_between_sources(source_a, source_b), + ) for source_b_name, source_b in sources_map.items() if source_a_name != source_b_name } From 7d4fb1d9ed8e7857749f98a8538ec0c99fd497c6 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Wed, 17 Jul 2024 13:30:53 +1200 Subject: [PATCH 10/11] docstring formatting + linting checks --- source_modelling/sources.py | 200 +++++------------------------------- 1 file changed, 26 insertions(+), 174 deletions(-) diff --git a/source_modelling/sources.py b/source_modelling/sources.py index b98d1f9..850614b 100644 --- a/source_modelling/sources.py +++ b/source_modelling/sources.py @@ -76,56 +76,27 @@ def from_lat_lon_depth(point_coordinates: np.ndarray, **kwargs) -> "Point": @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. - """ + """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). - """ + """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). - """ + """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). - """ + """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). - """ + """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( @@ -145,7 +116,6 @@ def fault_coordinates_to_wgs_depth_coordinates( 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( @@ -225,105 +195,52 @@ def from_corners(corners: np.ndarray) -> "Plane": @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. - """ + """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). - """ + """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). - """ + """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). - """ + """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. - """ + """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). - """ + """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). - """ + """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). - """ + """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). - """ + """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) - """ - + """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] @@ -333,12 +250,7 @@ def strike(self) -> float: @property def dip_dir(self) -> float: - """ - Returns - ------- - float - The bearing of the dip direction (from north; in degrees). - """ + """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]) @@ -351,12 +263,7 @@ def dip_dir(self) -> float: @property def dip(self) -> float: - """ - Returns - ------- - float - The dip angle of the fault. - """ + """float: The dip angle of the fault.""" return np.degrees(np.arcsin(np.abs(self.bottom_m - self.top_m) / self.width_m)) @staticmethod @@ -369,7 +276,7 @@ def from_centroid_strike_dip( length: float, width: float, ) -> "Plane": - """Create a fault plane from the centroid, strike, dip_dir, top, bottom, length, and width + """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 @@ -378,7 +285,7 @@ def from_centroid_strike_dip( Parameters ---------- centroid : np.ndarray - The centre of the fault plane in lat, lon coordinate.s + The centre of the fault plane in lat, lon coordinates. strike : float The strike of the fault (in degrees). dip_dir : Optional[float] @@ -410,21 +317,6 @@ def from_centroid_strike_dip( ) 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: @@ -522,9 +414,7 @@ def wgs_depth_coordinates_in_plane(self, global_coordinates: np.ndarray) -> bool """ try: - plane_coordinates = self.wgs_depth_coordinates_to_fault_coordinates( - global_coordinates - ) + self.wgs_depth_coordinates_to_fault_coordinates(global_coordinates) return True except ValueError: return False @@ -592,25 +482,12 @@ def area(self) -> float: @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). - """ + """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. - """ - + """float: The total length of each fault plane.""" return self.lengths.sum() @property @@ -627,47 +504,22 @@ def width(self) -> float: @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). - """ + """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. - """ + """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. - """ + """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. - """ + """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( From f6782dc070b4074531ae2d0da6b2736cf59ff233 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Mon, 29 Jul 2024 16:29:28 +1200 Subject: [PATCH 11/11] add type signatures to d0 and delta --- source_modelling/rupture_propagation.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/source_modelling/rupture_propagation.py b/source_modelling/rupture_propagation.py index 1516797..95097c8 100644 --- a/source_modelling/rupture_propagation.py +++ b/source_modelling/rupture_propagation.py @@ -20,7 +20,7 @@ """ from collections import defaultdict, namedtuple -from typing import Optional, Tuple +from typing import Optional import numpy as np @@ -34,7 +34,7 @@ JumpPair = namedtuple("JumpPair", ["from_point", "to_point"]) -def shaw_dieterich_distance_model(distance: float, d0, delta) -> float: +def shaw_dieterich_distance_model(distance: float, d0: float, delta: float) -> float: """ Compute fault jump probabilities using the Shaw-Dieterich distance model[0].