From 291ec90dfbb4e676de902d1e869caeb4ac199997 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Sebastian=20M=C3=BCller?= <mueller.seb@posteo.de>
Date: Sun, 14 Dec 2025 00:24:56 +0100
Subject: [PATCH 01/15] JBessel: better cor calculation

---
 src/gstools/covmodel/models.py | 7 +------
 1 file changed, 1 insertion(+), 6 deletions(-)

diff --git a/src/gstools/covmodel/models.py b/src/gstools/covmodel/models.py
index 3c47b351..5e5a5f52 100644
--- a/src/gstools/covmodel/models.py
+++ b/src/gstools/covmodel/models.py
@@ -1004,12 +1004,7 @@ def check_opt_arg(self):
     def cor(self, h):
         """J-Bessel correlation."""
         h = np.asarray(h, dtype=np.double)
-        h_gz = np.logical_not(np.isclose(h, 0))
-        hh = h[h_gz]
-        res = np.ones_like(h)
-        nu = self.nu
-        res[h_gz] = sps.gamma(nu + 1) * sps.jv(nu, hh) / (hh / 2.0) ** nu
-        return res
+        return sps.hyp0f1(self.nu + 1, -0.25 * h**2)
 
     def spectral_density(self, k):  # noqa: D102
         k = np.asarray(k, dtype=np.double)

From 1bdafc2d86f6556574c2479f253c84bf351e0260 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Sebastian=20M=C3=BCller?= <mueller.seb@posteo.de>
Date: Sun, 14 Dec 2025 00:25:51 +0100
Subject: [PATCH 02/15] tools: move derivative to tools.misc

---
 src/gstools/normalizer/base.py |  7 ++-----
 src/gstools/tools/misc.py      | 23 ++++++++++++++++++++++-
 2 files changed, 24 insertions(+), 6 deletions(-)

diff --git a/src/gstools/normalizer/base.py b/src/gstools/normalizer/base.py
index 2aedd4d9..0fa6f0b4 100644
--- a/src/gstools/normalizer/base.py
+++ b/src/gstools/normalizer/base.py
@@ -14,10 +14,7 @@
 import numpy as np
 import scipy.optimize as spo
 
-
-def _derivative(f, x, dx=1e-6):
-    """Central difference formula."""
-    return (f(x + dx) - f(x - dx)) / (2 * dx)
+from gstools.tools.misc import derivative
 
 
 class Normalizer:
@@ -60,7 +57,7 @@ def _normalize(self, data):
         return data
 
     def _derivative(self, data):
-        return _derivative(self._normalize, data, dx=self._dx)
+        return derivative(self._normalize, data, dx=self._dx)
 
     def _loglikelihood(self, data):
         add = -0.5 * np.size(data) * (np.log(2 * np.pi) + 1)
diff --git a/src/gstools/tools/misc.py b/src/gstools/tools/misc.py
index 3fe44678..90f6bad4 100755
--- a/src/gstools/tools/misc.py
+++ b/src/gstools/tools/misc.py
@@ -15,7 +15,28 @@
 
 from gstools.tools.geometric import format_struct_pos_dim, generate_grid
 
-__all__ = ["get_fig_ax", "list_format", "eval_func"]
+__all__ = ["derivative", "get_fig_ax", "list_format", "eval_func"]
+
+
+def derivative(f, x, dx=1e-6):
+    """Central difference formula.
+
+    Parameters
+    ----------
+    f : :any:`callable`
+        Function to differentiate.
+    x : array_like
+        Point(s) where to evaluate the derivative.
+    dx : :class:`float`, optional
+        Step size for the central difference. The default is 1e-6.
+
+    Returns
+    -------
+    :class:`numpy.ndarray`
+        Derivative of f at x.
+    """
+    x = np.asarray(x, dtype=np.double)
+    return (f(x + dx) - f(x - dx)) / (2 * dx)
 
 
 def get_fig_ax(fig=None, ax=None, ax_name="rectilinear"):  # pragma: no cover

From a37403084b63295c8cff8f2b6de1b7e0dc03e501 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Sebastian=20M=C3=BCller?= <mueller.seb@posteo.de>
Date: Sun, 14 Dec 2025 00:26:31 +0100
Subject: [PATCH 03/15] CovModel: add calc_roughness method and roughness
 property

