Add NCM profile fetching and loading functionality

- Implemented `fetch_ncm_profile` to retrieve geophysical profiles from the USGS National Crustal Model and save them as JSON (gzip-compressed option available).
- Added `load_ncm_profile` to load and parse the NCM JSON files into a pandas DataFrame, with optional simplification of adjacent layers.
- Introduced `_simplify_profile` to merge layers with similar shear-wave velocities.
- Enhanced `Layer` class with a method to compute layer-adjusted damping based on strain.
- Updated `adjusted_damping_curves` property to reflect layer-specific minimum damping.
- Added mean effective stress calculation in `calc_mean_eff_stress` to support new profile features.

Author: arkottke

.pre-commit-config.yaml

@@ -1,7 +1,7 @@
 repos:
     - repo: https://github.com/astral-sh/ruff-pre-commit
       # Ruff version.
-      rev: v0.14.14
+      rev: v0.15.4
       hooks:
           # Run the linter.
           - id: ruff
@@ -14,14 +14,10 @@ repos:
       hooks:
           - id: pyupgrade
             args: [--py3-plus, --py310-plus]
-    - repo: https://github.com/datarootsio/databooks
-      rev: 1.3.10
-      hooks:
-          - id: databooks-meta
-            #   - repo: https://github.com/pycqa/pydocstyle
-            #     rev: 6.1.1  # pick a git hash / tag to point to
-            #     hooks:
-            #       - id: pydocstyle
+    #   - repo: https://github.com/pycqa/pydocstyle
+    #     rev: 6.1.1  # pick a git hash / tag to point to
+    #     hooks:
+    #       - id: pydocstyle
     - repo: https://github.com/pre-commit/pre-commit-hooks
       rev: v6.0.0
       hooks:

examples/README.md

@@ -22,6 +22,8 @@ The examples included here are used to demonstrate the capabilities of
 14. Example with multiple input motions and simulated soil profiles.
 15. List of provided nonlinear curves.
 16. Example of using the logic tree.
+17. Building a site profile from the USGS National Crustal Model with
+    Darendeli nonlinear curves and κ₀-adjusted damping.
 
 To be added:
 

examples/data/ncm_example.json.gz

@@ -1,10 +1,16 @@
 """Generic velocity profiles."""
 
+import gzip as _gzip
+import json
+import os
+import urllib.request
 import warnings
 
 import numpy as np
 import pandas as pd
 
