add roughness information to CovModel

.github/workflows/main.yml

@@ -55,7 +55,7 @@ jobs:
       fail-fast: false
       matrix:
         os:
-          [ubuntu-latest, ubuntu-24.04-arm, windows-latest, macos-13, macos-14]
+          [ubuntu-latest, ubuntu-24.04-arm, windows-latest, macos-15-intel, macos-latest]
         # https://github.com/scipy/oldest-supported-numpy/blob/main/setup.cfg
         ver:
           - { py: "3.8", np: "==1.20.0", sp: "==1.5.4" }
@@ -64,15 +64,16 @@ jobs:
           - { py: "3.11", np: "==1.23.2", sp: "==1.9.2" }
           - { py: "3.12", np: "==1.26.2", sp: "==1.11.2" }
           - { py: "3.13", np: "==2.1.0", sp: "==1.14.1" }
-          - { py: "3.13", np: ">=2.1.0", sp: ">=1.14.1" }
+          - { py: "3.14", np: "==2.3.2", sp: "==1.16.1" }
+          - { py: "3.14", np: ">=2.3.2", sp: ">=1.16.1" }
         exclude:
           - os: ubuntu-24.04-arm
             ver: { py: "3.8", np: "==1.20.0", sp: "==1.5.4" }
-          - os: macos-14
+          - os: macos-latest
             ver: { py: "3.8", np: "==1.20.0", sp: "==1.5.4" }
-          - os: macos-14
+          - os: macos-latest
             ver: { py: "3.9", np: "==1.20.0", sp: "==1.5.4" }
-          - os: macos-14
+          - os: macos-latest
             ver: { py: "3.10", np: "==1.21.6", sp: "==1.7.2" }
     steps:
       - uses: actions/checkout@v4

examples/02_cov_model/07_roughness.py