---
 src/gstools/covmodel/base.py | 48 ++++++++++++++++++++++++++++++++++++
 1 file changed, 48 insertions(+)

diff --git a/src/gstools/covmodel/base.py b/src/gstools/covmodel/base.py
index 8e3b2ae3..0e459431 100644
--- a/src/gstools/covmodel/base.py
+++ b/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-3, dx=1e-6):
+        """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-3``
+        dx : :class:`float`, optional
+            Step size for the derivative calculation.
+            Default: ``1e-6``
+
+        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)
+
     def __eq__(self, other):
         """Compare CovModels."""
         if not isinstance(other, CovModel):

From a4c1b35156eec1efafa6c5e1f8b048b2de481d97 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Sebastian=20M=C3=BCller?= <mueller.seb@posteo.de>
Date: Sun, 14 Dec 2025 00:43:02 +0100
Subject: [PATCH 04/15] CovModel: add _roughness for all implemented models

---
 src/gstools/covmodel/models.py     | 38 ++++++++++++++++++++++++++++++
 src/gstools/covmodel/tpl_models.py | 12 ++++++++++
 2 files changed, 50 insertions(+)

diff --git a/src/gstools/covmodel/models.py b/src/gstools/covmodel/models.py
index 5e5a5f52..fcbd1a95 100644
--- a/src/gstools/covmodel/models.py
+++ b/src/gstools/covmodel/models.py
@@ -188,6 +188,8 @@ def _has_ppf(self):
     def calc_integral_scale(self):  # noqa: D102
         return self.len_rescaled * np.sqrt(np.pi) / 2.0
 
+    def _roughness(self):
+        return 2.0
 
 class Exponential(CovModel):
     r"""The Exponential covariance model.
@@ -272,6 +274,9 @@ def _has_ppf(self):
     def calc_integral_scale(self):  # noqa: D102
         return self.len_rescaled
 
+    def _roughness(self):
+        return 1.0
+
 
 class Stable(CovModel):
     r"""The stable covariance model.
@@ -347,6 +352,9 @@ def cor(self, h):
     def calc_integral_scale(self):  # noqa: D102
         return self.len_rescaled * sps.gamma(1.0 + 1.0 / self.alpha)
 
+    def _roughness(self):
+        return self.alpha
+
 
 class Matern(CovModel):
     r"""The Matérn covariance model.
@@ -457,6 +465,9 @@ def calc_integral_scale(self):  # noqa: D102
             / sps.beta(self.nu, 0.5)
         )
 
+    def _roughness(self):
+        return min(2.0, 2 * self.nu)
+
 
 class Integral(CovModel):
     r"""The Exponential Integral covariance model.
@@ -548,6 +559,9 @@ def calc_integral_scale(self):  # noqa: D102
             self.len_rescaled * self.nu * np.sqrt(np.pi) / (2 * self.nu + 2.0)
         )
 
+    def _roughness(self):
+        return min(2.0, self.nu)
+
 
 class Rational(CovModel):
     r"""The rational quadratic covariance model.
@@ -620,6 +634,9 @@ def calc_integral_scale(self):  # noqa: D102
             / 2.0
         )
 
+    def _roughness(self):
+        return 2.0
+
 
 class Cubic(CovModel):
     r"""The Cubic covariance model.
@@ -656,6 +673,9 @@ def cor(self, h):
         h = np.minimum(np.abs(h, dtype=np.double), 1.0)
         return 1.0 - 7 * h**2 + 8.75 * h**3 - 3.5 * h**5 + 0.75 * h**7
 
+    def _roughness(self):
+        return 2.0
+
 
 class Linear(CovModel):
     r"""The bounded linear covariance model.
@@ -692,6 +712,9 @@ def check_dim(self, dim):
         """Linear model is only valid in 1D."""
         return dim < 2
 
+    def _roughness(self):
+        return 1.0
+
 
 class Circular(CovModel):
     r"""The circular covariance model.
@@ -741,6 +764,9 @@ def check_dim(self, dim):
         """Circular model is only valid in 1D and 2D."""
         return dim < 3
 
+    def _roughness(self):
+        return 1.0
+
 
 class Spherical(CovModel):
     r"""The Spherical covariance model.