+from .tools import calc_mean_eff_stress
+
 
 def kea16_profile(depth: float, vs30: float, region: str) -> pd.DataFrame:
     """
@@ -95,3 +101,191 @@ def kea16_profile(depth: float, vs30: float, region: str) -> pd.DataFrame:
     )
 
     return df
+
+
+NCM_BASE_URL = "https://earthquake.usgs.gov/ws/nshmp/ncm/geophysical"
+
+
+def fetch_ncm_profile(
+    latitude: float,
+    longitude: float,
+    fpath: str | os.PathLike,
+    depth_start: float = 0,
+    depth_step: float = 5,
+    depth_end: float = 1e4,
+    gzip_output: bool = True,
+) -> None:
+    """Fetch a geophysical profile from the USGS National Crustal Model.
+
+    Queries the NCM geophysical endpoint for a single location and saves the
+    JSON response to *fpath*.
+
+    Parameters
+    ----------
+    latitude : float
+        Latitude of the site [degrees].
+    longitude : float
+        Longitude of the site [degrees].
+    fpath : str or path-like
+        Output file path.  If *gzip_output* is True the file is gzip-compressed.
+    depth_start : float, optional
+        Minimum depth [m] (default 0).
+    depth_step : float, optional
+        Depth increment [m] (default 5).
+    depth_end : float, optional
+        Maximum depth [m] (default 10 000).
+    gzip_output : bool, optional
+        If True (default), write a gzip-compressed JSON file.
+
+    Raises
+    ------
+    urllib.error.HTTPError
+        If the request to the NCM service fails.
+    RuntimeError
+        If the response status is not ``"sucess"`` (sic).
+    """
+    url = (
+        f"{NCM_BASE_URL}"
+        f"?location={latitude},{longitude}"
+        f"&depths={depth_start},{depth_step},{depth_end}"
+    )
+    with urllib.request.urlopen(url) as resp:
+        data = json.loads(resp.read().decode())
+
+    if data.get("status") not in ("sucess", "success"):
+        raise RuntimeError(f"NCM request failed with status: {data.get('status')}")
+
+    fpath = os.fspath(fpath)
+    if gzip_output:
+        with _gzip.open(fpath, "wt") as f:
+            json.dump(data, f)
+    else:
+        with open(fpath, "w") as f:
+            json.dump(data, f)
+
+
+def load_ncm_profile(
+    fpath: str | os.PathLike,
+    simplify: bool = False,
+    simplify_tol: float = 0.02,
+) -> tuple[pd.DataFrame, float]:
+    """Load an NCM geophysical profile and return a DataFrame.
+
+    Parameters
+    ----------
+    fpath : str or path-like
+        Path to the NCM JSON file (plain or gzip-compressed).
+    simplify : bool, optional
+        If True, merge adjacent layers whose shear-wave velocities differ by
+        less than *simplify_tol* (default False).
+    simplify_tol : float, optional
+        Relative tolerance on ``vel_shear`` for merging adjacent layers
+        (default 0.02, i.e. 2 %).
+
+    Returns
+    -------
+    df : pandas.DataFrame
+        DataFrame with columns ``depth`` [m], ``vel_shear`` [m/s], and
+        ``density`` [kg/m³].
+    water_table_depth : float
+        Depth to the water table [m].
+
+    Raises
+    ------
+    ValueError
+        If the file cannot be parsed or contains no results.
+    """
+    fpath = os.fspath(fpath)
+
+    # Try gzip first, fall back to plain JSON
+    try:
+        with _gzip.open(fpath, "rt") as f:
+            data = json.load(f)
+    except _gzip.BadGzipFile:
+        with open(fpath) as f:
+            data = json.load(f)
+
+    results = data["response"]["results"]
+    if not results:
+        raise ValueError("NCM response contains no results")
+
+    result = results[0]
+    water_table_depth = result["location"]["water_table_depth"]
+
+    depths = np.array(data["request"]["depths"]["depth_vector"], dtype=float)
+    vs = np.array(result["profile"]["vs"], dtype=float)
+    density = np.array(result["profile"]["density"], dtype=float)
+    unit_wt = density * 9.81 / 1000
+
+    df = pd.DataFrame({"depth": depths, "vel_shear": vs, "unit_wt": unit_wt})
+
+    if simplify:
+        df = _simplify_profile(df, simplify_tol)
+
+    # Compute mean effective stress
+    df["mean_eff_stress"] = calc_mean_eff_stress(
+        df["depth"].values, df["unit_wt"].values, water_table_depth
+    )
+    df["thickness"] = np.diff(df["depth"], append=df["depth"].iloc[-1])
+
+    return df, water_table_depth
+
+
+def _simplify_profile(
+    df: pd.DataFrame,
+    tol: float,
+) -> pd.DataFrame:
+    """Merge adjacent layers with similar shear-wave velocity.
+
+    Within each group of consecutive layers whose ``vel_shear`` values are
+    within *tol* of the first layer in the group, the time-averaged shear-wave
+    velocity and the thickness-weighted average unit weight are computed.
+
+    Parameters
+    ----------
+    df : pandas.DataFrame
+        Must contain ``depth``, ``vel_shear``, and ``unit_wt`` columns.
+    tol : float
+        Relative tolerance for grouping (e.g. 0.02 = 2 %).
+
+    Returns
+    -------
+    pandas.DataFrame
+        Simplified profile with the same columns.
+    """
+    depths = df["depth"].values
+    vs = df["vel_shear"].values
+    unit_wt = df["unit_wt"].values
+
+    # Compute layer thicknesses (last layer gets 0)
+    thicknesses = np.diff(depths, append=depths[-1])
+
+    groups: list[list[int]] = []
+    current_group: list[int] = [0]
+
+    for i in range(1, len(vs)):
+        ref_vs = vs[current_group[0]]
+        if abs(vs[i] - ref_vs) / ref_vs <= tol:
+            current_group.append(i)
+        else:
+            groups.append(current_group)
+            current_group = [i]
+    groups.append(current_group)
+
+    rows = []
+    for group in groups:
+        depth_top = depths[group[0]]
+        thick = thicknesses[group]
+        total_thick = thick.sum()
+
+        if total_thick > 0:
+            # Time-average velocity: Vs_avg = total_thickness / sum(thickness_i / Vs_i)
+            vs_avg = total_thick / np.sum(thick / vs[group])
+        else:
+            # Halfspace (last layer)
+            vs_avg = vs[group[0]]
+
+        unit_wt_avg = np.average(unit_wt[group], weights=np.maximum(thick, 1e-10))
+        rows.append({"depth": depth_top, "vel_shear": vs_avg, "unit_wt": unit_wt_avg})
+
+    return pd.DataFrame(rows)

examples/example-17.ipynb

@@ -1661,6 +1661,20 @@ def strain(self):
 
         return value
 
+    def _compute_damping(self, strain) -> float:
+        """Compute layer-adjusted damping at the given strain.
+
+        The soil type's minimum damping is replaced with the layer-specific
+        minimum damping.
+        """
+        try:
+            damping = self.soil_type.damping(strain)
+            damping -= self.soil_type.damping_min
+        except TypeError:
+            damping = 0.0
+
+        return damping + self.damping_min
+
     @strain.setter
     def strain(self, strain):
         if self.soil_type.is_nonlinear:
@@ -1676,19 +1690,22 @@ def strain(self, strain):
 
         self._shear_mod.value = self.initial_shear_mod * mod_reduc
 
-        try:
-            # Interpolate the damping at the strain, and then reduce by the
-            # minimum damping
-            damping = self.soil_type.damping(strain)
-            damping -= self.soil_type.damping_min
-        except TypeError:
-            # No iteration provided by damping
-            damping = 0
+        # Update the damping value
+        self._damping.value = self._compute_damping(strain)
 
-        # Add the layer-specific minimum damping
-        damping += self.damping_min
+    @property
+    def adjusted_damping_curves(self) -> np.recarray:
+        """Return the damping curve adjusted by the layer-specific minimum damping."""
+
+        if isinstance(self.soil_type.damping, (float, int)):
+            # No iteration provided by damping
+            strains = np.asarray([np.nan])
+            values = np.asarray([self.damping_min])
+        else:
+            strains = np.asarray(self.soil_type.damping.strains)
+            values = np.array([self._compute_damping(s) for s in strains])
 
-        self._damping.value = damping
+        return np.rec.array((strains, values), names=["strain", "damping"])
 
     @property
     def soil_type(self):

src/pystrata/generic.py

@@ -468,3 +468,70 @@ def adjust_damping_values(
     profile.reset_layers()
 
     return profile, site_atten_scatter
+
+
+def calc_mean_eff_stress(
+    depths: npt.ArrayLike,
+    unit_wt: npt.ArrayLike,
+    water_table_depth: float,
+    k0: float = 0.5,
+) -> np.ndarray:
+    """Compute mean effective stress at layer midpoints.
+
+    Calculates the total vertical stress from the overburden, subtracts pore
+    water pressure below the water table, and converts to mean effective stress
+    using the at-rest earth pressure coefficient.
+
+    Parameters
+    ----------
+    depths : array_like
+        Depth to the top of each layer [m].
+    unit_wt : array_like
+        Unit weight of each layer [kN/m³]. Values outside the typical range of
+        12–30 kN/m³ trigger a warning.
+    water_table_depth : float
+        Depth to the water table from the surface [m].
+    k0 : float, optional
+        At-rest earth pressure coefficient (default 0.5).
+
+    Returns
+    -------
+    np.ndarray
+        Mean effective stress at the midpoint of each layer [kN/m²].
+
+    Warns
+    -----
+    UserWarning
+        If any ``unit_wt`` values fall outside 12–30 kN/m³, which may indicate
+        incorrect units (e.g., kg/m³ instead of kN/m³).
+    """
+    import warnings
+
+    depths = np.asarray(depths, dtype=float)
+    unit_wt = np.asarray(unit_wt, dtype=float)
+
+    if np.any((unit_wt < 10) | (unit_wt > 30)):
+        warnings.warn(
+            "Some unit_wt values are outside the typical range of 10–30 kN/m³. "
+            "Ensure values are in kN/m³ (not kg/m³ or pcf).",
+            stacklevel=2,
+        )
+
+    UNIT_WT_WATER = 9.81  # kN/m³
+
+    # Layer thicknesses (last layer = halfspace with 0 thickness)
+    thicknesses = np.diff(depths, append=depths[-1])
+
+    # Total vertical stress at layer midpoints
+    depth_mids = depths + thicknesses / 2
+    stress_increments = unit_wt * thicknesses
+    stress_vert_total = np.cumsum(stress_increments) - stress_increments / 2
+
+    # Pore water pressure (zero above water table)
+    pore_pressure = np.maximum(0.0, UNIT_WT_WATER * (depth_mids - water_table_depth))
+
+    # Effective vertical stress -> mean effective stress
+    stress_vert_eff = stress_vert_total - pore_pressure
+    stress_mean = stress_vert_eff * (1 + 2 * k0) / 3
+
+    return stress_mean