source_modelling.community_fault_model
This module provides classes and functions for the New Zealand Community Fault Model.
The New Zealand Community Fault Model (CFM) is a collection of fault data for New
Zealand. The CFM is stored in a GeoPackage file and can be accessed using the
functions in this module. The function most_likely_nodal_plane can be used to
find the most likely nodal plane for a given centroid and two nodal planes, by
querying nearby faults in the CFM.
Examples
>>> from pathlib import Path
>>> model = get_community_fault_model()
>>> centroid = np.array([-45.1929,166.83])
>>> nodal_plane_1 = NodalPlane(strike=20, dip=35, rake=79)
>>> nodal_plane_2 = NodalPlane(strike=213, dip=56, rake=98)
>>> most_likely_plane = most_likely_nodal_plane(model, centroid, nodal_plane_1, nodal_plane_2)
>>> print(most_likely_plane)
NodalPlane(strike=20, dip=35, rake=79)
1"""This module provides classes and functions for the New Zealand Community Fault Model. 2 3The New Zealand Community Fault Model (CFM) is a collection of fault data for New 4Zealand. The CFM is stored in a GeoPackage file and can be accessed using the 5functions in this module. The function `most_likely_nodal_plane` can be used to 6find the most likely nodal plane for a given centroid and two nodal planes, by 7querying nearby faults in the CFM. 8 9Examples 10-------- 11>>> from pathlib import Path 12>>> model = get_community_fault_model() 13>>> centroid = np.array([-45.1929,166.83]) 14>>> nodal_plane_1 = NodalPlane(strike=20, dip=35, rake=79) 15>>> nodal_plane_2 = NodalPlane(strike=213, dip=56, rake=98) 16>>> most_likely_plane = most_likely_nodal_plane(model, centroid, nodal_plane_1, nodal_plane_2) 17>>> print(most_likely_plane) 18NodalPlane(strike=20, dip=35, rake=79) 19""" 20 21import re 22from dataclasses import dataclass 23from enum import Enum, Flag, auto 24from importlib import resources 25from importlib.resources.abc import Traversable 26from pathlib import Path 27from typing import NamedTuple, Optional 28 29import fiona 30import geopandas as gpd 31import numpy as np 32import numpy.typing as npt 33import shapely 34 35import source_modelling 36from qcore import coordinates, geo 37 38 39class NodalPlane(NamedTuple): 40 """Represents a nodal plane with strike, dip, and rake values.""" 41 42 strike: float 43 """The strike of the nodal plane.""" 44 45 dip: float 46 """The dip of the nodal plane.""" 47 48 rake: float 49 """The rake of the nodal plane.""" 50 51 52class CompassDirection(Enum): 53 """Enum for compass directions with their corresponding degrees.""" 54 55 N = 0 56 """North (0 degrees).""" 57 58 NW = 315 59 """Northwest (315 degrees).""" 60 61 W = 270 62 """West (270 degrees).""" 63 64 SW = 225 65 """Southwest (225 degrees).""" 66 67 S = 180 68 """South (180 degrees).""" 69 70 SE = 135 71 """Southeast (135 degrees).""" 72 73 E = 90 74 """East (90 degrees).""" 75 76 NE = 45 77 """Northeast (45 degrees).""" 78 79 80class MovementSense(Flag): 81 """Flag for movement sense types.""" 82 83 DEXTRAL = auto() 84 """Dextral movement.""" 85 86 REVERSE = auto() 87 """Reverse movement.""" 88 89 SINISTRAL = auto() 90 """Sinistral movement.""" 91 92 NORMAL = auto() 93 """Normal movement.""" 94 95 96class NameStatus(Enum): 97 """Enum for name status.""" 98 99 PUBLISHED = auto() 100 """Published status.""" 101 102 INFORMAL = auto() 103 """Informal status.""" 104 105 106class Lineage(Enum): 107 """Enum for lineage status.""" 108 109 UNMODIFIED = auto() 110 """Unmodified status.""" 111 112 MODIFIED = auto() 113 """Modified status.""" 114 115 NEW = auto() 116 """New status.""" 117 118 119class DipMethod(Enum): 120 """Enum for dip methods.""" 121 122 ELLIS2021 = auto() 123 """Ellis 2021 method.""" 124 125 OTHER = auto() 126 """Other method.""" 127 128 DOWN_DIP_INTERSECTION = auto() 129 """Down dip intersection method.""" 130 131 132class FaultStatus(Enum): 133 """Enum for fault status.""" 134 135 ACTIVE_SEISOGENIC = auto() 136 """Active seismogenic fault.""" 137 138 ACTIVE_NON_SEISOGENIC = auto() 139 """Active non-seismogenic fault.""" 140 141 INACTIVE = auto() 142 """Inactive fault.""" 143 144 145class NeotectonicDomain(Enum): 146 """Enum for neotectonic domains with their corresponding IDs.""" 147 148 HIKURANGI_SUBDUCTION_FRONT = 10 149 """Hikurangi Subduction Front (10).""" 150 ALPINE_FAULT = 15 151 """Alpine Fault (15).""" 152 PUYSEGUR_SUBDUCTION_FRONT = 26 153 """Puysegur Subduction Front (26).""" 154 155 NORTH_ISLAND_DEXTRAL_FAULT_BELT = 5 156 """North Island Dextral Fault Belt (5).""" 157 MARLBOROUGH_FAULT_SYSTEM = 14 158 """Marlborough Fault System (14).""" 159 PUYSEGUR_RIDGE_BANK = 25 160 """Puysegur Ridge Bank (25).""" 161 162 NORTH_WAIKATO_SOUTH_AUCKLAND = 2 163 """North Waikato South Auckland (2).""" 164 WESTERN_NORTH_ISLAND = 3 165 """Western North Island (3).""" 166 HAVRE_TROUGH_TAUPO_RIFT = 4 167 """Havre Trough Taupo Rift (4).""" 168 HIKURANGI_OUTER_RISE = 9 169 """Hikurangi Outer Rise (9).""" 170 NORTH_MERNOO_FRAC_ZONE = 16 171 """North Mernoo Fracture Zone (16).""" 172 PUYSEGUR_CASWELL_HIGH_OUTER_RISE = 27 173 """Puysegur Caswell High Outer Rise (27).""" 174 175 HIKURANGI_ACC_MARGIN = 7 176 """Hikurangi Accretionary Margin (7).""" 177 HIKURANGI_ACC_MARGIN_EAST_ZONE = 8 178 """Hikurangi Accretionary Margin East Zone (8).""" 179 KAPITI_MANAWATU = 12 180 """Kapiti Manawatu (12).""" 181 NW_SOUTH_ISLAND = 13 182 """Northwest South Island (13).""" 183 NE_CANTERBURY = 17 184 """Northeast Canterbury (17).""" 185 CENTRAL_CANTERBURY = 18 186 """Central Canterbury (18).""" 187 SOUTHERN_ALPS = 19 188 """Southern Alps (19).""" 189 OTAGO = 20 190 """Otago (20).""" 191 SOUTHERN_SOUTH_ISLAND = 21 192 """Southern South Island (21).""" 193 W_FIORDLAND_MARGIN_CASWELL_HIGH = 24 194 """West Fiordland Margin Caswell High (24).""" 195 196 EAST_CAPE_BLOCK = 6 197 """East Cape Block (6).""" 198 FIORDLAND_BLOCK = 23 199 """Fiordland Block (23).""" 200 201 NW_ZEALANDIA_PLATFORM = 1 202 """Northwest Zealandia Platform (1).""" 203 HIKURANGI_PLATEAU = 11 204 """Hikurangi Plateau (11).""" 205 SE_ZEALANDIA_PLATFORM = 22 206 """Southeast Zealandia Platform (22).""" 207 TASMAN_SEA_BASIN = 28 208 """Tasman Sea Basin (28).""" 209 210 211class Range(NamedTuple): 212 """Represents a range with minimum, preferred, and maximum values.""" 213 214 min: float 215 """The minimum value of the range.""" 216 pref: float 217 """The preferred value of the range.""" 218 max: float 219 """The maximum value of the range.""" 220 221 222@dataclass 223class CommunityFault: 224 """Represents a community fault with various attributes.""" 225 226 trace: shapely.LineString 227 """The trace of the fault.""" 228 fault_status: FaultStatus 229 """The status of the fault.""" 230 name: str 231 """The name of the fault.""" 232 domain: NeotectonicDomain 233 """The neotectonic domain of the fault.""" 234 dip_range: Range 235 """The dip range of the fault.""" 236 dip_dir: Optional[CompassDirection] 237 """The dip direction of the fault.""" 238 sense: MovementSense 239 """The primary movement sense of the fault.""" 240 secondary_sense: Optional[MovementSense] 241 """The secondary movement sense of the fault.""" 242 rake_range: Range 243 """The rake range of the fault.""" 244 slip_rate: Range 245 """The slip rate of the fault.""" 246 slip_rate_timeframe_pref: int 247 """The preferred slip rate timeframe.""" 248 slip_rate_timeframe_unit: Optional[int] 249 """The unit of the slip rate timeframe.""" 250 down_dip_depth90: float 251 """Not sure!""" 252 down_dip_depth_dfc: float 253 """Not sure!""" 254 down_dip_depth90_method: DipMethod 255 """Not sure!""" 256 down_dip_dfc_method: DipMethod 257 """Not sure!""" 258 up_dip_depth: Range 259 """The up dip depth range.""" 260 quality_code: int 261 """The quality code of the fault.""" 262 reference: str 263 """The reference for the fault data.""" 264 comments: Optional[str] 265 """Additional comments about the fault.""" 266 267 268def load_community_fault_model( 269 community_fault_model_shp_ffp: Path | Traversable, 270) -> list[CommunityFault]: 271 """Load the community fault model from a shapefile. 272 273 Parameters 274 ---------- 275 community_fault_model_shp_ffp : Path | Traversable 276 The file path to the shapefile. 277 278 Returns 279 ------- 280 list[CommunityFault] 281 A list of CommunityFault objects. 282 """ 283 faults = [] 284 with fiona.open(community_fault_model_shp_ffp) as fault_model_reader: 285 fault_status_map = { 286 "A-LS": FaultStatus.ACTIVE_SEISOGENIC, 287 "A-US": FaultStatus.ACTIVE_NON_SEISOGENIC, 288 "N-PS": FaultStatus.INACTIVE, 289 } 290 291 for feature in fault_model_reader: 292 faults.append( 293 CommunityFault( 294 trace=feature_trace(feature), 295 fault_status=fault_status_map[feature["properties"]["Fault_stat"]], 296 name=feature["properties"]["Name"], 297 domain=NeotectonicDomain(feature["properties"]["Domain_No"]), 298 dip_range=Range( 299 feature["properties"]["Dip_min"], 300 feature["properties"]["Dip_pref"], 301 feature["properties"]["Dip_max"], 302 ), 303 dip_dir=CompassDirection[feature["properties"]["Dip_dir"]] 304 if feature["properties"]["Dip_dir"].lower() 305 != "subvertical and variable" 306 else None, 307 sense=fault_sense_parse(feature["properties"]["Dom_sense"]), 308 secondary_sense=fault_sense_parse( 309 feature["properties"]["Sub_sense"] 310 ) 311 if feature["properties"]["Sub_sense"] 312 else None, 313 rake_range=Range( 314 feature["properties"]["Rake_minus"], 315 feature["properties"]["Rake_pref"], 316 feature["properties"]["Rake_plus"], 317 ), 318 slip_rate=Range( 319 feature["properties"]["SR_min"], 320 feature["properties"]["SR_pref"], 321 feature["properties"]["SR_max"], 322 ), 323 slip_rate_timeframe_pref=feature["properties"]["SRT_pref"], 324 slip_rate_timeframe_unit=parse_timeframe_unit( 325 feature["properties"]["SRT_gen"] 326 ), 327 up_dip_depth=Range( 328 feature["properties"]["UpdDth_min"], 329 feature["properties"]["UpdDth_prf"], 330 feature["properties"]["UpdDth_max"], 331 ), 332 down_dip_depth90=feature["properties"]["Depth_D90"], 333 down_dip_depth90_method=DipMethod( 334 feature["properties"]["Method_D90"] 335 ), 336 down_dip_depth_dfc=feature["properties"]["Depth_Dfc"], 337 down_dip_dfc_method=DipMethod(feature["properties"]["Method_Dfc"]), 338 quality_code=feature["properties"]["QualCode"], 339 reference=feature["properties"]["References"], 340 comments=feature["properties"]["Comments"], 341 ) 342 ) 343 return faults 344 345 346def get_community_fault_model() -> list[CommunityFault]: 347 """Get the community fault model. 348 349 Returns 350 ------- 351 list[CommunityFault] 352 A list of CommunityFault objects. 353 """ 354 return load_community_fault_model( 355 resources.files(source_modelling) / "NZ_CFM" / "NZ_CFM_v1_0.shp" 356 ) 357 358 359def community_fault_model_as_geodataframe() -> gpd.GeoDataFrame: 360 """ 361 Convert the community fault model to a GeoDataFrame. 362 363 Returns 364 ------- 365 gpd.GeoDataFrame 366 A GeoDataFrame containing the community fault model data, with the 367 coordinate reference system set to EPSG:2193 and the geometry column 368 set to 'trace'. The index is set to the fault names. 369 """ 370 model = get_community_fault_model() 371 372 # Transform the trace to WGS84 373 for fault in model: 374 fault.trace = shapely.transform(fault.trace, lambda coord: coord[:, ::-1]) 375 return gpd.GeoDataFrame( 376 [vars(fault) for fault in model], geometry="trace", crs="EPSG:2193" 377 ).set_index("name") 378 379 380def parse_timeframe_unit(unit: str) -> Optional[int]: 381 """Parse the timeframe unit from a string. 382 383 Parameters 384 ---------- 385 unit : str 386 The unit string. 387 388 Returns 389 ------- 390 Optional[int] 391 The parsed unit as an integer, or None if not applicable. 392 """ 393 if not unit: 394 return None 395 if "Ma" in unit: 396 return int(1e6) 397 398 match = re.match(r"(([\d,]+)s? )?yrs", unit) 399 if match: 400 return int(match.group(2).replace(",", "")) if match.group(2) else 1 401 402 raise ValueError(f"Unknown time unit {unit}") 403 404 405def fault_sense_parse(sense: str) -> MovementSense: 406 """Parse the fault sense from a string. 407 408 Parameters 409 ---------- 410 sense : str 411 The sense string. 412 413 Returns 414 ------- 415 MovementSense 416 The parsed MovementSense. 417 """ 418 parsed_sense = MovementSense(0) 419 if "dextral" in sense: 420 parsed_sense |= MovementSense.DEXTRAL 421 if "reverse" in sense: 422 parsed_sense |= MovementSense.REVERSE 423 if "sinistral" in sense: 424 parsed_sense |= MovementSense.SINISTRAL 425 if "normal" in sense: 426 parsed_sense |= MovementSense.NORMAL 427 return parsed_sense 428 429 430def most_likely_nodal_plane( 431 faults: list[CommunityFault], 432 centroid: npt.NDArray, 433 nodal_plane_1: NodalPlane, 434 nodal_plane_2: NodalPlane, 435 k_neighbours: int = 5, 436) -> NodalPlane: 437 """Find the most likely nodal plane by a nearest neighbour vote on 438 the strike values of the faults near the centroid. 439 440 Parameters 441 ---------- 442 faults : list[CommunityFault] 443 A list of CommunityFault objects. 444 centroid : npt.NDArray 445 The centroid of the fault. 446 nodal_plane_1 : NodalPlane 447 The first nodal plane. 448 nodal_plane_2 : NodalPlane 449 The second nodal plane. 450 k_neighbours : int, optional 451 The number of nearest neighbours to consider (default is 5). 452 453 Returns 454 ------- 455 NodalPlane 456 The most likely nodal plane. 457 """ 458 point = shapely.Point(coordinates.wgs_depth_to_nztm(centroid)) 459 line_segments = { 460 shapely.LineString(fault.trace.coords[i : i + 2]): line_segment_strike( 461 *fault.trace.coords[i : i + 2] 462 ) 463 for fault in faults 464 for i in range(len(fault.trace.coords) - 1) 465 } 466 closest_segments = sorted( 467 line_segments, key=lambda segment: segment.distance(point) 468 )[:k_neighbours] 469 nodal_plane_1_votes = sum( 470 1 471 / ( 472 segment.distance(point) 473 * abs(line_segments[segment] - nodal_plane_1.strike + 1e-5) 474 ) 475 for segment in closest_segments 476 if abs(line_segments[segment] - nodal_plane_1.strike) 477 < abs(line_segments[segment] - nodal_plane_2.strike) 478 ) 479 nodal_plane_2_votes = sum( 480 1 481 / ( 482 segment.distance(point) 483 * abs(line_segments[segment] - nodal_plane_1.strike + 1e-5) 484 ) 485 for segment in closest_segments 486 if abs(line_segments[segment] - nodal_plane_1.strike) 487 >= abs(line_segments[segment] - nodal_plane_2.strike) 488 ) 489 490 if nodal_plane_1_votes >= nodal_plane_2_votes: 491 return nodal_plane_1 492 493 return nodal_plane_2 494 495 496def line_segment_strike(point_a: npt.ArrayLike, point_b: npt.ArrayLike) -> float: 497 """Calculate the strike of a line segment defined by two points. 498 499 Parameters 500 ---------- 501 point_a : npt.ArrayLike 502 The first point of the line segment. 503 point_b : npt.ArrayLike 504 The second point of the line segment. 505 506 Returns 507 ------- 508 float 509 The strike of the line segment. 510 """ 511 point_a = np.asarray(point_a) 512 point_b = np.asarray(point_b) 513 514 return float( 515 coordinates.nztm_bearing_to_great_circle_bearing( 516 coordinates.nztm_to_wgs_depth(point_a), 517 np.linalg.norm(point_b - point_a) / 1000, 518 geo.oriented_bearing_wrt_normal( 519 np.array([1, 0, 0]), 520 np.append( 521 point_b - point_a, 522 0, 523 ), 524 np.array([0, 0, 1]), 525 ), 526 ) 527 ) 528 529 530def feature_trace(feature: fiona.Feature) -> shapely.LineString: 531 """Extract the trace of a fault feature as a LineString. 532 533 Parameters 534 ---------- 535 feature : fiona.Feature 536 The feature from which to extract the trace. 537 538 Returns 539 ------- 540 shapely.LineString 541 The extracted trace as a LineString. 542 """ 543 points = np.array(feature.geometry.coordinates)[:, ::-1] 544 strike = line_segment_strike(points[0], points[1]) 545 try: 546 compass_direction = CompassDirection[feature.properties["Dip_dir"]] 547 if strike > compass_direction.value: 548 points = points[::-1] 549 except KeyError: 550 pass 551 return shapely.LineString(points)
class
NodalPlane(typing.NamedTuple):
40class NodalPlane(NamedTuple): 41 """Represents a nodal plane with strike, dip, and rake values.""" 42 43 strike: float 44 """The strike of the nodal plane.""" 45 46 dip: float 47 """The dip of the nodal plane.""" 48 49 rake: float 50 """The rake of the nodal plane."""
Represents a nodal plane with strike, dip, and rake values.
class
CompassDirection(enum.Enum):
53class CompassDirection(Enum): 54 """Enum for compass directions with their corresponding degrees.""" 55 56 N = 0 57 """North (0 degrees).""" 58 59 NW = 315 60 """Northwest (315 degrees).""" 61 62 W = 270 63 """West (270 degrees).""" 64 65 SW = 225 66 """Southwest (225 degrees).""" 67 68 S = 180 69 """South (180 degrees).""" 70 71 SE = 135 72 """Southeast (135 degrees).""" 73 74 E = 90 75 """East (90 degrees).""" 76 77 NE = 45 78 """Northeast (45 degrees)."""
Enum for compass directions with their corresponding degrees.
class
MovementSense(enum.Flag):
81class MovementSense(Flag): 82 """Flag for movement sense types.""" 83 84 DEXTRAL = auto() 85 """Dextral movement.""" 86 87 REVERSE = auto() 88 """Reverse movement.""" 89 90 SINISTRAL = auto() 91 """Sinistral movement.""" 92 93 NORMAL = auto() 94 """Normal movement."""
Flag for movement sense types.
class
NameStatus(enum.Enum):
97class NameStatus(Enum): 98 """Enum for name status.""" 99 100 PUBLISHED = auto() 101 """Published status.""" 102 103 INFORMAL = auto() 104 """Informal status."""
Enum for name status.
class
Lineage(enum.Enum):
107class Lineage(Enum): 108 """Enum for lineage status.""" 109 110 UNMODIFIED = auto() 111 """Unmodified status.""" 112 113 MODIFIED = auto() 114 """Modified status.""" 115 116 NEW = auto() 117 """New status."""
Enum for lineage status.
class
DipMethod(enum.Enum):
120class DipMethod(Enum): 121 """Enum for dip methods.""" 122 123 ELLIS2021 = auto() 124 """Ellis 2021 method.""" 125 126 OTHER = auto() 127 """Other method.""" 128 129 DOWN_DIP_INTERSECTION = auto() 130 """Down dip intersection method."""
Enum for dip methods.
class
FaultStatus(enum.Enum):
133class FaultStatus(Enum): 134 """Enum for fault status.""" 135 136 ACTIVE_SEISOGENIC = auto() 137 """Active seismogenic fault.""" 138 139 ACTIVE_NON_SEISOGENIC = auto() 140 """Active non-seismogenic fault.""" 141 142 INACTIVE = auto() 143 """Inactive fault."""
Enum for fault status.
class
NeotectonicDomain(enum.Enum):
146class NeotectonicDomain(Enum): 147 """Enum for neotectonic domains with their corresponding IDs.""" 148 149 HIKURANGI_SUBDUCTION_FRONT = 10 150 """Hikurangi Subduction Front (10).""" 151 ALPINE_FAULT = 15 152 """Alpine Fault (15).""" 153 PUYSEGUR_SUBDUCTION_FRONT = 26 154 """Puysegur Subduction Front (26).""" 155 156 NORTH_ISLAND_DEXTRAL_FAULT_BELT = 5 157 """North Island Dextral Fault Belt (5).""" 158 MARLBOROUGH_FAULT_SYSTEM = 14 159 """Marlborough Fault System (14).""" 160 PUYSEGUR_RIDGE_BANK = 25 161 """Puysegur Ridge Bank (25).""" 162 163 NORTH_WAIKATO_SOUTH_AUCKLAND = 2 164 """North Waikato South Auckland (2).""" 165 WESTERN_NORTH_ISLAND = 3 166 """Western North Island (3).""" 167 HAVRE_TROUGH_TAUPO_RIFT = 4 168 """Havre Trough Taupo Rift (4).""" 169 HIKURANGI_OUTER_RISE = 9 170 """Hikurangi Outer Rise (9).""" 171 NORTH_MERNOO_FRAC_ZONE = 16 172 """North Mernoo Fracture Zone (16).""" 173 PUYSEGUR_CASWELL_HIGH_OUTER_RISE = 27 174 """Puysegur Caswell High Outer Rise (27).""" 175 176 HIKURANGI_ACC_MARGIN = 7 177 """Hikurangi Accretionary Margin (7).""" 178 HIKURANGI_ACC_MARGIN_EAST_ZONE = 8 179 """Hikurangi Accretionary Margin East Zone (8).""" 180 KAPITI_MANAWATU = 12 181 """Kapiti Manawatu (12).""" 182 NW_SOUTH_ISLAND = 13 183 """Northwest South Island (13).""" 184 NE_CANTERBURY = 17 185 """Northeast Canterbury (17).""" 186 CENTRAL_CANTERBURY = 18 187 """Central Canterbury (18).""" 188 SOUTHERN_ALPS = 19 189 """Southern Alps (19).""" 190 OTAGO = 20 191 """Otago (20).""" 192 SOUTHERN_SOUTH_ISLAND = 21 193 """Southern South Island (21).""" 194 W_FIORDLAND_MARGIN_CASWELL_HIGH = 24 195 """West Fiordland Margin Caswell High (24).""" 196 197 EAST_CAPE_BLOCK = 6 198 """East Cape Block (6).""" 199 FIORDLAND_BLOCK = 23 200 """Fiordland Block (23).""" 201 202 NW_ZEALANDIA_PLATFORM = 1 203 """Northwest Zealandia Platform (1).""" 204 HIKURANGI_PLATEAU = 11 205 """Hikurangi Plateau (11).""" 206 SE_ZEALANDIA_PLATFORM = 22 207 """Southeast Zealandia Platform (22).""" 208 TASMAN_SEA_BASIN = 28 209 """Tasman Sea Basin (28)."""
Enum for neotectonic domains with their corresponding IDs.
HIKURANGI_SUBDUCTION_FRONT =
<NeotectonicDomain.HIKURANGI_SUBDUCTION_FRONT: 10>
Hikurangi Subduction Front (10).
PUYSEGUR_SUBDUCTION_FRONT =
<NeotectonicDomain.PUYSEGUR_SUBDUCTION_FRONT: 26>
Puysegur Subduction Front (26).
NORTH_ISLAND_DEXTRAL_FAULT_BELT =
<NeotectonicDomain.NORTH_ISLAND_DEXTRAL_FAULT_BELT: 5>
North Island Dextral Fault Belt (5).
MARLBOROUGH_FAULT_SYSTEM =
<NeotectonicDomain.MARLBOROUGH_FAULT_SYSTEM: 14>
Marlborough Fault System (14).
NORTH_WAIKATO_SOUTH_AUCKLAND =
<NeotectonicDomain.NORTH_WAIKATO_SOUTH_AUCKLAND: 2>
North Waikato South Auckland (2).
HAVRE_TROUGH_TAUPO_RIFT =
<NeotectonicDomain.HAVRE_TROUGH_TAUPO_RIFT: 4>
Havre Trough Taupo Rift (4).
NORTH_MERNOO_FRAC_ZONE =
<NeotectonicDomain.NORTH_MERNOO_FRAC_ZONE: 16>
North Mernoo Fracture Zone (16).
PUYSEGUR_CASWELL_HIGH_OUTER_RISE =
<NeotectonicDomain.PUYSEGUR_CASWELL_HIGH_OUTER_RISE: 27>
Puysegur Caswell High Outer Rise (27).
HIKURANGI_ACC_MARGIN =
<NeotectonicDomain.HIKURANGI_ACC_MARGIN: 7>
Hikurangi Accretionary Margin (7).
HIKURANGI_ACC_MARGIN_EAST_ZONE =
<NeotectonicDomain.HIKURANGI_ACC_MARGIN_EAST_ZONE: 8>
Hikurangi Accretionary Margin East Zone (8).
W_FIORDLAND_MARGIN_CASWELL_HIGH =
<NeotectonicDomain.W_FIORDLAND_MARGIN_CASWELL_HIGH: 24>
West Fiordland Margin Caswell High (24).
NW_ZEALANDIA_PLATFORM =
<NeotectonicDomain.NW_ZEALANDIA_PLATFORM: 1>
Northwest Zealandia Platform (1).
SE_ZEALANDIA_PLATFORM =
<NeotectonicDomain.SE_ZEALANDIA_PLATFORM: 22>
Southeast Zealandia Platform (22).
class
Range(typing.NamedTuple):
212class Range(NamedTuple): 213 """Represents a range with minimum, preferred, and maximum values.""" 214 215 min: float 216 """The minimum value of the range.""" 217 pref: float 218 """The preferred value of the range.""" 219 max: float 220 """The maximum value of the range."""
Represents a range with minimum, preferred, and maximum values.
@dataclass
class
CommunityFault:
223@dataclass 224class CommunityFault: 225 """Represents a community fault with various attributes.""" 226 227 trace: shapely.LineString 228 """The trace of the fault.""" 229 fault_status: FaultStatus 230 """The status of the fault.""" 231 name: str 232 """The name of the fault.""" 233 domain: NeotectonicDomain 234 """The neotectonic domain of the fault.""" 235 dip_range: Range 236 """The dip range of the fault.""" 237 dip_dir: Optional[CompassDirection] 238 """The dip direction of the fault.""" 239 sense: MovementSense 240 """The primary movement sense of the fault.""" 241 secondary_sense: Optional[MovementSense] 242 """The secondary movement sense of the fault.""" 243 rake_range: Range 244 """The rake range of the fault.""" 245 slip_rate: Range 246 """The slip rate of the fault.""" 247 slip_rate_timeframe_pref: int 248 """The preferred slip rate timeframe.""" 249 slip_rate_timeframe_unit: Optional[int] 250 """The unit of the slip rate timeframe.""" 251 down_dip_depth90: float 252 """Not sure!""" 253 down_dip_depth_dfc: float 254 """Not sure!""" 255 down_dip_depth90_method: DipMethod 256 """Not sure!""" 257 down_dip_dfc_method: DipMethod 258 """Not sure!""" 259 up_dip_depth: Range 260 """The up dip depth range.""" 261 quality_code: int 262 """The quality code of the fault.""" 263 reference: str 264 """The reference for the fault data.""" 265 comments: Optional[str] 266 """Additional comments about the fault."""
Represents a community fault with various attributes.
CommunityFault( trace: shapely.geometry.linestring.LineString, fault_status: FaultStatus, name: str, domain: NeotectonicDomain, dip_range: Range, dip_dir: Optional[CompassDirection], sense: MovementSense, secondary_sense: Optional[MovementSense], rake_range: Range, slip_rate: Range, slip_rate_timeframe_pref: int, slip_rate_timeframe_unit: Optional[int], down_dip_depth90: float, down_dip_depth_dfc: float, down_dip_depth90_method: DipMethod, down_dip_dfc_method: DipMethod, up_dip_depth: Range, quality_code: int, reference: str, comments: Optional[str])
def
load_community_fault_model( community_fault_model_shp_ffp: pathlib.Path | importlib.resources.abc.Traversable) -> list[CommunityFault]:
269def load_community_fault_model( 270 community_fault_model_shp_ffp: Path | Traversable, 271) -> list[CommunityFault]: 272 """Load the community fault model from a shapefile. 273 274 Parameters 275 ---------- 276 community_fault_model_shp_ffp : Path | Traversable 277 The file path to the shapefile. 278 279 Returns 280 ------- 281 list[CommunityFault] 282 A list of CommunityFault objects. 283 """ 284 faults = [] 285 with fiona.open(community_fault_model_shp_ffp) as fault_model_reader: 286 fault_status_map = { 287 "A-LS": FaultStatus.ACTIVE_SEISOGENIC, 288 "A-US": FaultStatus.ACTIVE_NON_SEISOGENIC, 289 "N-PS": FaultStatus.INACTIVE, 290 } 291 292 for feature in fault_model_reader: 293 faults.append( 294 CommunityFault( 295 trace=feature_trace(feature), 296 fault_status=fault_status_map[feature["properties"]["Fault_stat"]], 297 name=feature["properties"]["Name"], 298 domain=NeotectonicDomain(feature["properties"]["Domain_No"]), 299 dip_range=Range( 300 feature["properties"]["Dip_min"], 301 feature["properties"]["Dip_pref"], 302 feature["properties"]["Dip_max"], 303 ), 304 dip_dir=CompassDirection[feature["properties"]["Dip_dir"]] 305 if feature["properties"]["Dip_dir"].lower() 306 != "subvertical and variable" 307 else None, 308 sense=fault_sense_parse(feature["properties"]["Dom_sense"]), 309 secondary_sense=fault_sense_parse( 310 feature["properties"]["Sub_sense"] 311 ) 312 if feature["properties"]["Sub_sense"] 313 else None, 314 rake_range=Range( 315 feature["properties"]["Rake_minus"], 316 feature["properties"]["Rake_pref"], 317 feature["properties"]["Rake_plus"], 318 ), 319 slip_rate=Range( 320 feature["properties"]["SR_min"], 321 feature["properties"]["SR_pref"], 322 feature["properties"]["SR_max"], 323 ), 324 slip_rate_timeframe_pref=feature["properties"]["SRT_pref"], 325 slip_rate_timeframe_unit=parse_timeframe_unit( 326 feature["properties"]["SRT_gen"] 327 ), 328 up_dip_depth=Range( 329 feature["properties"]["UpdDth_min"], 330 feature["properties"]["UpdDth_prf"], 331 feature["properties"]["UpdDth_max"], 332 ), 333 down_dip_depth90=feature["properties"]["Depth_D90"], 334 down_dip_depth90_method=DipMethod( 335 feature["properties"]["Method_D90"] 336 ), 337 down_dip_depth_dfc=feature["properties"]["Depth_Dfc"], 338 down_dip_dfc_method=DipMethod(feature["properties"]["Method_Dfc"]), 339 quality_code=feature["properties"]["QualCode"], 340 reference=feature["properties"]["References"], 341 comments=feature["properties"]["Comments"], 342 ) 343 ) 344 return faults
Load the community fault model from a shapefile.
Parameters
- community_fault_model_shp_ffp (Path | Traversable): The file path to the shapefile.
Returns
- list[CommunityFault]: A list of CommunityFault objects.
347def get_community_fault_model() -> list[CommunityFault]: 348 """Get the community fault model. 349 350 Returns 351 ------- 352 list[CommunityFault] 353 A list of CommunityFault objects. 354 """ 355 return load_community_fault_model( 356 resources.files(source_modelling) / "NZ_CFM" / "NZ_CFM_v1_0.shp" 357 )
Get the community fault model.
Returns
- list[CommunityFault]: A list of CommunityFault objects.
def
community_fault_model_as_geodataframe() -> geopandas.geodataframe.GeoDataFrame:
360def community_fault_model_as_geodataframe() -> gpd.GeoDataFrame: 361 """ 362 Convert the community fault model to a GeoDataFrame. 363 364 Returns 365 ------- 366 gpd.GeoDataFrame 367 A GeoDataFrame containing the community fault model data, with the 368 coordinate reference system set to EPSG:2193 and the geometry column 369 set to 'trace'. The index is set to the fault names. 370 """ 371 model = get_community_fault_model() 372 373 # Transform the trace to WGS84 374 for fault in model: 375 fault.trace = shapely.transform(fault.trace, lambda coord: coord[:, ::-1]) 376 return gpd.GeoDataFrame( 377 [vars(fault) for fault in model], geometry="trace", crs="EPSG:2193" 378 ).set_index("name")
Convert the community fault model to a GeoDataFrame.
Returns
- gpd.GeoDataFrame: A GeoDataFrame containing the community fault model data, with the coordinate reference system set to EPSG:2193 and the geometry column set to 'trace'. The index is set to the fault names.
def
parse_timeframe_unit(unit: str) -> Optional[int]:
381def parse_timeframe_unit(unit: str) -> Optional[int]: 382 """Parse the timeframe unit from a string. 383 384 Parameters 385 ---------- 386 unit : str 387 The unit string. 388 389 Returns 390 ------- 391 Optional[int] 392 The parsed unit as an integer, or None if not applicable. 393 """ 394 if not unit: 395 return None 396 if "Ma" in unit: 397 return int(1e6) 398 399 match = re.match(r"(([\d,]+)s? )?yrs", unit) 400 if match: 401 return int(match.group(2).replace(",", "")) if match.group(2) else 1 402 403 raise ValueError(f"Unknown time unit {unit}")
Parse the timeframe unit from a string.
Parameters
- unit (str): The unit string.
Returns
- Optional[int]: The parsed unit as an integer, or None if not applicable.
406def fault_sense_parse(sense: str) -> MovementSense: 407 """Parse the fault sense from a string. 408 409 Parameters 410 ---------- 411 sense : str 412 The sense string. 413 414 Returns 415 ------- 416 MovementSense 417 The parsed MovementSense. 418 """ 419 parsed_sense = MovementSense(0) 420 if "dextral" in sense: 421 parsed_sense |= MovementSense.DEXTRAL 422 if "reverse" in sense: 423 parsed_sense |= MovementSense.REVERSE 424 if "sinistral" in sense: 425 parsed_sense |= MovementSense.SINISTRAL 426 if "normal" in sense: 427 parsed_sense |= MovementSense.NORMAL 428 return parsed_sense
Parse the fault sense from a string.
Parameters
- sense (str): The sense string.
Returns
- MovementSense: The parsed MovementSense.
def
most_likely_nodal_plane( faults: list[CommunityFault], centroid: numpy.ndarray[typing.Any, numpy.dtype[+_ScalarType_co]], nodal_plane_1: NodalPlane, nodal_plane_2: NodalPlane, k_neighbours: int = 5) -> NodalPlane:
431def most_likely_nodal_plane( 432 faults: list[CommunityFault], 433 centroid: npt.NDArray, 434 nodal_plane_1: NodalPlane, 435 nodal_plane_2: NodalPlane, 436 k_neighbours: int = 5, 437) -> NodalPlane: 438 """Find the most likely nodal plane by a nearest neighbour vote on 439 the strike values of the faults near the centroid. 440 441 Parameters 442 ---------- 443 faults : list[CommunityFault] 444 A list of CommunityFault objects. 445 centroid : npt.NDArray 446 The centroid of the fault. 447 nodal_plane_1 : NodalPlane 448 The first nodal plane. 449 nodal_plane_2 : NodalPlane 450 The second nodal plane. 451 k_neighbours : int, optional 452 The number of nearest neighbours to consider (default is 5). 453 454 Returns 455 ------- 456 NodalPlane 457 The most likely nodal plane. 458 """ 459 point = shapely.Point(coordinates.wgs_depth_to_nztm(centroid)) 460 line_segments = { 461 shapely.LineString(fault.trace.coords[i : i + 2]): line_segment_strike( 462 *fault.trace.coords[i : i + 2] 463 ) 464 for fault in faults 465 for i in range(len(fault.trace.coords) - 1) 466 } 467 closest_segments = sorted( 468 line_segments, key=lambda segment: segment.distance(point) 469 )[:k_neighbours] 470 nodal_plane_1_votes = sum( 471 1 472 / ( 473 segment.distance(point) 474 * abs(line_segments[segment] - nodal_plane_1.strike + 1e-5) 475 ) 476 for segment in closest_segments 477 if abs(line_segments[segment] - nodal_plane_1.strike) 478 < abs(line_segments[segment] - nodal_plane_2.strike) 479 ) 480 nodal_plane_2_votes = sum( 481 1 482 / ( 483 segment.distance(point) 484 * abs(line_segments[segment] - nodal_plane_1.strike + 1e-5) 485 ) 486 for segment in closest_segments 487 if abs(line_segments[segment] - nodal_plane_1.strike) 488 >= abs(line_segments[segment] - nodal_plane_2.strike) 489 ) 490 491 if nodal_plane_1_votes >= nodal_plane_2_votes: 492 return nodal_plane_1 493 494 return nodal_plane_2
Find the most likely nodal plane by a nearest neighbour vote on the strike values of the faults near the centroid.
Parameters
- faults (list[CommunityFault]): A list of CommunityFault objects.
- centroid (npt.NDArray): The centroid of the fault.
- nodal_plane_1 (NodalPlane): The first nodal plane.
- nodal_plane_2 (NodalPlane): The second nodal plane.
- k_neighbours (int, optional): The number of nearest neighbours to consider (default is 5).
Returns
- NodalPlane: The most likely nodal plane.
def
line_segment_strike( point_a: Union[Buffer, numpy._typing._array_like._SupportsArray[numpy.dtype[Any]], numpy._typing._nested_sequence._NestedSequence[numpy._typing._array_like._SupportsArray[numpy.dtype[Any]]], bool, int, float, complex, str, bytes, numpy._typing._nested_sequence._NestedSequence[Union[bool, int, float, complex, str, bytes]]], point_b: Union[Buffer, numpy._typing._array_like._SupportsArray[numpy.dtype[Any]], numpy._typing._nested_sequence._NestedSequence[numpy._typing._array_like._SupportsArray[numpy.dtype[Any]]], bool, int, float, complex, str, bytes, numpy._typing._nested_sequence._NestedSequence[Union[bool, int, float, complex, str, bytes]]]) -> float:
497def line_segment_strike(point_a: npt.ArrayLike, point_b: npt.ArrayLike) -> float: 498 """Calculate the strike of a line segment defined by two points. 499 500 Parameters 501 ---------- 502 point_a : npt.ArrayLike 503 The first point of the line segment. 504 point_b : npt.ArrayLike 505 The second point of the line segment. 506 507 Returns 508 ------- 509 float 510 The strike of the line segment. 511 """ 512 point_a = np.asarray(point_a) 513 point_b = np.asarray(point_b) 514 515 return float( 516 coordinates.nztm_bearing_to_great_circle_bearing( 517 coordinates.nztm_to_wgs_depth(point_a), 518 np.linalg.norm(point_b - point_a) / 1000, 519 geo.oriented_bearing_wrt_normal( 520 np.array([1, 0, 0]), 521 np.append( 522 point_b - point_a, 523 0, 524 ), 525 np.array([0, 0, 1]), 526 ), 527 ) 528 )
Calculate the strike of a line segment defined by two points.
Parameters
- point_a (npt.ArrayLike): The first point of the line segment.
- point_b (npt.ArrayLike): The second point of the line segment.
Returns
- float: The strike of the line segment.
def
feature_trace(feature: fiona.model.Feature) -> shapely.geometry.linestring.LineString:
531def feature_trace(feature: fiona.Feature) -> shapely.LineString: 532 """Extract the trace of a fault feature as a LineString. 533 534 Parameters 535 ---------- 536 feature : fiona.Feature 537 The feature from which to extract the trace. 538 539 Returns 540 ------- 541 shapely.LineString 542 The extracted trace as a LineString. 543 """ 544 points = np.array(feature.geometry.coordinates)[:, ::-1] 545 strike = line_segment_strike(points[0], points[1]) 546 try: 547 compass_direction = CompassDirection[feature.properties["Dip_dir"]] 548 if strike > compass_direction.value: 549 points = points[::-1] 550 except KeyError: 551 pass 552 return shapely.LineString(points)
Extract the trace of a fault feature as a LineString.
Parameters
- feature (fiona.Feature): The feature from which to extract the trace.
Returns
- shapely.LineString: The extracted trace as a LineString.