source_modelling.srf
Module for handling SRF (Standard Rupture Format) files.
This module provides classes and functions for reading and writing SRF files, as well as representing their contents. See https://wiki.canterbury.ac.nz/display/QuakeCore/File+Formats+Used+On+GM for details on the SRF format.
Why Not qcore.srf?
You might use this module instead of the qcore.srf module because:
The
qcore.srfmodule does not support writing SRF files.Exposing SRF points as a pandas dataframe allows manipulation of the points using efficient vectorised operations. We use this in rupture propagation to delay ruptures by adding to the
tinitcolumn.There is better documentation for the new module than the old one.
You should use qcore.srf if you do not eventually intend to read all
points of the SRF file (it is memory efficient), or you are working
with code that already uses qcore.srf.
Classes
- SrfFile: Representation of an SRF file.
Exceptions
- SrfParseError: Exception raised for errors in parsing SRF files.
Functions
- read_srf: Read an SRF file into memory.
- write_srf: Write an SRF object to a filepath.
Example
>>> srf_file = srf.read_srf('/path/to/srf')
>>> srf_file.points['tinit'].max() # get the last time any point in the SRF ruptures
>>> srf_file.points['tinit'] += 1 # delay all points by one second
>>> coordinates.wgs_depth_to_nztm(srf_file.header[['elat', 'elon']].to_numpy())
<h1 id="get-the-coordinates-all-the-fault-plane-centres-in-the-rupture-in-nztm-format">^ get the coordinates all the fault plane centres in the rupture in NZTM format</h1>
etc...
>>> srf.write_srf('/path/to/srf', srf_file)
1"""Module for handling SRF (Standard Rupture Format) files. 2 3This module provides classes and functions for reading and writing SRF files, 4as well as representing their contents. 5See https://wiki.canterbury.ac.nz/display/QuakeCore/File+Formats+Used+On+GM 6for details on the SRF format. 7 8Why Not qcore.srf? 9------------------ 10You might use this module instead of the `qcore.srf` module because: 11 121. The `qcore.srf` module does not support writing SRF files. 13 142. Exposing SRF points as a pandas dataframe allows manipulation of 15 the points using efficient vectorised operations. We use this in 16 rupture propagation to delay ruptures by adding to the `tinit` column. 17 183. There is better documentation for the new module than the old one. 19 20You should use `qcore.srf` if you do not eventually intend to read all 21points of the SRF file (it is memory efficient), or you are working 22with code that already uses `qcore.srf`. 23 24Classes 25------- 26- SrfFile: Representation of an SRF file. 27 28Exceptions 29---------- 30- SrfParseError: Exception raised for errors in parsing SRF files. 31 32Functions 33--------- 34- read_srf: Read an SRF file into memory. 35- write_srf: Write an SRF object to a filepath. 36 37Example 38------- 39>>> srf_file = srf.read_srf('/path/to/srf') 40>>> srf_file.points['tinit'].max() # get the last time any point in the SRF ruptures 41>>> srf_file.points['tinit'] += 1 # delay all points by one second 42>>> coordinates.wgs_depth_to_nztm(srf_file.header[['elat', 'elon']].to_numpy()) 43# ^ get the coordinates all the fault plane centres in the rupture in NZTM format 44# etc... 45>>> srf.write_srf('/path/to/srf', srf_file) 46""" 47 48import dataclasses 49import functools 50import re 51from collections.abc import Sequence 52from pathlib import Path 53from typing import Optional, TextIO 54 55import numpy as np 56import pandas as pd 57import scipy as sp 58import shapely 59 60from qcore import coordinates 61from source_modelling import srf_reader 62 63PLANE_COUNT_RE = r"PLANE (\d+)" 64POINT_COUNT_RE = r"POINTS (\d+)" 65 66 67class Segments(Sequence): 68 """A read-only view for SRF segments.""" 69 70 def __init__(self, header: pd.DataFrame, points: pd.DataFrame): 71 self._header = header 72 self._points = points 73 74 def __getitem__(self, index: int) -> pd.DataFrame: 75 """Get the nth segment in the SRF. 76 77 Parameters 78 ---------- 79 index : int 80 The index of the segment. 81 82 Returns 83 ------- 84 int 85 The nth segment in the SRF. 86 """ 87 if not isinstance(index, int): 88 # NOTE: We are not covering this in test coverage because 89 # we intend to support slicing in the future. 90 raise TypeError( 91 "Segment index must an integer, not slice or tuple" 92 ) # pragma: no cover 93 points_offset = (self._header["nstk"] * self._header["ndip"]).cumsum() 94 if index == 0: 95 return self._points.iloc[: points_offset.iloc[index]] 96 return self._points.iloc[ 97 points_offset.iloc[index - 1] : points_offset.iloc[index] 98 ] 99 100 def __len__(self) -> int: 101 """int: The number of segments in the SRF.""" 102 return len(self._header) 103 104 105@dataclasses.dataclass 106class SrfFile: 107 """ 108 Representation of an SRF file. 109 110 Attributes 111 ---------- 112 version : str 113 The version of this SrfFile 114 header : pd.DataFrame 115 A list of SrfSegment objects representing the header of the SRF file. 116 The columns of the header are: 117 118 - elon: The centre longitude of the plane. 119 - elot: The centre latitude of the plane. 120 - nstk: The number of patches along strike for the plane. 121 - ndip: The number of patches along dip for the plane. 122 - len: The length of the plane (in km). 123 - wid: The width of the plane (in km). 124 - stk: The plane strike. 125 - dip: The plane dip. 126 - dtop: The top of the plane. 127 - shyp: The hypocentre location in strike coordinates. 128 - dhyp: The hypocentre location in dip coordinates. 129 130 131 points : pd.DataFrame 132 A list of SrfPoint objects representing the points in the SRF 133 file. The columns of the points dataframe are: 134 135 - lon: longitude of the patch. 136 - lat: latitude of the patch. 137 - dep: depth of the patch (in kilometres). 138 - stk: local strike. 139 - dip: local dip. 140 - area: area of the patch (in cm^2). 141 - tinit: initial rupture time for this patch (in seconds). 142 - dt: the timestep for all slipt* columns (in seconds). 143 - rake: local rake. 144 - slip1: total slip in the first component. 145 - slip2: total slip in the second component. 146 - slip3: total slip in the third component. 147 - slip: total slip. 148 149 The final two columns are computed from the SRF and are not saved to 150 disk. See the linked documentation on the SRF format for more details. 151 152 slipt{i}_array : csr_array 153 A sparse array containing the ith component of slip for each point and at each timestep, where 154 slipt{i}_array[i, j] is the slip for the ith patch at time t = j * dt. See also: SRFFile.slip. 155 156 References 157 ---------- 158 SRF File Format Doc: https://wiki.canterbury.ac.nz/display/QuakeCore/File+Formats+Used+On+GM 159 """ 160 161 version: str 162 header: pd.DataFrame 163 points: pd.DataFrame 164 slipt1_array: sp.sparse.csr_array 165 slipt2_array: sp.sparse.csr_array 166 slipt3_array: sp.sparse.csr_array 167 168 @property 169 def slip(self): 170 """csr_array: sparse array representing slip in all components.""" 171 slip_array = self.slipt1_array.power(2) 172 if self.slipt2_array: 173 slip_array += self.slipt2_array.power(2) 174 if self.slipt3_array: 175 slip_array += self.slipt3_array.power(2) 176 return slip_array.sqrt() 177 178 @property 179 def geometry(self) -> shapely.Geometry: 180 """shapely.Geometry: The shapely geometry of all segments in the SRF.""" 181 polygons = [] 182 for i, segment in enumerate(self.segments): 183 header = self.header.iloc[i] 184 nstk = header["nstk"] 185 ndip = header["ndip"] 186 corners = ( 187 segment[["lat", "lon"]] 188 .iloc[[0, nstk - 1, nstk * (ndip - 1), nstk * ndip - 1]] 189 .values 190 ) 191 if header["dip"] == 90: 192 polygons.append( 193 shapely.LineString(coordinates.wgs_depth_to_nztm(corners[:2])) 194 ) 195 else: 196 polygons.append( 197 shapely.convex_hull( 198 shapely.MultiPoint(coordinates.wgs_depth_to_nztm(corners)) 199 ) 200 ) 201 return shapely.union_all(polygons).normalize() 202 203 @property 204 def nt(self): 205 """int: The number of timeslices in the SRF.""" 206 return self.slipt1_array.shape[1] 207 208 @property 209 def dt(self): 210 """float: time resolution of SRF.""" 211 return self.points["dt"].iloc[0] 212 213 @property 214 def segments(self) -> Segments: 215 """Segments: A sequence of segments in the SRF.""" 216 return Segments(self.header, self.points) 217 218 219class SrfParseError(Exception): 220 """Exception raised for errors in parsing SRF files.""" 221 222 pass 223 224 225def read_float(srf_file: TextIO, label: Optional[str] = None) -> float: 226 """Read a float from an SRF file. 227 228 Parameters 229 ---------- 230 srf_file : TextIO 231 The SRF file to read from. 232 label : Optional[str] 233 A human friendly label for the floating point (for debugging 234 purposes), or None for no label. Defaults to None. 235 236 Raises 237 ------ 238 SrfParseError 239 If there is an error reading the float value from the SRF file. 240 241 Returns 242 ------- 243 float 244 The float read from the SRF file. 245 """ 246 while (cur := srf_file.read(1)).isspace(): 247 pass 248 float_str = cur 249 while not (cur := srf_file.read(1)).isspace(): 250 float_str += cur 251 try: 252 return float(float_str) 253 except ValueError: 254 if label: 255 raise SrfParseError(f'Expecting float ({label}), got: "{float_str}"') 256 else: 257 raise SrfParseError(f'Expecting float, got: "{float_str}"') 258 259 260def read_int(srf_file: TextIO, label: Optional[str] = None) -> int: 261 """Read a int from an SRF file. 262 263 Parameters 264 ---------- 265 srf_file : TextIO 266 The SRF file to read from. 267 268 Raises 269 ------ 270 SrfParseError 271 If there is an error reading the int value from the SRF file. 272 273 Returns 274 ------- 275 int 276 The int read from the SRF file. 277 """ 278 while (cur := srf_file.read(1)).isspace(): 279 pass 280 int_str = cur 281 while not (cur := srf_file.read(1)).isspace(): 282 int_str += cur 283 try: 284 return int(int_str) 285 except ValueError: 286 if label: 287 raise SrfParseError(f'Expecting int ({label}), got: "{int_str}"') 288 else: 289 raise SrfParseError(f'Expecting int, got: "{int_str}"') 290 291 292def read_srf(srf_ffp: Path) -> SrfFile: 293 """Read an SRF file into an SrfFile object. 294 295 Parameters 296 ---------- 297 srf_ffp : Path 298 The filepath of the SRF file. 299 300 Returns 301 ------- 302 SrfFile 303 The filepath of the SRF file. 304 """ 305 with open(srf_ffp, mode="r", encoding="utf-8") as srf_file_handle: 306 version = srf_file_handle.readline().strip() 307 308 plane_count_line = srf_file_handle.readline().strip() 309 plane_count_match = re.match(PLANE_COUNT_RE, plane_count_line) 310 if not plane_count_match: 311 raise SrfParseError( 312 f'Expecting PLANE header line, got: "{plane_count_line}"' 313 ) 314 plane_count = int(plane_count_match.group(1)) 315 segments = [] 316 317 for _ in range(plane_count): 318 segments.append( 319 { 320 "elon": read_float(srf_file_handle), 321 "elat": read_float(srf_file_handle), 322 "nstk": read_int(srf_file_handle), 323 "ndip": read_int(srf_file_handle), 324 "len": read_float(srf_file_handle), 325 "wid": read_float(srf_file_handle), 326 "stk": read_float(srf_file_handle), 327 "dip": read_float(srf_file_handle), 328 "dtop": read_float(srf_file_handle), 329 "shyp": read_float(srf_file_handle), 330 "dhyp": read_float(srf_file_handle), 331 } 332 ) 333 headers = pd.DataFrame(segments) 334 headers["nstk"] = headers["nstk"].astype(int) 335 headers["ndip"] = headers["ndip"].astype(int) 336 337 points_count_line = srf_file_handle.readline().strip() 338 points_count_match = re.match(POINT_COUNT_RE, points_count_line) 339 if not points_count_match: 340 raise SrfParseError( 341 f'Expecting POINTS header line, got: "{points_count_line}"' 342 ) 343 point_count = int(points_count_match.group(1)) 344 345 points_metadata, slipt1_array, slipt2_array, slipt3_array = ( 346 srf_reader.read_srf_points(srf_file_handle, point_count) 347 ) 348 points_df = pd.DataFrame( 349 data=points_metadata, 350 columns=[ 351 "lon", 352 "lat", 353 "dep", 354 "stk", 355 "dip", 356 "area", 357 "tinit", 358 "dt", 359 "rake", 360 "slip1", 361 "slip2", 362 "slip3", 363 ], 364 ) 365 366 points_df["slip"] = np.sqrt( 367 points_df["slip1"] ** 2 + points_df["slip2"] ** 2 + points_df["slip3"] ** 2 368 ) 369 return SrfFile( 370 version, headers, points_df, slipt1_array, slipt2_array, slipt3_array 371 ) 372 373 374def write_slip(srf_file: TextIO, slips: np.ndarray) -> None: 375 """Write out slip values to an SRF file. 376 377 Parameters 378 ---------- 379 srf_file : TextIO 380 The SRF file to write to. 381 slips : np.ndarray 382 The slip values to write. 383 """ 384 for i in range(0, len(slips), 6): 385 srf_file.write( 386 " " 387 + " ".join(f"{slips[j]:.5E}" for j in range(i, min(len(slips), i + 6))) 388 + "\n" 389 ) 390 391 392def write_srf_point(srf_file: TextIO, srf: SrfFile, point: pd.Series) -> None: 393 """Write out a single SRF point. 394 395 Parameters 396 ---------- 397 srf_file : TextIO 398 The SRF file to write to. 399 srf : SrfFile 400 The SRF file object to write. 401 point : pd.Series 402 The point to write. 403 """ 404 index = int(point["point_index"]) 405 # We need to get the raw slip values for the slip arrays to write out to the SRF. 406 # The slip values for the ith point in the SRF are stored in the ith row of the slip array. 407 # The indptr array contains the indices for each row in the slip array, so: 408 # 409 # - indptr[i] is the start index in the data array for the ith row and 410 # - indptr[i + 1] is the start index in the data array for the (i + 1)th row. 411 # 412 # Hence, data[indptr[i]:indptr[i+1]] collects the slip values for the ith row. 413 # Note we cannot use standard multi-dimensional array indexing arr[i, :] 414 # because scipy sparse arrays do not support this kind of indexing. 415 # 416 # Visually: 417 # 418 # +----+-----+----+---+---+--+ 419 # | | | | | | | indptr 420 # +----+---\-+----+---+-\-+--+ 421 # \ \ 422 # \ \ 423 # +----------+\-----------+\----------+-------------+ 424 # | ... | row i | row i + 1 | ... | data 425 # +----------+------------+-----------+-------------+ 426 # row_i = data[indptr[i]:indptr[i+1]] 427 slipt1 = ( 428 srf.slipt1_array.data[ 429 srf.slipt1_array.indptr[index] : srf.slipt1_array.indptr[index + 1] 430 ] 431 if srf.slipt1_array is not None 432 else None 433 ) 434 slipt2 = ( 435 srf.slipt2_array.data[ 436 srf.slipt2_array.indptr[index] : srf.slipt2_array.indptr[index + 1] 437 ] 438 if srf.slipt2_array is not None 439 else None 440 ) 441 slipt3 = ( 442 srf.slipt3_array.data[ 443 srf.slipt3_array.indptr[index] : srf.slipt3_array.indptr[index + 1] 444 ] 445 if srf.slipt3_array is not None 446 else None 447 ) 448 srf_file.write( 449 f"{point['lon']:.6f} {point['lat']:.6f} {point['dep']:g} {point['stk']:g} {point['dip']:g} {point['area']:.4E} {point['tinit']:.4f} {point['dt']:.6E}\n" 450 ) 451 srf_file.write( 452 f"{point['rake']:g} {point['slip1']:.4f} {len(slipt1) if slipt1 is not None else 0} {point['slip2']:.4f} {len(slipt2) if slipt2 is not None else 0} {point['slip3']:.4f} {len(slipt3) if slipt3 is not None else 0}\n" 453 ) 454 if slipt1 is not None: 455 write_slip(srf_file, slipt1) 456 if slipt2 is not None: 457 write_slip(srf_file, slipt2) 458 if slipt3 is not None: 459 write_slip(srf_file, slipt3) 460 461 462def write_srf(srf_ffp: Path, srf: SrfFile) -> None: 463 """Write an SRF object to a filepath. 464 465 Parameters 466 ---------- 467 srf_ffp : Path 468 The filepath to write the srf object to. 469 srf : SrfFile 470 The SRF object. 471 """ 472 with open(srf_ffp, mode="w", encoding="utf-8") as srf_file_handle: 473 srf_file_handle.write("1.0\n") 474 srf_file_handle.write(f"PLANE {len(srf.header)}\n") 475 # Cannot use srf.header.to_string because the newline separating headers is significant! 476 # This is ok because the number of headers is typically very small (< 100) 477 for _, plane in srf.header.iterrows(): 478 srf_file_handle.write( 479 "\n".join( 480 [ 481 f"{plane['elon']:.6f} {plane['elat']:.6f} {int(plane['nstk'])} {int(plane['ndip'])} {plane['len']:.4f} {plane['wid']:.4f}", 482 f"{plane['stk']:.4f} {plane['dip']:.4f} {plane['dtop']:.4f} {plane['shyp']:.4f} {plane['dhyp']:.4f}", 483 "", 484 ] 485 ) 486 ) 487 488 srf_file_handle.write(f"POINTS {len(srf.points)}\n") 489 srf.points["point_index"] = np.arange(len(srf.points)) 490 491 srf.points.apply( 492 functools.partial(write_srf_point, srf_file_handle, srf), axis=1 493 ) 494 srf.points = srf.points.drop(columns="point_index")
PLANE_COUNT_RE =
'PLANE (\\d+)'
POINT_COUNT_RE =
'POINTS (\\d+)'
class
Segments(collections.abc.Sequence):
68class Segments(Sequence): 69 """A read-only view for SRF segments.""" 70 71 def __init__(self, header: pd.DataFrame, points: pd.DataFrame): 72 self._header = header 73 self._points = points 74 75 def __getitem__(self, index: int) -> pd.DataFrame: 76 """Get the nth segment in the SRF. 77 78 Parameters 79 ---------- 80 index : int 81 The index of the segment. 82 83 Returns 84 ------- 85 int 86 The nth segment in the SRF. 87 """ 88 if not isinstance(index, int): 89 # NOTE: We are not covering this in test coverage because 90 # we intend to support slicing in the future. 91 raise TypeError( 92 "Segment index must an integer, not slice or tuple" 93 ) # pragma: no cover 94 points_offset = (self._header["nstk"] * self._header["ndip"]).cumsum() 95 if index == 0: 96 return self._points.iloc[: points_offset.iloc[index]] 97 return self._points.iloc[ 98 points_offset.iloc[index - 1] : points_offset.iloc[index] 99 ] 100 101 def __len__(self) -> int: 102 """int: The number of segments in the SRF.""" 103 return len(self._header)
A read-only view for SRF segments.
@dataclasses.dataclass
class
SrfFile:
106@dataclasses.dataclass 107class SrfFile: 108 """ 109 Representation of an SRF file. 110 111 Attributes 112 ---------- 113 version : str 114 The version of this SrfFile 115 header : pd.DataFrame 116 A list of SrfSegment objects representing the header of the SRF file. 117 The columns of the header are: 118 119 - elon: The centre longitude of the plane. 120 - elot: The centre latitude of the plane. 121 - nstk: The number of patches along strike for the plane. 122 - ndip: The number of patches along dip for the plane. 123 - len: The length of the plane (in km). 124 - wid: The width of the plane (in km). 125 - stk: The plane strike. 126 - dip: The plane dip. 127 - dtop: The top of the plane. 128 - shyp: The hypocentre location in strike coordinates. 129 - dhyp: The hypocentre location in dip coordinates. 130 131 132 points : pd.DataFrame 133 A list of SrfPoint objects representing the points in the SRF 134 file. The columns of the points dataframe are: 135 136 - lon: longitude of the patch. 137 - lat: latitude of the patch. 138 - dep: depth of the patch (in kilometres). 139 - stk: local strike. 140 - dip: local dip. 141 - area: area of the patch (in cm^2). 142 - tinit: initial rupture time for this patch (in seconds). 143 - dt: the timestep for all slipt* columns (in seconds). 144 - rake: local rake. 145 - slip1: total slip in the first component. 146 - slip2: total slip in the second component. 147 - slip3: total slip in the third component. 148 - slip: total slip. 149 150 The final two columns are computed from the SRF and are not saved to 151 disk. See the linked documentation on the SRF format for more details. 152 153 slipt{i}_array : csr_array 154 A sparse array containing the ith component of slip for each point and at each timestep, where 155 slipt{i}_array[i, j] is the slip for the ith patch at time t = j * dt. See also: SRFFile.slip. 156 157 References 158 ---------- 159 SRF File Format Doc: https://wiki.canterbury.ac.nz/display/QuakeCore/File+Formats+Used+On+GM 160 """ 161 162 version: str 163 header: pd.DataFrame 164 points: pd.DataFrame 165 slipt1_array: sp.sparse.csr_array 166 slipt2_array: sp.sparse.csr_array 167 slipt3_array: sp.sparse.csr_array 168 169 @property 170 def slip(self): 171 """csr_array: sparse array representing slip in all components.""" 172 slip_array = self.slipt1_array.power(2) 173 if self.slipt2_array: 174 slip_array += self.slipt2_array.power(2) 175 if self.slipt3_array: 176 slip_array += self.slipt3_array.power(2) 177 return slip_array.sqrt() 178 179 @property 180 def geometry(self) -> shapely.Geometry: 181 """shapely.Geometry: The shapely geometry of all segments in the SRF.""" 182 polygons = [] 183 for i, segment in enumerate(self.segments): 184 header = self.header.iloc[i] 185 nstk = header["nstk"] 186 ndip = header["ndip"] 187 corners = ( 188 segment[["lat", "lon"]] 189 .iloc[[0, nstk - 1, nstk * (ndip - 1), nstk * ndip - 1]] 190 .values 191 ) 192 if header["dip"] == 90: 193 polygons.append( 194 shapely.LineString(coordinates.wgs_depth_to_nztm(corners[:2])) 195 ) 196 else: 197 polygons.append( 198 shapely.convex_hull( 199 shapely.MultiPoint(coordinates.wgs_depth_to_nztm(corners)) 200 ) 201 ) 202 return shapely.union_all(polygons).normalize() 203 204 @property 205 def nt(self): 206 """int: The number of timeslices in the SRF.""" 207 return self.slipt1_array.shape[1] 208 209 @property 210 def dt(self): 211 """float: time resolution of SRF.""" 212 return self.points["dt"].iloc[0] 213 214 @property 215 def segments(self) -> Segments: 216 """Segments: A sequence of segments in the SRF.""" 217 return Segments(self.header, self.points)
Representation of an SRF file.
Attributes
- version (str): The version of this SrfFile
header (pd.DataFrame): A list of SrfSegment objects representing the header of the SRF file. The columns of the header are:
- elon: The centre longitude of the plane.
- elot: The centre latitude of the plane.
- nstk: The number of patches along strike for the plane.
- ndip: The number of patches along dip for the plane.
- len: The length of the plane (in km).
- wid: The width of the plane (in km).
- stk: The plane strike.
- dip: The plane dip.
- dtop: The top of the plane.
- shyp: The hypocentre location in strike coordinates.
- dhyp: The hypocentre location in dip coordinates.
points (pd.DataFrame): A list of SrfPoint objects representing the points in the SRF file. The columns of the points dataframe are:
- lon: longitude of the patch.
- lat: latitude of the patch.
- dep: depth of the patch (in kilometres).
- stk: local strike.
- dip: local dip.
- area: area of the patch (in cm^2).
- tinit: initial rupture time for this patch (in seconds).
- dt: the timestep for all slipt* columns (in seconds).
- rake: local rake.
- slip1: total slip in the first component.
- slip2: total slip in the second component.
- slip3: total slip in the third component.
- slip: total slip.
The final two columns are computed from the SRF and are not saved to disk. See the linked documentation on the SRF format for more details.
- slipt{i}_array (csr_array): A sparse array containing the ith component of slip for each point and at each timestep, where slipt{i}_array[i, j] is the slip for the ith patch at time t = j * dt. See also: SRFFile.slip.
References
SRF File Format Doc: https://wiki.canterbury.ac.nz/display/QuakeCore/File+Formats+Used+On+GM
SrfFile( version: str, header: pandas.core.frame.DataFrame, points: pandas.core.frame.DataFrame, slipt1_array: scipy.sparse._csr.csr_array, slipt2_array: scipy.sparse._csr.csr_array, slipt3_array: scipy.sparse._csr.csr_array)
slip
169 @property 170 def slip(self): 171 """csr_array: sparse array representing slip in all components.""" 172 slip_array = self.slipt1_array.power(2) 173 if self.slipt2_array: 174 slip_array += self.slipt2_array.power(2) 175 if self.slipt3_array: 176 slip_array += self.slipt3_array.power(2) 177 return slip_array.sqrt()
csr_array: sparse array representing slip in all components.
geometry: shapely.lib.Geometry
179 @property 180 def geometry(self) -> shapely.Geometry: 181 """shapely.Geometry: The shapely geometry of all segments in the SRF.""" 182 polygons = [] 183 for i, segment in enumerate(self.segments): 184 header = self.header.iloc[i] 185 nstk = header["nstk"] 186 ndip = header["ndip"] 187 corners = ( 188 segment[["lat", "lon"]] 189 .iloc[[0, nstk - 1, nstk * (ndip - 1), nstk * ndip - 1]] 190 .values 191 ) 192 if header["dip"] == 90: 193 polygons.append( 194 shapely.LineString(coordinates.wgs_depth_to_nztm(corners[:2])) 195 ) 196 else: 197 polygons.append( 198 shapely.convex_hull( 199 shapely.MultiPoint(coordinates.wgs_depth_to_nztm(corners)) 200 ) 201 ) 202 return shapely.union_all(polygons).normalize()
shapely.Geometry: The shapely geometry of all segments in the SRF.
nt
204 @property 205 def nt(self): 206 """int: The number of timeslices in the SRF.""" 207 return self.slipt1_array.shape[1]
int: The number of timeslices in the SRF.
class
SrfParseError(builtins.Exception):
220class SrfParseError(Exception): 221 """Exception raised for errors in parsing SRF files.""" 222 223 pass
Exception raised for errors in parsing SRF files.
def
read_float(srf_file: <class 'TextIO'>, label: Optional[str] = None) -> float:
226def read_float(srf_file: TextIO, label: Optional[str] = None) -> float: 227 """Read a float from an SRF file. 228 229 Parameters 230 ---------- 231 srf_file : TextIO 232 The SRF file to read from. 233 label : Optional[str] 234 A human friendly label for the floating point (for debugging 235 purposes), or None for no label. Defaults to None. 236 237 Raises 238 ------ 239 SrfParseError 240 If there is an error reading the float value from the SRF file. 241 242 Returns 243 ------- 244 float 245 The float read from the SRF file. 246 """ 247 while (cur := srf_file.read(1)).isspace(): 248 pass 249 float_str = cur 250 while not (cur := srf_file.read(1)).isspace(): 251 float_str += cur 252 try: 253 return float(float_str) 254 except ValueError: 255 if label: 256 raise SrfParseError(f'Expecting float ({label}), got: "{float_str}"') 257 else: 258 raise SrfParseError(f'Expecting float, got: "{float_str}"')
Read a float from an SRF file.
Parameters
- srf_file (TextIO): The SRF file to read from.
- label (Optional[str]): A human friendly label for the floating point (for debugging purposes), or None for no label. Defaults to None.
Raises
- SrfParseError: If there is an error reading the float value from the SRF file.
Returns
- float: The float read from the SRF file.
def
read_int(srf_file: <class 'TextIO'>, label: Optional[str] = None) -> int:
261def read_int(srf_file: TextIO, label: Optional[str] = None) -> int: 262 """Read a int from an SRF file. 263 264 Parameters 265 ---------- 266 srf_file : TextIO 267 The SRF file to read from. 268 269 Raises 270 ------ 271 SrfParseError 272 If there is an error reading the int value from the SRF file. 273 274 Returns 275 ------- 276 int 277 The int read from the SRF file. 278 """ 279 while (cur := srf_file.read(1)).isspace(): 280 pass 281 int_str = cur 282 while not (cur := srf_file.read(1)).isspace(): 283 int_str += cur 284 try: 285 return int(int_str) 286 except ValueError: 287 if label: 288 raise SrfParseError(f'Expecting int ({label}), got: "{int_str}"') 289 else: 290 raise SrfParseError(f'Expecting int, got: "{int_str}"')
Read a int from an SRF file.
Parameters
- srf_file (TextIO): The SRF file to read from.
Raises
- SrfParseError: If there is an error reading the int value from the SRF file.
Returns
- int: The int read from the SRF file.
293def read_srf(srf_ffp: Path) -> SrfFile: 294 """Read an SRF file into an SrfFile object. 295 296 Parameters 297 ---------- 298 srf_ffp : Path 299 The filepath of the SRF file. 300 301 Returns 302 ------- 303 SrfFile 304 The filepath of the SRF file. 305 """ 306 with open(srf_ffp, mode="r", encoding="utf-8") as srf_file_handle: 307 version = srf_file_handle.readline().strip() 308 309 plane_count_line = srf_file_handle.readline().strip() 310 plane_count_match = re.match(PLANE_COUNT_RE, plane_count_line) 311 if not plane_count_match: 312 raise SrfParseError( 313 f'Expecting PLANE header line, got: "{plane_count_line}"' 314 ) 315 plane_count = int(plane_count_match.group(1)) 316 segments = [] 317 318 for _ in range(plane_count): 319 segments.append( 320 { 321 "elon": read_float(srf_file_handle), 322 "elat": read_float(srf_file_handle), 323 "nstk": read_int(srf_file_handle), 324 "ndip": read_int(srf_file_handle), 325 "len": read_float(srf_file_handle), 326 "wid": read_float(srf_file_handle), 327 "stk": read_float(srf_file_handle), 328 "dip": read_float(srf_file_handle), 329 "dtop": read_float(srf_file_handle), 330 "shyp": read_float(srf_file_handle), 331 "dhyp": read_float(srf_file_handle), 332 } 333 ) 334 headers = pd.DataFrame(segments) 335 headers["nstk"] = headers["nstk"].astype(int) 336 headers["ndip"] = headers["ndip"].astype(int) 337 338 points_count_line = srf_file_handle.readline().strip() 339 points_count_match = re.match(POINT_COUNT_RE, points_count_line) 340 if not points_count_match: 341 raise SrfParseError( 342 f'Expecting POINTS header line, got: "{points_count_line}"' 343 ) 344 point_count = int(points_count_match.group(1)) 345 346 points_metadata, slipt1_array, slipt2_array, slipt3_array = ( 347 srf_reader.read_srf_points(srf_file_handle, point_count) 348 ) 349 points_df = pd.DataFrame( 350 data=points_metadata, 351 columns=[ 352 "lon", 353 "lat", 354 "dep", 355 "stk", 356 "dip", 357 "area", 358 "tinit", 359 "dt", 360 "rake", 361 "slip1", 362 "slip2", 363 "slip3", 364 ], 365 ) 366 367 points_df["slip"] = np.sqrt( 368 points_df["slip1"] ** 2 + points_df["slip2"] ** 2 + points_df["slip3"] ** 2 369 ) 370 return SrfFile( 371 version, headers, points_df, slipt1_array, slipt2_array, slipt3_array 372 )
Read an SRF file into an SrfFile object.
Parameters
- srf_ffp (Path): The filepath of the SRF file.
Returns
- SrfFile: The filepath of the SRF file.
def
write_slip(srf_file: <class 'TextIO'>, slips: numpy.ndarray) -> None:
375def write_slip(srf_file: TextIO, slips: np.ndarray) -> None: 376 """Write out slip values to an SRF file. 377 378 Parameters 379 ---------- 380 srf_file : TextIO 381 The SRF file to write to. 382 slips : np.ndarray 383 The slip values to write. 384 """ 385 for i in range(0, len(slips), 6): 386 srf_file.write( 387 " " 388 + " ".join(f"{slips[j]:.5E}" for j in range(i, min(len(slips), i + 6))) 389 + "\n" 390 )
Write out slip values to an SRF file.
Parameters
- srf_file (TextIO): The SRF file to write to.
- slips (np.ndarray): The slip values to write.
def
write_srf_point( srf_file: <class 'TextIO'>, srf: SrfFile, point: pandas.core.series.Series) -> None:
393def write_srf_point(srf_file: TextIO, srf: SrfFile, point: pd.Series) -> None: 394 """Write out a single SRF point. 395 396 Parameters 397 ---------- 398 srf_file : TextIO 399 The SRF file to write to. 400 srf : SrfFile 401 The SRF file object to write. 402 point : pd.Series 403 The point to write. 404 """ 405 index = int(point["point_index"]) 406 # We need to get the raw slip values for the slip arrays to write out to the SRF. 407 # The slip values for the ith point in the SRF are stored in the ith row of the slip array. 408 # The indptr array contains the indices for each row in the slip array, so: 409 # 410 # - indptr[i] is the start index in the data array for the ith row and 411 # - indptr[i + 1] is the start index in the data array for the (i + 1)th row. 412 # 413 # Hence, data[indptr[i]:indptr[i+1]] collects the slip values for the ith row. 414 # Note we cannot use standard multi-dimensional array indexing arr[i, :] 415 # because scipy sparse arrays do not support this kind of indexing. 416 # 417 # Visually: 418 # 419 # +----+-----+----+---+---+--+ 420 # | | | | | | | indptr 421 # +----+---\-+----+---+-\-+--+ 422 # \ \ 423 # \ \ 424 # +----------+\-----------+\----------+-------------+ 425 # | ... | row i | row i + 1 | ... | data 426 # +----------+------------+-----------+-------------+ 427 # row_i = data[indptr[i]:indptr[i+1]] 428 slipt1 = ( 429 srf.slipt1_array.data[ 430 srf.slipt1_array.indptr[index] : srf.slipt1_array.indptr[index + 1] 431 ] 432 if srf.slipt1_array is not None 433 else None 434 ) 435 slipt2 = ( 436 srf.slipt2_array.data[ 437 srf.slipt2_array.indptr[index] : srf.slipt2_array.indptr[index + 1] 438 ] 439 if srf.slipt2_array is not None 440 else None 441 ) 442 slipt3 = ( 443 srf.slipt3_array.data[ 444 srf.slipt3_array.indptr[index] : srf.slipt3_array.indptr[index + 1] 445 ] 446 if srf.slipt3_array is not None 447 else None 448 ) 449 srf_file.write( 450 f"{point['lon']:.6f} {point['lat']:.6f} {point['dep']:g} {point['stk']:g} {point['dip']:g} {point['area']:.4E} {point['tinit']:.4f} {point['dt']:.6E}\n" 451 ) 452 srf_file.write( 453 f"{point['rake']:g} {point['slip1']:.4f} {len(slipt1) if slipt1 is not None else 0} {point['slip2']:.4f} {len(slipt2) if slipt2 is not None else 0} {point['slip3']:.4f} {len(slipt3) if slipt3 is not None else 0}\n" 454 ) 455 if slipt1 is not None: 456 write_slip(srf_file, slipt1) 457 if slipt2 is not None: 458 write_slip(srf_file, slipt2) 459 if slipt3 is not None: 460 write_slip(srf_file, slipt3)
Write out a single SRF point.
Parameters
- srf_file (TextIO): The SRF file to write to.
- srf (SrfFile): The SRF file object to write.
- point (pd.Series): The point to write.
463def write_srf(srf_ffp: Path, srf: SrfFile) -> None: 464 """Write an SRF object to a filepath. 465 466 Parameters 467 ---------- 468 srf_ffp : Path 469 The filepath to write the srf object to. 470 srf : SrfFile 471 The SRF object. 472 """ 473 with open(srf_ffp, mode="w", encoding="utf-8") as srf_file_handle: 474 srf_file_handle.write("1.0\n") 475 srf_file_handle.write(f"PLANE {len(srf.header)}\n") 476 # Cannot use srf.header.to_string because the newline separating headers is significant! 477 # This is ok because the number of headers is typically very small (< 100) 478 for _, plane in srf.header.iterrows(): 479 srf_file_handle.write( 480 "\n".join( 481 [ 482 f"{plane['elon']:.6f} {plane['elat']:.6f} {int(plane['nstk'])} {int(plane['ndip'])} {plane['len']:.4f} {plane['wid']:.4f}", 483 f"{plane['stk']:.4f} {plane['dip']:.4f} {plane['dtop']:.4f} {plane['shyp']:.4f} {plane['dhyp']:.4f}", 484 "", 485 ] 486 ) 487 ) 488 489 srf_file_handle.write(f"POINTS {len(srf.points)}\n") 490 srf.points["point_index"] = np.arange(len(srf.points)) 491 492 srf.points.apply( 493 functools.partial(write_srf_point, srf_file_handle, srf), axis=1 494 ) 495 srf.points = srf.points.drop(columns="point_index")
Write an SRF object to a filepath.
Parameters
- srf_ffp (Path): The filepath to write the srf object to.
- srf (SrfFile): The SRF object.