diff --git a/CHANGELOG.md b/CHANGELOG.md index 5027d6e2..6c9e7984 100755 --- a/CHANGELOG.md +++ b/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 diff --git a/examples/07_transformations/02_discrete.py b/examples/07_transformations/02_discrete.py new file mode 100755 index 00000000..29384da2 --- /dev/null +++ b/examples/07_transformations/02_discrete.py @@ -0,0 +1,37 @@ +""" +discrete fields +------------- + +Here we transform a field to a discrete field with values. +If we do not give thresholds, the pairwise means of the given +values are taken as thresholds. +If thresholds are given, arbitrary values can be applied to the field. +""" +import numpy as np +import gstools as gs + +# structured field with a size of 100x100 and a grid-size of 0.5x0.5 +x = y = np.arange(200) * 0.5 +model = gs.Gaussian(dim=2, var=1, len_scale=5) +srf = gs.SRF(model, seed=20170519) + +# create 5 equidistanly spaced values, thresholds are the arithmetic means +discrete_values = np.linspace(np.min(srf.field), np.max(srf.field), 5) +srf.structured([x, y]) +gs.transform.discrete(srf, discrete_values) +srf.plot() + +# calculate thresholds for equal shares +# but apply different values to the separated classes +discrete_values2 = [0, -1, 2, -3, 4] +srf.structured([x, y]) +gs.transform.discrete(srf, discrete_values2, thresholds="equal") +srf.plot() + +# user defined thresholds +thresholds = [-1, 1] +# apply different values to the separated classes +discrete_values3 = [0, 1, 10] +srf.structured([x, y]) +gs.transform.discrete(srf, discrete_values3, thresholds=thresholds) +srf.plot() diff --git a/examples/07_transformations/02_zinn_harvey.py b/examples/07_transformations/03_zinn_harvey.py similarity index 100% rename from examples/07_transformations/02_zinn_harvey.py rename to examples/07_transformations/03_zinn_harvey.py diff --git a/examples/07_transformations/03_bimodal.py b/examples/07_transformations/04_bimodal.py similarity index 100% rename from examples/07_transformations/03_bimodal.py rename to examples/07_transformations/04_bimodal.py diff --git a/examples/07_transformations/04_combinations.py b/examples/07_transformations/05_combinations.py similarity index 100% rename from examples/07_transformations/04_combinations.py rename to examples/07_transformations/05_combinations.py diff --git a/examples/07_transformations/README.rst b/examples/07_transformations/README.rst index 4d210244..d3a2ba06 100644 --- a/examples/07_transformations/README.rst +++ b/examples/07_transformations/README.rst @@ -12,6 +12,7 @@ common transformations: .. autosummary:: binary + discrete boxcox zinnharvey normal_force_moments diff --git a/gstools/transform/__init__.py b/gstools/transform/__init__.py index 8a5a84c8..1937e6a7 100644 --- a/gstools/transform/__init__.py +++ b/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", diff --git a/gstools/transform/field.py b/gstools/transform/field.py index c3a90e43..63b49e4b 100644 --- a/gstools/transform/field.py +++ b/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", @@ -63,8 +65,72 @@ def binary(fld, divide=None, upper=None, lower=None): divide = fld.mean if divide is None else divide upper = fld.mean + np.sqrt(fld.model.sill) if upper is None else upper lower = fld.mean - np.sqrt(fld.model.sill) if lower is None else lower - fld.field[fld.field > divide] = upper - fld.field[fld.field <= divide] = lower + discrete(fld, [lower, upper], thresholds=[divide]) + + +def discrete(fld, values, thresholds="arithmetic"): + """ + 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:`str` or :any:`np.ndarray`, optional + the thresholds, where the value classes are separated + possible values are: + * "arithmetic": the mean of the 2 neighbouring values + * "equal": devide the field into equal parts + * an array of explicitly given thresholds + Default: "arithmetic" + """ + if fld.field is None: + print("discrete: no field stored in SRF class.") + else: + if thresholds == "arithmetic": + # just in case, sort the values + values = np.sort(values) + thresholds = (values[1:] + values[:-1]) / 2 + elif thresholds == "equal": + values = np.array(values) + n = len(values) + p = np.arange(1, n) / n # n-1 equal subdivisions of [0, 1] + rescale = np.sqrt(fld.model.sill * 2) + # use quantile of the normal distribution to get equal ratios + thresholds = fld.mean + rescale * erfinv(2 * p - 1) + else: + if len(values) != len(thresholds) + 1: + raise ValueError( + "discrete transformation: " + + "len(values) != len(thresholds) + 1" + ) + values = np.array(values) + thresholds = np.array(thresholds) + # check thresholds + if not np.all(thresholds[:-1] < thresholds[1:]): + raise ValueError( + "discrete transformation: " + + "thresholds need to be ascending." + ) + # use a separate result so the intermediate results are not affected + result = np.empty_like(fld.field) + # handle edge cases + result[fld.field <= thresholds[0]] = values[0] + result[fld.field > thresholds[-1]] = values[-1] + for i, value in enumerate(values[1:-1]): + result[ + np.logical_and( + thresholds[i] < fld.field, fld.field <= thresholds[i + 1] + ) + ] = value + # overwrite the field + fld.field = result def boxcox(fld, lmbda=1, shift=0): @@ -269,12 +335,12 @@ def _zinnharvey(field, conn="high", mean=None, var=None): mean = np.mean(field) if var is None: var = np.var(field) - field = np.abs((field - mean) / var) + field = np.abs((field - mean) / np.sqrt(var)) field = 2 * erf(field / np.sqrt(2)) - 1 field = np.sqrt(2) * erfinv(field) if conn == "high": field = -field - return field * var + mean + return field * np.sqrt(var) + mean def _normal_force_moments(field, mean=0, var=1): diff --git a/tests/test_srf.py b/tests/test_srf.py index 81661541..082a7ec3 100644 --- a/tests/test_srf.py +++ b/tests/test_srf.py @@ -250,6 +250,25 @@ 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) + np.testing.assert_array_equal(np.unique(srf.field), [-1, 0, 1]) + + srf((self.x_grid, self.y_grid), seed=self.seed, mesh_type="structured") + values = [-1, 0, 1] + tf.discrete(srf, values, thresholds="arithmetic") + np.testing.assert_array_equal(np.unique(srf.field), [-1.0, 0.0, 1.0]) + + srf((self.x_grid, self.y_grid), seed=self.seed, mesh_type="structured") + values = [-1, 0, 0.5, 1] + tf.discrete(srf, values, thresholds="equal") + np.testing.assert_array_equal(np.unique(srf.field), values) def test_incomprrandmeth(self): self.cov_model = Gaussian(dim=2, var=0.5, len_scale=1.0)