@@ -780,6 +806,9 @@ def check_dim(self, dim):
         """Spherical model is only valid in 1D, 2D and 3D."""
         return dim < 4
 
+    def _roughness(self):
+        return 1.0
+
 
 class HyperSpherical(CovModel):
     r"""The Hyper-Spherical covariance model.
@@ -841,6 +870,9 @@ def spectral_density(self, k):  # noqa: D102
         )
         return res
 
+    def _roughness(self):
+        return 1.0
+
 
 class SuperSpherical(CovModel):
     r"""The Super-Spherical covariance model.
@@ -917,6 +949,9 @@ def cor(self, h):
         )
         return res
 
+    def _roughness(self):
+        return 1.0
+
 
 class JBessel(CovModel):
     r"""The J-Bessel hole model.
@@ -1020,3 +1055,6 @@ def spectral_density(self, k):  # noqa: D102
             * (1.0 - (kk * self.len_rescaled) ** 2) ** (self.nu - self.dim / 2)
         )
         return res
+
+    def _roughness(self):
+        return 2.0
diff --git a/src/gstools/covmodel/tpl_models.py b/src/gstools/covmodel/tpl_models.py
index 46472d1b..7b6e7d34 100644
--- a/src/gstools/covmodel/tpl_models.py
+++ b/src/gstools/covmodel/tpl_models.py
@@ -211,6 +211,9 @@ def spectral_density(self, k):  # noqa: D102
             k, self.dim, self.len_rescaled, self.hurst, self.len_low_rescaled
         )
 
+    def _roughness(self):
+        return 2 * self.hurst
+
 
 class TPLExponential(TPLCovModel):
     r"""Truncated-Power-Law with Exponential modes.
@@ -344,6 +347,9 @@ def spectral_density(self, k):  # noqa: D102
             k, self.dim, self.len_rescaled, self.hurst, self.len_low_rescaled
         )
 
+    def _roughness(self):
+        return 2 * self.hurst
+
 
 class TPLStable(TPLCovModel):
     r"""Truncated-Power-Law with Stable modes.
@@ -499,6 +505,9 @@ def correlation(self, r):
             - self.len_low_rescaled ** (2 * self.hurst)
         )
 
+    def _roughness(self):
+        return 2 * self.hurst
+
 
 class TPLSimple(CovModel):
     r"""The simply truncated power law model.
@@ -573,3 +582,6 @@ def default_opt_arg_bounds(self):
     def cor(self, h):
         """TPL Simple - normalized correlation function."""
         return np.maximum(1 - np.abs(h, dtype=np.double), 0.0) ** self.nu
+
+    def _roughness(self):
+        return 1.0

From 3ecc53eeae96d5517259cd007725b11bb0697d3f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Sebastian=20M=C3=BCller?= <mueller.seb@posteo.de>
Date: Mon, 15 Dec 2025 09:55:56 +0100
Subject: [PATCH 05/15] roughness: use forward difference near origin for
 roughness estimate

---
 src/gstools/covmodel/base.py   | 8 ++++----
 src/gstools/covmodel/models.py | 1 +
 src/gstools/tools/misc.py      | 7 ++++++-
 3 files changed, 11 insertions(+), 5 deletions(-)

diff --git a/src/gstools/covmodel/base.py b/src/gstools/covmodel/base.py
index 0e459431..4e28bcd7 100644
--- a/src/gstools/covmodel/base.py
+++ b/src/gstools/covmodel/base.py
@@ -1179,7 +1179,7 @@ def roughness(self):
             return self._roughness()
         return self.calc_roughness()
 
