From 09b4c10c6716cc28acd4f779292f4d9a6088ae93 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Wed, 17 Jul 2024 12:55:18 +1200 Subject: [PATCH 1/6] add srf module --- requirements.txt | 2 + source_modelling/srf.py | 380 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 382 insertions(+) create mode 100644 source_modelling/srf.py diff --git a/requirements.txt b/requirements.txt index e69de29..5da331c 100644 --- a/requirements.txt +++ b/requirements.txt @@ -0,0 +1,2 @@ +numpy +pandas diff --git a/source_modelling/srf.py b/source_modelling/srf.py new file mode 100644 index 0000000..a8e6051 --- /dev/null +++ b/source_modelling/srf.py @@ -0,0 +1,380 @@ +""" +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. + + +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. +""" + +import dataclasses +import re +from pathlib import Path +from typing import Optional, TextIO + +import numpy as np +import pandas as pd + +PLANE_COUNT_RE = r"PLANE (\d+)" +POINT_COUNT_RE = r"POINTS (\d+)" + + +@dataclasses.dataclass +class SrfFile: + """ + 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. + points : pd.DataFrame + A list of SrfPoint objects representing the points in the SRF file. + """ + + version: str + header: pd.DataFrame + points: pd.DataFrame + + +class SrfParseError(Exception): + """Exception raised for errors in parsing SRF files.""" + + pass + + +def read_version(srf_file: TextIO): + """Read the version value from an srf file. + + Parameters + ---------- + srf_file : TextIO + The SRF file object to read from. + """ + return float(srf_file.readline()) + + +def write_version(srf_file: TextIO): + """Write version value to an srf file. + + Parameters + ---------- + srf_file : TextIO + The SRF file object to write to. + """ + srf_file.write("1.0\n") + + +def read_srf_headers(srf_file: TextIO) -> pd.DataFrame: + """Read the header section of an SRF file. + + Parameters + ---------- + srf_file : TextIO + The SRF file object to read the headers from. + + Raises + ------ + SrfParseError + An SrfParseError is raised if an error occurs parsing the SRF file headers. + + Returns + ------- + list[SrfSegment] + A list of all the planes listed in the header of the SRF file. + """ + plane_count_line = srf_file.readline().strip() + plane_count_match = re.match(PLANE_COUNT_RE, plane_count_line) + if not plane_count_match: + raise SrfParseError(f'Expecting PLANE header line, got: "{plane_count_line}"') + plane_count = int(plane_count_match.group(1)) + segments = [] + + for _ in range(plane_count): + segments.append( + { + "elon": read_float(srf_file), + "elat": read_float(srf_file), + "nstk": read_int(srf_file), + "ndip": read_int(srf_file), + "len": read_float(srf_file), + "wid": read_float(srf_file), + "stk": read_float(srf_file), + "dip": read_float(srf_file), + "dtop": read_float(srf_file), + "shyp": read_float(srf_file), + "dhyp": read_float(srf_file), + } + ) + return pd.DataFrame(segments) + + +def read_float(srf_file: TextIO, label: Optional[str] = None) -> float: + """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. + """ + while (cur := srf_file.read(1)).isspace(): + pass + float_str = cur + while not (cur := srf_file.read(1)).isspace(): + float_str += cur + try: + return float(float_str) + except ValueError: + if label: + raise SrfParseError(f'Expecting float ({label}), got: "{float_str}"') + else: + raise SrfParseError(f'Expecting float, got: "{float_str}"') + + +def read_int(srf_file: TextIO, label: Optional[str] = None) -> int: + """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. + """ + while (cur := srf_file.read(1)).isspace(): + pass + int_str = cur + while not (cur := srf_file.read(1)).isspace(): + int_str += cur + try: + return int(int_str) + except ValueError: + if label: + raise SrfParseError(f'Expecting int ({label}), got: "{int_str}"') + else: + raise SrfParseError(f'Expecting int, got: "{int_str}"') + + +def read_points_count(srf_file: TextIO) -> int: + """Read the number of points contained in an srf file. + + Parameters + ---------- + srf_file : TextIO + The SRF file object to read from. + + Raises + ------ + SrfParseError + If the current line does not match the expected syntax for the number + of points in an SRF file (see POINT_COUNT_RE). + + Returns + ------- + int + The number of points contained in the SRF file. + """ + points_count_line = srf_file.readline().strip() + points_count_match = re.match(POINT_COUNT_RE, points_count_line) + if not points_count_match: + raise SrfParseError(f'Expecting POINTS header line, got: "{points_count_line}"') + return int(points_count_match.group(1)) + + +def read_srf_point(srf_file: TextIO) -> dict[str, int | float]: + """Read a single SRF point from a file handle. + + Parameters + ---------- + srf_file : TextIO + The file handle of the SRF file, pointing to the start of the + line. + + Returns + ------- + dict[str, int | float] + A dictionary containing the values for the point. + """ + # skip optional first spaces + row = { + "lon": read_float(srf_file, label="lon"), + "lat": read_float(srf_file, label="lat"), + "dep": read_float(srf_file, label="dep"), + "stk": read_float(srf_file, label="stk"), + "dip": read_float(srf_file, label="dip"), + "area": read_float(srf_file, label="area"), + "tinit": read_float(srf_file, label="tinit"), + "dt": read_float(srf_file, label="dt"), + "rake": read_float(srf_file, label="rake"), + "slip1": read_float(srf_file, label="slip1"), + } + nt1 = read_int(srf_file, label="nt1") + row["slip2"] = read_float(srf_file, label="slip2") + nt2 = read_int(srf_file, label="nt2") + row["slip3"] = read_float(srf_file, label="slip3") + nt3 = read_int(srf_file, label="nt3") + row["slipt1"] = np.fromiter( + (read_float(srf_file, label="slipt1") for _ in range(nt1)), float + ) + row["slipt2"] = np.fromiter( + (read_float(srf_file, label="slipt2") for _ in range(nt2)), float + ) + row["slipt3"] = np.fromiter( + (read_float(srf_file, label="slipt2") for _ in range(nt3)), float + ) + return row + + +def read_srf(srf_ffp: Path) -> SrfFile: + """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. + """ + with open(srf_ffp, mode="r", encoding="utf-8") as srf_file_handle: + version = read_version(srf_file_handle) + headers = read_srf_headers(srf_file_handle) + point_count = read_points_count(srf_file_handle) + points = pd.DataFrame( + (read_srf_point(srf_file_handle) for _ in range(point_count)) + ) + return SrfFile(version, headers, points) + + +def write_srf_header(srf_file: TextIO, header: pd.DataFrame) -> None: + """Write a list of SRF plane segments as a header. + + Parameters + ---------- + srf_file : TextIO + The SRF file to write to. + header : pd.DataFrame + The list of SrfSegments to write. + """ + srf_file.write(f"PLANE {len(header)}\n") + srf_file.write( + header.to_string( + index=False, + header=None, + formatters={ + "elon": lambda elon: f"{elon:.6f}", + "elat": lambda elat: f"{elat:.6f}", + "len": lambda len: f"{len:.4f}", + "wid": lambda wid: f"{wid:.4f}", + }, + ) + + "\n" + ) + + +def write_point_count(srf_file: TextIO, point_count: int) -> None: + """Write out a POINTS declaration for the number of points in an SRF file.. + + Parameters + ---------- + srf_file : TextIO + The SRF file to write to. + point_count : int + The number of points to write. + """ + srf_file.write(f"POINTS {point_count}\n") + + +def write_slip(srf_file: TextIO, slips: np.ndarray) -> None: + """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. + """ + for i in range(0, len(slips), 6): + srf_file.write( + " " + + " ".join(f"{slips[j]:.5E}" for j in range(i, min(len(slips), i + 6))) + + "\n" + ) + + +def write_srf_point(srf_file: TextIO, point: pd.Series) -> None: + """Write out a single SRF point. + + Parameters + ---------- + srf_file : TextIO + The SRF file to write to. + point : pd.Series + The point to write. + """ + srf_file.write( + 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" + ) + srf_file.write( + f"{point['rake']:g} {point['slip1']:.4f} {len(point['slipt1'])} {point['slip2']:.4f} {len(point['slipt2'])} {point['slip3']:.4f} {len(point['slipt3'])}\n" + ) + if len(point["slipt1"]): + write_slip(srf_file, point["slipt1"]) + if len(point["slipt2"]): + write_slip(srf_file, point["slipt2"]) + if len(point["slipt3"]): + write_slip(srf_file, point["slipt3"]) + + +def write_srf(srf_ffp: Path, srf: SrfFile) -> None: + """Write an SRF object to a filepath. + + Parameters + ---------- + srf_ffp : Path + The filepath to write the srf object to. + srf : SrfFile + The SRF object. + """ + with open(srf_ffp, mode="w", encoding="utf-8") as srf_file_handle: + write_version(srf_file_handle) + write_srf_header(srf_file_handle, srf.header) + write_point_count(srf_file_handle, len(srf.points)) + srf.points.apply(lambda point: write_srf_point(srf_file_handle, point), axis=1) From bb09fda0a36dd5656976bf30f3909d46080d7798 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Wed, 17 Jul 2024 13:04:46 +1200 Subject: [PATCH 2/6] improve module docstring --- source_modelling/srf.py | 19 ++++++++++++++++--- 1 file changed, 16 insertions(+), 3 deletions(-) diff --git a/source_modelling/srf.py b/source_modelling/srf.py index a8e6051..be225de 100644 --- a/source_modelling/srf.py +++ b/source_modelling/srf.py @@ -7,15 +7,28 @@ for details on the SRF format. -Classes: +Classes +------- - SrfFile: Representation of an SRF file. -Exceptions: +Exceptions +---------- - SrfParseError: Exception raised for errors in parsing SRF files. -Functions: +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()) +# ^ get the coordinates all the fault plane centres in the rupture in NZTM format +# etc... +>>> srf.write_srf('/path/to/srf', srf_file) """ import dataclasses From 3f147e046f5b949e87e3f9d149b6120314c6245d Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Wed, 17 Jul 2024 13:08:08 +1200 Subject: [PATCH 3/6] Add motivation to module docstring --- source_modelling/srf.py | 18 ++++++++++++++++-- 1 file changed, 16 insertions(+), 2 deletions(-) diff --git a/source_modelling/srf.py b/source_modelling/srf.py index be225de..db13d7c 100644 --- a/source_modelling/srf.py +++ b/source_modelling/srf.py @@ -1,11 +1,25 @@ -""" -Module for handling SRF (Standard Rupture Format) files. +"""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: + +1. The `qcore.srf` module does not support writing SRF files. + +2. 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 `tinit` column. + +3. 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 ------- From 0f90110bc2556bf36f9d20ea7ff2ab6e5462aed5 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Fri, 19 Jul 2024 13:40:12 +1200 Subject: [PATCH 4/6] cut down on the number of functions in the module --- source_modelling/srf.py | 82 ++++++++--------------------------------- 1 file changed, 16 insertions(+), 66 deletions(-) diff --git a/source_modelling/srf.py b/source_modelling/srf.py index db13d7c..d6a8095 100644 --- a/source_modelling/srf.py +++ b/source_modelling/srf.py @@ -46,6 +46,7 @@ """ import dataclasses +import functools import re from pathlib import Path from typing import Optional, TextIO @@ -83,28 +84,6 @@ class SrfParseError(Exception): pass -def read_version(srf_file: TextIO): - """Read the version value from an srf file. - - Parameters - ---------- - srf_file : TextIO - The SRF file object to read from. - """ - return float(srf_file.readline()) - - -def write_version(srf_file: TextIO): - """Write version value to an srf file. - - Parameters - ---------- - srf_file : TextIO - The SRF file object to write to. - """ - srf_file.write("1.0\n") - - def read_srf_headers(srf_file: TextIO) -> pd.DataFrame: """Read the header section of an SRF file. @@ -216,32 +195,6 @@ def read_int(srf_file: TextIO, label: Optional[str] = None) -> int: raise SrfParseError(f'Expecting int, got: "{int_str}"') -def read_points_count(srf_file: TextIO) -> int: - """Read the number of points contained in an srf file. - - Parameters - ---------- - srf_file : TextIO - The SRF file object to read from. - - Raises - ------ - SrfParseError - If the current line does not match the expected syntax for the number - of points in an SRF file (see POINT_COUNT_RE). - - Returns - ------- - int - The number of points contained in the SRF file. - """ - points_count_line = srf_file.readline().strip() - points_count_match = re.match(POINT_COUNT_RE, points_count_line) - if not points_count_match: - raise SrfParseError(f'Expecting POINTS header line, got: "{points_count_line}"') - return int(points_count_match.group(1)) - - def read_srf_point(srf_file: TextIO) -> dict[str, int | float]: """Read a single SRF point from a file handle. @@ -300,12 +253,22 @@ def read_srf(srf_ffp: Path) -> SrfFile: The filepath of the SRF file. """ with open(srf_ffp, mode="r", encoding="utf-8") as srf_file_handle: - version = read_version(srf_file_handle) + version = float(srf_file_handle.readline()) + headers = read_srf_headers(srf_file_handle) - point_count = read_points_count(srf_file_handle) + + points_count_line = srf_file_handle.readline().strip() + points_count_match = re.match(POINT_COUNT_RE, points_count_line) + if not points_count_match: + raise SrfParseError( + f'Expecting POINTS header line, got: "{points_count_line}"' + ) + point_count = int(points_count_match.group(1)) + points = pd.DataFrame( (read_srf_point(srf_file_handle) for _ in range(point_count)) ) + return SrfFile(version, headers, points) @@ -335,19 +298,6 @@ def write_srf_header(srf_file: TextIO, header: pd.DataFrame) -> None: ) -def write_point_count(srf_file: TextIO, point_count: int) -> None: - """Write out a POINTS declaration for the number of points in an SRF file.. - - Parameters - ---------- - srf_file : TextIO - The SRF file to write to. - point_count : int - The number of points to write. - """ - srf_file.write(f"POINTS {point_count}\n") - - def write_slip(srf_file: TextIO, slips: np.ndarray) -> None: """Write out slip values to an SRF file. @@ -401,7 +351,7 @@ def write_srf(srf_ffp: Path, srf: SrfFile) -> None: The SRF object. """ with open(srf_ffp, mode="w", encoding="utf-8") as srf_file_handle: - write_version(srf_file_handle) + srf_file_handle.write("1.0\n") write_srf_header(srf_file_handle, srf.header) - write_point_count(srf_file_handle, len(srf.points)) - srf.points.apply(lambda point: write_srf_point(srf_file_handle, point), axis=1) + srf_file_handle.write(f"POINTS {len(srf.points)}\n") + srf.points.apply(functools.partial(write_srf_point, srf_file_handle), axis=1) From c9c58cba0dabec98f0e6c11dda39d065f78d2274 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Fri, 19 Jul 2024 14:03:39 +1200 Subject: [PATCH 5/6] further reduce functions in srf module --- source_modelling/srf.py | 100 ++++++++++++++++------------------------ 1 file changed, 40 insertions(+), 60 deletions(-) diff --git a/source_modelling/srf.py b/source_modelling/srf.py index d6a8095..48b1f64 100644 --- a/source_modelling/srf.py +++ b/source_modelling/srf.py @@ -84,50 +84,6 @@ class SrfParseError(Exception): pass -def read_srf_headers(srf_file: TextIO) -> pd.DataFrame: - """Read the header section of an SRF file. - - Parameters - ---------- - srf_file : TextIO - The SRF file object to read the headers from. - - Raises - ------ - SrfParseError - An SrfParseError is raised if an error occurs parsing the SRF file headers. - - Returns - ------- - list[SrfSegment] - A list of all the planes listed in the header of the SRF file. - """ - plane_count_line = srf_file.readline().strip() - plane_count_match = re.match(PLANE_COUNT_RE, plane_count_line) - if not plane_count_match: - raise SrfParseError(f'Expecting PLANE header line, got: "{plane_count_line}"') - plane_count = int(plane_count_match.group(1)) - segments = [] - - for _ in range(plane_count): - segments.append( - { - "elon": read_float(srf_file), - "elat": read_float(srf_file), - "nstk": read_int(srf_file), - "ndip": read_int(srf_file), - "len": read_float(srf_file), - "wid": read_float(srf_file), - "stk": read_float(srf_file), - "dip": read_float(srf_file), - "dtop": read_float(srf_file), - "shyp": read_float(srf_file), - "dhyp": read_float(srf_file), - } - ) - return pd.DataFrame(segments) - - def read_float(srf_file: TextIO, label: Optional[str] = None) -> float: """Read a float from an SRF file. @@ -255,7 +211,32 @@ def read_srf(srf_ffp: Path) -> SrfFile: with open(srf_ffp, mode="r", encoding="utf-8") as srf_file_handle: version = float(srf_file_handle.readline()) - headers = read_srf_headers(srf_file_handle) + plane_count_line = srf_file_handle.readline().strip() + plane_count_match = re.match(PLANE_COUNT_RE, plane_count_line) + if not plane_count_match: + raise SrfParseError( + f'Expecting PLANE header line, got: "{plane_count_line}"' + ) + plane_count = int(plane_count_match.group(1)) + segments = [] + + for _ in range(plane_count): + segments.append( + { + "elon": read_float(srf_file_handle), + "elat": read_float(srf_file_handle), + "nstk": read_int(srf_file_handle), + "ndip": read_int(srf_file_handle), + "len": read_float(srf_file_handle), + "wid": read_float(srf_file_handle), + "stk": read_float(srf_file_handle), + "dip": read_float(srf_file_handle), + "dtop": read_float(srf_file_handle), + "shyp": read_float(srf_file_handle), + "dhyp": read_float(srf_file_handle), + } + ) + headers = pd.DataFrame(segments) points_count_line = srf_file_handle.readline().strip() points_count_match = re.match(POINT_COUNT_RE, points_count_line) @@ -282,20 +263,6 @@ def write_srf_header(srf_file: TextIO, header: pd.DataFrame) -> None: header : pd.DataFrame The list of SrfSegments to write. """ - srf_file.write(f"PLANE {len(header)}\n") - srf_file.write( - header.to_string( - index=False, - header=None, - formatters={ - "elon": lambda elon: f"{elon:.6f}", - "elat": lambda elat: f"{elat:.6f}", - "len": lambda len: f"{len:.4f}", - "wid": lambda wid: f"{wid:.4f}", - }, - ) - + "\n" - ) def write_slip(srf_file: TextIO, slips: np.ndarray) -> None: @@ -352,6 +319,19 @@ def write_srf(srf_ffp: Path, srf: SrfFile) -> None: """ with open(srf_ffp, mode="w", encoding="utf-8") as srf_file_handle: srf_file_handle.write("1.0\n") - write_srf_header(srf_file_handle, srf.header) + srf_file_handle.write(f"PLANE {len(srf.header)}\n") + srf_file_handle.write( + srf.header.to_string( + index=False, + header=None, + formatters={ + "elon": lambda elon: f"{elon:.6f}", + "elat": lambda elat: f"{elat:.6f}", + "len": lambda len: f"{len:.4f}", + "wid": lambda wid: f"{wid:.4f}", + }, + ) + + "\n" + ) srf_file_handle.write(f"POINTS {len(srf.points)}\n") srf.points.apply(functools.partial(write_srf_point, srf_file_handle), axis=1) From 2d9689872e4fa08fc67713a40ce53e31a03c6a1a Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Fri, 19 Jul 2024 14:04:35 +1200 Subject: [PATCH 6/6] remove empty function --- source_modelling/srf.py | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/source_modelling/srf.py b/source_modelling/srf.py index 48b1f64..c857b5b 100644 --- a/source_modelling/srf.py +++ b/source_modelling/srf.py @@ -253,18 +253,6 @@ def read_srf(srf_ffp: Path) -> SrfFile: return SrfFile(version, headers, points) -def write_srf_header(srf_file: TextIO, header: pd.DataFrame) -> None: - """Write a list of SRF plane segments as a header. - - Parameters - ---------- - srf_file : TextIO - The SRF file to write to. - header : pd.DataFrame - The list of SrfSegments to write. - """ - - def write_slip(srf_file: TextIO, slips: np.ndarray) -> None: """Write out slip values to an SRF file.