Discrete transform

CHANGELOG.md

@@ -15,6 +15,7 @@ All notable changes to **GSTools** will be documented in this file.
   - Universal
   - External Drift Kriging
   - Detrended Kriging
+- a new transformation function for discrete fields has been added #70
 
 ### Changes
 - Python versions 2.7 and 3.4 are no longer supported #40 #43

examples/07_transformations/02_discrete.py

@@ -0,0 +1,20 @@
+"""
+discrete fields
+-------------
+
+Here we transform a field to a discrete field with five values.
+If we do not give thresholds, the pairwise means of the given
+values are taken as thresholds.
+"""
+import numpy as np
+import gstools as gs
+
+# structured field with a size of 100x100 and a grid-size of 1x1
+x = y = range(100)
+model = gs.Gaussian(dim=2, var=1, len_scale=10)
+srf = gs.SRF(model, seed=20170519)
+srf.structured([x, y])
+# create 5 equidistanly spaced values
+discrete_values = np.linspace(np.min(srf.field), np.max(srf.field), 5)
+gs.transform.discrete(srf, discrete_values)
+srf.plot()

examples/07_transformations/README.rst

@@ -12,6 +12,7 @@ common transformations:
 
 .. autosummary::
    binary
+   discrete
    boxcox
    zinnharvey
    normal_force_moments

gstools/transform/__init__.py

@@ -9,6 +9,7 @@
 
 .. autosummary::
    binary
+   discrete
    boxcox
    zinnharvey
    normal_force_moments
@@ -22,6 +23,7 @@
 
 from gstools.transform.field import (
     binary,
+    discrete,
     boxcox,
     zinnharvey,
     normal_force_moments,
@@ -33,6 +35,7 @@
 
 __all__ = [
     "binary",
+    "discrete",
     "boxcox",
     "zinnharvey",
     "normal_force_moments",

gstools/transform/field.py

@@ -8,6 +8,7 @@
 
 .. autosummary::
    binary
+   discrete
    boxcox
    zinnharvey
    normal_force_moments
@@ -26,6 +27,7 @@
 
 __all__ = [
     "binary",
+    "discrete",
     "boxcox",
     "zinnharvey",
     "normal_force_moments",
@@ -67,6 +69,56 @@ def binary(fld, divide=None, upper=None, lower=None):
         fld.field[fld.field <= divide] = lower
 
 
+def discrete(fld, values, thresholds=None):
+    """
+    Discrete transformation.
+
+    After this transformation, the field has only `len(values)` discrete
+    values.
+
+    Parameters
+    ----------
+    fld : :any:`Field`
+        Spatial Random Field class containing a generated field.
+        Field will be transformed inplace.
+    values : :any:`np.ndarray`
+        The discrete values the field will take
+    thresholds : :class:`numpy.ndarray`, optional
+        the thresholds, where the value classes are separated
+        Default: mean of the neighbouring values
+    """
+    if fld.field is None:
+        print("discrete: no field stored in SRF class.")
+    else:
+        if thresholds is None:
+            # just in case, sort the values
+            values = np.sort(values)
+            thresholds = (values[1:] + values[:-1]) / 2
+        else:
+            if len(values) != len(thresholds) + 1:
+                raise ValueError(
+                    "discrete transformation: len(values) != len(thresholds) + 1"
+                )
+            values = np.array(values)
+            thresholds = np.array(thresholds)
+            for i in range(len(thresholds)):
+                if not (values[i] <= thresholds[i] < values[i + 1]):
+                    raise ValueError(
+                        "discrete transformation: thresholds must lie between values"
+                    )
+
+        # handle edge cases
+        fld.field[fld.field <= thresholds[0]] = values[0]
+        fld.field[fld.field > thresholds[-1]] = values[-1]
+
+        for i in range(len(values[:-2])):
+            fld.field[
+                np.logical_and(
+                    thresholds[i] <= fld.field, fld.field < thresholds[i + 1]
+                )
+            ] = values[i + 1]
+
+
 def boxcox(fld, lmbda=1, shift=0):
     """
     Box-Cox transformation.

tests/test_srf.py

@@ -250,6 +250,13 @@ def test_transform(self):
         tf.binary(srf)
         srf((self.x_grid, self.y_grid), seed=self.seed, mesh_type="structured")
         tf.boxcox(srf)
+        srf((self.x_grid, self.y_grid), seed=self.seed, mesh_type="structured")
+        values = np.linspace(np.min(srf.field), np.max(srf.field), 3)
+        tf.discrete(srf, values)
+        srf((self.x_grid, self.y_grid), seed=self.seed, mesh_type="structured")
+        values = [-1, 0, 1]
+        thresholds = [-0.9, 0.1]
+        tf.discrete(srf, values, thresholds)
 
     def test_incomprrandmeth(self):
         self.cov_model = Gaussian(dim=2, var=0.5, len_scale=1.0)