From 596a07b99682f0a7b32ce84ec88aaeb9d7bee9cf Mon Sep 17 00:00:00 2001 From: Denis Anikiev Date: Fri, 11 Sep 2026 02:26:09 +0300 Subject: [PATCH 01/14] Update citation guidelines in README to reflect DOI changes for versions 1.0.0 and 1.0.1 --- README.md | 44 +++++++++++++++++++++++++++++++++----------- 1 file changed, 33 insertions(+), 11 deletions(-) diff --git a/README.md b/README.md index 8fafe1e..3682aaf 100644 --- a/README.md +++ b/README.md @@ -198,31 +198,53 @@ All notable changes are documented in [CHANGELOG.md](CHANGELOG.md), following ## Citing -If you use PyNetDesign in your research, please cite the release you used. The DOI badge -at the top of this page always resolves to the most recently archived release, and the -Zenodo record of each release also carries its own version DOI. +If you use PyNetDesign in your research, please cite it as follows: -Version 1.0.0 was deposited manually, before the automated GitHub–Zenodo integration was -enabled, so it has a DOI of its own: +**Anikiev, D. (2026).** *PyNetDesign: a Python framework for microseismic monitoring network design.* Zenodo. + -**Anikiev, D. (2025).** *PyNetDesign: a Python framework for microseismic monitoring network design* (Version v1.0.0) [Computer software]. Zenodo. - +In BibTeX format: + +```bibtex +@software{Anikiev2026PyNetDesign, + author = {Anikiev, Denis}, + title = {{PyNetDesign}: a {P}ython framework for microseismic monitoring network design}, + year = {2026}, + publisher = {Zenodo}, + doi = {10.5281/zenodo.22695001}, + url = {https://github.com/danikiev/PyNetDesign} +} +``` + +The DOI above represents all versions and will always resolve to the latest one. +To cite a particular version of PyNetDesign, please use the following format, e.g. for +version 1.0.1: + +**Anikiev, D. (2026).** *PyNetDesign: a Python framework for microseismic monitoring network design (1.0.1).* Zenodo. + In BibTeX format: ```bibtex -@software{Anikiev2025PyNetDesign, +@software{Anikiev2026PyNetDesignVersion, author = {Anikiev, Denis}, title = {{PyNetDesign}: a {P}ython framework for microseismic monitoring network design}, - year = {2025}, + year = {2026}, publisher = {Zenodo}, - version = {v1.0.0}, + version = {1.0.1}, license = {MIT}, - doi = {10.5281/zenodo.14945917}, + doi = {10.5281/zenodo.22695002}, url = {https://github.com/danikiev/PyNetDesign} } ``` +### Earlier release + +Version 1.0.0 was deposited manually, before the automated GitHub–Zenodo integration was +enabled, and therefore belongs to a separate Zenodo lineage. It is cross-linked to the +current record and remains citable on its own: + + **Note:** magnitudes computed with version 1.0.0 are affected by an error in the seismic moment calculation that was corrected in version 1.0.1. If you are citing or reusing results obtained with 1.0.0, please see [CHANGELOG.md](CHANGELOG.md) and issue From d270bab08506d0aca3ced52222bfc3b73998ef53 Mon Sep 17 00:00:00 2001 From: Denis Anikiev Date: Fri, 11 Sep 2026 04:17:59 +0300 Subject: [PATCH 02/14] Add SV/SH wave modes, polarization projection and per-station free surface Ports the medium-independent parts of the reference implementation that the homogeneous public package was missing, so that published DAS sensitivity results are reproducible here. Verified against the reference on the manuscript homogeneous case (1 km vertical DAS with 1 m channel spacing, and the 51-geophone comparison array): P, S and PS agree bit-exactly over all 4860 grid points, with identical NaN masks, for both fs_mode='off' and the default fs_mode='auto'. Previously S differed by 0.077 magnitude units and the geophone array by up to 0.152. Wave modes and radiation patterns: - wave_mode accepts 'P', 'SV', 'SH', 'S' and 'PS'; 'S' is composite and keeps the more detectable of the SV and SH branches through np.fmin - each phase carries its root-mean-square radiation-pattern magnitude over the focal sphere of a double couple: sqrt(4/15) for P, sqrt(7/30) for SV and sqrt(1/6) for SH. The value 0.63 used previously is the combined S-wave magnitude, since + = 2/5, and does not apply to SV alone - rad_pattern_p, rad_pattern_s and rad_patterns override the defaults, so earlier conventions stay reproducible; rad_pattern_s=0.63 recovers 1.0.x Receiver projection: - get_phase_station_directionality projects the arriving polarization vectors onto the local cable tangent, replacing the 1/|cos| and 1/|sin| shortcut that was only valid for a vertical cable - stations with Components == 'Z' project onto the vertical axis; 3C stations record the full vector and are unaffected - get_ray_station_directionality is kept as a deprecated wrapper. Its rays_p and rays_s arguments are removed because nothing in the homogeneous package could ever populate them Free surface: - compute_free_surface_coefficients applies the factor of two per receiver, gated on depth or on an optional Surface column of weights, instead of amplifying every receiver in the geometry - DetectionParameters gains fs_mode, fs_level and fs_deviation; the deprecated free_surface flag maps onto fs_mode and validation stays idempotent Also fixes the station-count guard in mag_detectable, which compared against the grid axis rather than the station axis. Adds 105 tests across five modules covering wave-mode expansion, the radiation constants and their combination to sqrt(2/5), polarization projection against closed forms for a vertical cable, SH being unobservable on a vertical cable, the free-surface modes and weights, network detectability, and the end-to-end phase behaviour. All ten examples still run. --- CHANGELOG.md | 27 + pynetdesign/modelling/classes.py | 85 +++- pynetdesign/modelling/core.py | 17 +- pynetdesign/modelling/homo.py | 177 ++++--- pynetdesign/modelling/utils.py | 765 ++++++++++++++++++++++------- pytests/test_detectability.py | 107 ++++ pytests/test_directionality.py | 192 ++++++++ pytests/test_free_surface.py | 125 +++++ pytests/test_sensitivity_phases.py | 172 +++++++ pytests/test_wave_modes.py | 125 +++++ 10 files changed, 1537 insertions(+), 255 deletions(-) create mode 100644 pytests/test_detectability.py create mode 100644 pytests/test_directionality.py create mode 100644 pytests/test_free_surface.py create mode 100644 pytests/test_sensitivity_phases.py create mode 100644 pytests/test_wave_modes.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 0e82d54..830eef3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,33 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +### Added + +- phase-resolved wave modes: `wave_mode` accepts `'P'`, `'SV'`, `'SH'`, `'S'` and `'PS'`, where `'S'` is a composite mode keeping the more detectable of the SV and SH branches +- root-mean-square radiation-pattern magnitudes per phase, `R_P = sqrt(4/15) ~ 0.52`, `R_SV = sqrt(7/30) ~ 0.48` and `R_SH = sqrt(1/6) ~ 0.41`, averaged over the focal sphere of a double-couple source after Boore and Boatwright (1984) and Hallo and Eisner (2013), exposed through `phase_radiation_pattern` +- `rad_pattern_p`, `rad_pattern_s` and `rad_patterns` overrides on `mag_sensitivity_grid`, so that results computed with another convention stay reproducible; pass `rad_pattern_s=0.63` to use the combined S-wave value +- receiver projection from the arriving polarization vectors, through `get_phase_station_directionality`, which is correct for inclined and curved cables rather than only for a vertical one +- projection of single vertical-component stations onto the vertical axis, driven by a `Components` column; three-component stations record the full vector and are unaffected +- per-station free surface amplification through `compute_free_surface_coefficients`, with `fs_mode`, `fs_level` and `fs_deviation` parameters and an optional `Surface` column carrying per-station weights +- helpers `normalize_wave_mode`, `wave_mode_phases`, `wave_mode_has_p`, `wave_mode_has_s`, `wave_mode_s_phases`, `resolve_radiation_pattern`, `check_min_stations` and `geometry_requires_receiver_projection` +- test modules covering wave modes and radiation patterns, receiver projection, the free surface correction, network detectability, and phase-resolved sensitivity end to end + +### Changed + +- **S-wave results differ from 1.0.x.** The S branch now uses the phase-resolved radiation patterns instead of the combined value 0.63, which raises S-wave thresholds by `2/3*log10(0.63/sqrt(7/30))`, about 0.077 magnitude units. Pass `rad_pattern_s=0.63` to recover the previous behaviour +- **the free surface correction is applied per receiver.** In 1.0.x `free_surface=True` amplified every receiver irrespective of its depth; it now applies only to receivers at the free surface, so geometries mixing surface and downhole receivers are treated correctly. The default mode is `'auto'`, which derives the indicator from the receiver depths or from the `Surface` column when present +- `DetectionParameters` takes `fs_mode`, `fs_level` and `fs_deviation`; the `free_surface` flag is deprecated but still accepted, mapping to `fs_mode='on'` or `'off'` +- `mag_detectable` accepts the phase-resolved wave modes +- receiver projection is decided from the geometry itself, so `use_station_directionality` is only needed to force tangent projection for a station geometry + +### Removed + +- the `rays_p` and `rays_s` arguments of `get_ray_station_directionality`, which could never be populated because the homogeneous package has no ray tracer. The function is deprecated in favour of `get_phase_station_directionality` and now emits a `DeprecationWarning` + +### Fixed + +- the station-count guard in `mag_detectable` compared against the grid axis instead of the station axis, so a request for more receivers than the geometry has was not reported clearly + ## [v1.0.1] - 2026-09-10 Bugfix release correcting the seismic moment calculation. diff --git a/pynetdesign/modelling/classes.py b/pynetdesign/modelling/classes.py index 3d68509..aa6db1b 100644 --- a/pynetdesign/modelling/classes.py +++ b/pynetdesign/modelling/classes.py @@ -1,6 +1,55 @@ import numpy as np +# Detection modes for the free surface correction +VALID_FS_MODES = ('auto', 'on', 'off') + + class DetectionParameters: + r""" + Container for the parameters controlling event detection. + + Parameters + ---------- + min_SNR_p : :obj:`float` + Minimum signal-to-noise ratio required to detect a P wave on one receiver. + min_SNR_s : :obj:`float` + Minimum signal-to-noise ratio required to detect an S wave on one receiver. + min_stations_p : :obj:`int` + Minimum number of receivers on which the event must be detected with P waves. + min_stations_s : :obj:`int` + Minimum number of receivers on which the event must be detected with S waves. + f_p : :obj:`float`, optional + Representative frequency of the P wave (Hz). + If ``None``, the peak frequency is computed from the model. + f_s : :obj:`float`, optional + Representative frequency of the S wave (Hz). + If ``None``, the peak frequency is computed from the model. + f_p_corner : :obj:`float`, optional, default: ``None`` (100 Hz) + Corner frequency of the P wave (Hz), limiting the peak frequency. + f_s_corner : :obj:`float`, optional, default: ``None`` (100 Hz) + Corner frequency of the S wave (Hz), limiting the peak frequency. + fs_mode : :obj:`str`, optional, default: ``None`` (``'auto'``) + Free surface correction mode, one of ``'auto'``, ``'on'`` or ``'off'``. + In ``'auto'`` mode the correction is applied to receivers that are within + ``fs_deviation`` of the free surface level, or according to the ``Surface`` + column of the geometry when present. + fs_level : :obj:`float`, optional, default: ``None`` (0.0) + Level of the free surface in metres above the reference level. + fs_deviation : :obj:`float`, optional, default: ``None`` (1.0) + Permitted deviation from the free surface level in metres. + free_surface : :obj:`bool`, optional + Deprecated override kept for backward compatibility. ``True`` maps to + ``fs_mode='on'`` and ``False`` to ``fs_mode='off'``. If ``None``, ``fs_mode`` + is used. + + Notes + ----- + The free surface correction is applied per receiver, so a geometry combining + surface and downhole receivers is handled correctly. This differs from + PyNetDesign 1.0.x, where ``free_surface=True`` amplified every receiver + irrespective of its depth. + """ + def __init__(self, min_SNR_p: float, min_SNR_s: float, @@ -10,7 +59,10 @@ def __init__(self, f_s: float=None, f_p_corner: float=None, f_s_corner: float=None, - free_surface: bool=False): + fs_mode: str=None, + fs_level: float=None, + fs_deviation: float=None, + free_surface: bool=None): self.min_SNR_p = min_SNR_p self.min_SNR_s = min_SNR_s @@ -20,6 +72,9 @@ def __init__(self, self.f_s = f_s self.f_p_corner = f_p_corner self.f_s_corner = f_s_corner + self.fs_mode = fs_mode + self.fs_level = fs_level + self.fs_deviation = fs_deviation self.free_surface = free_surface def validate(self): @@ -55,3 +110,31 @@ def validate(self): raise ValueError("S-wave corner frequency must be positive") else: self.f_s_corner = 100 + + # The deprecated free_surface flag, when given, decides the mode. + # Written back to fs_mode so that repeated validation is idempotent. + if self.free_surface is not None: + mapped = 'on' if self.free_surface else 'off' + if self.fs_mode is not None and self.fs_mode != mapped: + raise ValueError("Provide either free_surface or fs_mode, not both") + self.fs_mode = mapped + + if self.fs_mode is None: + self.fs_mode = 'auto' + elif not isinstance(self.fs_mode, str): + raise ValueError("Free surface correction mode must be a string") + elif self.fs_mode not in VALID_FS_MODES: + raise ValueError( + "Free surface correction mode must be one of the following: " + + ", ".join(VALID_FS_MODES) + ) + + if self.fs_level is None: + self.fs_level = 0.0 + elif np.ndim(self.fs_level) != 0: + raise ValueError("Free surface level must be a scalar") + + if self.fs_deviation is None: + self.fs_deviation = 1.0 + elif self.fs_deviation < 0: + raise ValueError("Free surface level deviation must be non-negative") diff --git a/pynetdesign/modelling/core.py b/pynetdesign/modelling/core.py index f39fffa..4676d50 100644 --- a/pynetdesign/modelling/core.py +++ b/pynetdesign/modelling/core.py @@ -30,7 +30,7 @@ def get_mag_sensitivity_local(grid_coords: np.ndarray, params : :obj:`pynetdesign.modelling.classes.DetectionParameters` Detection parameters wave_mode: :obj:`str`, optional, default: 'PS' - Wave mode to use, can be 'P','S' or 'PS' + Wave mode to use, can be 'P', 'SV', 'SH', 'S' or 'PS' strict_nan_check: :obj:`bool`, optional, default: False flag to check if NaN values appear in the result @@ -50,11 +50,9 @@ def get_mag_sensitivity_local(grid_coords: np.ndarray, # Get gauge length if exists gauge_length = geometry_df.attrs.get('gauge_length', None) - # Decide on station directionality - if gauge_length is None: - use_station_directionality = False - else: - use_station_directionality = True + # Receiver projection is decided from the geometry itself: a DAS geometry is + # projected onto the local cable tangent and a single-component station onto + # the vertical axis, so no explicit flag is needed here. # Get minimum detectable amplitudes min_amps_p, min_amps_s = retrieve_min_amps(geometry_df=geometry_df, @@ -86,8 +84,9 @@ def get_mag_sensitivity_local(grid_coords: np.ndarray, f_p_corner=params.f_p_corner, f_s_corner=params.f_s_corner, wave_mode=wave_mode, - use_free_surface=params.free_surface, - use_station_directionality=use_station_directionality, + fs_mode=params.fs_mode, + fs_level=params.fs_level, + fs_deviation=params.fs_deviation, return_stations=True, strict_nan_check=strict_nan_check) @@ -120,7 +119,7 @@ def get_mag_sensitivity(grid_coords: np.ndarray, params : :obj:`pynetdesign.modelling.classes.DetectionParameters` Detection parameters wave_mode: :obj:`str`, optional, default: 'PS' - Wave mode to use, can be 'P','S' or 'PS' + Wave mode to use, can be 'P', 'SV', 'SH', 'S' or 'PS' strict_nan_check: :obj:`bool`, optional, default: False flag to check if NaN values appear in the result diff --git a/pynetdesign/modelling/homo.py b/pynetdesign/modelling/homo.py index a1a377c..4687cc7 100644 --- a/pynetdesign/modelling/homo.py +++ b/pynetdesign/modelling/homo.py @@ -15,14 +15,17 @@ def mag_sensitivity_grid(grid_coords: np.ndarray, min_amps_s: np.ndarray = None, min_stations_p: int = None, min_stations_s: int = None, - rad_pattern_p: float = 0.52, - rad_pattern_s: float = 0.63, + rad_pattern_p: float = None, + rad_pattern_s: float = None, + rad_patterns: dict = None, f_p: float = None, f_s: float = None, f_p_corner: float = None, f_s_corner: float = None, wave_mode: str = 'PS', - use_free_surface: bool = False, + fs_mode: str = 'auto', + fs_level: float = None, + fs_deviation: float = None, use_station_directionality: bool = False, return_stations: bool = False, strict_nan_check: bool = False, @@ -53,10 +56,17 @@ def mag_sensitivity_grid(grid_coords: np.ndarray, Minimum number of stations on which event must be detected with P waves min_stations_s : :obj:`int`, optional, default: None (3) Minimum number of stations on which event must be detected with S waves - rad_pattern_p : :obj:`float`, optional, default: 0.52 - Radiation pattern factor for P waves - rad_pattern_s : :obj:`float`, optional, default: 0.63 - Radiation pattern factor for S waves + rad_pattern_p : :obj:`float`, optional + Override of the radiation pattern factor for P waves. + If ``None``, the root-mean-square value :math:`\sqrt{4/15} \approx 0.52` is used. + rad_pattern_s : :obj:`float`, optional + Override of the radiation pattern factor applied to both S phases. + If ``None``, the root-mean-square values :math:`\sqrt{7/30} \approx 0.48` for SV + and :math:`\sqrt{1/6} \approx 0.41` for SH are used. Pass 0.63 to reproduce + results computed with the combined S-wave value. + rad_patterns : :obj:`dict`, optional + Per-phase overrides of the radiation pattern factor, e.g. ``{'SV': 0.63}``. + Takes precedence over ``rad_pattern_p`` and ``rad_pattern_s``. f_p : :obj:`float`, optional Frequency of the P wave (Hz) f_s : :obj:`float`, optional @@ -66,14 +76,22 @@ def mag_sensitivity_grid(grid_coords: np.ndarray, f_s_corner : :obj:`float`, optional Corner frequency of the S wave (Hz) wave_mode: :obj:`str`, optional, default: 'PS' - Wave mode to use, can be 'P','S' or 'PS' - use_free_surface : :obj:`bool`, optional, default: False - Use free surface correction. - If True, this will consider amplification of the amplitudes when recording on the daily surface - due to the free surface boundary condition. + Wave mode to use, can be 'P', 'SV', 'SH', 'S' or 'PS' + fs_mode : :obj:`str`, optional, default: 'auto' + Free surface correction mode, can be 'auto', 'on' or 'off'. + The correction accounts for the amplification of amplitudes recorded at the + daily surface due to the free surface boundary condition, and is applied per + station so that geometries mixing surface and downhole receivers are handled + correctly. See :func:`~pynetdesign.modelling.utils.compute_free_surface_coefficients`. + fs_level : :obj:`float`, optional, default: None (0.0) + Level of the free surface in metres above the reference level. + fs_deviation : :obj:`float`, optional, default: None (1.0) + Permitted deviation from the free surface level in metres. use_station_directionality : :obj:`bool`, optional, default: False - Use station directionality correction. - This will consider sensitivity of the measurement only in direction parallel to the adjacent stations + Force projection of the arriving polarization onto the local tangent defined by + the adjacent stations. Projection is applied automatically for DAS geometries + and for stations recording a single vertical component, so this is only needed + to force tangent projection for a station geometry. return_stations : :obj:`bool`, optional, default: False Return magnitude sensitivity for each station, don't take into account detectability strict_nan_check: :obj:`bool`, optional, default: False @@ -119,6 +137,12 @@ def mag_sensitivity_grid(grid_coords: np.ndarray, if station_coords is None: station_coords = geometry_df[['X', 'Y', 'Z']].to_numpy() + # Get the free surface information from geometry, if provided + if geometry_df is not None and 'Surface' in geometry_df.columns: + station_coords_on_surface = geometry_df['Surface'].to_numpy() + else: + station_coords_on_surface = None + # Handle the case if station_coords has only one point (3 elements) if station_coords.ndim == 1 and station_coords.size == 3: station_coords = station_coords[np.newaxis, :] # reshape to (1, 3) @@ -138,6 +162,21 @@ def mag_sensitivity_grid(grid_coords: np.ndarray, if len(velocity_df) != 1: raise ValueError("Velocity model is not homogeneous!") + # Expand the wave mode into the phases that have to be evaluated + wave_mode = normalize_wave_mode(wave_mode) + phases = wave_mode_phases(wave_mode) + s_phases = wave_mode_s_phases(wave_mode) + + # Check the free surface parameters + if fs_level is None: + fs_level = 0.0 + if fs_deviation is None: + fs_deviation = 1.0 + elif fs_deviation < 0: + raise ValueError("Free surface level deviation must be non-negative") + if fs_mode not in ('auto', 'on', 'off'): + raise ValueError("Free surface mode must be 'auto', 'on' or 'off'") + # Retrieve medium parameters density = velocity_df.at[0,'Rho'] v_p = velocity_df.at[0,'Vp'] @@ -153,11 +192,16 @@ def mag_sensitivity_grid(grid_coords: np.ndarray, Q_p=Q_p, min_amps_p=min_amps_p, min_amps_s=min_amps_s, - rad_pattern_p=rad_pattern_p, - rad_pattern_s=rad_pattern_s, wave_mode=wave_mode, station_coords=station_coords) + # Check the minimum number of stations required for detection + if not return_stations: + check_min_stations(min_stations_p=min_stations_p, + min_stations_s=min_stations_s, + n_stations=len(station_coords), + wave_mode=wave_mode) + # Precompute distances between all grid points and all stations # Shape: (n_grid_points, n_stations) distances = np.linalg.norm(grid_coords[:, np.newaxis] - station_coords, axis=2) @@ -165,33 +209,39 @@ def mag_sensitivity_grid(grid_coords: np.ndarray, if np.any(distances == 0): raise ValueError("Zero distance: source and receiver must have different coordinates") - # Calculate free surface correction coefficient if needed - if use_free_surface: - if v_p is None or v_s is None: - raise ValueError("Both P- and S-wave velocities must be declared for free surface correction!") - # Simple approximation - fs_coef_p = fs_coef_s = 2*np.ones_like(distances) + # Calculate free surface correction coefficients, applied per station + fs_coef_p, fs_coef_s = compute_free_surface_coefficients(station_coords=station_coords, + grid_coords=grid_coords, + station_coords_on_surface=station_coords_on_surface, + fs_mode=fs_mode, + fs_level=fs_level, + fs_deviation=fs_deviation) + + # Calculate receiver projection coefficients if needed. Projection applies + # automatically to DAS geometries and to single-component stations. + geometry_projection_required = geometry_requires_receiver_projection(geometry_df) + if use_station_directionality or geometry_projection_required: + directionality = get_phase_station_directionality(station_coords=station_coords, + grid_coords=grid_coords, + phases=phases, + geometry_df=geometry_df, + force_tangent=use_station_directionality and not geometry_projection_required, + strict_nan_check=strict_nan_check) else: - # If not using free surface correction, use ones (no effect on calculations) - fs_coef_p = fs_coef_s = np.ones_like(distances) - - # Calculate station directionality if needed - if use_station_directionality: - stations_dir_coef_p, stations_dir_coef_s = get_ray_station_directionality(station_coords=station_coords, - grid_coords=grid_coords, - strict_nan_check=strict_nan_check) - else: - # If not using directionality, use ones (no effect on calculations) - stations_dir_coef_p = stations_dir_coef_s = np.ones_like(distances) + # If not using projection, use ones (no effect on calculations) + directionality = {phase: np.ones_like(distances) for phase in phases} - # Helper function to calculate Mw_min for a given wave type - def calculate_Mw_min(f, fcorner, density, v, Q, rad_pattern, min_amps, fs_coef, stations_dir_coef, amps_type): + # Helper function to calculate Mw_min for a given phase + def calculate_Mw_min(f, fcorner, v, Q, phase, min_amps, fs_coef): M0 = calculate_M0(density=density, v=v, Q=Q, - rad_pattern=rad_pattern, + rad_pattern=resolve_radiation_pattern(phase, + rad_pattern_p=rad_pattern_p, + rad_pattern_s=rad_pattern_s, + rad_patterns=rad_patterns), r=distances, - amps = min_amps * stations_dir_coef / fs_coef, + amps = min_amps * directionality[phase] / fs_coef, f=f, fcorner=fcorner, amps_type=amps_type) @@ -202,30 +252,27 @@ def calculate_Mw_min(f, fcorner, density, v, Q, rad_pattern, min_amps, fs_coef, Mw_min_stations = np.full((2, *distances.shape), np.nan) # Calculate Mw_min for P waves if applicable - if wave_mode in ('P', 'PS'): + if wave_mode_has_p(wave_mode): Mw_min_stations[0] = calculate_Mw_min(f=f_p, fcorner=f_p_corner, - density=density, v=v_p, Q=Q_p, - rad_pattern=rad_pattern_p, + phase='P', min_amps=min_amps_p, - fs_coef=fs_coef_p, - stations_dir_coef=stations_dir_coef_p, - amps_type=amps_type) - - # Calculate Mw_min for S waves if applicable - if wave_mode in ('S', 'PS'): - Mw_min_stations[1] = calculate_Mw_min(f=f_s, - fcorner=f_s_corner, - density=density, - v=v_s, - Q=Q_s, - rad_pattern=rad_pattern_s, - min_amps=min_amps_s, - fs_coef=fs_coef_s, - stations_dir_coef=stations_dir_coef_s, - amps_type=amps_type) + fs_coef=fs_coef_p) + + # Calculate Mw_min for S waves if applicable. For the composite S mode the + # more detectable of the SV and SH branches is kept, i.e. the smaller magnitude. + if s_phases: + mw_s = [calculate_Mw_min(f=f_s, + fcorner=f_s_corner, + v=v_s, + Q=Q_s, + phase=phase, + min_amps=min_amps_s, + fs_coef=fs_coef_s) + for phase in s_phases] + Mw_min_stations[1] = mw_s[0] if len(mw_s) == 1 else np.fmin(mw_s[0], mw_s[1]) # Return depending on the request if return_stations: @@ -412,7 +459,7 @@ def calculate_M0(density: float, return M0 -def validate_parameters(density, v_p, v_s, Q_p, Q_s, min_amps_p, min_amps_s, rad_pattern_p, rad_pattern_s, wave_mode, station_coords): +def validate_parameters(density, v_p, v_s, Q_p, Q_s, min_amps_p, min_amps_s, wave_mode, station_coords): r""" Validates the input parameters. @@ -432,12 +479,8 @@ def validate_parameters(density, v_p, v_s, Q_p, Q_s, min_amps_p, min_amps_s, rad Minimal measurable displacement amplitudes of on stations [a1, a2, ...] for P waves min_amps_s : :obj:`numpy.ndarray` Minimal measurable displacement amplitudes on stations [a1, a2, ...] for S waves - rad_pattern_p : :obj:`float` - Radiation pattern factor for P waves - rad_pattern_s : :obj:`float` - Radiation pattern factor for S waves wave_mode: :obj:`str` - Wave mode to use, can be 'P','S' or 'PS' + Wave mode to use, can be 'P', 'SV', 'SH', 'S' or 'PS' station_coords : :obj:`numpy.ndarray` Array of station coordinates [[xr1, yr1, zr1], [xr2, yr2, zr2], ...] @@ -445,8 +488,6 @@ def validate_parameters(density, v_p, v_s, Q_p, Q_s, min_amps_p, min_amps_s, rad ------ ValueError : if density, wave velocity or attenuation factor is not positive - ValueError : - if radiation pattern factor for P or S wave is less or equal to 0 or higher than 1 ValueError : if number of minimal measurable amplitudes does not match number of stations ValueError : @@ -462,7 +503,7 @@ def validate_parameters(density, v_p, v_s, Q_p, Q_s, min_amps_p, min_amps_s, rad raise ValueError("Both P wave and S wave attenuation factors can not be None.") # Validate P wave parameters if applicable - if wave_mode in ('P', 'PS'): + if wave_mode_has_p(wave_mode): if any(x <= 0 for x in (v_p, Q_p)): raise ValueError("P wave parameters must be positive.") if min_amps_p is None: @@ -472,11 +513,9 @@ def validate_parameters(density, v_p, v_s, Q_p, Q_s, min_amps_p, min_amps_s, rad raise ValueError("Invalid minimal measurable amplitudes for P waves.") elif len(min_amps_p) != len(station_coords): raise ValueError("Invalid minimal measurable amplitudes for P waves.") - if not 0 < rad_pattern_p <= 1: - raise ValueError("Invalid radiation pattern factor for P waves.") # Validate S wave parameters if applicable - if wave_mode in ('S', 'PS'): + if wave_mode_has_s(wave_mode): if any(x <= 0 for x in (v_s, Q_s)): raise ValueError("S wave parameters must be positive.") if min_amps_s is None: @@ -486,5 +525,3 @@ def validate_parameters(density, v_p, v_s, Q_p, Q_s, min_amps_p, min_amps_s, rad raise ValueError("Invalid minimal measurable amplitudes for S waves.") elif len(min_amps_s) != len(station_coords): raise ValueError("Invalid minimal measurable amplitudes for S waves.") - if not 0 < rad_pattern_s <= 1: - raise ValueError("Invalid radiation pattern factor for S waves.") \ No newline at end of file diff --git a/pynetdesign/modelling/utils.py b/pynetdesign/modelling/utils.py index 04703c4..347b9c1 100644 --- a/pynetdesign/modelling/utils.py +++ b/pynetdesign/modelling/utils.py @@ -6,6 +6,185 @@ from typing import Union, Optional import warnings +VALID_WAVE_MODES = ("P", "SV", "SH", "S", "PS") +S_WAVE_PHASES = ("SV", "SH") + +# Root-mean-square radiation-pattern magnitudes over the focal sphere of a +# double-couple source, after Boore and Boatwright (1984) and Hallo and Eisner +# (2013). Note that the frequently quoted value sqrt(2/5) ~ 0.63 is the combined +# S-wave magnitude, since + = 2/5, and therefore does not apply +# to SV or SH individually. +RAD_PATTERN_P_DEFAULT = float(np.sqrt(4.0 / 15.0)) # approx 0.516 +RAD_PATTERN_SV_DEFAULT = float(np.sqrt(7.0 / 30.0)) # approx 0.483 +RAD_PATTERN_SH_DEFAULT = float(np.sqrt(1.0 / 6.0)) # approx 0.408 + + +def normalize_wave_mode(wave_mode: str) -> str: + r""" + Return the canonical form of a wave mode. + + ``'S'`` is a composite mode representing the better recorded of the two + isotropic S polarizations, SV and SH. + + Parameters + ---------- + wave_mode : :obj:`str` + Wave mode, one of ``'P'``, ``'SV'``, ``'SH'``, ``'S'`` or ``'PS'``, + in any letter case and with surrounding whitespace ignored. + + Returns + ------- + mode : :obj:`str` + Canonical upper-case wave mode. + + Raises + ------ + ValueError + If ``wave_mode`` is not one of the supported modes. + """ + mode = str(wave_mode).strip().upper() + if mode not in VALID_WAVE_MODES: + raise ValueError( + f"Invalid wave_mode '{wave_mode}'. Must be one of {VALID_WAVE_MODES}." + ) + return mode + + +def wave_mode_phases(wave_mode: str) -> tuple: + r""" + Expand a wave mode into the physical phases that must be evaluated. + + Parameters + ---------- + wave_mode : :obj:`str` + Wave mode, see :func:`normalize_wave_mode`. + + Returns + ------- + phases : :obj:`tuple` of :obj:`str` + Phases to evaluate, drawn from ``'P'``, ``'SV'`` and ``'SH'``. + """ + mode = normalize_wave_mode(wave_mode) + if mode == "P": + return ("P",) + if mode in S_WAVE_PHASES: + return (mode,) + if mode == "S": + return S_WAVE_PHASES + return ("P", *S_WAVE_PHASES) + + +def wave_mode_has_p(wave_mode: str) -> bool: + r"""Return whether ``wave_mode`` includes the P branch.""" + return "P" in wave_mode_phases(wave_mode) + + +def wave_mode_s_phases(wave_mode: str) -> tuple: + r"""Return the S-family phases included in ``wave_mode``.""" + return tuple(phase for phase in wave_mode_phases(wave_mode) if phase in S_WAVE_PHASES) + + +def wave_mode_has_s(wave_mode: str) -> bool: + r"""Return whether ``wave_mode`` includes any S-family branch.""" + return bool(wave_mode_s_phases(wave_mode)) + + +def phase_radiation_pattern(phase: str) -> float: + r""" + Return the root-mean-square radiation-pattern magnitude of a phase. + + The values are averages over the focal sphere of a double-couple source + :cite:p:`BooreBoatwright1984`, as used for microseismicity with random source + orientations :cite:p:`HalloEisner2013`: + + .. math:: + R_P = \sqrt{4/15} \approx 0.52, \quad + R_{SV} = \sqrt{7/30} \approx 0.48, \quad + R_{SH} = \sqrt{1/6} \approx 0.41. + + Parameters + ---------- + phase : :obj:`str` + Phase, one of ``'P'``, ``'SV'`` or ``'SH'``. + + Returns + ------- + rad_pattern : :obj:`float` + Radiation-pattern magnitude of the phase. + + Raises + ------ + ValueError + If ``phase`` is a composite mode rather than a single phase. + + Notes + ----- + The value :math:`\sqrt{2/5} \approx 0.63` often quoted for S-waves is the + magnitude of the *combined* S-wave amplitude, since + :math:`\langle R_{SV}^2 \rangle + \langle R_{SH}^2 \rangle = 2/5`, and is + therefore not applicable to SV or SH individually. + """ + phase = normalize_wave_mode(phase) + if phase == "P": + return RAD_PATTERN_P_DEFAULT + if phase == "SV": + return RAD_PATTERN_SV_DEFAULT + if phase == "SH": + return RAD_PATTERN_SH_DEFAULT + raise ValueError("Radiation pattern is only defined for P, SV, and SH.") + + +def resolve_radiation_pattern(phase: str, + rad_pattern_p: float = None, + rad_pattern_s: float = None, + rad_patterns: dict = None) -> float: + r""" + Resolve the radiation-pattern magnitude to use for a phase. + + Overrides take precedence over the root-mean-square defaults of + :func:`phase_radiation_pattern`, which makes it possible to reproduce results + published with other conventions, for instance the combined S-wave value + :math:`R_S = 0.63`. + + Parameters + ---------- + phase : :obj:`str` + Phase, one of ``'P'``, ``'SV'`` or ``'SH'``. + rad_pattern_p : :obj:`float`, optional + Override for the P phase. If ``None``, the default is used. + rad_pattern_s : :obj:`float`, optional + Override applied to both S phases. If ``None``, the defaults are used. + rad_patterns : :obj:`dict`, optional + Per-phase overrides, e.g. ``{'SV': 0.63}``. Takes precedence over + ``rad_pattern_p`` and ``rad_pattern_s``. + + Returns + ------- + rad_pattern : :obj:`float` + Radiation-pattern magnitude to use. + + Raises + ------ + ValueError + If a resolved value is not in the interval (0, 1]. + """ + phase = normalize_wave_mode(phase) + value = None + if rad_patterns: + normalized = {normalize_wave_mode(key): val for key, val in rad_patterns.items()} + value = normalized.get(phase) + if value is None: + if phase == 'P': + value = rad_pattern_p + else: + value = rad_pattern_s + if value is None: + return phase_radiation_pattern(phase) + value = float(value) + if not 0 < value <= 1: + raise ValueError(f"Invalid radiation pattern factor for {phase}: {value}.") + return value + def calculate_Mw(M0): r""" Calculates the moment magnitude :math:`M_w` according to :ref:`Kanamori (1977) ` using formula @@ -608,7 +787,7 @@ def mag_detectable(Mw_min_stations: np.ndarray, min_stations_s : :obj:`int`, optional, default: 3 Minimum number of stations on which event must be detected with S waves wave_mode: :obj:`str`, optional, default: 'PS' - Wave mode to use, can be 'P','S' or 'PS' + Wave mode to use, can be 'P', 'SV', 'SH', 'S' or 'PS' Returns ------- @@ -620,21 +799,24 @@ def mag_detectable(Mw_min_stations: np.ndarray, ValueError : if minimum number of stations on which event must be detected with P or S waves is less than 1 or exceeds number of stations """ + wave_mode = normalize_wave_mode(wave_mode) + if min_stations_p is None: min_stations_p = 3 if min_stations_s is None: min_stations_s = 3 - if not 0 < min_stations_p <= len(Mw_min_stations[0]): - raise ValueError("Invalid minimum number of stations for P waves.") - if not 0 < min_stations_s <= len(Mw_min_stations[1]): - raise ValueError("Invalid minimum number of stations for S waves.") + # Number of stations is the last axis of the (2, n_grid_points, n_stations) array + check_min_stations(min_stations_p=min_stations_p, + min_stations_s=min_stations_s, + n_stations=Mw_min_stations.shape[2], + wave_mode=wave_mode) # Compute final Mw_min_grid based on wave mode if wave_mode == 'P': # For P waves, find the nth smallest value (where n = min_stations_p) Mw_min_grid = np.partition(Mw_min_stations[0], min_stations_p - 1, axis=1)[:, min_stations_p - 1] - elif wave_mode == 'S': + elif wave_mode in ('S', 'SV', 'SH'): # For S waves, find the nth smallest value (where n = min_stations_s) Mw_min_grid = np.partition(Mw_min_stations[1], min_stations_s - 1, axis=1)[:, min_stations_s - 1] else: # wave_mode == 'PS' @@ -645,204 +827,437 @@ def mag_detectable(Mw_min_stations: np.ndarray, return Mw_min_grid -def get_ray_station_directionality(station_coords: np.ndarray, - grid_coords: np.ndarray, - rays_p: list=None, - rays_s: list=None, - strict_nan_check: bool=False): +def check_min_stations(min_stations_p: int, + min_stations_s: int, + n_stations: int, + wave_mode: str = 'PS'): r""" - Derive directionality coefficients for P and S waves - for each station using rays to this station. - Takes into account the sequence of stations. + Check that the minimum numbers of stations required for detection are valid. + + Parameters + ---------- + min_stations_p : :obj:`int` + Minimum number of stations on which event must be detected with P waves. + min_stations_s : :obj:`int` + Minimum number of stations on which event must be detected with S waves. + n_stations : :obj:`int` + Total number of stations in the geometry. + wave_mode : :obj:`str`, optional, default: ``'PS'`` + Wave mode to use, can be 'P', 'SV', 'SH', 'S' or 'PS'. + + Raises + ------ + ValueError + If a minimum number of stations is not positive or exceeds the number of stations. + + Notes + ----- + Each of ``min_stations_p`` and ``min_stations_s`` is checked only if it is not + ``None`` and only if the corresponding branch is part of ``wave_mode``. + """ + wave_mode = normalize_wave_mode(wave_mode) + if wave_mode_has_p(wave_mode) and min_stations_p is not None: + if min_stations_p <= 0: + raise ValueError("Minimum number of stations on which event must be detected with P waves must be positive.") + if min_stations_p > n_stations: + raise ValueError("Minimum number of stations on which event must be detected with P waves must not exceed number of stations.") + if wave_mode_has_s(wave_mode) and min_stations_s is not None: + if min_stations_s <= 0: + raise ValueError("Minimum number of stations on which event must be detected with S waves must be positive.") + if min_stations_s > n_stations: + raise ValueError("Minimum number of stations on which event must be detected with S waves must not exceed number of stations.") + + +def compute_free_surface_coefficients(station_coords: np.ndarray, + grid_coords: np.ndarray, + station_coords_on_surface: np.ndarray = None, + fs_mode: str = 'auto', + fs_level: float = 0.0, + fs_deviation: float = 1.0): + r""" + Compute free surface amplification coefficients for each station and grid point. + + A receiver recording at the free surface sees an amplitude approximately twice + that of the incident wave, which is the exact amplification for normal or + near-normal incidence. The coefficient is applied per station, so that a + geometry combining surface receivers with receivers at depth is treated + correctly. Parameters ---------- station_coords : :obj:`numpy.ndarray` Array of station coordinates [[xr1, yr1, zr1], [xr2, yr2, zr2], ...] - grid_coords : :obj:`numpy.ndarray`, optional + grid_coords : :obj:`numpy.ndarray` Array of grid coordinates [[xs1, ys1, zs1], [xs2, ys2, zs2], ...] - rays_p: :obj:`list` of :obj:`numpy.ndarray` - If provided, must contain exactly one ray path from each (station, grid) pair, - sorted from that station (first row) to that grid point (last row). - The path's shape is ``(m,3)``. - Used for P-wave directionality. - rays_s: :obj:`list` of :obj:`numpy.ndarray` - Same as rays_p but for S-wave directionality. - strict_nan_check: :obj:`bool`, optional, default: False - If True, raise a RuntimeError if NaNs appear in the resulting - directionality coefficients. + station_coords_on_surface : :obj:`numpy.ndarray`, optional + Per-station indicator of whether the station is on the free surface, as a + 1D array of the same length as the number of stations, with values in + [0, 1]. Values between 0 and 1 act as a weighting factor between no + amplification and full amplification. If ``None``, the indicator is + derived from the station depths. + fs_mode : :obj:`str`, optional, default: ``'auto'`` + Free surface correction mode, can be ``'auto'``, ``'on'`` or ``'off'``. + If ``'on'``, the correction is applied to stations close to the free + surface irrespective of ``station_coords_on_surface``. If ``'off'``, the + correction is never applied. If ``'auto'``, ``station_coords_on_surface`` + is used when available, otherwise stations are considered to be on the + free surface when their Z coordinate deviates from the free surface level + by no more than ``fs_deviation``. + fs_level : :obj:`float`, optional, default: 0.0 + Level of the free surface in metres above the reference level. + fs_deviation : :obj:`float`, optional, default: 1.0 + Permitted deviation from the free surface level in metres, must be + non-negative. Returns ------- - dir_coef_p : :obj:`numpy.ndarray` - Array of correction coefficients for P wave, for each grid point and station - dir_coef_s : :obj:`numpy.ndarray` - Array of correction coefficients for S wave, for each grid point and station + fs_coef_p : :obj:`numpy.ndarray` + Free surface coefficients for P waves, of shape (n_grid_points, n_stations) + fs_coef_s : :obj:`numpy.ndarray` + Free surface coefficients for S waves, of shape (n_grid_points, n_stations) Raises ------ - ValueError : - if len(station_coords) is less than 2 - RuntimeError : - if strict_nan_check is True and NaN appears in the result + ValueError + If ``fs_level`` is not a scalar, ``fs_deviation`` is negative, + ``station_coords_on_surface`` has the wrong shape or values outside + [0, 1], or ``fs_mode`` is not one of the supported modes. Notes ----- - For each station it computes the correction coefficient for P and S wave. - If no `rays_p` and `rays_s` are provided, we treat the medium as homogeneous. - For P wave the correction coefficient is equal to the inverse of the absolute value of cosine of the angle between - the line connecting the current station with each grid point (CG line) - and the line connecting the current station with the next station (CN line). - If there is no next station (the current station is the last one), - the line connecting the current station with the previous station (CP line) is used instead of the CN line. - If the angle between the CN and CP lines is larger than 30 degrees - or if the distance between the current station and the next station (CN distance) - is 10 times larger than the distance between the current station and the previous station (CP distance), - the CP line is used for calculation of the angle with CG line, whereas for the next station the CN line is always used. - For S wave the correction coefficient is equal to the inverse of the absolute value of the sine of the same angle as described for the P wave. - This is a vectorized version of the function. It computes the directionality - coefficients for all grid points and stations simultaneously, which should be much - more efficient for large datasets. - - If `rays_p` is provided, we do NOT assume a simple straight ray (line) - from grid point to station; instead, we use the ray path for that pair, - extract only the first segment near the station, and treat that as CG. - Ray is considered to start from station. - The same logic applies if `rays_s` is provided for the S-wave. + The free surface level is negated internally to match the convention that Z + increases downwards. """ - # Check number of stations - if len(station_coords) < 2: - raise ValueError("There must be at least two stations.") - - # Determine which wave(s) we are actually computing - # per user specification: - # - if both rays_p and rays_s are None => compute both P & S, assume straight ray - # - if rays_p is not None and rays_s is None => only P - # - if rays_p is None and rays_s is not None => only S - # - if both are provided => both P & S - straight_rays = False - if rays_p is not None and rays_s is not None: - compute_p, compute_s = True, True - elif rays_p is not None: # only rays_p is provided - compute_p, compute_s = True, False - elif rays_s is not None: # only rays_s is provided - compute_p, compute_s = False, True - else: - straight_rays = True - compute_p, compute_s = True, True # neither is provided => both + if station_coords.ndim == 1 and station_coords.size == 3: + station_coords = station_coords[np.newaxis, :] + if grid_coords.ndim == 1 and grid_coords.size == 3: + grid_coords = grid_coords[np.newaxis, :] - # Get lengths n_stations = station_coords.shape[0] n_grid_points = grid_coords.shape[0] - # Calculate vectors between stations - station_vectors = np.diff(station_coords, axis=0) + # Negate the free surface level to match the downward-positive Z convention + if np.ndim(fs_level) != 0: + raise ValueError("Free surface level must be a scalar.") + fs_depth = -fs_level + + if fs_deviation < 0: + raise ValueError("Free surface level deviation must be non-negative.") + + # Check and normalize the optional per-station weights + surface_weights = None + if station_coords_on_surface is not None: + station_coords_on_surface = np.asarray(station_coords_on_surface) + if station_coords_on_surface.ndim != 1 or station_coords_on_surface.shape[0] != n_stations: + raise ValueError("station_coords_on_surface must be a 1D array with the same length as the number of stations.") + surface_weights = station_coords_on_surface.astype(float, copy=False) + if np.any((surface_weights < 0) | (surface_weights > 1)): + raise ValueError("station_coords_on_surface must contain only values in the range [0, 1].") + + # Simple free-surface approximation: amplitudes are doubled at the surface + full_fs_coef = 2.0 + + # Decide how much of the full amplification applies to each station + if fs_mode == 'auto': + if surface_weights is None: + surface_weights = (np.abs(station_coords[:, 2] - fs_depth) <= fs_deviation).astype(float) + elif fs_mode == 'on': + surface_weights = (np.abs(station_coords[:, 2] - fs_depth) <= fs_deviation).astype(float) + elif fs_mode == 'off': + surface_weights = np.zeros(n_stations, dtype=float) + else: + raise ValueError("Invalid free surface mode. Must be 'auto', 'on', or 'off'.") + + # Interpolate between no amplification (1.0) and full amplification + fs_coef = np.broadcast_to( + 1.0 + (full_fs_coef - 1.0) * surface_weights[np.newaxis, :], + (n_grid_points, n_stations) + ).astype(float, copy=True) + + # The same coefficients are used for P and S waves + return fs_coef, fs_coef + - # Pad the first and last rows to handle edge cases +def _unit_vectors(vectors: np.ndarray) -> np.ndarray: + r"""Normalize vectors along the last axis, leaving zero-length rows as NaN.""" + norms = np.linalg.norm(vectors, axis=-1, keepdims=True) + with np.errstate(invalid='ignore', divide='ignore'): + return vectors / norms + + +def _station_reference_vectors(station_coords: np.ndarray) -> np.ndarray: + r""" + Estimate the local cable tangent at each channel from its neighbours. + + The forward difference to the next channel is used by default. Where the cable + turns by more than 30 degrees, or where the forward step exceeds ten times the + backward step, the backward difference is used instead, so that a channel at a + kink or at the end of a segment still receives a sensible tangent. + """ + if len(station_coords) < 2: + raise ValueError("There must be at least two stations.") + + station_vectors = np.diff(station_coords, axis=0) cn_vectors = np.vstack([station_vectors, station_vectors[-1]]) cp_vectors = np.vstack([station_vectors[0], station_vectors]) - # Calculate distances for distance checks - cn_distances = np.linalg.norm(cn_vectors, axis=1) - cp_distances = np.linalg.norm(cp_vectors, axis=1) + eps = 1e-12 + cn_distances = np.linalg.norm(cn_vectors, axis=1) + eps + cp_distances = np.linalg.norm(cp_vectors, axis=1) + eps - # Calculate angles between CN and CP vectors - angles_cn_cp = np.arccos(np.sum(cn_vectors * cp_vectors, axis=1) / - (cn_distances * cp_distances)) - - # Determine whether to use CP vector based on angle and distance conditions + cos_vals = np.sum(cn_vectors * cp_vectors, axis=1) / (cn_distances * cp_distances) + cos_vals = np.clip(cos_vals, -1.0, 1.0) + angles_cn_cp = np.arccos(cos_vals) use_cp = (angles_cn_cp > np.radians(30)) | (cn_distances > 10 * cp_distances) + return np.where(use_cp[:, np.newaxis], cp_vectors, cn_vectors) + - # Prepare vectors for angle calculation - vectors_for_angle = np.where(use_cp[:, np.newaxis], cp_vectors, cn_vectors) - - # Function to build cg vector using P and S rays - def build_cg_vectors(rays): - r""" - Given ray paths in ``rays``, build an array of CG vectors - of shape ``(n_grid_points, n_stations, 3)``. - Use the first segment of the provided ray path near the station. - """ - - cg = np.zeros((n_grid_points, n_stations, 3), dtype=float) - - for i_st in range(n_stations): - for j_gr in range(n_grid_points): - iray = i_st*n_grid_points + j_gr - path = rays[iray] - if path is None: - raise ValueError(f"No ray found for station #{i_st} and grid #{j_gr}") - if path.shape[0] < 2: - raise ValueError( - f"Ray path for station #{i_st}, grid #{j_gr} " - f"has too few points ({path.shape[0]})." - ) - stc = station_coords[i_st] - if not np.all(np.isclose(path[0], stc, atol=1e-03)): - #raise ValueError(f"Ray {iray} ({path[0]}) does not start at the station {i_st} ({stc})") - warnings.warn(f"Ray {iray} start point ({path[0]}) differ from station {i_st} coordinates ({stc})") - cg[j_gr, i_st, :] = path[1] - path[0] - return cg - - if straight_rays: - cg_vectors = grid_coords[:, np.newaxis, :] - station_coords[np.newaxis, :, :] +def _polarization_vectors(ray_vectors: np.ndarray, phase: str) -> np.ndarray: + r""" + Arriving polarization unit vectors for a phase. + + P is polarized along the ray, SH horizontally and perpendicular to the vertical + ray plane, and SV perpendicular to both. For a vertical ray the SH direction is + degenerate, so a stable horizontal fallback basis is used instead. + """ + phase = normalize_wave_mode(phase) + n = _unit_vectors(ray_vectors) + if phase == "P": + return n + + vertical = np.array([0.0, 0.0, 1.0]) + sh = np.cross(vertical, n) + sh_norm = np.linalg.norm(sh, axis=-1) + + fallback_x = np.cross(np.array([1.0, 0.0, 0.0]), n) + fallback_y = np.cross(np.array([0.0, 1.0, 0.0]), n) + fallback_x_norm = np.linalg.norm(fallback_x, axis=-1) + fallback = np.where((fallback_x_norm > 1e-12)[..., np.newaxis], fallback_x, fallback_y) + sh = np.where((sh_norm > 1e-12)[..., np.newaxis], sh, fallback) + sh = _unit_vectors(sh) + + if phase == "SH": + return sh + if phase == "SV": + return _unit_vectors(np.cross(sh, n)) + + raise ValueError("Polarization vectors are defined only for P, SV, and SH.") + + +def _validate_station_components(components) -> pd.Series: + r"""Validate a ``Components`` column, filling missing entries with ``'3C'``.""" + component_series = pd.Series(components).dropna() + invalid_components = sorted(set(component_series).difference({'3C', 'Z'})) + if invalid_components: + raise ValueError( + "Components values must be either '3C' or 'Z'. " + f"Invalid values: {', '.join(invalid_components)}" + ) + return pd.Series(components).fillna('3C') + + +def geometry_requires_receiver_projection(geometry_df: Optional[pd.DataFrame]) -> bool: + r""" + Return whether a geometry needs receiver-component projection. + + A DAS geometry, identified by the presence of a gauge length, always requires + projection onto the local cable tangent. A station geometry requires projection + only if at least one station records a single vertical component. + + Parameters + ---------- + geometry_df : :obj:`pandas.DataFrame` or ``None`` + Geometry DataFrame. + + Returns + ------- + required : :obj:`bool` + Whether projection is required. + """ + if geometry_df is None: + return False + if geometry_df.attrs.get('gauge_length') is not None: + return True + if 'Components' not in geometry_df.columns: + return False + components = _validate_station_components(geometry_df['Components']) + return bool((components == 'Z').any()) + + +def get_phase_station_directionality(station_coords: np.ndarray, + grid_coords: np.ndarray, + phases=("P", "SV", "SH"), + geometry_df: Optional[pd.DataFrame] = None, + force_tangent: bool = False, + strict_nan_check: bool = False) -> dict: + r""" + Derive receiver projection coefficients for P, SV and SH polarizations. + + The recorded amplitude of a phase is the projection of its arriving + polarization onto the direction the receiver is sensitive to. For a phase + :math:`\phi` the correction applied to the minimum detectable amplitude is + + .. math:: + C_{\phi} = \frac{1}{\left| \mathbf{d} \cdot \mathbf{e}_{\phi} \right|}, + + where :math:`\mathbf{e}_{\phi}` is the arriving polarization and + :math:`\mathbf{d}` is the local cable tangent for DAS, or the vertical axis for + a single-component station. Three-component stations record the full vector and + therefore use a coefficient of one. + + Parameters + ---------- + station_coords : :obj:`numpy.ndarray` + Array of station coordinates [[xr1, yr1, zr1], [xr2, yr2, zr2], ...] + grid_coords : :obj:`numpy.ndarray` + Array of grid coordinates [[xs1, ys1, zs1], [xs2, ys2, zs2], ...] + phases : sequence of :obj:`str`, optional, default: ``("P", "SV", "SH")`` + Phases to compute coefficients for. Must be drawn from 'P', 'SV' and 'SH'. + geometry_df : :obj:`pandas.DataFrame`, optional + Geometry DataFrame, used to detect a DAS geometry from its gauge length and + to read a ``Components`` column for station geometries. If ``None``, + projection onto the local tangent is used. + force_tangent : :obj:`bool`, optional, default: ``False`` + Force projection onto the local tangent even when ``geometry_df`` describes + a station geometry. + strict_nan_check : :obj:`bool`, optional, default: ``False`` + Flag to check whether NaN values appear in the result. + + Returns + ------- + coefficients : :obj:`dict` + Mapping from phase name to an array of coefficients of shape + (n_grid_points, n_stations). + + Raises + ------ + ValueError + If a requested phase is not P, SV or SH, if the geometry and station + coordinates disagree in length, or if a DAS geometry also carries a + ``Components`` column. + RuntimeError + If ``strict_nan_check`` is True and NaN appears in the result. + + Notes + ----- + A vanishing projection means the phase cannot be recorded at all, for example + SH on a vertical cable, and yields NaN so that the phase is excluded rather + than reported as infinitely detectable. + + Straight rays between each station and each grid point are assumed, which is + exact in a homogeneous medium. + """ + station_coords = np.asarray(station_coords, dtype=float) + grid_coords = np.asarray(grid_coords, dtype=float) + if station_coords.ndim == 1 and station_coords.size == 3: + station_coords = station_coords[np.newaxis, :] + if grid_coords.ndim == 1 and grid_coords.size == 3: + grid_coords = grid_coords[np.newaxis, :] + + phases = tuple(dict.fromkeys(normalize_wave_mode(phase) for phase in phases)) + if any(phase not in ("P", "SV", "SH") for phase in phases): + raise ValueError("Directionality phases must be P, SV, or SH.") + + use_tangent_projection = force_tangent or geometry_df is None + components = None + if geometry_df is not None: + if len(geometry_df) != len(station_coords): + raise ValueError("Geometry and station coordinates must have the same number of stations.") + has_gauge_length = geometry_df.attrs.get('gauge_length') is not None + has_components = 'Components' in geometry_df.columns + if has_gauge_length: + if has_components and geometry_df['Components'].notna().any(): + raise ValueError("Components column is only allowed for station geometries, not DAS geometries.") + use_tangent_projection = True + elif has_components: + components = _validate_station_components(geometry_df['Components']).to_numpy() + else: + components = np.full(len(station_coords), '3C', dtype=object) + + if use_tangent_projection: + reference_unit = _unit_vectors(_station_reference_vectors(station_coords)) else: - cg_vectors_p = build_cg_vectors(rays_p) if compute_p else None - cg_vectors_s = build_cg_vectors(rays_s) if compute_s else None - - # Function to compute directionality from a set of CG vectors - # cg_vectors shape: (n_grid_points, n_stations, 3) - def compute_dir_coeffs_for_cg_vectors(cg_vectors): - r""" - Given ``cg_vectors`` of shape ``(n_grid_points, n_stations, 3)``, compute - ``dir_coef_p`` and ``dir_coef_s`` using the angle with ``vectors_for_angle``. - Returns ``(dir_coef_p, dir_coef_s)`` each shaped as ``(n_grid_points, n_stations)``. - """ - # Normalize vectors_for_angle per station - # shape: (n_stations, 3) - station_ref_norm = vectors_for_angle / np.linalg.norm(vectors_for_angle, axis=1, keepdims=True) - - # Now we want the dot-product with each cg_vector - # cg_vectors shape: (n_grid_points, n_stations, 3) - # We'll normalize cg_vectors along axis=2 - cg_norm = cg_vectors / np.linalg.norm(cg_vectors, axis=2, keepdims=True) - - # cos(angle) = abs( dot( cg_norm, station_ref_norm ) ) - # We'll expand station_ref_norm so we can broadcast: shape -> (1, n_stations, 3) - dot_vals = np.sum(cg_norm * station_ref_norm[np.newaxis, :, :], axis=2) - cos_vals = np.abs(dot_vals) - # sin(angle) = sqrt(1 - cos^2(angle)) - sin_vals = np.sqrt(1.0 - cos_vals**2) - - # Zero or near-zero cos/sin => we'd get inf => we can turn them into NaN - cos_near_zero = np.isclose(cos_vals, 0.0) - sin_near_zero = np.isclose(sin_vals, 0.0) - - cos_vals[cos_near_zero] = np.nan - sin_vals[sin_near_zero] = np.nan - - dir_coef_p = 1.0 / cos_vals # shape (n_grid_points, n_stations) - dir_coef_s = 1.0 / sin_vals # shape (n_grid_points, n_stations) - - if strict_nan_check: - if np.any(np.isnan(dir_coef_p)) or np.any(np.isnan(dir_coef_s)): - current_function = inspect.currentframe().f_code.co_name - raise RuntimeError( - f"NaN values detected while computing directionality " - f"coefficients in {current_function}." + vertical_unit = np.array([0.0, 0.0, 1.0]) + + # Straight rays from each station to each grid point + ray_vectors = grid_coords[:, np.newaxis, :] - station_coords[np.newaxis, :, :] + + result = {} + for phase in phases: + polarization = _polarization_vectors(ray_vectors, phase) + if use_tangent_projection: + projection = np.abs( + np.sum(polarization * reference_unit[np.newaxis, :, :], axis=2) + ) + projection[np.isclose(projection, 0.0)] = np.nan + result[phase] = 1.0 / projection + else: + coefficients = np.ones(polarization.shape[:2], dtype=float) + z_mask = components == 'Z' + if np.any(z_mask): + projection = np.abs( + np.sum(polarization[:, z_mask, :] * vertical_unit, axis=2) ) + projection[np.isclose(projection, 0.0)] = np.nan + coefficients[:, z_mask] = 1.0 / projection + result[phase] = coefficients - return dir_coef_p, dir_coef_s + if strict_nan_check: + if any(np.any(np.isnan(values)) for values in result.values()): + current_function = inspect.currentframe().f_code.co_name + raise RuntimeError( + f"NaN values detected while computing directionality coefficients in {current_function}." + ) - # Compute directionality coefficients for whichever waves we need - if straight_rays: - dir_coef_p, dir_coef_s = compute_dir_coeffs_for_cg_vectors(cg_vectors) - else: - dir_coef_p = None - dir_coef_s = None - if compute_p: - dir_coef_p, _ = compute_dir_coeffs_for_cg_vectors(cg_vectors_p) - if compute_s: - _, dir_coef_s = compute_dir_coeffs_for_cg_vectors(cg_vectors_s) - - return dir_coef_p, dir_coef_s + return result + + +def get_ray_station_directionality(station_coords: np.ndarray, + grid_coords: np.ndarray, + strict_nan_check: bool = False): + r""" + Derive directionality coefficients for P and S waves along a cable. + + .. deprecated:: + Use :func:`get_phase_station_directionality`, which resolves SV and SH + separately and projects the actual polarization vectors. + + Parameters + ---------- + station_coords : :obj:`numpy.ndarray` + Array of station coordinates [[xr1, yr1, zr1], [xr2, yr2, zr2], ...] + grid_coords : :obj:`numpy.ndarray` + Array of grid coordinates [[xs1, ys1, zs1], [xs2, ys2, zs2], ...] + strict_nan_check : :obj:`bool`, optional, default: ``False`` + Flag to check whether NaN values appear in the result. + + Returns + ------- + dir_coef_p : :obj:`numpy.ndarray` + Directionality coefficients for P waves, of shape (n_grid_points, n_stations) + dir_coef_s : :obj:`numpy.ndarray` + Directionality coefficients for SV waves, of shape (n_grid_points, n_stations) + + Notes + ----- + Returns the P and SV coefficients of + :func:`get_phase_station_directionality` with tangent projection forced. For a + vertical cable these are the familiar :math:`1/|\cos\theta|` and + :math:`1/|\sin\theta|`; for an inclined or curved cable they are the correct + polarization projections rather than that approximation. + """ + warnings.warn( + "get_ray_station_directionality is deprecated; use " + "get_phase_station_directionality, which resolves SV and SH separately.", + DeprecationWarning, + stacklevel=2, + ) + coefficients = get_phase_station_directionality( + station_coords=station_coords, + grid_coords=grid_coords, + phases=("P", "SV"), + geometry_df=None, + force_tangent=True, + strict_nan_check=strict_nan_check, + ) + return coefficients["P"], coefficients["SV"] diff --git a/pytests/test_detectability.py b/pytests/test_detectability.py new file mode 100644 index 0000000..eb5995d --- /dev/null +++ b/pytests/test_detectability.py @@ -0,0 +1,107 @@ +"""Tests for the network detectability criterion. + +The network threshold at a grid point is the N-th smallest per-station threshold, +so that the event is required to be detected on at least N receivers. For joint +P and S detection the stricter of the two thresholds applies. +""" +import numpy as np +import pytest + +from pynetdesign.modelling.utils import check_min_stations, mag_detectable + + +def stations_grid(p_values, s_values): + """Build a (2, n_grid, n_stations) array from per-station values.""" + p = np.asarray(p_values, dtype=float) + s = np.asarray(s_values, dtype=float) + if p.ndim == 1: + p = p[np.newaxis, :] + if s.ndim == 1: + s = s[np.newaxis, :] + return np.stack([p, s]) + + +class TestNthSmallest: + def test_p_mode_picks_nth_smallest(self): + mw = stations_grid([0.5, -1.0, 0.2, -0.4], [9.0, 9.0, 9.0, 9.0]) + # sorted P: -1.0, -0.4, 0.2, 0.5 + assert mag_detectable(mw, min_stations_p=1, wave_mode='P')[0] == pytest.approx(-1.0) + assert mag_detectable(mw, min_stations_p=2, wave_mode='P')[0] == pytest.approx(-0.4) + assert mag_detectable(mw, min_stations_p=3, wave_mode='P')[0] == pytest.approx(0.2) + assert mag_detectable(mw, min_stations_p=4, wave_mode='P')[0] == pytest.approx(0.5) + + def test_more_stations_required_means_higher_threshold(self): + mw = stations_grid(np.linspace(-2.0, 1.0, 8), np.full(8, 9.0)) + thresholds = [mag_detectable(mw, min_stations_p=n, wave_mode='P')[0] + for n in range(1, 9)] + assert np.all(np.diff(thresholds) > 0) + + @pytest.mark.parametrize("mode", ["S", "SV", "SH"]) + def test_s_family_modes_use_the_s_row(self, mode): + mw = stations_grid([9.0, 9.0, 9.0], [0.3, -0.7, 0.1]) + assert mag_detectable(mw, min_stations_s=2, wave_mode=mode)[0] == pytest.approx(0.1) + + def test_default_is_three_stations(self): + mw = stations_grid([0.5, -1.0, 0.2, -0.4], [0.5, -1.0, 0.2, -0.4]) + assert mag_detectable(mw, wave_mode='P')[0] == pytest.approx(0.2) + + +class TestJointDetection: + def test_ps_takes_the_stricter_branch(self): + mw = stations_grid([0.0, 0.0, 0.0], [-1.0, -1.0, -1.0]) + # P needs the larger magnitude, so P governs + assert mag_detectable(mw, min_stations_p=3, min_stations_s=3, + wave_mode='PS')[0] == pytest.approx(0.0) + + def test_ps_is_never_easier_than_either_branch(self): + rng = np.random.default_rng(0) + p = rng.normal(0.0, 1.0, (20, 6)) + s = rng.normal(-0.5, 1.0, (20, 6)) + mw = np.stack([p, s]) + ps = mag_detectable(mw, min_stations_p=3, min_stations_s=3, wave_mode='PS') + only_p = mag_detectable(mw, min_stations_p=3, wave_mode='P') + only_s = mag_detectable(mw, min_stations_s=3, wave_mode='S') + assert np.all(ps >= only_p - 1e-12) + assert np.all(ps >= only_s - 1e-12) + + def test_independent_station_counts(self): + mw = stations_grid([0.5, -1.0, 0.2, -0.4], [0.9, -1.2, 0.4, -0.2]) + # P with 2 stations -> -0.4 ; S with 4 stations -> 0.9 ; max -> 0.9 + assert mag_detectable(mw, min_stations_p=2, min_stations_s=4, + wave_mode='PS')[0] == pytest.approx(0.9) + + +class TestGuards: + def test_rejects_more_stations_than_available(self): + """Regression guard: the check must look at the station axis, not the grid + axis. This array has 2 grid points and 3 stations.""" + mw = np.full((2, 2, 3), 0.0) + with pytest.raises(ValueError): + mag_detectable(mw, min_stations_p=4, wave_mode='P') + + def test_accepts_exactly_all_stations(self): + mw = np.full((2, 2, 3), 0.0) + assert mag_detectable(mw, min_stations_p=3, wave_mode='P').shape == (2,) + + def test_rejects_non_positive(self): + mw = np.full((2, 2, 3), 0.0) + with pytest.raises(ValueError): + mag_detectable(mw, min_stations_p=0, wave_mode='P') + + def test_only_the_requested_branch_is_checked(self): + """An impossible S count must not break a P-only request.""" + mw = np.full((2, 2, 3), 0.0) + assert mag_detectable(mw, min_stations_p=2, min_stations_s=99, + wave_mode='P').shape == (2,) + + def test_rejects_invalid_wave_mode(self): + mw = np.full((2, 2, 3), 0.0) + with pytest.raises(ValueError): + mag_detectable(mw, wave_mode='Q') + + def test_check_min_stations_directly(self): + check_min_stations(min_stations_p=3, min_stations_s=3, n_stations=3, wave_mode='PS') + with pytest.raises(ValueError): + check_min_stations(min_stations_p=4, min_stations_s=1, n_stations=3, wave_mode='PS') + # None means "unset" and is skipped + check_min_stations(min_stations_p=None, min_stations_s=None, n_stations=1, wave_mode='PS') diff --git a/pytests/test_directionality.py b/pytests/test_directionality.py new file mode 100644 index 0000000..eec7468 --- /dev/null +++ b/pytests/test_directionality.py @@ -0,0 +1,192 @@ +"""Tests for receiver projection of the arriving polarization. + +For a vertical cable the projection factors reduce to the familiar +|cos(theta)| for P and |sin(theta)| for SV, where theta is the angle between the +propagation direction and the cable, and SH cannot be recorded at all. These +closed forms are used as the reference here. +""" +import numpy as np +import pandas as pd +import pytest + +from pynetdesign.modelling.utils import ( + geometry_requires_receiver_projection, + get_phase_station_directionality, + get_ray_station_directionality, +) + + +def vertical_cable(n=5, spacing=100.0): + """Vertical cable at the origin, Z increasing downwards.""" + return np.column_stack([np.zeros(n), np.zeros(n), np.arange(n) * spacing]) + + +def das_geometry(station_coords, gauge_length=10.0): + df = pd.DataFrame({ + 'Name': [f"C{i+1}" for i in range(len(station_coords))], + 'X': station_coords[:, 0], + 'Y': station_coords[:, 1], + 'Z': station_coords[:, 2], + 'NoiseLevel': np.full(len(station_coords), 1e-9), + }) + df.attrs['noise_type'] = 'strain rate' + df.attrs['gauge_length'] = gauge_length + return df + + +class TestVerticalCable: + @pytest.mark.parametrize("theta_deg", [15.0, 30.0, 45.0, 60.0, 75.0]) + def test_p_and_sv_match_closed_forms(self, theta_deg): + theta = np.radians(theta_deg) + stations = vertical_cable() + # One grid point at angle theta from the cable, seen from the first channel + r = 1000.0 + grid = np.array([[r * np.sin(theta), 0.0, stations[0, 2] + r * np.cos(theta)]]) + + coefficients = get_phase_station_directionality( + station_coords=stations, grid_coords=grid, force_tangent=True) + + assert coefficients["P"][0, 0] == pytest.approx(1.0 / abs(np.cos(theta)), rel=1e-12) + assert coefficients["SV"][0, 0] == pytest.approx(1.0 / abs(np.sin(theta)), rel=1e-12) + + def test_sh_is_unobservable_on_a_vertical_cable(self): + """SH is horizontal and perpendicular to the ray plane, so a vertical + cable has no sensitivity to it at all.""" + stations = vertical_cable() + grid = np.array([[500.0, 0.0, 250.0], [1500.0, 200.0, 800.0]]) + coefficients = get_phase_station_directionality( + station_coords=stations, grid_coords=grid, force_tangent=True) + assert np.all(np.isnan(coefficients["SH"])) + + def test_axial_ray_is_blind_to_sv(self): + """A ray travelling along the cable gives sin(theta) = 0, so SV cannot be + recorded. This is the low-sensitivity zone directly below a vertical array.""" + stations = vertical_cable() + grid = np.array([[0.0, 0.0, 2000.0]]) # straight below the cable + coefficients = get_phase_station_directionality( + station_coords=stations, grid_coords=grid, force_tangent=True) + assert np.all(np.isnan(coefficients["SV"])) + assert np.allclose(coefficients["P"], 1.0) # fully axial, perfect for P + + def test_broadside_ray_is_blind_to_p(self): + """A ray arriving perpendicular to the cable gives cos(theta) = 0.""" + stations = vertical_cable() + grid = np.array([[1000.0, 0.0, stations[2, 2]]]) # level with channel 3 + coefficients = get_phase_station_directionality( + station_coords=stations, grid_coords=grid, force_tangent=True) + assert np.isnan(coefficients["P"][0, 2]) + assert coefficients["SV"][0, 2] == pytest.approx(1.0, rel=1e-12) + + def test_coefficients_are_at_least_one(self): + """Projection can only reduce the recorded amplitude, so the correction + applied to the threshold is never smaller than one.""" + stations = vertical_cable() + grid = np.array([[300.0, 100.0, 150.0], [1200.0, -400.0, 900.0]]) + coefficients = get_phase_station_directionality( + station_coords=stations, grid_coords=grid, force_tangent=True) + for phase in ("P", "SV"): + finite = coefficients[phase][np.isfinite(coefficients[phase])] + assert np.all(finite >= 1.0 - 1e-12) + + +class TestGeometryDrivenProjection: + def test_das_geometry_requires_projection(self): + stations = vertical_cable() + assert geometry_requires_receiver_projection(das_geometry(stations)) is True + + def test_plain_station_geometry_does_not(self): + df = pd.DataFrame({'Name': ['A', 'B'], 'X': [0.0, 1.0], 'Y': [0.0, 0.0], + 'Z': [0.0, 0.0], 'NoiseLevel': [1e-9, 1e-9]}) + df.attrs['noise_type'] = 'velocity' + assert geometry_requires_receiver_projection(df) is False + + def test_none_geometry_does_not(self): + assert geometry_requires_receiver_projection(None) is False + + def test_three_component_stations_are_not_projected(self): + """A 3C station records the full vector, so no amplitude is lost.""" + stations = vertical_cable(n=3) + df = pd.DataFrame({'Name': ['A', 'B', 'C'], + 'X': stations[:, 0], 'Y': stations[:, 1], 'Z': stations[:, 2], + 'NoiseLevel': np.full(3, 1e-9), + 'Components': ['3C', '3C', '3C']}) + df.attrs['noise_type'] = 'velocity' + grid = np.array([[500.0, 0.0, 250.0]]) + coefficients = get_phase_station_directionality( + station_coords=stations, grid_coords=grid, geometry_df=df) + for phase in ("P", "SV", "SH"): + assert np.allclose(coefficients[phase], 1.0) + + def test_vertical_component_station_is_projected(self): + """A Z-only station sees the vertical projection of the polarization.""" + stations = vertical_cable(n=2) + df = pd.DataFrame({'Name': ['A', 'B'], + 'X': stations[:, 0], 'Y': stations[:, 1], 'Z': stations[:, 2], + 'NoiseLevel': np.full(2, 1e-9), + 'Components': ['Z', '3C']}) + df.attrs['noise_type'] = 'velocity' + theta = np.radians(35.0) + r = 1000.0 + grid = np.array([[r * np.sin(theta), 0.0, r * np.cos(theta)]]) + coefficients = get_phase_station_directionality( + station_coords=stations, grid_coords=grid, geometry_df=df) + # Station 0 records only Z, so P is scaled by 1/|cos(theta)| + assert coefficients["P"][0, 0] == pytest.approx(1.0 / abs(np.cos(theta)), rel=1e-12) + # Station 1 is 3C and therefore unaffected + assert coefficients["P"][0, 1] == pytest.approx(1.0) + + def test_rejects_invalid_components(self): + stations = vertical_cable(n=2) + df = pd.DataFrame({'Name': ['A', 'B'], + 'X': stations[:, 0], 'Y': stations[:, 1], 'Z': stations[:, 2], + 'NoiseLevel': np.full(2, 1e-9), + 'Components': ['Z', 'N']}) + df.attrs['noise_type'] = 'velocity' + with pytest.raises(ValueError): + get_phase_station_directionality(station_coords=stations, + grid_coords=np.array([[100.0, 0.0, 50.0]]), + geometry_df=df) + + def test_rejects_components_on_a_das_geometry(self): + stations = vertical_cable(n=3) + df = das_geometry(stations) + df['Components'] = ['Z', 'Z', 'Z'] + with pytest.raises(ValueError): + get_phase_station_directionality(station_coords=stations, + grid_coords=np.array([[100.0, 0.0, 50.0]]), + geometry_df=df) + + +class TestValidation: + def test_rejects_unknown_phase(self): + with pytest.raises(ValueError): + get_phase_station_directionality(station_coords=vertical_cable(), + grid_coords=np.array([[100.0, 0.0, 50.0]]), + phases=("P", "S")) + + def test_requires_two_stations_for_tangent(self): + with pytest.raises(ValueError): + get_phase_station_directionality( + station_coords=np.array([[0.0, 0.0, 0.0]]), + grid_coords=np.array([[100.0, 0.0, 50.0]]), + force_tangent=True) + + def test_strict_nan_check_raises(self): + stations = vertical_cable() + grid = np.array([[0.0, 0.0, 2000.0]]) # axial ray, SV unobservable + with pytest.raises(RuntimeError): + get_phase_station_directionality(station_coords=stations, grid_coords=grid, + force_tangent=True, strict_nan_check=True) + + +class TestDeprecatedWrapper: + def test_warns_and_returns_p_and_sv(self): + stations = vertical_cable() + theta = np.radians(40.0) + r = 1000.0 + grid = np.array([[r * np.sin(theta), 0.0, r * np.cos(theta)]]) + with pytest.warns(DeprecationWarning): + dir_p, dir_s = get_ray_station_directionality(station_coords=stations, + grid_coords=grid) + assert dir_p[0, 0] == pytest.approx(1.0 / abs(np.cos(theta)), rel=1e-12) + assert dir_s[0, 0] == pytest.approx(1.0 / abs(np.sin(theta)), rel=1e-12) diff --git a/pytests/test_free_surface.py b/pytests/test_free_surface.py new file mode 100644 index 0000000..28ebeb3 --- /dev/null +++ b/pytests/test_free_surface.py @@ -0,0 +1,125 @@ +"""Tests for the per-station free surface amplification. + +A receiver at the free surface records approximately twice the amplitude of the +incident wave, which lowers its detection threshold by 2/3*log10(2) ~ 0.201 +magnitude units. The correction must apply only to receivers that are actually at +the surface, so that a geometry mixing surface and downhole receivers is handled +correctly. +""" +import numpy as np +import pytest + +from pynetdesign.modelling.utils import compute_free_surface_coefficients + +GRID = np.array([[500.0, 0.0, 250.0], [1500.0, 200.0, 900.0]]) + + +def mixed_geometry(): + """Two receivers at the surface and two at depth.""" + return np.array([ + [0.0, 0.0, 0.0], + [100.0, 0.0, 0.5], # within the default 1 m deviation + [0.0, 0.0, 500.0], + [100.0, 0.0, 1000.0], + ]) + + +class TestModes: + def test_auto_gates_on_depth(self): + fs_p, fs_s = compute_free_surface_coefficients( + station_coords=mixed_geometry(), grid_coords=GRID, fs_mode='auto') + assert fs_p.shape == (len(GRID), 4) + np.testing.assert_allclose(fs_p[0], [2.0, 2.0, 1.0, 1.0]) + np.testing.assert_allclose(fs_p, fs_s) + + def test_off_disables_everything(self): + fs_p, _ = compute_free_surface_coefficients( + station_coords=mixed_geometry(), grid_coords=GRID, fs_mode='off') + assert np.all(fs_p == 1.0) + + def test_on_ignores_supplied_weights(self): + """'on' re-derives the indicator from depth rather than trusting weights.""" + weights = np.zeros(4) + fs_p, _ = compute_free_surface_coefficients( + station_coords=mixed_geometry(), grid_coords=GRID, + station_coords_on_surface=weights, fs_mode='on') + np.testing.assert_allclose(fs_p[0], [2.0, 2.0, 1.0, 1.0]) + + def test_auto_honours_supplied_weights(self): + weights = np.array([0.0, 1.0, 1.0, 0.0]) + fs_p, _ = compute_free_surface_coefficients( + station_coords=mixed_geometry(), grid_coords=GRID, + station_coords_on_surface=weights, fs_mode='auto') + np.testing.assert_allclose(fs_p[0], [1.0, 2.0, 2.0, 1.0]) + + def test_fractional_weights_blend(self): + weights = np.array([0.5, 0.25, 0.0, 1.0]) + fs_p, _ = compute_free_surface_coefficients( + station_coords=mixed_geometry(), grid_coords=GRID, + station_coords_on_surface=weights, fs_mode='auto') + np.testing.assert_allclose(fs_p[0], [1.5, 1.25, 1.0, 2.0]) + + def test_rejects_unknown_mode(self): + with pytest.raises(ValueError): + compute_free_surface_coefficients(station_coords=mixed_geometry(), + grid_coords=GRID, fs_mode='sometimes') + + +class TestLevelAndDeviation: + def test_level_shifts_the_surface(self): + """fs_level is measured upwards, Z downwards, so a level of 100 m puts the + surface at Z = -100 m.""" + stations = np.array([[0.0, 0.0, -100.0], [0.0, 0.0, 0.0]]) + fs_p, _ = compute_free_surface_coefficients( + station_coords=stations, grid_coords=GRID, fs_mode='auto', fs_level=100.0) + np.testing.assert_allclose(fs_p[0], [2.0, 1.0]) + + def test_deviation_widens_the_band(self): + stations = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 5.0], [0.0, 0.0, 50.0]]) + tight, _ = compute_free_surface_coefficients( + station_coords=stations, grid_coords=GRID, fs_mode='auto', fs_deviation=1.0) + loose, _ = compute_free_surface_coefficients( + station_coords=stations, grid_coords=GRID, fs_mode='auto', fs_deviation=10.0) + np.testing.assert_allclose(tight[0], [2.0, 1.0, 1.0]) + np.testing.assert_allclose(loose[0], [2.0, 2.0, 1.0]) + + def test_rejects_negative_deviation(self): + with pytest.raises(ValueError): + compute_free_surface_coefficients(station_coords=mixed_geometry(), + grid_coords=GRID, fs_deviation=-1.0) + + def test_rejects_non_scalar_level(self): + with pytest.raises(ValueError): + compute_free_surface_coefficients(station_coords=mixed_geometry(), + grid_coords=GRID, + fs_level=np.array([0.0, 1.0])) + + +class TestWeightValidation: + def test_rejects_wrong_length(self): + with pytest.raises(ValueError): + compute_free_surface_coefficients( + station_coords=mixed_geometry(), grid_coords=GRID, + station_coords_on_surface=np.zeros(3)) + + @pytest.mark.parametrize("bad", [-0.1, 1.1]) + def test_rejects_out_of_range(self, bad): + with pytest.raises(ValueError): + compute_free_surface_coefficients( + station_coords=mixed_geometry(), grid_coords=GRID, + station_coords_on_surface=np.full(4, bad)) + + +class TestMagnitudeEffect: + def test_amplification_lowers_threshold_by_expected_amount(self): + """Doubling the recorded amplitude halves the required moment, which is + 2/3*log10(2) ~ 0.2007 magnitude units.""" + expected = 2 / 3 * np.log10(2.0) + assert expected == pytest.approx(0.200686, abs=1e-6) + + def test_single_station_input_is_accepted(self): + fs_p, _ = compute_free_surface_coefficients( + station_coords=np.array([0.0, 0.0, 0.0]), + grid_coords=np.array([100.0, 0.0, 50.0]), fs_mode='auto') + assert fs_p.shape == (1, 1) + assert fs_p[0, 0] == pytest.approx(2.0) diff --git a/pytests/test_sensitivity_phases.py b/pytests/test_sensitivity_phases.py new file mode 100644 index 0000000..b9238de --- /dev/null +++ b/pytests/test_sensitivity_phases.py @@ -0,0 +1,172 @@ +"""End-to-end tests of phase-resolved magnitude sensitivity on a small grid. + +These pin the behaviour introduced with SV/SH support: the composite S mode keeps +the more detectable branch, joint detection takes the stricter branch, and the +radiation-pattern and free-surface options shift the result by known amounts. +""" +import numpy as np +import pandas as pd +import pytest + +from pynetdesign.modelling.homo import mag_sensitivity_grid +from pynetdesign.modelling.utils import ( + RAD_PATTERN_SV_DEFAULT, + generate_grid, +) + +VELOCITY = pd.DataFrame({ + 'Depth': [0], 'Vp': [2370.0], 'VpVsRatio': [2.0], 'Qp': [100.0], 'Qs': [100.0], + 'Rho': [2500.0], 'Ep': [0.0], 'Dt': [0.0], 'Gm': [0.0], +}) + + +def das_geometry(n=20, spacing=50.0, noise=1.95e-9, gauge_length=10.0): + """Vertical DAS cable at the origin, Z increasing downwards.""" + z = np.arange(n) * spacing + df = pd.DataFrame({ + 'Name': [f"C{i+1}" for i in range(n)], + 'X': np.zeros(n), 'Y': np.zeros(n), 'Z': z, + 'NoiseLevel': np.full(n, noise), + }) + df.attrs['noise_type'] = 'strain rate' + df.attrs['gauge_length'] = gauge_length + df.attrs['_is_combined'] = False + return df + + +def surface_geometry(noise=1e-10): + """Three 3C receivers on the free surface.""" + df = pd.DataFrame({ + 'Name': ['A', 'B', 'C'], + 'X': [0.0, 1000.0, 500.0], 'Y': [0.0, 0.0, 866.0], 'Z': [0.0, 0.0, 0.0], + 'NoiseLevel': np.full(3, noise), + }) + df.attrs['noise_type'] = 'velocity' + df.attrs['_is_combined'] = False + return df + + +def small_grid(): + return generate_grid(x=np.arange(200.0, 1200.0, 200.0), + y=np.array([0.0]), + z=np.arange(100.0, 900.0, 200.0)) + + +def run(geometry_df, wave_mode, grid=None, **kwargs): + grid = small_grid() if grid is None else grid + n = len(geometry_df) + amps = geometry_df['NoiseLevel'].to_numpy() * 2.0 + gauge = geometry_df.attrs.get('gauge_length') + if gauge is not None: + amps = amps * gauge + kwargs.setdefault('min_stations_p', 3) + kwargs.setdefault('min_stations_s', 3) + return mag_sensitivity_grid(grid_coords=grid, velocity_df=VELOCITY, + geometry_df=geometry_df, + min_amps_p=amps, min_amps_s=amps, + wave_mode=wave_mode, **kwargs) + + +class TestCompositeS: + def test_s_equals_sv_on_a_vertical_cable(self): + """A vertical cable cannot record SH, so the composite S mode must fall + back to SV everywhere.""" + geom = das_geometry() + s = run(geom, 'S', fs_mode='off') + sv = run(geom, 'SV', fs_mode='off') + np.testing.assert_allclose(s, sv, rtol=1e-12) + + def test_sh_alone_is_undetectable_on_a_vertical_cable(self): + sh = run(das_geometry(), 'SH', fs_mode='off') + assert np.all(np.isnan(sh)) + + def test_composite_s_is_never_worse_than_either_branch(self): + geom = das_geometry() + s = run(geom, 'S', fs_mode='off') + sv = run(geom, 'SV', fs_mode='off') + finite = np.isfinite(s) & np.isfinite(sv) + assert np.all(s[finite] <= sv[finite] + 1e-12) + + +class TestJointMode: + def test_ps_is_the_stricter_of_p_and_s(self): + geom = das_geometry() + p = run(geom, 'P', fs_mode='off') + s = run(geom, 'S', fs_mode='off') + ps = run(geom, 'PS', fs_mode='off') + np.testing.assert_allclose(ps, np.maximum(p, s), rtol=1e-12) + + +class TestRadiationPatternOverride: + def test_combined_s_override_shifts_by_a_known_amount(self): + """Passing 0.63, the combined S value, reproduces results computed with + that convention. Since M0 is inversely proportional to R, the shift is + 2/3*log10(0.63 / sqrt(7/30)).""" + geom = das_geometry() + default = run(geom, 'SV', fs_mode='off') + overridden = run(geom, 'SV', fs_mode='off', rad_pattern_s=0.63) + expected = -2 / 3 * np.log10(0.63 / RAD_PATTERN_SV_DEFAULT) + finite = np.isfinite(default) & np.isfinite(overridden) + np.testing.assert_allclose(overridden[finite] - default[finite], + expected, rtol=1e-10) + assert expected == pytest.approx(-0.0769, abs=1e-4) + + def test_p_override_shifts_p_only(self): + geom = das_geometry() + default = run(geom, 'P', fs_mode='off') + overridden = run(geom, 'P', fs_mode='off', rad_pattern_p=0.52) + finite = np.isfinite(default) + expected = -2 / 3 * np.log10(0.52 / np.sqrt(4 / 15)) + np.testing.assert_allclose(overridden[finite] - default[finite], + expected, rtol=1e-10) + assert expected == pytest.approx(-0.00201, abs=1e-5) + + def test_dict_override(self): + geom = das_geometry() + a = run(geom, 'SV', fs_mode='off', rad_pattern_s=0.63) + b = run(geom, 'SV', fs_mode='off', rad_patterns={'SV': 0.63}) + np.testing.assert_allclose(a, b, rtol=1e-12) + + +class TestFreeSurfaceEffect: + def test_surface_receivers_gain_exactly_two_thirds_log_two(self): + geom = surface_geometry() + off = run(geom, 'P', fs_mode='off') + on = run(geom, 'P', fs_mode='on') + expected = -2 / 3 * np.log10(2.0) + np.testing.assert_allclose(on - off, expected, rtol=1e-12) + + def test_auto_matches_on_for_an_all_surface_geometry(self): + geom = surface_geometry() + auto = run(geom, 'P', fs_mode='auto') + on = run(geom, 'P', fs_mode='on') + np.testing.assert_allclose(auto, on, rtol=1e-12) + + def test_auto_matches_off_for_a_downhole_geometry(self): + """No channel of this cable is at the surface, so 'auto' changes nothing.""" + geom = das_geometry(n=20, spacing=50.0) + geom['Z'] = geom['Z'] + 100.0 # move the whole cable below the surface + auto = run(geom, 'P', fs_mode='auto') + off = run(geom, 'P', fs_mode='off') + np.testing.assert_allclose(auto, off, rtol=1e-12) + + def test_surface_column_overrides_depth(self): + geom = surface_geometry() + geom['Surface'] = [0.0, 0.0, 0.0] + auto = run(geom, 'P', fs_mode='auto') + off = run(geom, 'P', fs_mode='off') + np.testing.assert_allclose(auto, off, rtol=1e-12) + + +class TestWaveModeValidation: + def test_rejects_unknown_mode(self): + with pytest.raises(ValueError): + run(das_geometry(), 'Q') + + def test_rejects_bad_fs_mode(self): + with pytest.raises(ValueError): + run(das_geometry(), 'P', fs_mode='maybe') + + def test_rejects_too_many_required_stations(self): + with pytest.raises(ValueError): + run(das_geometry(n=5), 'P', min_stations_p=99) diff --git a/pytests/test_wave_modes.py b/pytests/test_wave_modes.py new file mode 100644 index 0000000..eec927f --- /dev/null +++ b/pytests/test_wave_modes.py @@ -0,0 +1,125 @@ +"""Tests for wave-mode handling and radiation-pattern constants.""" +import numpy as np +import pytest + +from pynetdesign.modelling.utils import ( + RAD_PATTERN_P_DEFAULT, + RAD_PATTERN_SH_DEFAULT, + RAD_PATTERN_SV_DEFAULT, + S_WAVE_PHASES, + VALID_WAVE_MODES, + normalize_wave_mode, + phase_radiation_pattern, + resolve_radiation_pattern, + wave_mode_has_p, + wave_mode_has_s, + wave_mode_phases, + wave_mode_s_phases, +) + + +class TestNormalizeWaveMode: + @pytest.mark.parametrize("mode", VALID_WAVE_MODES) + def test_accepts_all_valid_modes(self, mode): + assert normalize_wave_mode(mode) == mode + + @pytest.mark.parametrize("given,expected", [ + ("p", "P"), (" ps ", "PS"), ("sv", "SV"), ("Sh", "SH"), ("s", "S"), + ]) + def test_case_and_whitespace_insensitive(self, given, expected): + assert normalize_wave_mode(given) == expected + + @pytest.mark.parametrize("mode", ["", "X", "PSV", "P S", None]) + def test_rejects_invalid(self, mode): + with pytest.raises(ValueError): + normalize_wave_mode(mode) + + +class TestWaveModeExpansion: + @pytest.mark.parametrize("mode,phases", [ + ("P", ("P",)), + ("SV", ("SV",)), + ("SH", ("SH",)), + ("S", ("SV", "SH")), + ("PS", ("P", "SV", "SH")), + ]) + def test_phases(self, mode, phases): + assert wave_mode_phases(mode) == phases + + @pytest.mark.parametrize("mode,has_p,has_s", [ + ("P", True, False), + ("SV", False, True), + ("SH", False, True), + ("S", False, True), + ("PS", True, True), + ]) + def test_branch_predicates(self, mode, has_p, has_s): + assert wave_mode_has_p(mode) is has_p + assert wave_mode_has_s(mode) is has_s + + def test_s_phases_subset(self): + assert wave_mode_s_phases("PS") == S_WAVE_PHASES + assert wave_mode_s_phases("P") == () + assert wave_mode_s_phases("SV") == ("SV",) + + +class TestRadiationPattern: + def test_rms_values(self): + """RMS radiation-pattern magnitudes over the focal sphere of a double couple.""" + assert phase_radiation_pattern("P") == pytest.approx(np.sqrt(4 / 15), rel=1e-12) + assert phase_radiation_pattern("SV") == pytest.approx(np.sqrt(7 / 30), rel=1e-12) + assert phase_radiation_pattern("SH") == pytest.approx(np.sqrt(1 / 6), rel=1e-12) + + def test_rounded_values(self): + assert phase_radiation_pattern("P") == pytest.approx(0.52, abs=5e-3) + assert phase_radiation_pattern("SV") == pytest.approx(0.48, abs=5e-3) + assert phase_radiation_pattern("SH") == pytest.approx(0.41, abs=5e-3) + + def test_sv_and_sh_combine_to_the_combined_s_value(self): + """ + = 2/5, so the combined S value is sqrt(2/5) ~ 0.63.""" + combined = np.sqrt(RAD_PATTERN_SV_DEFAULT**2 + RAD_PATTERN_SH_DEFAULT**2) + assert combined == pytest.approx(np.sqrt(2 / 5), rel=1e-12) + assert combined == pytest.approx(0.63, abs=3e-3) + + def test_sv_is_more_efficient_than_sh(self): + assert RAD_PATTERN_SV_DEFAULT > RAD_PATTERN_SH_DEFAULT + + @pytest.mark.parametrize("mode", ["S", "PS"]) + def test_rejects_composite_modes(self, mode): + with pytest.raises(ValueError): + phase_radiation_pattern(mode) + + +class TestResolveRadiationPattern: + def test_defaults_when_no_override(self): + for phase in ("P", "SV", "SH"): + assert resolve_radiation_pattern(phase) == phase_radiation_pattern(phase) + + def test_p_override(self): + assert resolve_radiation_pattern("P", rad_pattern_p=0.52) == pytest.approx(0.52) + # An S override must not leak into P + assert resolve_radiation_pattern("P", rad_pattern_s=0.63) == RAD_PATTERN_P_DEFAULT + + def test_s_override_applies_to_both_s_phases(self): + """Passing 0.63 reproduces results computed with the combined S value.""" + assert resolve_radiation_pattern("SV", rad_pattern_s=0.63) == pytest.approx(0.63) + assert resolve_radiation_pattern("SH", rad_pattern_s=0.63) == pytest.approx(0.63) + + def test_dict_takes_precedence(self): + value = resolve_radiation_pattern("SV", rad_pattern_s=0.63, + rad_patterns={"SV": 0.5}) + assert value == pytest.approx(0.5) + # Phases absent from the dict fall back to the scalar override + assert resolve_radiation_pattern("SH", rad_pattern_s=0.63, + rad_patterns={"SV": 0.5}) == pytest.approx(0.63) + + def test_dict_keys_are_normalized(self): + assert resolve_radiation_pattern("SV", rad_patterns={"sv": 0.4}) == pytest.approx(0.4) + + @pytest.mark.parametrize("bad", [0.0, -0.1, 1.5]) + def test_rejects_out_of_range(self, bad): + with pytest.raises(ValueError): + resolve_radiation_pattern("P", rad_pattern_p=bad) + + def test_unity_is_allowed(self): + assert resolve_radiation_pattern("P", rad_pattern_p=1.0) == pytest.approx(1.0) From da3dabbca2afa50f2c6e6543e54b3110ad13c46b Mon Sep 17 00:00:00 2001 From: Denis Anikiev Date: Fri, 11 Sep 2026 11:30:26 +0300 Subject: [PATCH 03/14] Add synthetic geometry builders and geometry writers Ports the geometry tooling from the reference implementation, so that the layouts used in the examples can be produced from code instead of shipped as hand-written tables, and so that the Surface and Components columns the modelling code already understands can actually be expressed in a file. Geometry builders, in the new pynetdesign.modelling.geometry module: - generate_geometry dispatches on mode to generate_borehole_geometry for vertical, inclined and deviated boreholes, generate_surface_geometry for line, zigzag, L, square and double-line layouts, and generate_darkfiber_geometry for a continuous path with prescribed turns - each returns geometry text lines that read_geometry accepts, with an optional depth-dependent noise profile for boreholes - generated output is byte-identical to the reference across all nine layout cases, checked by diffing 2234 lines of generated geometry Geometry I/O: - save_geometry writes either generated lines or a geometry DataFrame back in the read_geometry format, restoring Z to an elevation so the round trip is lossless - decimate_geometry keeps every n-th channel while preserving attrs, which is how channel spacing studies thin a dense array - read_geometry reads the optional Surface and Components columns and validates them: Surface must be a weight in [0, 1], Components must be '3C' or 'Z', and Components is rejected on a DAS geometry - read_geometry records the path it read in attrs['file_path'] Fixes the default dark-fibre turn distances, which were fixed at 6, 12 and 18 km and therefore exceeded the default 5 km cable, making the default call raise. They are now one quarter, one half and three quarters of the cable, so the default works at any length and is unchanged for a 24 km cable. Adds 39 tests covering the builders, the borehole geometry against closed forms for vertical and inclined cables, save and read round trips, the optional columns and their validation, decimation, and a generated array driving a sensitivity computation end to end. The reference parity established for the magnitude core is unaffected, and all ten examples still run. --- CHANGELOG.md | 6 + pynetdesign/modelling/__init__.py | 2 + pynetdesign/modelling/geometry.py | 548 ++++++++++++++++++++++++++++++ pynetdesign/modelling/io.py | 260 +++++++++++++- pytests/test_geometry_io.py | 264 ++++++++++++++ 5 files changed, 1074 insertions(+), 6 deletions(-) create mode 100644 pynetdesign/modelling/geometry.py create mode 100644 pytests/test_geometry_io.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 830eef3..9ca0b27 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -16,6 +16,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - projection of single vertical-component stations onto the vertical axis, driven by a `Components` column; three-component stations record the full vector and are unaffected - per-station free surface amplification through `compute_free_surface_coefficients`, with `fs_mode`, `fs_level` and `fs_deviation` parameters and an optional `Surface` column carrying per-station weights - helpers `normalize_wave_mode`, `wave_mode_phases`, `wave_mode_has_p`, `wave_mode_has_s`, `wave_mode_s_phases`, `resolve_radiation_pattern`, `check_min_stations` and `geometry_requires_receiver_projection` +- `pynetdesign.modelling.geometry`, a module of synthetic geometry builders: `generate_geometry` dispatches to `generate_borehole_geometry` for vertical, inclined and deviated boreholes, to `generate_surface_geometry` for line, zigzag, L, square and double-line layouts, and to `generate_darkfiber_geometry` for a continuous dark-fibre path with turns. Each returns geometry text lines that `read_geometry` can read back, including a depth-dependent noise profile for boreholes +- `save_geometry`, writing either generated text lines or a geometry DataFrame back to the `read_geometry` format, with `Z` restored to an elevation so that a round trip is lossless +- `decimate_geometry`, keeping every n-th channel while preserving the metadata, for studying the effect of channel spacing on a dense array +- `read_geometry` reads the optional `Surface` and `Components` columns, so the per-station free surface weights and the single vertical component stations introduced above can be described in a geometry file +- `read_geometry` records the path it read in `attrs['file_path']` - test modules covering wave modes and radiation patterns, receiver projection, the free surface correction, network detectability, and phase-resolved sensitivity end to end ### Changed @@ -33,6 +38,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Fixed - the station-count guard in `mag_detectable` compared against the grid axis instead of the station axis, so a request for more receivers than the geometry has was not reported clearly +- the default dark-fibre turn distances are now fractions of the cable length, at one quarter, one half and three quarters, so that `generate_geometry(mode='darkfiber')` works for any cable. The previous fixed distances of 6, 12 and 18 km exceeded the default cable length and made the default call fail; for a 24 km cable the layout is unchanged ## [v1.0.1] - 2026-09-10 diff --git a/pynetdesign/modelling/__init__.py b/pynetdesign/modelling/__init__.py index 7ca235c..baac255 100644 --- a/pynetdesign/modelling/__init__.py +++ b/pynetdesign/modelling/__init__.py @@ -1,4 +1,6 @@ from .core import * +from . import geometry +from .geometry import * from .homo import * from .io import * from .utils import * diff --git a/pynetdesign/modelling/geometry.py b/pynetdesign/modelling/geometry.py new file mode 100644 index 0000000..157a3b2 --- /dev/null +++ b/pynetdesign/modelling/geometry.py @@ -0,0 +1,548 @@ +"""Synthetic DAS geometry builders used by examples and tests. + +These helpers generate simple borehole, surface, and dark-fiber layouts in +the tab-delimited geometry format consumed by :mod:`pynetdesign.modelling.io`. +""" + +import numpy as np + +__all__ = [ + 'generate_geometry', + 'generate_borehole_geometry', + 'generate_surface_geometry', + 'generate_darkfiber_geometry', +] + + +def _format_header_value(value: float) -> str: + """Format numeric header values without forcing integer rounding.""" + value = float(value) + if value.is_integer(): + return str(int(value)) + return f"{value:g}" + + +def generate_geometry( + mode: str, + # Common parameters + center_x: float = 0.0, + center_y: float = 0.0, + spacing: float = 10.0, + noise_level: float = 5.0e-5, + gauge_length: float = 10.0, + cable_length: float = 5000.0, + # Borehole parameters + start_depth: float = 0.0, + azimuth: float = 0.0, + dip: float = 90.0, + surface_noise: float = None, + noise_decay_depth: float = 150.0, + noise_low_depth: float = 200.0, + # Surface / dark-fiber common + depth: float = 0.0, + # Surface zigzag + surface_shape: str = 'zigzag', + zigzag_amplitude: float = 500.0, + zigzag_segment_length: float = 1000.0, + # Surface L-shape + l_arm_x: float = 3000.0, + l_arm_y: float = 3000.0, + # Surface square / rectangle + rect_width: float = 3000.0, + rect_height: float = 3000.0, + # Surface double lines + double_line_separation: float = 500.0, + # Dark-fiber parameters + start: dict = None, + turns: dict = None, + behavior: dict = None, + seed: int = 42, +): + """Generate a DAS geometry file in PyNetDesign tab-delimited format.""" + if mode == 'borehole': + coords = generate_borehole_geometry( + center_x, + center_y, + spacing, + start_depth, + cable_length, + azimuth, + dip, + ) + elif mode == 'surface': + coords = generate_surface_geometry( + center_x, + center_y, + spacing, + depth, + cable_length, + surface_shape, + zigzag_amplitude, + zigzag_segment_length, + l_arm_x, + l_arm_y, + rect_width, + rect_height, + double_line_separation, + ) + elif mode == 'darkfiber': + coords = generate_darkfiber_geometry( + center_x, + center_y, + spacing, + depth, + cable_length, + start, + turns, + behavior, + seed, + ) + else: + raise ValueError( + f"Unknown mode: {mode}. Use 'borehole', 'surface', or 'darkfiber'." + ) + + header = ( + f"Northing(m)\tEasting(m)\tElevation(m)\t" + f"NoiseLevel(1/s)\tGaugeLength={_format_header_value(gauge_length)}(m)" + ) + lines = [header] + + for x, y, z in coords: + current_noise = noise_level + if mode == 'borehole' and surface_noise is not None: + depth_m = -z + if depth_m < noise_decay_depth: + current_noise = surface_noise + elif depth_m < noise_low_depth: + log_n0 = np.log(surface_noise) + log_nf = np.log(noise_level) + t = ( + (depth_m - noise_decay_depth) + / (noise_low_depth - noise_decay_depth) + ) + log_n = log_n0 + (log_nf - log_n0) * t + current_noise = np.exp(log_n) + lines.append(f"{y:.2f}\t{x:.2f}\t{z:.2f}\t{current_noise:.6e}") + + return lines + + +def generate_borehole_geometry( + center_x: float, + center_y: float, + spacing: float, + start_depth: float, + cable_length: float, + azimuth: float, + dip: float, +): + """Generate borehole DAS channel coordinates.""" + measured_depths = start_depth + _build_sample_distances(cable_length, spacing) + + dip_rad = np.radians(dip) + az_rad = np.radians(azimuth) + + coords = [] + for measured_depth in measured_depths: + true_vertical_depth = measured_depth * np.sin(dip_rad) + horizontal_offset = measured_depth * np.cos(dip_rad) + dx = horizontal_offset * np.sin(az_rad) + dy = horizontal_offset * np.cos(az_rad) + x = center_x + dx + y = center_y + dy + coords.append((x, y, -true_vertical_depth)) + + return coords + + +def generate_surface_geometry( + center_x: float, + center_y: float, + spacing: float, + depth: float, + cable_length: float, + shape: str, + zigzag_amplitude: float, + zigzag_segment_length: float, + l_arm_x: float, + l_arm_y: float, + rect_width: float, + rect_height: float, + double_line_separation: float, +): + """Generate surface DAS channel coordinates.""" + if shape == 'line': + waypoints = _surface_line(center_x, center_y, cable_length) + elif shape == 'zigzag': + waypoints = _surface_zigzag( + center_x, + center_y, + cable_length, + zigzag_amplitude, + zigzag_segment_length, + ) + elif shape == 'L': + waypoints = _surface_l(center_x, center_y, l_arm_x, l_arm_y) + elif shape == 'square': + waypoints = _surface_square(center_x, center_y, rect_width, rect_height) + elif shape == 'double_line': + waypoints = _surface_double_line( + center_x, + center_y, + cable_length, + double_line_separation, + ) + else: + raise ValueError(f"Unknown surface_shape: {shape}") + + return _interpolate_polyline(waypoints, spacing, depth) + + +def generate_darkfiber_geometry( + center_x: float, + center_y: float, + spacing: float, + depth: float, + cable_length: float, + start: dict = None, + turns: dict = None, + behavior: dict = None, + seed: int = 42, +): + """Generate a continuous dark-fiber cable path.""" + rng = np.random.RandomState(seed) + start_segment = _normalize_darkfiber_start(start) + + if turns is None: + # Turn at one quarter, one half and three quarters of the cable, so that the + # default layout works for any cable length. For a 24 km cable these are the + # 6, 12 and 18 km turns used historically. + turns = {0.25 * cable_length: 0.0, + 0.50 * cable_length: 90.0, + 0.75 * cable_length: 0.0} + if not isinstance(turns, dict): + raise ValueError("turns must be a dictionary for darkfiber mode.") + + turn_distances = sorted(turns.keys()) + current_dist = 0.0 + for turn_distance in turn_distances: + if turn_distance <= current_dist or turn_distance >= cable_length: + raise ValueError( + "turn distances must be strictly increasing and lie within " + "cable_length." + ) + current_dist = turn_distance + + segment_behavior = _normalize_darkfiber_behavior( + behavior, + turn_distances, + start_segment, + ) + + seg_lengths = [] + seg_azimuths = [] + seg_defs = [start_segment] + + current_dist = 0.0 + current_azimuth = start_segment['azimuth'] + for turn_distance in turn_distances: + seg_lengths.append(turn_distance - current_dist) + seg_azimuths.append(current_azimuth) + current_dist = turn_distance + current_azimuth = turns[turn_distance] + seg_defs.append(segment_behavior[turn_distance]) + + if current_dist < cable_length: + seg_lengths.append(cable_length - current_dist) + seg_azimuths.append(current_azimuth) + + all_coords = [] + current_xy = np.array([0.0, 0.0], dtype=float) + for index, (segment_length, segment_azimuth, segment_def) in enumerate( + zip(seg_lengths, seg_azimuths, seg_defs) + ): + seg_coords = _generate_darkfiber_segment( + start_xy=current_xy, + length=segment_length, + azimuth=segment_azimuth, + depth=depth, + spacing=spacing, + segment=segment_def, + rng=rng, + ) + + if index > 0: + seg_coords = seg_coords[1:] + all_coords.extend(seg_coords) + + if seg_coords: + current_xy = np.array(seg_coords[-1][:2], dtype=float) + + if not all_coords: + return [] + + all_coords = np.asarray(all_coords, dtype=float) + if start_segment['start_x'] is not None: + all_coords[:, 0] += start_segment['start_x'] + all_coords[:, 1] += start_segment['start_y'] + else: + min_x, max_x = np.min(all_coords[:, 0]), np.max(all_coords[:, 0]) + min_y, max_y = np.min(all_coords[:, 1]), np.max(all_coords[:, 1]) + all_coords[:, 0] += center_x - (min_x + max_x) / 2 + all_coords[:, 1] += center_y - (min_y + max_y) / 2 + + return _resample_polyline_coords(all_coords, spacing) + + +def _build_sample_distances(total_length, spacing): + """Return cumulative distances sampled at the requested spacing.""" + if spacing <= 0: + raise ValueError("spacing must be positive") + + total_length = float(total_length) + n_full_steps = int(np.floor(total_length / spacing)) + sample_s = spacing * np.arange(n_full_steps + 1, dtype=float) + + if sample_s.size == 0: + sample_s = np.array([0.0], dtype=float) + + if not np.isclose(sample_s[-1], total_length): + sample_s = np.append(sample_s, total_length) + + return sample_s + + +def _resample_polyline_coords(coords, spacing): + """Resample a 3D polyline so channel spacing follows the cable path.""" + coords = np.asarray(coords, dtype=float) + + if coords.ndim != 2 or coords.shape[1] != 3: + raise ValueError("coords must have shape (n_points, 3)") + if coords.shape[0] == 0: + return [] + if coords.shape[0] == 1: + return [tuple(coords[0])] + + seg_vecs = np.diff(coords, axis=0) + seg_lens = np.linalg.norm(seg_vecs, axis=1) + keep = np.concatenate(([True], seg_lens > 0.0)) + coords = coords[keep] + + if coords.shape[0] == 1: + return [tuple(coords[0])] + + seg_lens = np.linalg.norm(np.diff(coords, axis=0), axis=1) + cum_len = np.concatenate(([0.0], np.cumsum(seg_lens))) + sample_s = _build_sample_distances(cum_len[-1], spacing) + + x = np.interp(sample_s, cum_len, coords[:, 0]) + y = np.interp(sample_s, cum_len, coords[:, 1]) + z = np.interp(sample_s, cum_len, coords[:, 2]) + + return list(zip(x, y, z)) + + +def _surface_line(center_x, center_y, length): + """Straight line along the X axis.""" + return [ + (center_x - length / 2, center_y), + (center_x + length / 2, center_y), + ] + + +def _surface_zigzag(center_x, center_y, length, amplitude, segment_length): + """Zigzag path along the X axis with alternating Y offsets.""" + x_start = center_x - length / 2 + n_seg = int(length / segment_length) + points = [] + for idx in range(n_seg + 1): + x = x_start + idx * segment_length + y = center_y + amplitude * (1 if idx % 2 == 0 else -1) + points.append((x, y)) + return points + + +def _surface_l(center_x, center_y, arm_x, arm_y): + """L-shaped path: horizontal arm followed by vertical arm.""" + return [ + (center_x - arm_x / 2, center_y), + (center_x + arm_x / 2, center_y), + (center_x + arm_x / 2, center_y + arm_y), + ] + + +def _surface_square(center_x, center_y, width, height): + """Rectangular perimeter.""" + return [ + (center_x - width / 2, center_y - height / 2), + (center_x + width / 2, center_y - height / 2), + (center_x + width / 2, center_y + height / 2), + (center_x - width / 2, center_y + height / 2), + (center_x - width / 2, center_y - height / 2), + ] + + +def _surface_double_line(center_x, center_y, length, separation): + """Two parallel lines connected at one end.""" + y1 = center_y - separation / 2 + y2 = center_y + separation / 2 + return [ + (center_x - length / 2, y1), + (center_x + length / 2, y1), + (center_x + length / 2, y2), + (center_x - length / 2, y2), + ] + + +def _interpolate_polyline(waypoints, spacing, depth): + """Interpolate (x, y) waypoints at uniform spacing along the path.""" + elevation = -depth + coords = [(x, y, elevation) for x, y in waypoints] + return _resample_polyline_coords(coords, spacing) + + +def _normalize_darkfiber_start(start): + """Validate and normalize the first dark-fiber segment definition.""" + defaults = { + 'start_x': None, + 'start_y': None, + 'azimuth': 90.0, + 'behavior': 'sinusoid', + 'sinusoid_amplitude': 500.0, + 'sinusoid_period': 2000.0, + 'amp_variation_amplitude': 200.0, + 'amp_variation_period': 3000.0, + } + + if start is None: + start = {} + if not isinstance(start, dict): + raise ValueError("start must be a dictionary for darkfiber mode.") + + segment = defaults.copy() + segment.update(start) + + if segment['behavior'] not in ('straight', 'sinusoid'): + raise ValueError( + "start['behavior'] must be either 'straight' or 'sinusoid'." + ) + + if (segment['start_x'] is None) != (segment['start_y'] is None): + raise ValueError( + "start['start_x'] and start['start_y'] must be provided together." + ) + + if segment['behavior'] == 'sinusoid': + for key in ( + 'sinusoid_amplitude', + 'sinusoid_period', + 'amp_variation_amplitude', + 'amp_variation_period', + ): + if segment[key] is None: + raise ValueError( + f"start['{key}'] is required for sinusoid behavior." + ) + + return segment + + +def _normalize_darkfiber_behavior(behavior, turn_distances, start_segment): + """Validate and normalize per-segment behavior definitions.""" + if behavior is None: + behavior = {td: start_segment['behavior'] for td in turn_distances} + if not isinstance(behavior, dict): + raise ValueError("behavior must be a dictionary for darkfiber mode.") + + if set(behavior.keys()) != set(turn_distances): + raise ValueError( + "behavior keys must match turns keys exactly for darkfiber mode." + ) + + normalized = {} + last_sinusoid = { + 'sinusoid_amplitude': start_segment['sinusoid_amplitude'], + 'sinusoid_period': start_segment['sinusoid_period'], + 'amp_variation_amplitude': start_segment['amp_variation_amplitude'], + 'amp_variation_period': start_segment['amp_variation_period'], + } + + for turn_distance in turn_distances: + spec = behavior[turn_distance] + if isinstance(spec, str): + spec = {'behavior': spec} + elif isinstance(spec, dict): + spec = spec.copy() + else: + raise ValueError( + "Each behavior entry must be either a string or a dictionary." + ) + + segment_behavior = spec.get('behavior') + if segment_behavior not in ('straight', 'sinusoid'): + raise ValueError( + f"behavior[{turn_distance}] must define 'straight' or " + f"'sinusoid' behavior." + ) + + segment = {'behavior': segment_behavior} + if segment_behavior == 'sinusoid': + for key in last_sinusoid: + segment[key] = spec.get(key, last_sinusoid[key]) + last_sinusoid = {key: segment[key] for key in last_sinusoid} + + normalized[turn_distance] = segment + + return normalized + + +def _generate_darkfiber_segment(start_xy, length, azimuth, depth, spacing, + segment, rng): + """Generate one dark-fiber segment before final spacing resampling.""" + dense_spacing = min(max(spacing / 5.0, 0.5), 2.0) + sample_s = _build_sample_distances(length, dense_spacing) + + az_rad = np.radians(azimuth) + tangent = np.array([np.sin(az_rad), np.cos(az_rad)]) + perpendicular = np.array([-tangent[1], tangent[0]]) + + base_x = start_xy[0] + sample_s * tangent[0] + base_y = start_xy[1] + sample_s * tangent[1] + wobble = np.zeros_like(sample_s) + + if segment['behavior'] == 'sinusoid': + amp_phase = rng.uniform(0, 2 * np.pi) + local_amp = ( + segment['sinusoid_amplitude'] + + segment['amp_variation_amplitude'] + * np.sin( + 2 * np.pi * sample_s / segment['amp_variation_period'] + + amp_phase + ) + ) + wobble = local_amp * np.sin( + 2 * np.pi * sample_s / segment['sinusoid_period'] + ) + + taper_radius = min(segment['sinusoid_period'] / 2.0, length / 2.0) + if taper_radius > 0: + taper = np.ones_like(sample_s) + + start_mask = sample_s < taper_radius + taper[start_mask] *= 0.5 * ( + 1 - np.cos(np.pi * sample_s[start_mask] / taper_radius) + ) + + end_dist = length - sample_s + end_mask = end_dist < taper_radius + taper[end_mask] *= 0.5 * ( + 1 - np.cos(np.pi * end_dist[end_mask] / taper_radius) + ) + + wobble *= taper + + x = base_x + wobble * perpendicular[0] + y = base_y + wobble * perpendicular[1] + z = np.full_like(sample_s, -depth) + return list(zip(x, y, z)) diff --git a/pynetdesign/modelling/io.py b/pynetdesign/modelling/io.py index 8d99525..8340ae5 100644 --- a/pynetdesign/modelling/io.py +++ b/pynetdesign/modelling/io.py @@ -2,6 +2,7 @@ import pandas as pd import re import copy +from os import PathLike from typing import List, Dict, Union def combine_geometry(*dataframes: pd.DataFrame, names: Union[List[str], Dict[int, str]] = None) -> pd.DataFrame: @@ -207,7 +208,8 @@ def read_geometry(file_path: str) -> pd.DataFrame: This function reads a CSV (ASCII) file containing geometry information, ensuring that coordinate columns (Latitude/Northing, Longitude/Easting, Elevation) and - NoiseLevel column are present. The 'Name' column is optional. + NoiseLevel column are present. The 'Name', 'Surface' and 'Components' columns + are optional. The resulting DataFrame has columns renamed and units converted as necessary. If header has 'GaugeLength=', the appropriate gauge length is read and stored in the attributes. @@ -226,9 +228,12 @@ def read_geometry(file_path: str) -> pd.DataFrame: - 'X' for Longitude or Easting (converted to meters) - 'Z' for Depth (negative of Elevation, converted to meters) - 'NoiseLevel' for Noise Level + - 'Surface' for the free surface indicator (if present in the input file) + - 'Components' for the recorded station components (if present in the input file) The DataFrame also includes metadata: + - ``df.attrs['file_path']``: Path of the file that was read - ``df.attrs['original_coord_units']``: Original units of coordinate columns - ``df.attrs['noise_type']``: Type of noise measurement - ``df.attrs['noise_unit']``: Original unit of noise measurement @@ -250,8 +255,14 @@ def read_geometry(file_path: str) -> pd.DataFrame: m/s^2 (acceleration), and 1/s (strain). - The 'Z' column is converted to represent depth (negative of elevation). - The 'Name' column, if present, is interpreted as strings. - - Only known columns ('Name', 'Latitude', 'Northing', 'Longitude', 'Easting', - 'Elevation', 'NoiseLevel') or combinations like 'Latitude/Northing' or + - The 'Surface' column, if present, is a weight in [0, 1] marking receivers at the + free surface, where a fractional value blends between no amplification and full + free surface amplification. + - The 'Components' column, if present, must be '3C' or 'Z' and is only allowed for + station geometries, not for DAS geometries. + - Only known columns ('Name', 'Latitude', 'Northing', 'Longitude', 'Easting', + 'Elevation', 'NoiseLevel', 'Surface', 'Components') or combinations like + 'Latitude/Northing' or 'Latitude(WGS84)/Northing(m)' are read from the file. - For Latitude and Longitude conversion to m is not implemented, inputs are treated as meters - The function will raise an error if any unknown columns are present in the input file. @@ -322,7 +333,7 @@ def extract_unit(col_name): # Check for required columns, known and unknown columns required_columns = {'Y': False, 'X': False, 'Z': False, 'NoiseLevel': False} # Define known columns - known_columns = ['Name', 'Latitude', 'Northing', 'Longitude', 'Easting', 'Elevation', 'NoiseLevel'] + known_columns = ['Name', 'Latitude', 'Northing', 'Longitude', 'Easting', 'Elevation', 'NoiseLevel', 'Surface', 'Components'] # Define unknown columns unknown_columns = [] # Default lat/lon flags @@ -370,7 +381,7 @@ def extract_unit(col_name): raise ValueError(f"Unknown columns found: {', '.join(unknown_columns)}") # Define data types for each column (float for all except 'Name' which is str) - dtype_mapping = {col: str if 'Name' in col else float for col in columns_to_read} + dtype_mapping = {col: str if ('Name' in col or 'Components' in col) else float for col in columns_to_read} # Read the CSV file, using only the known columns df = pd.read_csv(file_path, delimiter=r"\s+", dtype=dtype_mapping, usecols=columns_to_read) @@ -388,6 +399,10 @@ def extract_unit(col_name): column_mapping[col] = 'NoiseLevel' elif 'Name' in base_name: column_mapping[col] = 'Name' + elif 'Surface' in base_name: + column_mapping[col] = 'Surface' + elif 'Components' in base_name: + column_mapping[col] = 'Components' # Rename the columns df.rename(columns=column_mapping, inplace=True) @@ -448,8 +463,21 @@ def extract_unit(col_name): noise_unit = 'ε' else: raise ValueError(f"Unrecognized unit for NoiseLevel: {noise_unit}") - + + # Check that the surface indicator is a weight between 0 and 1 + if 'Surface' in df.columns: + if not df['Surface'].between(0, 1).all(): + raise ValueError("Surface values must be in the range of 0 to 1") + + # Check the recorded components of station geometries + if 'Components' in df.columns: + if gauge_length is not None and df['Components'].notna().any(): + raise ValueError("Components column is only allowed for station geometries, not DAS geometries.") + _validate_geometry_components(df['Components']) + # Edit DataFrame attributes: + # Store the path of the file that was read + df.attrs['file_path'] = file_path # Store original coordinate units in DataFrame attributes df.attrs['original_coord_units'] = list(coord_units)[0] if coord_units else None @@ -472,6 +500,226 @@ def extract_unit(col_name): return df +def _format_geometry_header_value(value: float) -> str: + r"""Format header values without introducing unnecessary decimals.""" + value = float(value) + if value.is_integer(): + return str(int(value)) + return f"{value:g}" + + +def _get_noise_header_unit(geometry_df: pd.DataFrame) -> str: + r""" + Infer a writable noise unit for a geometry DataFrame. + + Parameters + ---------- + geometry_df : :obj:`pandas.DataFrame` + Geometry DataFrame carrying ``attrs['noise_type']`` or ``attrs['noise_unit']``. + + Returns + ------- + noise_unit : :obj:`str` + Unit string to write into the header. + + Raises + ------ + ValueError + If neither the noise type nor the noise unit is available. + """ + noise_type = geometry_df.attrs.get('noise_type') + if noise_type == 'displacement': + return 'm' + if noise_type == 'velocity': + return 'm/s' + if noise_type == 'acceleration': + return 'm/s^2' + if noise_type == 'strain rate': + return '1/s' + if noise_type == 'strain': + return 'strain' + + noise_unit = geometry_df.attrs.get('noise_unit') + if noise_unit is None: + raise ValueError( + "Geometry DataFrame must define either attrs['noise_type'] or " + "attrs['noise_unit'] to be written." + ) + return noise_unit + + +def _validate_geometry_components(components): + r"""Validate that recorded components are either ``'3C'`` or ``'Z'``.""" + component_values = pd.Series(components).dropna() + invalid_components = sorted(set(component_values).difference({'3C', 'Z'})) + if invalid_components: + invalid_text = ', '.join(str(value) for value in invalid_components) + raise ValueError( + "Components values must be either '3C' or 'Z'. " + f"Invalid values: {invalid_text}" + ) + + +def _geometry_dataframe_to_lines(geometry_df: pd.DataFrame) -> List[str]: + r""" + Convert a geometry DataFrame into lines that :func:`read_geometry` can read back. + + Parameters + ---------- + geometry_df : :obj:`pandas.DataFrame` + Non-combined geometry DataFrame. + + Returns + ------- + lines : :obj:`list` of :obj:`str` + Header line followed by one line per station. + + Raises + ------ + ValueError + If the geometry is combined, required columns are missing, or a + ``Components`` column is present on a DAS geometry. + + Notes + ----- + ``Z`` is written back as an elevation, that is with the opposite sign, so that a + round trip through :func:`read_geometry` reproduces the original coordinates. + """ + if is_combined_geometry(geometry_df): + raise ValueError("Geometry must not be combined to save it to a single file.") + + required_columns = {'Y', 'X', 'Z', 'NoiseLevel'} + missing_columns = required_columns.difference(geometry_df.columns) + if missing_columns: + missing = ', '.join(sorted(missing_columns)) + raise ValueError(f"Geometry DataFrame is missing required columns: {missing}") + + export_columns = {} + if 'Name' in geometry_df.columns and not geometry_df['Name'].isna().all(): + export_columns['Name'] = geometry_df['Name'] + + export_columns['Northing(m)'] = geometry_df['Y'] + export_columns['Easting(m)'] = geometry_df['X'] + export_columns['Elevation(m)'] = -geometry_df['Z'] + export_columns[f"NoiseLevel({_get_noise_header_unit(geometry_df)})"] = ( + geometry_df['NoiseLevel'] + ) + + if 'Surface' in geometry_df.columns and not geometry_df['Surface'].isna().all(): + export_columns['Surface'] = geometry_df['Surface'] + + if 'Components' in geometry_df.columns and not geometry_df['Components'].isna().all(): + if geometry_df.attrs.get('gauge_length') is not None: + raise ValueError("Components column is only allowed for station geometries, not DAS geometries.") + _validate_geometry_components(geometry_df['Components']) + export_columns['Components'] = geometry_df['Components'] + + export_df = pd.DataFrame(export_columns, index=geometry_df.index) + + header = '\t'.join(export_df.columns.tolist()) + gauge_length = geometry_df.attrs.get('gauge_length') + if gauge_length is not None: + header += f"\tGaugeLength={_format_geometry_header_value(gauge_length)}(m)" + + if export_df.empty: + return [header] + + data_text = export_df.to_csv( + sep='\t', + index=False, + header=False, + lineterminator='\n', + float_format='%.12g', + ) + return [header, *data_text.rstrip('\n').splitlines()] + + +def save_geometry(geometry: Union[List[str], pd.DataFrame], file_path: Union[str, PathLike]): + r""" + Save geometry to disk in a format compatible with :func:`read_geometry`. + + Accepts either the geometry text lines returned by the generators in + :mod:`pynetdesign.modelling.geometry` or a geometry :obj:`pandas.DataFrame` + previously read by :func:`read_geometry`. + + Parameters + ---------- + geometry : :obj:`list` of :obj:`str` or :obj:`pandas.DataFrame` + Geometry content to write. + file_path : :obj:`str` or path-like + Output file path. + + Returns + ------- + file_path : :obj:`str` or path-like + The output path, returned for convenience. + + Raises + ------ + ValueError + If the geometry content is empty or invalid for writing. + TypeError + If ``geometry`` is neither a list of strings nor a geometry DataFrame. + """ + if isinstance(geometry, pd.DataFrame): + geometry_lines = _geometry_dataframe_to_lines(geometry) + elif isinstance(geometry, list): + geometry_lines = geometry + if not geometry_lines: + raise ValueError("geometry lines must not be empty.") + if not all(isinstance(line, str) for line in geometry_lines): + raise TypeError("geometry lines must contain only strings.") + else: + raise TypeError( + "geometry must be either a list of strings or a pandas DataFrame." + ) + + with open(file_path, 'w') as f: + f.write('\n'.join(geometry_lines)) + + return file_path + + +def decimate_geometry(geometry_df: pd.DataFrame, step: int) -> pd.DataFrame: + r""" + Decimate a non-combined geometry by keeping every ``step``-th station. + + This is useful for checking how the channel spacing of a dense DAS array + affects sensitivity, by thinning an existing geometry rather than regenerating it. + + Parameters + ---------- + geometry_df : :obj:`pandas.DataFrame` + The input geometry, with one row per station and metadata in ``attrs``. + step : :obj:`int` + Decimation factor, so that ``step=5`` keeps rows 0, 5, 10 and so on. + + Returns + ------- + decimated_geometry_df : :obj:`pandas.DataFrame` + New geometry with every ``step``-th station, the index reset and all + original ``attrs`` preserved. + + Raises + ------ + ValueError + If ``step`` is not a positive integer, or if the geometry is combined. + """ + if not isinstance(step, int) or isinstance(step, bool) or step < 1: + raise ValueError(f"Decimation step must be a positive integer, got {step!r}") + + if is_combined_geometry(geometry_df): + raise ValueError("The input DataFrame is a combined one. Decimation is not supported for combined DataFrames.") + + # Select every step-th row + decimated_geometry_df = geometry_df.iloc[::step].reset_index(drop=True).copy() + + # Preserve metadata + decimated_geometry_df.attrs = geometry_df.attrs.copy() + + return decimated_geometry_df + + def read_velocity(file_path: str): r""" Read velocity model data from a whitespace-delimited CSV (ASCII) file and return a DataFrame. diff --git a/pytests/test_geometry_io.py b/pytests/test_geometry_io.py new file mode 100644 index 0000000..8836650 --- /dev/null +++ b/pytests/test_geometry_io.py @@ -0,0 +1,264 @@ +"""Tests for the synthetic geometry generators and the geometry writers.""" +import numpy as np +import pandas as pd +import pytest + +import pynetdesign.modelling as pndmod +from pynetdesign.modelling.geometry import ( + generate_borehole_geometry, + generate_geometry, +) +from pynetdesign.modelling.io import ( + combine_geometry, + decimate_geometry, + read_geometry, + save_geometry, +) + +STATION_HEADER = "Name\tNorthing(m)\tEasting(m)\tElevation(m)\tNoiseLevel(m/s)" + + +@pytest.fixture +def written(tmp_path): + """Write geometry lines to a temporary file and read them back.""" + def _written(lines, name="geometry.txt"): + path = tmp_path / name + save_geometry(lines, path) + return read_geometry(str(path)) + return _written + + +class TestGenerators: + @pytest.mark.parametrize("shape", ["line", "zigzag", "L", "square", "double_line"]) + def test_surface_shapes_are_readable(self, shape, written): + lines = generate_geometry(mode="surface", cable_length=3000.0, spacing=50.0, + surface_shape=shape) + df = written(lines, f"surface_{shape}.txt") + assert len(df) > 1 + assert {'X', 'Y', 'Z', 'NoiseLevel'}.issubset(df.columns) + + def test_borehole_channel_count_and_spacing(self, written): + lines = generate_geometry(mode="borehole", cable_length=999.0, spacing=1.0, + noise_level=1.95e-9, gauge_length=10.0) + df = written(lines, "borehole.txt") + assert len(df) == 1000 + assert df.attrs['gauge_length'] == pytest.approx(10.0) + assert df.attrs['noise_type'] == 'strain rate' + np.testing.assert_allclose(np.diff(df['Z']), 1.0) + + def test_vertical_borehole_is_a_straight_column(self): + coords = np.array(generate_borehole_geometry( + center_x=100.0, center_y=-50.0, spacing=25.0, start_depth=0.0, + cable_length=500.0, azimuth=0.0, dip=90.0)) + assert np.allclose(coords[:, 0], 100.0) + assert np.allclose(coords[:, 1], -50.0) + assert coords[0, 2] == pytest.approx(0.0) + assert coords[-1, 2] == pytest.approx(-500.0) + + def test_inclined_borehole_steps_out(self): + """A 60 degree dip advances horizontally by cos(60) per unit of cable.""" + coords = np.array(generate_borehole_geometry( + center_x=0.0, center_y=0.0, spacing=100.0, start_depth=0.0, + cable_length=1000.0, azimuth=90.0, dip=60.0)) + measured = 1000.0 + assert coords[-1, 0] == pytest.approx(measured * np.cos(np.radians(60.0))) + assert coords[-1, 1] == pytest.approx(0.0, abs=1e-9) + assert coords[-1, 2] == pytest.approx(-measured * np.sin(np.radians(60.0))) + + def test_darkfiber_defaults_work_for_any_cable_length(self): + """The default turns are fractions of the cable, so the default layout is + valid whatever length is requested.""" + for length in (4000.0, 5000.0, 24000.0): + lines = generate_geometry(mode="darkfiber", cable_length=length, + spacing=100.0) + assert len(lines) > 2 + + def test_darkfiber_rejects_turns_beyond_the_cable(self): + with pytest.raises(ValueError): + generate_geometry(mode="darkfiber", cable_length=1000.0, spacing=50.0, + turns={5000.0: 90.0}) + + def test_rejects_unknown_mode(self): + with pytest.raises(ValueError): + generate_geometry(mode="spiral", cable_length=1000.0) + + def test_rejects_unknown_surface_shape(self): + with pytest.raises(ValueError): + generate_geometry(mode="surface", cable_length=1000.0, + surface_shape="triangle") + + +class TestSaveGeometry: + def test_writes_generated_lines(self, tmp_path): + lines = generate_geometry(mode="borehole", cable_length=200.0, spacing=50.0) + path = tmp_path / "out.txt" + save_geometry(lines, path) + assert path.read_text().splitlines() == lines + + def test_dataframe_round_trip(self, tmp_path): + lines = generate_geometry(mode="borehole", cable_length=500.0, spacing=25.0, + noise_level=1.95e-9, gauge_length=10.0) + first = tmp_path / "a.txt" + save_geometry(lines, first) + df = read_geometry(str(first)) + + second = tmp_path / "b.txt" + save_geometry(df, second) + back = read_geometry(str(second)) + + for col in ('X', 'Y', 'Z', 'NoiseLevel'): + np.testing.assert_allclose(df[col], back[col]) + assert back.attrs['gauge_length'] == pytest.approx(df.attrs['gauge_length']) + assert back.attrs['noise_type'] == df.attrs['noise_type'] + + def test_elevation_sign_survives_the_round_trip(self, tmp_path, written): + df = written([STATION_HEADER, "A\t0\t0\t-250\t1e-10"], "one.txt") + assert df['Z'][0] == pytest.approx(250.0) # Z is depth, positive down + out = tmp_path / "one_out.txt" + save_geometry(df, out) + assert "-250" in out.read_text() # written back as elevation + assert read_geometry(str(out))['Z'][0] == pytest.approx(250.0) + + @pytest.mark.parametrize("noise_type,unit", [ + ('velocity', 'm/s'), ('acceleration', 'm/s^2'), + ('displacement', 'm'), ('strain rate', '1/s'), + ]) + def test_noise_unit_is_written_from_the_noise_type(self, tmp_path, noise_type, unit): + df = pd.DataFrame({'Name': ['A'], 'X': [0.0], 'Y': [0.0], 'Z': [0.0], + 'NoiseLevel': [1e-9]}) + df.attrs['noise_type'] = noise_type + df.attrs['_is_combined'] = False + path = tmp_path / "n.txt" + save_geometry(df, path) + assert f"NoiseLevel({unit})" in path.read_text().splitlines()[0] + + def test_rejects_geometry_without_noise_metadata(self, tmp_path): + df = pd.DataFrame({'X': [0.0], 'Y': [0.0], 'Z': [0.0], 'NoiseLevel': [1e-9]}) + df.attrs['_is_combined'] = False + with pytest.raises(ValueError): + save_geometry(df, tmp_path / "bad.txt") + + def test_rejects_missing_columns(self, tmp_path): + df = pd.DataFrame({'X': [0.0], 'Y': [0.0]}) + df.attrs['noise_type'] = 'velocity' + df.attrs['_is_combined'] = False + with pytest.raises(ValueError): + save_geometry(df, tmp_path / "bad.txt") + + def test_rejects_combined_geometry(self, tmp_path, written): + a = written([STATION_HEADER, "A\t0\t0\t0\t1e-10"], "a.txt") + b = written([STATION_HEADER, "B\t100\t0\t0\t1e-10"], "b.txt") + with pytest.raises(ValueError): + save_geometry(combine_geometry(a, b), tmp_path / "combined.txt") + + def test_rejects_wrong_type(self, tmp_path): + with pytest.raises(TypeError): + save_geometry("not a list", tmp_path / "bad.txt") + with pytest.raises(TypeError): + save_geometry([1, 2, 3], tmp_path / "bad.txt") + + def test_rejects_empty_lines(self, tmp_path): + with pytest.raises(ValueError): + save_geometry([], tmp_path / "bad.txt") + + +class TestOptionalColumns: + def test_surface_and_components_round_trip(self, tmp_path, written): + lines = [ + STATION_HEADER + "\tSurface\tComponents", + "A\t0\t0\t0\t1e-10\t1.0\t3C", + "B\t1000\t0\t-500\t1e-10\t0.0\tZ", + "C\t500\t866\t0\t1e-10\t0.5\t3C", + ] + df = written(lines, "stations.txt") + assert list(df['Surface']) == [1.0, 0.0, 0.5] + assert list(df['Components']) == ['3C', 'Z', '3C'] + + out = tmp_path / "stations_out.txt" + save_geometry(df, out) + back = read_geometry(str(out)) + assert list(back['Surface']) == [1.0, 0.0, 0.5] + assert list(back['Components']) == ['3C', 'Z', '3C'] + + def test_surface_must_be_a_weight(self, written): + with pytest.raises(ValueError): + written([STATION_HEADER + "\tSurface", "A\t0\t0\t0\t1e-10\t1.5"], "bad.txt") + + def test_components_must_be_3c_or_z(self, written): + with pytest.raises(ValueError): + written([STATION_HEADER + "\tComponents", "A\t0\t0\t0\t1e-10\tN"], "bad.txt") + + def test_components_not_allowed_on_das(self, written): + lines = [ + "Northing(m)\tEasting(m)\tElevation(m)\tNoiseLevel(1/s)\tComponents\tGaugeLength=10(m)", + "0\t0\t0\t1e-9\tZ", + ] + with pytest.raises(ValueError): + written(lines, "das_components.txt") + + def test_save_rejects_components_on_das(self, tmp_path, written): + df = written([ + "Northing(m)\tEasting(m)\tElevation(m)\tNoiseLevel(1/s)\tGaugeLength=10(m)", + "0\t0\t0\t1e-9", "0\t0\t-10\t1e-9", + ], "das.txt") + df['Components'] = ['Z', 'Z'] + with pytest.raises(ValueError): + save_geometry(df, tmp_path / "bad.txt") + + +class TestDecimate: + def test_keeps_every_nth_station(self, written): + lines = generate_geometry(mode="borehole", cable_length=990.0, spacing=10.0, + gauge_length=10.0, noise_level=1e-9) + df = written(lines, "das.txt") + dec = decimate_geometry(df, 5) + assert len(dec) == int(np.ceil(len(df) / 5)) + np.testing.assert_allclose(dec['Z'].to_numpy(), df['Z'].to_numpy()[::5]) + + def test_preserves_attributes(self, written): + lines = generate_geometry(mode="borehole", cable_length=200.0, spacing=10.0, + gauge_length=10.0, noise_level=1e-9) + df = written(lines, "das.txt") + dec = decimate_geometry(df, 2) + assert dec.attrs['gauge_length'] == df.attrs['gauge_length'] + assert dec.attrs['noise_type'] == df.attrs['noise_type'] + + def test_step_of_one_is_a_copy(self, written): + df = written([STATION_HEADER, "A\t0\t0\t0\t1e-10", "B\t10\t0\t0\t1e-10"], "s.txt") + dec = decimate_geometry(df, 1) + assert len(dec) == len(df) + assert dec is not df + + @pytest.mark.parametrize("bad", [0, -1, 2.5, "5", True]) + def test_rejects_invalid_step(self, bad, written): + df = written([STATION_HEADER, "A\t0\t0\t0\t1e-10"], "s.txt") + with pytest.raises(ValueError): + decimate_geometry(df, bad) + + def test_rejects_combined_geometry(self, written): + a = written([STATION_HEADER, "A\t0\t0\t0\t1e-10"], "a.txt") + b = written([STATION_HEADER, "B\t100\t0\t0\t1e-10"], "b.txt") + with pytest.raises(ValueError): + decimate_geometry(combine_geometry(a, b), 2) + + +class TestGeneratorsAreUsable: + def test_generated_das_drives_a_sensitivity_computation(self, written): + """The generators feed straight into the modelling path.""" + lines = generate_geometry(mode="borehole", cable_length=500.0, spacing=50.0, + noise_level=1.95e-9, gauge_length=10.0) + geometry_df = written(lines, "das.txt") + velocity_df = pd.DataFrame({ + 'Depth': [0], 'Vp': [2370.0], 'VpVsRatio': [2.0], 'Qp': [100.0], + 'Qs': [100.0], 'Rho': [2500.0], 'Ep': [0.0], 'Dt': [0.0], 'Gm': [0.0]}) + velocity_df.attrs['model_type'] = 'homo' + params = pndmod.classes.DetectionParameters( + min_SNR_p=2, min_SNR_s=2, min_stations_p=3, min_stations_s=3) + grid = pndmod.utils.generate_grid(x=np.array([500.0, 1000.0]), + y=np.array([0.0]), + z=np.array([250.0, 500.0])) + mw = pndmod.core.get_mag_sensitivity(grid_coords=grid, geometry_df=geometry_df, + velocity_df=velocity_df, params=params, + wave_mode='PS') + assert mw.shape == (4,) + assert np.all(np.isfinite(mw)) From 84fa4f2cfaace805d65afe9c073bc618ca4e2fee Mon Sep 17 00:00:00 2001 From: Denis Anikiev Date: Fri, 11 Sep 2026 19:12:27 +0300 Subject: [PATCH 04/14] Add Sphinx documentation in HTML and PDF, Pages deployment and one environment Gives the package a published documentation site and consolidates the dependency declaration, so that the environment file pins only the interpreter and everything else is a module declared in pyproject.toml. Single environment: - environment.yml keeps only python and pip; runtime dependencies move to pyproject.toml and the tooling into the dev and docs optional groups. Install with conda env create -f environment.yml -n pnd followed by pip install -e . or pip install -e ".[dev,docs]" - environment-dev.yml is removed. The two files duplicated the dependency list and had drifted: both advertised scipy and plotly, neither of which the package imports, while cmcrameri, which it does import, was missing - pyproject.toml gains requires-python >=3.9,<3.13, the MIT license and its file, project URLs, and fuller keywords and classifiers. Python 3.8 is no longer advertised because it is untested and the package needs 3.9. The build requirement moves to setuptools >= 77 so the PEP 639 license form can be used, which the installed setuptools resolves to License-Expression: MIT - install.bat creates one environment, pnd, from the single file and asks whether to add the development tools, instead of choosing between two conda files and two environment names. Dry-run checked for yes, no and a bare Enter, which defaults to the runtime install, and the resulting pip command was checked with pip install --dry-run in a clean environment. The variable holding the install specification is deliberately not named PIP_something, because pip reads PIP_