source_modelling.fsp

  1"""Module for handling FSP files.
  2
  3This module provides classes and functions for reading FSP files,
  4as well as representing their contents.
  5See http://equake-rc.info/SRCMOD/fileformats/fsp/
  6for details on the FSP format.
  7
  8Classes
  9-------
 10- FSPFile: Representation of an FSP file.
 11
 12Exceptions
 13----------
 14- FSPParseError: Exception raised for errors in parsing FSP files.
 15
 16Example
 17-------
 18>>> fsp_file = FSPFile.read_from_file(fsp_ffp)
 19>>> (fsp_file.data['trup'] + fsp_file.data['rise']).max() # Get time of final rise for subfaults.
 20"""
 21
 22import dataclasses
 23import re
 24from pathlib import Path
 25from typing import Callable, Optional
 26
 27import pandas as pd
 28import parse
 29
 30HEADER_PATTERN = """EventTAG: {event_tag}
 31Loc : LAT = {latitude:f} LON = {longitude:f} DEP = {depth:f}
 32Size : LEN = {length:f} km WID = {width:f} km Mw = {magnitude:f} Mo = {moment:e} Nm
 33Mech : STRK = {strike:f} DIP = {dip:f} RAKE = {rake:f} Htop = {htop:f} km
 34Rupt : HypX = {hypx:f} km Hypz = {hypz:f} km avTr = {average_rise_time:f} s avVr = {average_rupture_speed:f} km/s
 35Invs : Nx = {nx:d} Nz = {nz:d} Fmin = {fmin:f} Hz Fmax = {fmax:f} Hz
 36Invs : Dx = {dx:f} km Dz = {dz:f} km
 37Invs : Ntw = {time_window_count:d} Nsg = {segment_count:d} (# of time-windows,# of fault segments)
 38Invs : LEN = {time_window_length:f} s SHF = {time_shift:f} s (time-window length and time-shift)
 39SVF : {slip_velocity_function} (type of slip-velocity function used)
 40"""
 41
 42
 43def _normalise_value(value: float) -> Optional[float]:
 44    """Normalise an FSP parameter value.
 45
 46    Parameters
 47    ----------
 48    value : float
 49        The value to normalise.
 50
 51    Returns
 52    -------
 53    Optional[float]
 54        None if value is 999 or less than zero.
 55    """
 56    return None if value == 999 or value < 0 else value
 57
 58
 59class FSPParseError(Exception):
 60    """Exception raised for errors in parsing FSP files."""
 61
 62    pass
 63
 64
 65@dataclasses.dataclass
 66class FSPFile:
 67    """A representation of the FSP file format used to represent source models in SRCMOD.
 68
 69    Attributes
 70    ----------
 71    event_tag : str
 72        An identifier for the FSPFile, e.g. "s2004IRIANx01WEIx".
 73    latitude : float
 74        Latitude of hypocentre.
 75    longitude : float
 76        Longitude of hypocentre.
 77    depth : float
 78        Depth of hypocentre (in km).
 79    hypx : float or None
 80        The x-coordinate of the hypocentre (along strike, in km).
 81    hypz : float or None
 82        The z-coordinate of the hypocentre (along dip, in km).
 83    length : float
 84        The length of the rupture (in km).
 85    width : float
 86        The width of the rupture (in km).
 87    strike : float
 88        Rupture plane strike.
 89    dip : float
 90        Rupture plane dip.
 91    rake : float
 92        Rupture plane rake.
 93    htop : float
 94        Top depth of rupture plane (in km).
 95    magnitude : float
 96        Rupture magnitude.
 97    moment : float
 98        Rupture moment (in Nm).
 99    average_rise_time : float or None
100        Average rise time of subfaults.
101    average_rupture_speed : float or None
102        Average speed of rupture front.
103    slip_velocity_function : str
104        Description of the slip velocity function used
105    nx : float or None
106        The number of subfaults along strike.
107    nz : float or None
108        The number of subfaults along dip.
109    dx : float
110        The length of each subfault.
111    dz : float
112        The width of each subfault.
113    fmin : float or None
114        The minimum frequency in the inversion.
115    fmax : float or None
116        The maximum frequency in the inversion.
117    time_window_count : float or None
118        The slip (or rake) time windows for each subfault.
119    time_window_length : float or None
120        The time window length for each subfault.
121    time_shift : float or None
122        The time shift for each time window.
123    segment_count : int
124        Number of segments (subfaults).
125    data : pd.DataFrame
126        The data frame of info on each subfault.
127    """
128
129    event_tag: str
130
131    # Hypocentre parameters
132    latitude: float
133    longitude: float
134    depth: float
135    hypx: Optional[float]
136    hypz: Optional[float]
137
138    # Fault and parameters
139    length: float
140    width: float
141    strike: float
142    dip: float
143    rake: float
144    htop: float
145
146    # Rupture parameters
147    magnitude: float
148    moment: float
149    average_rise_time: Optional[float]
150    average_rupture_speed: Optional[float]
151    slip_velocity_function: str
152
153    # Model subfault sizes
154    nx: Optional[int]
155    nz: Optional[int]
156    dx: float
157    dz: float
158
159    # Inversion parameters
160    fmin: Optional[float]
161    fmax: Optional[float]
162
163    # Time window parameters
164    time_window_count: Optional[float]
165    time_window_length: Optional[float]
166    time_shift: Optional[float]
167
168    # Number of segments (== len(data))
169    segment_count: int
170    data: pd.DataFrame
171
172    @classmethod
173    def read_from_file(cls: Callable, fsp_ffp: Path) -> "FSPFile":
174        """Parse an FSPFile.
175
176        Parameters
177        ----------
178        fsp_ffp : Path
179            Path to the FSP file.
180
181        Returns
182        -------
183        FSPFile
184            The parsed FSPFile.
185        """
186        with open(fsp_ffp, "r") as fsp_file_handle:
187            header = ""
188            # Collect and normalise all the lines that make up the header
189            # of the file.
190            for line in fsp_file_handle:
191                if line.startswith("% Data :"):
192                    break
193                if (
194                    line.startswith("% -")
195                    or line.strip() == "%"
196                    or line.startswith("% Event :")
197                ):
198                    continue
199                # Strip the leading "% ", and deduplicate the spaces in the line.
200                # This is required to normalise the string so that the parse
201                # module can handle the rest of the parsing
202                header += " ".join(line.lstrip("% ").split()) + "\n"
203
204            # Now we can parse the header according to the header pattern.
205            parse_result = parse.parse(
206                HEADER_PATTERN,
207                header,
208            )
209            if not parse_result:
210                raise FSPParseError("Failed to parse FSP file header")
211            metadata = parse_result.named
212
213            # now sniff ahead for the columns
214            for line in fsp_file_handle:
215                if re.match(r"%\s+LAT\s+LON", line):
216                    break
217            else:
218                raise FSPParseError("Cannot find columns for FSP file!")
219            columns = line.lower().lstrip("% ").split()
220
221            data = pd.read_csv(
222                fsp_file_handle,
223                delimiter=r"\s+",
224                header=None,
225                names=columns,
226                comment="%",
227            )
228            data = data.rename(columns={"x==ew": "x", "y==ns": "y"})
229            fsp_file = cls(data=data, **metadata)
230
231            # A lot of parameters can be 999 or -999 to indicate "no known
232            # value", so we can detect that and set to None
233            fsp_file.hypx = _normalise_value(fsp_file.hypx)
234            fsp_file.hypz = _normalise_value(fsp_file.hypz)
235
236            fsp_file.average_rise_time = _normalise_value(fsp_file.average_rise_time)
237            fsp_file.average_rupture_speed = _normalise_value(
238                fsp_file.average_rupture_speed
239            )
240
241            fsp_file.nx = _normalise_value(fsp_file.nx)
242            fsp_file.nz = _normalise_value(fsp_file.nz)
243
244            fsp_file.fmin = _normalise_value(fsp_file.fmin)
245            fsp_file.fmax = _normalise_value(fsp_file.fmax)
246
247            fsp_file.time_window_count = _normalise_value(fsp_file.time_window_count)
248            fsp_file.time_window_length = _normalise_value(fsp_file.time_window_length)
249            fsp_file.time_shift = _normalise_value(fsp_file.time_shift)
250            return fsp_file