-    def calc_roughness(self, x=1e-3, dx=1e-6):
+    def calc_roughness(self, x=1e-4, dx=1e-4):
         """Calculate the roughness of the model.
 
         This ignores the nugget of the model.
@@ -1188,10 +1188,10 @@ def calc_roughness(self, x=1e-3, dx=1e-6):
         ----------
         x : :class:`float`, optional
             Point at which the derivative is calculated.
-            Default: ``1e-3``
+            Default: ``1e-4``
         dx : :class:`float`, optional
             Step size for the derivative calculation.
-            Default: ``1e-6``
+            Default: ``1e-4``
 
         Returns
         -------
@@ -1213,7 +1213,7 @@ def f(h):
             """Function for derivative calculation."""
             return np.log(1 - self.cor(np.exp(h)))
 
-        return derivative(f, np.log(x), dx=dx)
+        return derivative(f, np.log(x), dx=dx, order=1)
 
     def __eq__(self, other):
         """Compare CovModels."""
diff --git a/src/gstools/covmodel/models.py b/src/gstools/covmodel/models.py
index fcbd1a95..319f6a1e 100644
--- a/src/gstools/covmodel/models.py
+++ b/src/gstools/covmodel/models.py
@@ -191,6 +191,7 @@ def calc_integral_scale(self):  # noqa: D102
     def _roughness(self):
         return 2.0
 
+
 class Exponential(CovModel):
     r"""The Exponential covariance model.
 
diff --git a/src/gstools/tools/misc.py b/src/gstools/tools/misc.py
index 90f6bad4..8bc74eb7 100755
--- a/src/gstools/tools/misc.py
+++ b/src/gstools/tools/misc.py
@@ -18,7 +18,7 @@
 __all__ = ["derivative", "get_fig_ax", "list_format", "eval_func"]
 
 
-def derivative(f, x, dx=1e-6):
+def derivative(f, x, dx=1e-6, order=2):
     """Central difference formula.
 
     Parameters
@@ -29,6 +29,9 @@ def derivative(f, x, dx=1e-6):
         Point(s) where to evaluate the derivative.
     dx : :class:`float`, optional
         Step size for the central difference. The default is 1e-6.
+    order : :class:`int`, optional
+        Order of the derivative approximation. Either 1 (forward difference) or
+        2 (central difference). The default is 2.
 
     Returns
     -------
@@ -36,6 +39,8 @@ def derivative(f, x, dx=1e-6):
         Derivative of f at x.
     """
     x = np.asarray(x, dtype=np.double)
+    if order == 1:
+        return (f(x + dx) - f(x)) / dx
     return (f(x + dx) - f(x - dx)) / (2 * dx)
 
 

From 1a12d5daabfcf0baaae136f8ecd97d0ccfddb44e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Sebastian=20M=C3=BCller?= <mueller.seb@posteo.de>
Date: Fri, 26 Dec 2025 13:28:18 +0100
Subject: [PATCH 06/15] TPL: set correct roughness when len_low > 0

---
 src/gstools/covmodel/tpl_models.py | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/src/gstools/covmodel/tpl_models.py b/src/gstools/covmodel/tpl_models.py
index 7b6e7d34..bc894c77 100644
--- a/src/gstools/covmodel/tpl_models.py
+++ b/src/gstools/covmodel/tpl_models.py
@@ -212,6 +212,8 @@ def spectral_density(self, k):  # noqa: D102
         )
 
     def _roughness(self):
+        if self.len_low > 0:
+            return 2.0  # roughness of gaussian model
         return 2 * self.hurst
 
 
@@ -348,6 +350,8 @@ def spectral_density(self, k):  # noqa: D102
         )
 
     def _roughness(self):
+        if self.len_low > 0:
+            return 1.0  # roughness of exponential model
         return 2 * self.hurst
 
 
@@ -506,6 +510,8 @@ def correlation(self, r):
         )
 
     def _roughness(self):
+        if self.len_low > 0:
+            return self.alpha  # roughness of stable model
         return 2 * self.hurst
 
 

From e742117421021e713a2fae8a6998958e4bd7914b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Sebastian=20M=C3=BCller?= <mueller.seb@posteo.de>
Date: Fri, 26 Dec 2025 13:30:34 +0100
Subject: [PATCH 07/15] CI: fix macos runner images

---
 .github/workflows/main.yml | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml
index 79993388..5805ff9b 100644
--- a/.github/workflows/main.yml
+++ b/.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" }
@@ -68,11 +68,11 @@ jobs:
         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

From 76983ce1b3ccb22f213a0b35430bd5f8a8f086a1 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Sebastian=20M=C3=BCller?= <mueller.seb@posteo.de>
Date: Sat, 27 Dec 2025 00:49:29 +0100
Subject: [PATCH 08/15] CI: test py314

