diff --git a/.gitignore b/.gitignore index 557aca80..c905252b 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,7 @@ build/ dist/ qcore.egg-info +qcore_utils.egg-info **/*.pyc **/*.coverage .cache/ diff --git a/qcore/test/test_xyts/test_xyts.py b/qcore/test/test_xyts/test_xyts.py index 08ac3d4e..1a87eb95 100644 --- a/qcore/test/test_xyts/test_xyts.py +++ b/qcore/test/test_xyts/test_xyts.py @@ -1,14 +1,71 @@ """Test module for XYTS file processing using pytest fixtures.""" +import struct +import tempfile from pathlib import Path from urllib import request import numpy as np import pytest +from hypothesis import HealthCheck, given, settings +from hypothesis import strategies as st +from hypothesis.extra.numpy import arrays from qcore import xyts +def _make_proc_local_header( + endian: str, + x0: int, y0: int, z0: int, t0: int, + local_nx: int, local_ny: int, local_nz: int, + nx: int, ny: int, nz: int, nt: int, + dx: float, dy: float, hh: float, dt: float, + mrot: float, mlat: float, mlon: float, +) -> bytes: + """Return a 72-byte proc-local (tsheader_procP3) header as bytes.""" + pfx = ">" if endian == ">" else "<" + ints = struct.pack( + f"{pfx}11i", + x0, y0, z0, t0, + local_nx, local_ny, local_nz, + nx, ny, nz, nt, + ) + floats = struct.pack( + f"{pfx}7f", + dx, dy, hh, dt, mrot, mlat, mlon, + ) + return ints + floats # 44 + 28 = 72 bytes + + +@st.composite +def xyzts_file_data(draw: st.DrawFn) -> tuple: + """Hypothesis composite strategy generating valid XYZTS file parameters and payload. + + Returns + ------- + tuple + (endian, ncomp, local_nx, local_ny, local_nz, nx, ny, nt, payload) + where payload has shape (nt, ncomp, local_nz, local_ny, local_nx). + """ + endian = draw(st.sampled_from([">", "<"])) + ncomp = draw(st.sampled_from([3, 6, 9])) + local_nx = draw(st.integers(min_value=1, max_value=8)) + local_ny = draw(st.integers(min_value=1, max_value=8)) + # local_nz > 1 distinguishes XYZTS from surface proc-local XYTS + local_nz = draw(st.integers(min_value=2, max_value=4)) + nx = draw(st.integers(min_value=local_nx, max_value=20)) + ny = draw(st.integers(min_value=local_ny, max_value=20)) + nt = draw(st.integers(min_value=1, max_value=4)) + payload = draw( + arrays( + dtype=np.float32, + shape=(nt, ncomp, local_nz, local_ny, local_nx), + elements=st.floats(allow_nan=False, allow_infinity=False, width=32), + ) + ) + return endian, ncomp, local_nx, local_ny, local_nz, nx, ny, nt, payload + + @pytest.fixture(scope="session") def xyts_file() -> xyts.XYTSFile: """Provide path to the XYTS test file.""" @@ -192,3 +249,99 @@ def test_tslice_get( test_output = xyts_file.tslice_get(step, comp=comp) sample_array = np.fromfile(sample_file, dtype="3 None: + """XYTSFile correctly detects and parses a synthetic XYZTS file (round-trip). + + Checks auto-detection, header field parsing, 5-D data shape, and payload + value fidelity across all valid ncomp values and both endiannesses. + """ + endian, ncomp, local_nx, local_ny, local_nz, nx, ny, nt, payload = data + header = _make_proc_local_header( + endian, + x0=0, y0=0, z0=0, t0=0, + local_nx=local_nx, local_ny=local_ny, local_nz=local_nz, + nx=nx, ny=ny, nz=local_nz, nt=nt, + dx=0.4, dy=0.4, hh=0.1, dt=0.02, + mrot=0.0, mlat=-43.5, mlon=172.0, + ) + file_bytes = header + payload.astype(f"{endian}f4").tobytes() + with tempfile.TemporaryDirectory() as tmpdir: + fpath = Path(tmpdir) / "test_xyzts-0" + fpath.write_bytes(file_bytes) + xf = xyts.XYTSFile(fpath) + assert int(xf.local_nz) == local_nz + assert int(xf.local_ny) == local_ny + assert int(xf.local_nx) == local_nx + assert xf.ncomp == ncomp + assert int(xf.nt) == nt + assert xf.data is not None + assert xf.data.ndim == 5 + assert xf.data.shape == (nt, ncomp, local_nz, local_ny, local_nx) + assert xf.data == pytest.approx(payload) + + +@given(data=xyzts_file_data()) +def test_xyzts_tslice_get_raises(data: tuple) -> None: + """tslice_get should raise ValueError for volumetric XYZTS files.""" + endian, ncomp, local_nx, local_ny, local_nz, nx, ny, nt, payload = data + header = _make_proc_local_header( + endian, + x0=0, y0=0, z0=0, t0=0, + local_nx=local_nx, local_ny=local_ny, local_nz=local_nz, + nx=nx, ny=ny, nz=local_nz, nt=nt, + dx=0.4, dy=0.4, hh=0.1, dt=0.02, + mrot=0.0, mlat=-43.5, mlon=172.0, + ) + file_bytes = header + payload.astype(f"{endian}f4").tobytes() + with tempfile.TemporaryDirectory() as tmpdir: + fpath = Path(tmpdir) / "test_xyzts-0" + fpath.write_bytes(file_bytes) + xf = xyts.XYTSFile(fpath) + with pytest.raises(ValueError, match="tslice_get"): + xf.tslice_get(0) + + +@given(data=xyzts_file_data()) +def test_xyzts_pgv_raises(data: tuple) -> None: + """pgv() should raise ValueError for volumetric XYZTS files.""" + endian, ncomp, local_nx, local_ny, local_nz, nx, ny, nt, payload = data + header = _make_proc_local_header( + endian, + x0=0, y0=0, z0=0, t0=0, + local_nx=local_nx, local_ny=local_ny, local_nz=local_nz, + nx=nx, ny=ny, nz=local_nz, nt=nt, + dx=0.4, dy=0.4, hh=0.1, dt=0.02, + mrot=0.0, mlat=-43.5, mlon=172.0, + ) + file_bytes = header + payload.astype(f"{endian}f4").tobytes() + with tempfile.TemporaryDirectory() as tmpdir: + fpath = Path(tmpdir) / "test_xyzts-0" + fpath.write_bytes(file_bytes) + xf = xyts.XYTSFile(fpath) + with pytest.raises(ValueError, match="pgv"): + xf.pgv() + + +@given(data=xyzts_file_data()) +def test_xyzts_meta_only(data: tuple) -> None: + """meta_only=True should work for XYZTS files and leave data=None.""" + endian, ncomp, local_nx, local_ny, local_nz, nx, ny, nt, payload = data + header = _make_proc_local_header( + endian, + x0=0, y0=0, z0=0, t0=0, + local_nx=local_nx, local_ny=local_ny, local_nz=local_nz, + nx=nx, ny=ny, nz=local_nz, nt=nt, + dx=0.4, dy=0.4, hh=0.1, dt=0.02, + mrot=0.0, mlat=-43.5, mlon=172.0, + ) + file_bytes = header + payload.astype(f"{endian}f4").tobytes() + with tempfile.TemporaryDirectory() as tmpdir: + fpath = Path(tmpdir) / "test_xyzts-0" + fpath.write_bytes(file_bytes) + xf = xyts.XYTSFile(fpath, meta_only=True) + assert xf.data is None + assert int(xf.local_nz) == local_nz + assert xf.ncomp == ncomp diff --git a/qcore/xyts.py b/qcore/xyts.py index dbd3d17a..1e475dce 100644 --- a/qcore/xyts.py +++ b/qcore/xyts.py @@ -1,15 +1,40 @@ """ -This module provides functionality to read xyts files. +This module provides functionality to read xyts and xyzts files. Extended Summary ---------------- -This module includes the XYTSFile class, which represents an XYTS file. It -allows users to load metadata, retrieve data, and calculate PGV (Peak Ground -Velocity) and MMI (Modified Mercalli Intensity) values from the XYTS file. +This module includes the XYTSFile class, which represents an XYTS or XYZTS +file. It allows users to load metadata, retrieve data, and calculate PGV +(Peak Ground Velocity) and MMI (Modified Mercalli Intensity) values from the +XYTS file. + +Two file layouts are supported: + +1. **Standard XYTS** (60-byte ``tsheader`` header, ``nz=1``): + Surface XY-plane timeslices with exactly 3 velocity components. + Data shape: ``(nt, 3, ny, nx)``. + +2. **Proc-local timeslice** (72-byte ``tsheader_procP`` / ``tsheader_procP3`` + header): + Written per MPI rank. Two sub-variants are auto-detected from the header: + + a. **Proc-local XYTS** – ``local_nz == 1``, ``ncomp == 3``. + Data shape: ``(nt, 3, local_ny, local_nx)``. + Requires ``proc_local_file=True`` (backwards-compatible flag). + + b. **Proc-local XYZTS** (EMOD3D ≥ v3.0.13, ``ts_xyz`` output) – + ``local_nz > 1``, ``ncomp ∈ {3, 6, 9}``. + Detected automatically without any user flag. + Data shape: ``(nt, ncomp, local_nz, local_ny, local_nx)``. + Filenames typically follow the ``*_xyzts-*`` pattern. + +For proc-local files ``ncomp`` is derived from the file size so that no +explicit ``ncomp`` field is required in the header. Classes ---------------- -- XYTSFile: Represents an XYTS file and provides methods to interact with it. +- XYTSFile: Represents an XYTS / XYZTS file and provides methods to interact + with it. Notes ----- @@ -25,8 +50,14 @@ Examples -------- -# Load an XYTS file -xyts_file = XYTSFile("example.x3d") +# Load a standard XYTS file +xyts_file = XYTSFile("example.e3d") + +# Load a proc-local XYTS file (backwards-compatible flag) +xyts_file = XYTSFile("example_xyts-000000", proc_local_file=True) + +# Load a volumetric XYZTS file (auto-detected – no flag needed) +xyzts_file = XYTSFile("example_xyzts-000000") # Retrieve corners of the simulation domain corners = xyts_file.corners() @@ -39,13 +70,19 @@ """ import dataclasses +from enum import Enum from math import cos, radians, sin from pathlib import Path import numpy as np from qcore import geo -from enum import Enum + +# Maximum plausible grid dimension. Values up to this threshold are treated +# as valid local_nx / local_ny during endianness detection; byte-swapped +# representations of small integers exceed this bound by several orders of +# magnitude (e.g. 100 byte-swapped → 1 677 721 600). +_MAX_GRID_DIM = 0xFFFF class Component(Enum): @@ -60,9 +97,9 @@ class Component(Enum): @dataclasses.dataclass class XYTSFile: """ - Represents an XYTS file containing time slices on the X-Y plane (z = 1, top level). - This class provides methods to read metadata, retrieve data, and calculate - PGV (Peak Ground Velocity) and MMI (Modified Mercalli Intensity) values. + Represents an XYTS or XYZTS file. + + Supports three file layouts – see module docstring for details. Assumptions: - dip = 0: Simulation domain is flat. @@ -73,9 +110,9 @@ class XYTSFile: y0: Starting y-coordinate. z0: Starting z-coordinate. t0: Starting time. - local_nx: Number of local x-coordinates (for proc-local files only). - local_ny: Number of local y-coordinates (for proc-local files only). - local_nz: Number of local z-coordinates (for proc-local files only). + local_nx: Number of local x-coordinates (proc-local files only). + local_ny: Number of local y-coordinates (proc-local files only). + local_nz: Number of local z-coordinates (proc-local files only). nx: Total number of x-coordinates. ny: Total number of y-coordinates. nz: Total number of z-coordinates. @@ -87,6 +124,9 @@ class XYTSFile: mrot: Rotation angle for model origin. mlat: Latitude of the model origin. mlon: Longitude of the model origin. + ncomp: Number of components per grid point. Always 3 for standard + XYTS and proc-local XYTS files. For XYZTS files this is 3, 6, + or 9 as set by ``ts_xyz_ncomp`` in the EMOD3D configuration. dxts: Original simulation grid spacing in the x-direction. dyts: Original simulation grid spacing in the y-direction. nx_sim: Original simulation size in the x-direction. @@ -99,6 +139,9 @@ class XYTSFile: sinP: Sine of the dip angle. rot_matrix: Rotation matrix for components. data: Memory-mapped array containing the data. + Shape is ``(nt, ncomp, ny, nx)`` for standard XYTS, + ``(nt, ncomp, local_ny, local_nx)`` for proc-local XYTS, and + ``(nt, ncomp, local_nz, local_ny, local_nx)`` for XYZTS. ll_map: Longitude-latitude map for data. Methods: @@ -113,10 +156,11 @@ class XYTSFile: Returns the simulation region as a tuple (x_min, x_max, y_min, y_max). tslice_get(step, comp=-1, outfile=None): - Retrieves timeslice data. + Retrieves timeslice data (standard and proc-local XYTS only). pgv(mmi=False, pgvout=None, mmiout=None): - Retrieves PGV map and optionally calculates MMI. + Retrieves PGV map and optionally calculates MMI (standard XYTS + only). """ # Header values @@ -137,6 +181,7 @@ class XYTSFile: mlat: float mlon: float # Derived values + ncomp: int dxts: int dyts: int nx_sim: int @@ -171,12 +216,16 @@ def __init__( Parameters ---------- xyts_path : Path | str - Path to the xyts file. + Path to the xyts / xyzts file. meta_only : bool If True, only loads metadata and doesn't prepare gridpoint datum locations (slower). proc_local_file : bool - If True, indicates a proc-local file. + If True, forces reading with a 72-byte proc-local header. This + is still required for proc-local XYTS files whose ``local_nz`` + equals 1, because those cannot be distinguished from standard XYTS + files by header inspection alone. XYZTS files (``local_nz > 1``) + are detected automatically and do not need this flag. round_dt : bool If True, round the dt value to 4dp (present only for backwards compatibility). @@ -184,29 +233,71 @@ def __init__( Raises ------ ValueError - ValueError: If the file is not an XY timeslice file. + If the file layout cannot be determined from the header bytes. """ xytf = open(xyts_path, "rb") self.xyts_path = xyts_path - # determine endianness, an x-y timeslice has 1 z value - nz = np.fromfile(xytf, dtype=">i4", count=7)[-1] - if nz == 0x00000001: + # Detect endianness and file variant by inspecting the first 7 i32 words. + # + # Standard XYTS (60-byte tsheader): + # offset 0: [x0, y0, z0, t0] (4 × i32) + # offset 16: [nx, ny, nz, nt] (4 × i32) ← word 7 = nz == 1 + # offset 32: [dx, dy, hh, dt, mrot, mlat, mlon] (7 × f32) + # + # Proc-local XYTS / XYZTS (72-byte tsheader_procP / tsheader_procP3): + # offset 0: [x0, y0, z0, t0] (4 × i32) + # offset 16: [local_nx, local_ny, local_nz] (3 × i32) ← word 7 = local_nz + # offset 28: [nx, ny, nz, nt] (4 × i32) + # offset 44: [dx, dy, hh, dt, mrot, mlat, mlon] (7 × f32) + # + # For standard XYTS, word 7 (nz) is always 1, giving a reliable endianness + # probe via the byte-swapped value 0x01000000. For XYZTS files local_nz > 1, + # so the probe fails and we fall through to the endianness tiebreaker below. + raw7_be = np.fromfile(xytf, dtype=">i4", count=7) + seventh_be = int(raw7_be[-1]) + + if seventh_be == 0x00000001: endian = ">" - elif nz == 0x01000000: + proc_local = proc_local_file + elif seventh_be == 0x01000000: endian = "<" + proc_local = proc_local_file else: - xytf.close() - raise ValueError("File is not an XY timeslice file: %s" % (xyts_path)) + # local_nz > 1: proc-local XYZTS file. Determine endianness from + # local_nx (word 5, index 4): byte-swapping a small integer produces + # a value orders of magnitude larger, so only one interpretation will + # be in the plausible range 1…_MAX_GRID_DIM. + local_nx_be = int(raw7_be[4]) + raw7_le = np.frombuffer(raw7_be.tobytes(), dtype=" 1) : (nt, ncomp, local_nz, local_ny, local_nx) + if proc_local: + if int(self.local_nz) > 1: + # Volumetric XYZTS: 5-D tensor + shape: tuple[int, ...] = ( + int(self.nt), + self.ncomp, + int(self.local_nz), + int(self.local_ny), + int(self.local_nx), + ) + else: + # Surface proc-local XYTS: 4-D tensor (backwards-compatible) + shape = ( + int(self.nt), + self.ncomp, + int(self.local_ny), + int(self.local_nx), + ) self.data = np.memmap( xyts_path, dtype="%sf4" % (endian), mode="r", offset=72, - shape=(self.nt, len(self.comps), self.local_ny, self.local_nx), + shape=shape, ) else: # memory map for data section @@ -270,7 +406,7 @@ def __init__( dtype="%sf4" % (endian), mode="r", offset=60, - shape=(self.nt, len(self.comps), self.ny, self.nx), + shape=(int(self.nt), self.ncomp, int(self.ny), int(self.nx)), ) # create longitude, latitude map for data @@ -360,6 +496,10 @@ def tslice_get( ) -> np.ndarray: """Retrieves timeslice data. + This method operates on the surface XY plane and is intended for + standard XYTS and proc-local XYTS files (``local_nz == 1``). For + volumetric XYZTS files access ``self.data`` directly. + Parameters ---------- step : int @@ -371,7 +511,17 @@ def tslice_get( ------- np.ndarray Retrieved timeslice data. + + Raises + ------ + ValueError + If called on a volumetric XYZTS file (``local_nz > 1``). """ + if self.local_nz is not None and int(self.local_nz) > 1: + raise ValueError( + "tslice_get() is not supported for volumetric XYZTS files " + "(local_nz > 1). Access self.data directly." + ) match comp: case Component.MAGNITUDE: return np.linalg.norm(self.data[step, :3, :, :], axis=0) @@ -396,6 +546,9 @@ def pgv( ) -> None | np.ndarray | tuple[np.ndarray, np.ndarray]: """Retrieves PGV and/or MMI map. + This method is intended for standard XYTS files. It is not applicable + to volumetric XYZTS files. + Parameters ---------- mmi : bool @@ -410,7 +563,17 @@ def pgv( None | np.ndarray | Tuple[np.ndarray, np.ndarray] PGV map or tuple of (PGV map, MMI map) or None (if both are written to a file). + + Raises + ------ + ValueError + If called on a volumetric XYZTS file (``local_nz > 1``). """ + if self.local_nz is not None and int(self.local_nz) > 1: + raise ValueError( + "pgv() is not supported for volumetric XYZTS files " + "(local_nz > 1)." + ) # PGV as timeslices reduced to maximum value at each point pgv = np.zeros(self.nx * self.ny) for ts in range(self.t0, self.nt):