Add a Python version of RConv

py-rconv/README.md

@@ -0,0 +1,116 @@
+# rconv-py
+
+A pure-Python rewrite of SeisSol's `rconv` tool.
+
+`rconv` converts kinematic rupture models from the [Standard Rupture
+Format](https://strike.scec.org/scecpedia/Standard_Rupture_Format) (SRF) to
+SeisSol's NetCDF Rupture Format (NRF). It can also emit an XDMF file for
+visualisation in ParaView.
+
+The original `rconv` is a C++ tool that depends on a deprecated version of
+PROJ.4 and on the SeisSol build system. This package replaces it with a small,
+self-contained Python module that uses modern `pyproj` and `netCDF4` — nothing
+else.
+
+## Installation
+
+```sh
+pip install .
+```
+
+Dependencies (installed automatically):
+
+* `numpy >= 1.24`
+* `pyproj >= 3.4`
+* `netCDF4 >= 1.6`
+
+## Quick start
+
+```sh
+# Convert SRF → NRF for SeisSol input.
+rconv to-nrf input.srf output.nrf \
+    --mcs "+proj=tmerc +datum=WGS84 +k=0.9996 +lon_0=-118.5150 +lat_0=34.3440 +axis=enu"
+
+# Produce an XDMF visualisation file for ParaView.
+rconv to-xdmf input.srf visualization.xmf
+```
+
+Equivalent invocation as a module (no entry-point installation needed):
+
+```sh
+python -m rconv to-nrf input.srf output.nrf --mcs "..."
+```
+
+## CLI
+
+Two subcommands cover the two distinct conversions:
+
+```
+rconv to-nrf INPUT OUTPUT --mcs PROJ_STRING [--normalize-onset] [--source-proj PROJ_STRING]
+rconv to-xdmf INPUT OUTPUT [--vcs PROJ_STRING] [--source-proj PROJ_STRING]
+```
+
+* `--mcs` — PROJ string for the *mesh* coordinate system. **Required** for NRF
+  output. Determines both how lon/lat/depth becomes (x, y, z) and how
+  slip-vector components are oriented (via the `+axis=` option).
+* `--vcs` — PROJ string for the *visualisation* coordinate system. Defaults to
+  `+proj=geocent +datum=WGS84 +units=m +no_defs` (a geocentric globe view).
+* `--normalize-onset` — subtract the minimum onset time across all sources.
+  Useful when the SRF's absolute clock is meaningless.
+* `--source-proj` — override the source CRS (default: WGS84 longlat). Useful
+  only when the SRF has already-projected coordinates.
+
+## Programmatic API
+
+```python
+from rconv import parse_srf, CoordinateMapping, write_nrf, write_xmf
+
+sources = parse_srf("northridge.srf")
+mapping = CoordinateMapping(
+    "+proj=tmerc +datum=WGS84 +k=0.9996 +lon_0=-118.515 +lat_0=34.344 +axis=ned"
+)
+write_nrf("northridge.nrf", sources, mapping)
+write_xmf("northridge.xmf", sources, mapping)
+```
+
+## What `+axis=` means
+
+The `+axis=` value in your PROJ string controls how a canonical
+(north, east, down) triple — the natural frame of SRF coordinates — is
+permuted and signed into the mesh CS. It is three characters from `enwsud`:
+
+* `+axis=enu` — x = east, y = north, z = up (PROJ's default if you omit
+  `+axis=` entirely).
+* `+axis=ned` — x = north, y = east, z = down.
+* `+axis=seu` — x = south, y = east, z = up.
+
+The same string also controls how slip vectors (along strike, in fault plane,
+along normal) are oriented in your mesh. If your moment tensors look wrong,
+double-check `+axis=` first.
+
+## Differences from the C++ tool
+
+* **CLI**: two focused subcommands instead of five short flags.
+* **PROJ**: uses modern `pyproj` (which wraps PROJ 8/9), not the legacy
+  PROJ.4 API the C++ tool depended on.
+* **Output**: NRF output is functionally identical and consumed by SeisSol's
+  NRF reader without modification. Two cosmetic differences:
+  - The `Subfault_units` compound type is replaced by a plain string
+    attribute on `subfaults`. SeisSol's reader never inspects either form,
+    so this is invisible to downstream use; only `ncdump` output differs.
+  - Slip-rate dimensions of size 0 are written as UNLIMITED rather than
+    fixed-size-0. Again, the reader doesn't care.
+* **Examples**: ships with a test suite that round-trips the Northridge SRF
+  from `SeisSol/Examples`.
+
+## Development
+
+Run the tests:
+
+```sh
+pip install -e .[dev]
+pytest
+```
+
+The Northridge integration tests are skipped automatically if the SRF file
+isn't available at `~/Examples/Northridge/`.

py-rconv/pyproject.toml

@@ -0,0 +1,43 @@
+[build-system]
+requires = ["setuptools>=64", "wheel"]
+build-backend = "setuptools.build_meta"
+
+[project]
+name = "rconv"
+version = "0.1.0"
+description = "Convert SRF kinematic rupture models to SeisSol's NetCDF Rupture Format (NRF)."
+readme = "README.md"
+requires-python = ">=3.10"
+license = { text = "BSD-3-Clause" }
+authors = [
+    { name = "SeisSol Group" },
+]
+keywords = ["seismology", "srf", "nrf", "seissol", "rupture", "earthquake"]
+classifiers = [
+    "Development Status :: 4 - Beta",
+    "Intended Audience :: Science/Research",
+    "License :: OSI Approved :: BSD License",
+    "Programming Language :: Python :: 3",
+    "Programming Language :: Python :: 3 :: Only",
+    "Topic :: Scientific/Engineering",
+]
+dependencies = [
+    "numpy>=1.24",
+    "pyproj>=3.4",
+    "netCDF4>=1.6",
+]
+
+[project.optional-dependencies]
+dev = [
+    "pytest>=7",
+]
+
+[project.scripts]
+rconv = "rconv.cli:main"
+
+[project.urls]
+Homepage = "https://github.com/SeisSol/SeisSol"
+Documentation = "https://seissol.readthedocs.io/en/latest/standard-rupture-format.html"
+
+[tool.setuptools.packages.find]
+where = ["src"]

py-rconv/src/rconv/__init__.py

@@ -0,0 +1,45 @@
+"""rconv-py — convert SRF kinematic rupture models to SeisSol's NRF.
+
+A pure-Python replacement for the C++ ``rconv`` tool that ships with SeisSol.
+This package is dependency-free with respect to SeisSol; it relies only on
+``numpy``, ``pyproj`` and ``netCDF4``.
+
+The public API is small:
+
+* :func:`rconv.srf.parse_srf` — read an SRF file into :class:`PointSource`\\ s.
+* :class:`rconv.mapping.CoordinateMapping` — transform geographic coords and
+  fault-local vectors into a mesh CS.
+* :func:`rconv.nrf.write_nrf` — emit a NetCDF Rupture Format file for SeisSol.
+* :func:`rconv.xmf.write_xmf` — emit an XDMF file for ParaView.
+
+Example
+-------
+
+>>> from rconv import parse_srf, CoordinateMapping, write_nrf
+>>> sources = parse_srf("northridge.srf")
+>>> mapping = CoordinateMapping(
+...     "+proj=tmerc +datum=WGS84 +k=0.9996 +lon_0=-118.515 +lat_0=34.344 "
+...     "+units=m +axis=ned"
+... )
+>>> write_nrf("northridge.nrf", sources, mapping)
+"""
+
+from __future__ import annotations
+
+from .mapping import CoordinateMapping, Vec3
+from .nrf import read_nrf, write_nrf
+from .srf import PointSource, SRFParseError, parse_srf
+from .xmf import write_xmf
+
+__all__ = [
+    "CoordinateMapping",
+    "PointSource",
+    "SRFParseError",
+    "Vec3",
+    "parse_srf",
+    "read_nrf",
+    "write_nrf",
+    "write_xmf",
+]
+
+__version__ = "0.1.0"

py-rconv/src/rconv/__main__.py

@@ -0,0 +1,10 @@
+"""Allow ``python -m rconv ...`` invocation."""
+
+from __future__ import annotations
+
+import sys
+
+from .cli import main
+
+if __name__ == "__main__":
+    sys.exit(main())

py-rconv/src/rconv/cli.py

@@ -0,0 +1,148 @@
+"""Modern command-line interface for ``rconv``.
+
+The classic C++ tool accepts five short flags that mix concerns (some required
+together, others mutually independent). Here we split into two focused
+sub-commands so the help is self-documenting and shell-completion is happy:
+
+* ``rconv to-nrf INPUT OUTPUT --mcs ...`` — produce a SeisSol input file.
+* ``rconv to-xdmf INPUT OUTPUT [--vcs ...]`` — produce a visualisation file.
+
+Both commands accept ``--source-proj`` to override the default input CRS
+(WGS84 longlat) on the rare occasion an SRF arrives in something else.
+"""
+
+from __future__ import annotations
+
+import argparse
+import logging
+import sys
+from pathlib import Path
+
+from .mapping import CoordinateMapping
+from .nrf import write_nrf
+from .srf import parse_srf
+from .xmf import write_xmf
+
+_DEFAULT_VCS = "+proj=geocent +datum=WGS84 +units=m +no_defs"
+
+
+def main(argv: list[str] | None = None) -> int:
+    """Entry point. Returns a process exit code."""
+    parser = _build_parser()
+    args = parser.parse_args(argv)
+
+    logging.basicConfig(
+        level=logging.DEBUG if args.verbose else logging.INFO,
+        format="%(levelname)s %(name)s: %(message)s",
+    )
+
+    if args.command == "to-nrf":
+        return _cmd_to_nrf(args)
+    if args.command == "to-xdmf":
+        return _cmd_to_xdmf(args)
+    parser.print_help()
+    return 2
+
+
+# ---------------------------------------------------------------------------
+# argparse plumbing
+# ---------------------------------------------------------------------------
+
+def _build_parser() -> argparse.ArgumentParser:
+    parser = argparse.ArgumentParser(
+        prog="rconv",
+        description=(
+            "Convert SRF (Standard Rupture Format) files to NRF (NetCDF "
+            "Rupture Format) for SeisSol, or to XDMF for visualisation."
+        ),
+    )
+    parser.add_argument(
+        "-v", "--verbose", action="store_true", help="Enable debug logging.",
+    )
+    sub = parser.add_subparsers(dest="command", required=True, metavar="COMMAND")
+
+    nrf = sub.add_parser(
+        "to-nrf",
+        help="Convert SRF → NRF.",
+        description="Convert an SRF file to a NetCDF Rupture Format file for SeisSol.",
+    )
+    nrf.add_argument("input", type=Path, help="Input SRF file (.srf).")
+    nrf.add_argument("output", type=Path, help="Output NRF file (.nrf).")
+    nrf.add_argument(
+        "--mcs", required=True, metavar="PROJ_STRING",
+        help=(
+            "PROJ string describing the mesh coordinate system, e.g. "
+            '"+proj=utm +zone=10 +datum=WGS84 +units=m +axis=ned +no_defs". '
+            "The +axis= option determines how slip vectors are oriented in the "
+            "MCS (default: enu)."
+        ),
+    )
+    nrf.add_argument(
+        "--normalize-onset", action="store_true",
+        help="Subtract the minimum onset time across sources from every source.",
+    )
+    nrf.add_argument(
+        "--source-proj", default=None, metavar="PROJ_STRING",
+        help=(
+            "Override the source CRS (defaults to WGS84 longlat). Useful only "
+            "if your SRF carries already-projected coordinates."
+        ),
+    )
+
+    xmf = sub.add_parser(
+        "to-xdmf",
+        help="Convert SRF → XDMF (for ParaView).",
+        description="Convert an SRF file to an XDMF visualisation file.",
+    )
+    xmf.add_argument("input", type=Path, help="Input SRF file (.srf).")
+    xmf.add_argument("output", type=Path, help="Output XDMF file (.xmf/.xdmf).")
+    xmf.add_argument(
+        "--vcs", default=_DEFAULT_VCS, metavar="PROJ_STRING",
+        help=(
+            "PROJ string for the visualisation coordinate system. "
+            f"Defaults to {_DEFAULT_VCS!r} (geocentric — gives a globe view)."
+        ),
+    )
+    xmf.add_argument(
+        "--source-proj", default=None, metavar="PROJ_STRING",
+        help="Override the source CRS (defaults to WGS84 longlat).",
+    )
+
+    return parser
+
+
+# ---------------------------------------------------------------------------
+# commands
+# ---------------------------------------------------------------------------
+
+def _cmd_to_nrf(args: argparse.Namespace) -> int:
+    sources = _read_srf(args.input)
+    if not sources:
+        print(f"error: {args.input} produced no usable point sources", file=sys.stderr)
+        return 1
+    mapping = CoordinateMapping(args.mcs, source_proj=args.source_proj)
+    write_nrf(args.output, sources, mapping, normalize_onset=args.normalize_onset)
+    print(f"Wrote {args.output} ({len(sources)} sources).")
+    return 0
+
+
+def _cmd_to_xdmf(args: argparse.Namespace) -> int:
+    sources = _read_srf(args.input)
+    if not sources:
+        print(f"error: {args.input} produced no usable point sources", file=sys.stderr)
+        return 1
+    mapping = CoordinateMapping(args.vcs, source_proj=args.source_proj)
+    write_xmf(args.output, sources, mapping)
+    print(f"Wrote {args.output} ({len(sources)} sources).")
+    return 0
+
+
+def _read_srf(path: Path) -> list:
+    print(f"Reading {path} ...", end=" ", flush=True)
+    sources = parse_srf(path)
+    print(f"{len(sources)} sources.")
+    return sources
+
+
+if __name__ == "__main__":
+    sys.exit(main())