---
 .github/workflows/main.yml | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml
index 5805ff9b..da78b645 100644
--- a/.github/workflows/main.yml
+++ b/.github/workflows/main.yml
@@ -64,7 +64,8 @@ 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" }

From 31994c8fd826275a92d8f65cc29cb9a9b9573433 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Sebastian=20M=C3=BCller?= <mueller.seb@posteo.de>
Date: Sat, 27 Dec 2025 00:49:47 +0100
Subject: [PATCH 09/15] pyproject: update specifiers

---
 pyproject.toml | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/pyproject.toml b/pyproject.toml
index 1dae8c0b..6edfb1bf 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -12,7 +12,7 @@ authors = [
     { name = "Sebastian Müller, Lennart Schüler", email = "info@geostat-framework.org" },
 ]
 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",

From 705d7bfd50c4ace677dee1fb83efd3009652b4c7 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Sebastian=20M=C3=BCller?= <mueller.seb@posteo.de>
Date: Sat, 27 Dec 2025 01:05:54 +0100
Subject: [PATCH 10/15] tests: add tests for roughness attribute

---
 tests/test_covmodel.py | 18 ++++++++++++++++++
 1 file changed, 18 insertions(+)

diff --git a/tests/test_covmodel.py b/tests/test_covmodel.py
index 62ea8a1c..bf079472 100644
--- a/tests/test_covmodel.py
+++ b/tests/test_covmodel.py
@@ -377,6 +377,24 @@ class Gau_err(CovModel):
         self.assertAlmostEqual(cor, model_cov.cor(2.5))
         self.assertAlmostEqual(cor, model_cor.cor(2.5))
 
+    def test_roughness(self):
+        self.assertAlmostEqual(Gaussian().roughness, 2.0)
+        self.assertAlmostEqual(Exponential().roughness, 1.0)
+        self.assertAlmostEqual(Stable(alpha=1.2).roughness, 1.2)
+        self.assertAlmostEqual(Matern(nu=0.4).roughness, 0.8)
+        self.assertAlmostEqual(Integral(nu=0.6).roughness, 0.6)
+        self.assertAlmostEqual(TPLGaussian(hurst=0.4).roughness, 0.8)
+        self.assertAlmostEqual(
+            TPLGaussian(hurst=0.4, len_low=1.0).roughness, 2.0
+        )
+        self.assertAlmostEqual(
+            TPLStable(hurst=0.4, alpha=1.2, len_low=1.0).roughness,
+            1.2,
+        )
+        self.assertAlmostEqual(Gaussian(nugget=0.5).roughness, 0.0)
+        self.assertTrue(np.isinf(Gaussian(var=0.0).roughness))
+        self.assertAlmostEqual(Gau_cor().roughness, 2.0, places=2)
+
     def test_rescale(self):
         model1 = Exponential()
         model2 = Exponential(rescale=2.1)

From e5edf27b6044475cb28640f85ecef189387a99d3 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Sebastian=20M=C3=BCller?= <mueller.seb@posteo.de>
Date: Sat, 27 Dec 2025 01:36:40 +0100
Subject: [PATCH 11/15] examples: add example for roughness information of a
 model

---
 examples/02_cov_model/07_roughness.py | 127 ++++++++++++++++++++++++++
 1 file changed, 127 insertions(+)
 create mode 100644 examples/02_cov_model/07_roughness.py