@@ -0,0 +1,162 @@
+r"""
+Roughness
+=========
+
+The roughness describes the power-law behavior of a variogram at the origin
+([Wu2016]_):
+
+.. math::
+   \gamma(r) \sim c \cdot r^\alpha \quad (r \to 0)
+
+The exponent :math:`\alpha` is the roughness information, bounded by
+:math:`0 \le \alpha \le 2`.
+A value of 0 corresponds to a nugget effect and 2 is the smooth limit for random fields.
+You can access it via :any:`CovModel.roughness`.
+On a log-log plot, the slope of the variogram near the origin approaches this value.
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+import gstools as gs
+
+###############################################################################
+# Variogram behavior near the origin
+# ----------------------------------
+#
+# Compare variograms near the origin for models with different roughness.
+
+models = [
+    gs.Exponential(len_scale=1.0),
+    gs.Gaussian(len_scale=1.0),
+    gs.Stable(len_scale=1.0, alpha=0.7),
+]
+
+###############################################################################
+# Use a small-lag grid and fit the slope on a log-log scale to estimate the
+# roughness numerically.
+
+r = np.logspace(-4, -1, 200)
+fig, ax = plt.subplots()
+
+for model in models:
+    gamma = model.variogram(r)
+    slope = np.polyfit(np.log(r[:20]), np.log(gamma[:20]), 1)[0]
+    ax.loglog(
+        r, gamma, label=rf"{model.name} ($\alpha={model.roughness:.2f}$)"
+    )
+    print(f"{model.name}: roughness={model.roughness:.2f}, fit={slope:.2f}")
+
+ax.set_xlabel("r")
+ax.set_ylabel(r"$\gamma(r)$")
+ax.set_title("Variogram near the origin (log-log)")
+ax.legend()
+
+###############################################################################
+# A nugget masks the power-law behavior at the origin, so roughness is 0.
+
+nugget_model = gs.Gaussian(nugget=1.0)
+print("Gaussian with nugget roughness:", nugget_model.roughness)
+
+###############################################################################
+# Roughness and random fields
+# ---------------------------
+#
+# From the theory of fractal stochastic processes (Mandelbrot and Van Ness
+# 1968) ([Mandelbrot1968]_), the roughness can be interpreted as:
+#
+# 1. Persistent (:math:`1 < \alpha \le 2`): smooth behavior, slowly increasing
+#    variogram, long memory (e.g. Gaussian-like).
+# 2. Antipersistent (:math:`0 < \alpha < 1`): rough behavior, steep variogram
+#    near the origin, short memory.
+# 3. No memory (:math:`\alpha = 1`): linear slope at the origin (Exponential).
+#
+# The Integral model lets us control the roughness with its parameter ``nu``.
+# For this model, the roughness is given by :math:`\alpha = \min(2, \nu)`.
+
+###############################################################################
+# Set up a common grid and plotting scales so the realizations are comparable.
+
+sep = 100  # separation point (mid field)
+ext = 10  # field extend
+imext = 2 * [0, ext]
+x = y = np.linspace(0, ext, 2 * sep + 1)
+xm = np.linspace(0, 5, 1001)
+vmin, vmax = -3, 3
+
+###############################################################################
+# Create Integral models with the same integral scale but different roughness.
+
+rough = [0.1, 1, 10]
+mod = [gs.Integral(dim=2, integral_scale=1, nu=nu) for nu in rough]
+
+###############################################################################
+# .. note::
+#
+#    Rough fields require a higher ``mode_no`` so the randomization method
+#    samples the high-frequency part of the spectrum sufficiently well.
+
+###############################################################################
+# Generate random field realizations with a shared seed for fair comparison.
+
+srf = []
+for m in mod:
+    mode_no = int(1000 / m.roughness)  # increase mode_no for rough fields
+    s = gs.SRF(m, seed=20110101, mode_no=mode_no)
+    s.structured((x, y))
+    srf.append(s)
+
+###############################################################################
+# Plot variograms, fields, and cross sections column-wise by roughness.
+
+fig, axes = plt.subplots(3, 3, figsize=(10, 10))
+for i in range(3):
+    label = rf"$\alpha(\gamma)={mod[i].roughness:.2f}$"
+    axes[0, i].plot(xm, mod[i].variogram(xm), label=label)
+    axes[0, i].legend(loc="lower right")
+    axes[0, i].set_ylim([-0.05, 1.05])
+    axes[0, i].set_xlim([-0.25, 5.25])
+    axes[0, i].grid()
+    axes[1, i].imshow(
+        srf[i].field.T,
+        origin="lower",
+        interpolation="bicubic",
+        vmin=vmin,
+        vmax=vmax,
+        extent=imext,
+    )
+    axes[1, i].axhline(y=5, color="k")
+    axes[2, i].plot(x, srf[i].field[:, sep])
+    axes[2, i].set_ylim([vmin, vmax])
+    axes[2, i].set_xlim([0, ext])
+    axes[2, i].grid()
+
+axes[0, 0].set_ylabel(r"$\gamma$")
+axes[1, 0].set_ylabel(r"$y$")
+axes[2, 0].set_ylabel(r"$f$")
+axes[2, 0].set_xlabel(r"$x$")
+axes[2, 1].set_xlabel(r"$x$")
+axes[2, 2].set_xlabel(r"$x$")
+fig.tight_layout()
+# sphinx_gallery_thumbnail_number = 2
+
+###############################################################################
+# Illustration of the impact of the roughness information on spatial random
+# fields. Each column shows the variogram on top, a single field realization in
+# the middle and the highlighted cross section of the field on the bottom. The
+# left column shows a very rough (antipersistent) field with a steep increase
+# of the variogram at the origin, the middle column shows an Exponential like
+# variogram with a linear slope at the origin resulting in a less rough (no
+# memory) field and the right column shows a Gaussian like variogram together
+# with a very smooth (persistent) field. All variograms have the same integral
+# scale and x- and y-axis are given in multiples of the integral scale.
+
+###############################################################################
+# References
+# ----------
+# .. [Wu2016] Wu, W.-Y., and C. Y. Lim. 2016. "ESTIMATION OF SMOOTHNESS OF A
+#    STATIONARY GAUSSIAN RANDOM FIELD." Statistica Sinica 26 (4): 1729-1745.
+#    https://doi.org/10.5705/ss.202014.0109.
+# .. [Mandelbrot1968] Mandelbrot, B. B., and J. W. Van Ness. 1968. "Fractional
+#    Brownian Motions, Fractional Noises and Applications." SIAM Review 10 (4):
+#    422-437. https://doi.org/10.1137/1010093.

pyproject.toml

@@ -12,7 +12,7 @@ authors = [
     { name = "Sebastian Müller, Lennart Schüler", email = "[email protected]" },
 ]
 readme = "README.md"
-license = "LGPL-3.0"
+license = "LGPL-3.0-or-later"
 dynamic = ["version"]
 classifiers = [
     "Development Status :: 5 - Production/Stable",
@@ -34,6 +34,7 @@ classifiers = [
     "Programming Language :: Python :: 3.11",
     "Programming Language :: Python :: 3.12",
     "Programming Language :: Python :: 3.13",
+    "Programming Language :: Python :: 3.14",
     "Topic :: Scientific/Engineering",
     "Topic :: Scientific/Engineering :: GIS",
     "Topic :: Scientific/Engineering :: Hydrology",

src/gstools/covmodel/base.py

@@ -52,6 +52,7 @@
     pos2latlon,
     rotated_main_axes,
 )
+from gstools.tools.misc import derivative
 
 __all__ = ["CovModel"]
 
@@ -1167,6 +1168,53 @@ def is_isotropic(self):
         """:any:`bool`: State if a model is isotropic."""
         return np.all(np.isclose(self.anis, 1.0))
 
+    @property
+    def roughness(self):
+        """:any:`float`: roughness information of the model. Zero for any present nugget."""
+        if self.nugget > 0:
+            return 0.0
+        if np.isclose(self.var, 0):
+            return np.inf
+        if hasattr(self, "_roughness"):
+            return self._roughness()
+        return self.calc_roughness()
+
+    def calc_roughness(self, x=1e-4, dx=1e-4):
+        """Calculate the roughness of the model.
+
+        This ignores the nugget of the model.
+
+        Parameters
+        ----------
+        x : :class:`float`, optional
+            Point at which the derivative is calculated.
+            Default: ``1e-4``
+        dx : :class:`float`, optional
+            Step size for the derivative calculation.
+            Default: ``1e-4``
+
+        Returns
+        -------
+        roughness : :class:`float`
+            Roughness of the model.
+
+        Notes
+        -----
+        The roughness is defined as the derivative of the log-log
+        correlation function at zero:
+
+            * ``roughness = d( log(1 - cor(r)) ) / d( log(r) ) | r->0``
+
+        This is calculated numerically by evaluating the derivative
+        at a small distance `x`.
+        """
+
+        def f(h):
+            """Function for derivative calculation."""
+            return np.log(1 - self.cor(np.exp(h)))
+
+        return derivative(f, np.log(x), dx=dx, order=1)
+
     def __eq__(self, other):
         """Compare CovModels."""
         if not isinstance(other, CovModel):