diff --git a/examples/02_cov_model/07_roughness.py b/examples/02_cov_model/07_roughness.py
new file mode 100644
index 00000000..2ba27433
--- /dev/null
+++ b/examples/02_cov_model/07_roughness.py
@@ -0,0 +1,127 @@
+r"""
+Roughness
+=========
+
+The roughness describes the power-law behavior of a variogram at the origin
+([Wu2016]_):
+
+.. math::
+   \gamma(r) \sim c \cdot r^\alpha
+
+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.
+
+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.
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+import gstools as gs
+
+models = [
+    gs.Exponential(len_scale=1.0),
+    gs.Gaussian(len_scale=1.0),
+    gs.Stable(len_scale=1.0, alpha=0.7),
+]
+
+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
+# ---------------------------
+#
+# 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)`.
+
+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
+
+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.
+
+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)
+
+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()
+
+###############################################################################
+# 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.

From 524b5f15ca02aa700710eeb62d735e33f4813405 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Sebastian=20M=C3=BCller?= <mueller.seb@posteo.de>
Date: Sat, 27 Dec 2025 11:24:13 +0100
Subject: [PATCH 12/15] update roughness example

---
 examples/02_cov_model/07_roughness.py | 57 +++++++++++++++++++++------
 1 file changed, 46 insertions(+), 11 deletions(-)

diff --git a/examples/02_cov_model/07_roughness.py b/examples/02_cov_model/07_roughness.py
index 2ba27433..ae7719ba 100644
--- a/examples/02_cov_model/07_roughness.py
+++ b/examples/02_cov_model/07_roughness.py
@@ -9,17 +9,10 @@
    \gamma(r) \sim c \cdot r^\alpha
 
 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.
-
-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.
+: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
@@ -27,12 +20,22 @@
 
 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()
 
@@ -59,9 +62,21 @@
 # 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]
@@ -69,6 +84,9 @@
 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]
 
@@ -78,6 +96,9 @@
 #    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
@@ -85,6 +106,9 @@
     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}$"
@@ -114,6 +138,7 @@
 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
@@ -125,3 +150,13 @@
 # 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.

From 4ea4c9aadc4725c211129b7cd2b58f9b81f24151 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Sebastian=20M=C3=BCller?= <mueller.seb@posteo.de>
Date: Sat, 27 Dec 2025 14:43:51 +0100
Subject: [PATCH 13/15] JBessel: add analytical integral scale

---
 src/gstools/covmodel/models.py | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/src/gstools/covmodel/models.py b/src/gstools/covmodel/models.py
index 319f6a1e..f2059b12 100644
--- a/src/gstools/covmodel/models.py
+++ b/src/gstools/covmodel/models.py
@@ -1057,5 +1057,11 @@ def spectral_density(self, k):  # noqa: D102
         )
         return res
 
+    def calc_integral_scale(self):  # noqa: D102
+        return (
+            self.len_rescaled
+            * np.sqrt(np.pi) * sps.gamma(self.nu + 1) / sps.gamma(self.nu + 0.5)
+        )
+
     def _roughness(self):
         return 2.0

From 99f4bffb8646463637ce2f636fb7a0ac608c3eb6 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Sebastian=20M=C3=BCller?= <mueller.seb@posteo.de>
Date: Sat, 27 Dec 2025 14:52:19 +0100
Subject: [PATCH 14/15] black

---
 src/gstools/covmodel/models.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/src/gstools/covmodel/models.py b/src/gstools/covmodel/models.py
index f2059b12..d0bdcbf0 100644
--- a/src/gstools/covmodel/models.py
+++ b/src/gstools/covmodel/models.py
@@ -1060,7 +1060,9 @@ def spectral_density(self, k):  # noqa: D102
     def calc_integral_scale(self):  # noqa: D102
         return (
             self.len_rescaled
-            * np.sqrt(np.pi) * sps.gamma(self.nu + 1) / sps.gamma(self.nu + 0.5)
+            * np.sqrt(np.pi)
+            * sps.gamma(self.nu + 1)
+            / sps.gamma(self.nu + 0.5)
         )
 
     def _roughness(self):

From d410fc1c1ce2ca8465a06c8c833d8bece8ddd169 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Sebastian=20M=C3=BCller?= <mueller.seb@posteo.de>
Date: Sat, 27 Dec 2025 14:52:44 +0100
Subject: [PATCH 15/15] more precise roughness def

---
 examples/02_cov_model/07_roughness.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/examples/02_cov_model/07_roughness.py b/examples/02_cov_model/07_roughness.py
index ae7719ba..e392a6b7 100644
--- a/examples/02_cov_model/07_roughness.py
+++ b/examples/02_cov_model/07_roughness.py
@@ -6,7 +6,7 @@
 ([Wu2016]_):
 
 .. math::
-   \gamma(r) \sim c \cdot r^\alpha
+   \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`.