From eb0671d2e0fe0dd1ba5cd6d78825506a3d73d12e Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 21 Dec 2020 14:56:33 +0100 Subject: [PATCH 01/97] normalize: add base class --- gstools/normalize/base.py | 174 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 174 insertions(+) create mode 100644 gstools/normalize/base.py diff --git a/gstools/normalize/base.py b/gstools/normalize/base.py new file mode 100644 index 000000000..5cd547076 --- /dev/null +++ b/gstools/normalize/base.py @@ -0,0 +1,174 @@ +# -*- coding: utf-8 -*- +""" +GStools subpackage providing the base class for normalizers. + +.. currentmodule:: gstools.normalize.base + +The following classes are provided + +.. autosummary:: + Normalizer +""" +import warnings +import numpy as np +import scipy.misc as spm +import scipy.optimize as spo + + +class Normalizer: + """Normalizer class. + + Parameters + ---------- + data : array_like, optional + Input data to fit the transformation in order to gain normality. + The default is None. + **parameter + Specified parameters given by name. If not given, default values + will be applied. + """ + + def __init__(self, data=None, **parameter): + # only use values, that have a provided default value + for key, value in self.default_parameter().items(): + setattr(self, key, parameter.get(key, value)) + if data is not None: + self.fit(data) + # precision for printing + self._prec = 3 + + def default_parameter(self): + """Get default parameters for the transformation. + + Returns + ------- + :class:`dict` + Default parameters. + """ + return {} + + def denormalize(self, values): + """Transform to input distribution. + + Parameters + ---------- + values : array_like + Input values (normal distributed). + + Returns + ------- + :class:`numpy.ndarray` + Denormalized values. + """ + return values + + def normalize(self, values): + """Transform to normal distribution. + + Parameters + ---------- + values : array_like + Input values (not normal distributed). + + Returns + ------- + :class:`numpy.ndarray` + Normalized values. + """ + return values + + def derivative(self, values): + """Factor for normal PDF to gain target PDF. + + Parameters + ---------- + values : array_like + Input values. + + Returns + ------- + :class:`numpy.ndarray` + Derivative of the normalization transformation function. + """ + return spm.derivative(self.normalize, np.asanyarray(values), dx=1e-6) + + def loglikelihood(self, data): + """Log-Likelihood for given data with current parameters. + + Parameters + ---------- + data : array_like + Input data to fit the transformation in order to gain normality. + + Returns + ------- + :class:`float` + Log-Likelihood of the given data. + + Notes + ----- + This loglikelihood function could be missing additive constants, + that are not needed for optimization. + """ + res = -0.5 * np.size(data) * np.log(np.var(self.normalize(data))) + res += np.sum(np.log(np.maximum(1e-16, self.derivative(data)))) + return res + + def fit(self, data, skip=None, **kwargs): + """Fitting the transformation to data by maximizing Log-Likelihood. + + Parameters + ---------- + data : array_like + Input data to fit the transformation in order to gain normality. + skip : :class:`list` of :class:`str` or :any:`None`, optional + Names of parameters to be skiped in fitting. + The default is None. + **kwargs + Keyword arguments passed to :any:`scipy.optimize.minimize_scalar` + when only one parameter present or :any:`scipy.optimize.minimize`. + + Returns + ------- + :class:`dict` + Optimal paramters given by names. + """ + skip = [] if skip is None else skip + all_names = list(self.default_parameter().keys()) + para_names = [p for p in all_names if p not in skip] + + def _neg_llf(par, dat): + for name, val in zip(para_names, np.atleast_1d(par)): + setattr(self, name, val) + return -self.loglikelihood(dat) + + if len(para_names) == 0: # transformations without para. (no opti.) + warnings.warn("Normalizer.fit: no parameters for fitting.") + return {} + if len(para_names) == 1: # one-para. transformations (simple opti.) + # default bracket like in scipy's boxcox (if not given) + kwargs.setdefault("bracket", (-2, 2)) + out = spo.minimize_scalar(_neg_llf, args=(data,), **kwargs) + else: # general case + # init guess from current values (if x0 not given) + kwargs.setdefault("x0", [getattr(self, p) for p in para_names]) + out = spo.minimize(_neg_llf, args=(data,), **kwargs) + # save optimization results + self._opti = out + for name, val in zip(para_names, np.atleast_1d(out.x)): + setattr(self, name, val) + return {name: getattr(self, name) for name in all_names} + + def __str__(self): + """Return String representation.""" + return self.__repr__() + + def __repr__(self): + """Return String representation.""" + out_str = self.__class__.__name__ + para_strs = [] + for p in self.default_parameter(): + para_strs.append( + "{0}={1:.{2}}".format(p, float(getattr(self, p)), self._prec) + ) + return out_str + "(" + ", ".join(para_strs) + ")" From c950ef8a1860559ab2b3cb8d632d62e62b286b97 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 21 Dec 2020 14:57:12 +0100 Subject: [PATCH 02/97] normalize: add normalize transformations --- gstools/normalize/__init__.py | 45 +++++ gstools/normalize/normalizer.py | 337 ++++++++++++++++++++++++++++++++ 2 files changed, 382 insertions(+) create mode 100644 gstools/normalize/__init__.py create mode 100644 gstools/normalize/normalizer.py diff --git a/gstools/normalize/__init__.py b/gstools/normalize/__init__.py new file mode 100644 index 000000000..02fb586ea --- /dev/null +++ b/gstools/normalize/__init__.py @@ -0,0 +1,45 @@ +# -*- coding: utf-8 -*- +""" +GStools subpackage providing normalization routines. + +.. currentmodule:: gstools.normalize + +Base-Normalizer +^^^^^^^^^^^^^^^ + +.. autosummary:: + Normalizer + +Field-Normalizer +^^^^^^^^^^^^^^^^ + +.. autosummary:: + LogNormal + BoxCox + BoxCoxShift + YeoJohnson + Modulus + Manly + +---- +""" + +from gstools.normalize.base import Normalizer +from gstools.normalize.normalizer import ( + LogNormal, + BoxCox, + BoxCoxShift, + YeoJohnson, + Modulus, + Manly, +) + +__all__ = [ + "Normalizer", + "LogNormal", + "BoxCox", + "BoxCoxShift", + "YeoJohnson", + "Modulus", + "Manly", +] diff --git a/gstools/normalize/normalizer.py b/gstools/normalize/normalizer.py new file mode 100644 index 000000000..bbaf8e4ac --- /dev/null +++ b/gstools/normalize/normalizer.py @@ -0,0 +1,337 @@ +# -*- coding: utf-8 -*- +""" +GStools subpackage providing different normalizer transformations. + +.. currentmodule:: gstools.normalize.normalizer + +The following classes are provided + +.. autosummary:: + LogNormal + BoxCox + BoxCoxShift + YeoJohnson + Modulus + Manly +""" +import numpy as np +from gstools.normalize.base import Normalizer + + +class LogNormal(Normalizer): + r"""Log-normal fields. + + Notes + ----- + This parameter-free transformation is given by: + + .. math:: + y=\log(x) + """ + + def denormalize(self, values): + """Transform to log-normal distribution.""" + return np.exp(values) + + def normalize(self, values): + """Transform to normal distribution.""" + return np.log(values) + + def derivative(self, values): + """Factor for normal PDF to gain target PDF.""" + return np.power(values, -1) + + +class BoxCox(Normalizer): + r"""Box-Cox (1964) transformed fields [1]_. + + Parameters + ---------- + data : array_like, optional + Input data to fit the transformation in order to gain normality. + The default is None. + lmbda : :class:`float`, optional + Shape parameter. Default: 1 + + Notes + ----- + This transformation is given by: + + .. math:: + y=\begin{cases} + \frac{x^{\lambda} - 1}{\lambda} & \lambda\neq 0 \\ + \log(x) & \lambda = 0 + \end{cases} + + References + ---------- + .. [1] G.E.P. Box and D.R. Cox, "An Analysis of Transformations", Journal + of the Royal Statistical Society B, 26, 211-252 (1964). + """ + + def default_parameter(self): + """Get default parameters.""" + return {"lmbda": 1} + + def denormalize(self, values): + """Transform to target distribution.""" + if np.isclose(self.lmbda, 0): + return np.exp(values) + return (1 + np.multiply(values, self.lmbda)) ** (1 / self.lmbda) + + def normalize(self, values): + """Transform to normal distribution.""" + if np.isclose(self.lmbda, 0): + return np.log(values) + return (np.power(values, self.lmbda) - 1) / self.lmbda + + def derivative(self, values): + """Factor for normal PDF to gain target PDF.""" + return np.power(values, self.lmbda - 1) + + +class BoxCoxShift(Normalizer): + r"""Box-Cox (1964) transformed fields including shifting [1]_. + + Parameters + ---------- + data : array_like, optional + Input data to fit the transformation in order to gain normality. + The default is None. + lmbda : :class:`float`, optional + Shape parameter. Default: 1 + shift : :class:`float`, optional + Shift parameter. Default: 0 + + Notes + ----- + This transformation is given by: + + .. math:: + y=\begin{cases} + \frac{(x+s)^{\lambda} - 1}{\lambda} & \lambda\neq 0 \\ + \log(x+s) & \lambda = 0 + \end{cases} + + Fitting the shift parameter is rather hard. You should consider skipping + "shift" during fitting: + + >>> data = range(5) + >>> norm = BoxCoxShift(shift=0.5) + >>> norm.fit(data, skip=["shift"]) + {'shift': 0.5, 'lmbda': 0.6747515267420799} + + References + ---------- + .. [1] G.E.P. Box and D.R. Cox, "An Analysis of Transformations", Journal + of the Royal Statistical Society B, 26, 211-252 (1964). + """ + + def default_parameter(self): + """Get default parameters.""" + return {"shift": 0, "lmbda": 1} + + def denormalize(self, values): + """Transform to target distribution.""" + if np.isclose(self.lmbda, 0): + return np.exp(values) - self.shift + return (1 + np.multiply(values, self.lmbda)) ** ( + 1 / self.lmbda + ) - self.shift + + def normalize(self, values): + """Transform to normal distribution.""" + if np.isclose(self.lmbda, 0): + return np.log(np.add(values, self.shift)) + return (np.add(values, self.shift) ** self.lmbda - 1) / self.lmbda + + def derivative(self, values): + """Factor for normal PDF to gain target PDF.""" + return np.power(np.add(values, self.shift), self.lmbda - 1) + + +class YeoJohnson(Normalizer): + r"""Yeo-Johnson (2000) transformed fields [1]_. + + Parameters + ---------- + data : array_like, optional + Input data to fit the transformation in order to gain normality. + The default is None. + lmbda : :class:`float`, optional + Shape parameter. Default: 1 + + Notes + ----- + This transformation is given by: + + .. math:: + y=\begin{cases} + \frac{(x+1)^{\lambda} - 1}{\lambda} + & x\geq 0,\, \lambda\neq 0 \\ + \log(x+1) + & x\geq 0,\, \lambda = 0 \\ + \frac{(|x|+1)^{2-\lambda} - 1}{2-\lambda} + & x<0,\, \lambda\neq 2 \\ + -\log(|x|+1) + & x<0,\, \lambda = 2 + \end{cases} + + + References + ---------- + .. [1] I.K. Yeo and R.A. Johnson, "A new family of power transformations to + improve normality or symmetry." Biometrika, 87(4), pp.954-959, + (2000). + """ + + def default_parameter(self): + """Get default parameters.""" + return {"lmbda": 1} + + def denormalize(self, values): + """Transform to target distribution.""" + values = np.asanyarray(values) + res = np.zeros_like(values, dtype=np.double) + pos = values >= 0 + # when values >= 0 + if np.isclose(self.lmbda, 0): + res[pos] = np.expm1(values[pos]) + else: # self.lmbda != 0 + res[pos] = ( + np.power(values[pos] * self.lmbda + 1, 1 / self.lmbda) - 1 + ) + # when values < 0 + if np.isclose(self.lmbda, 2): + res[~pos] = -np.expm1(-values[~pos]) + else: # self.lmbda != 2 + res[~pos] = 1 - np.power( + -(2 - self.lmbda) * values[~pos] + 1, 1 / (2 - self.lmbda) + ) + return res + + def normalize(self, values): + """Transform to normal distribution.""" + values = np.asanyarray(values) + res = np.zeros_like(values, dtype=np.double) + pos = values >= 0 + # when values >= 0 + if np.isclose(self.lmbda, 0): + res[pos] = np.log1p(values[pos]) + else: # self.lmbda != 0 + res[pos] = (np.power(values[pos] + 1, self.lmbda) - 1) / self.lmbda + # when values < 0 + if np.isclose(self.lmbda, 2): + res[~pos] = -np.log1p(-values[~pos]) + else: # self.lmbda != 2 + res[~pos] = -(np.power(-values[~pos] + 1, 2 - self.lmbda) - 1) / ( + 2 - self.lmbda + ) + return res + + def derivative(self, values): + """Factor for normal PDF to gain target PDF.""" + return (np.abs(values) + 1) ** (np.sign(values) * (self.lmbda - 1)) + + +class Modulus(Normalizer): + r"""Modulus or John-Draper (1980) transformed fields [1]_. + + Parameters + ---------- + data : array_like, optional + Input data to fit the transformation in order to gain normality. + The default is None. + lmbda : :class:`float`, optional + Shape parameter. Default: 1 + + Notes + ----- + This transformation is given by: + + .. math:: + y=\begin{cases} + \mathrm{sgn}(x)\frac{(|x|+1)^{\lambda} - 1}{\lambda} & \lambda\neq 0 \\ + \mathrm{sgn}(x)\log(|x|+1) & \lambda = 0 + \end{cases} + + References + ---------- + .. [1] J. A. John, and N. R. Draper, + “An Alternative Family of Transformations.” Journal + of the Royal Statistical Society C, 29.2, 190-197, (1980) + """ + + def default_parameter(self): + """Get default parameters.""" + return {"lmbda": 1} + + def denormalize(self, values): + """Transform to target distribution.""" + if np.isclose(self.lmbda, 0): + return np.sign(values) * np.expm1(np.abs(values)) + return np.sign(values) * ( + (1 + self.lmbda * np.abs(values)) ** (1 / self.lmbda) - 1 + ) + + def normalize(self, values): + """Transform to normal distribution.""" + if np.isclose(self.lmbda, 0): + return np.sign(values) * np.log1p(np.abs(values)) + return ( + np.sign(values) + * ((np.abs(values) + 1) ** self.lmbda - 1) + / self.lmbda + ) + + def derivative(self, values): + """Factor for normal PDF to gain target PDF.""" + return np.power(np.abs(values) + 1, self.lmbda - 1) + + +class Manly(Normalizer): + r"""Manly (1971) transformed fields [1]_. + + Parameters + ---------- + data : array_like, optional + Input data to fit the transformation in order to gain normality. + The default is None. + lmbda : :class:`float`, optional + Shape parameter. Default: 1 + + Notes + ----- + This transformation is given by: + + .. math:: + y=\begin{cases} + \frac{\exp(\lambda x) - 1}{\lambda} & \lambda\neq 0 \\ + x & \lambda = 0 + \end{cases} + + References + ---------- + .. [1] B. F. J. Manly, "Exponential data transformations.", + Journal of the Royal Statistical Society D, 25.1, 37-42 (1976). + """ + + def default_parameter(self): + """Get default parameters.""" + return {"lmbda": 1} + + def denormalize(self, values): + """Transform to target distribution.""" + if np.isclose(self.lmbda, 0): + return values + return np.log1p(np.multiply(values, self.lmbda)) / self.lmbda + + def normalize(self, values): + """Transform to normal distribution.""" + if np.isclose(self.lmbda, 0): + return values + return np.expm1(np.multiply(values, self.lmbda)) / self.lmbda + + def derivative(self, values): + """Factor for normal PDF to gain target PDF.""" + return np.exp(np.multiply(values, self.lmbda)) From 15efa3f7709358d2986048e8b5d7a45550c799c3 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 21 Dec 2020 15:12:19 +0100 Subject: [PATCH 03/97] normalize: add to gstools init --- gstools/__init__.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/gstools/__init__.py b/gstools/__init__.py index a23a8dfd2..1599bafad 100644 --- a/gstools/__init__.py +++ b/gstools/__init__.py @@ -21,6 +21,7 @@ random tools transform + normalize Classes ======= @@ -106,7 +107,9 @@ vario_estimate_axis """ # Hooray! -from gstools import field, variogram, random, covmodel, tools, krige, transform +from gstools import ( + field, variogram, random, covmodel, tools, krige, transform, normalize +) from gstools.field import SRF from gstools.tools import ( rotated_main_axes, @@ -151,7 +154,7 @@ __all__ = ["__version__"] __all__ += ["covmodel", "field", "variogram", "krige", "random", "tools"] -__all__ += ["transform"] +__all__ += ["transform", "normalize"] __all__ += [ "CovModel", "Gaussian", From 8d2bd32061f4aa5f504060bc81559dc9c5864ab5 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 21 Dec 2020 16:10:31 +0100 Subject: [PATCH 04/97] normalize: add to docs --- docs/source/normalize.rst | 8 ++++++++ docs/source/package.rst | 1 + gstools/__init__.py | 9 ++++++++- gstools/normalize/__init__.py | 6 ++++-- gstools/normalize/normalizer.py | 20 ++++++++++---------- 5 files changed, 31 insertions(+), 13 deletions(-) create mode 100644 docs/source/normalize.rst diff --git a/docs/source/normalize.rst b/docs/source/normalize.rst new file mode 100644 index 000000000..160356404 --- /dev/null +++ b/docs/source/normalize.rst @@ -0,0 +1,8 @@ +gstools.normalize +================= + +.. automodule:: gstools.normalize + +.. raw:: latex + + \clearpage diff --git a/docs/source/package.rst b/docs/source/package.rst index 73d3db961..487351a91 100644 --- a/docs/source/package.rst +++ b/docs/source/package.rst @@ -18,3 +18,4 @@ GSTools API random.rst tools.rst transform.rst + normalize.rst diff --git a/gstools/__init__.py b/gstools/__init__.py index 1599bafad..1073c531b 100644 --- a/gstools/__init__.py +++ b/gstools/__init__.py @@ -108,7 +108,14 @@ """ # Hooray! from gstools import ( - field, variogram, random, covmodel, tools, krige, transform, normalize + field, + variogram, + random, + covmodel, + tools, + krige, + transform, + normalize, ) from gstools.field import SRF from gstools.tools import ( diff --git a/gstools/normalize/__init__.py b/gstools/normalize/__init__.py index 02fb586ea..2dab088b2 100644 --- a/gstools/normalize/__init__.py +++ b/gstools/normalize/__init__.py @@ -8,20 +8,22 @@ ^^^^^^^^^^^^^^^ .. autosummary:: + :toctree: generated + Normalizer Field-Normalizer ^^^^^^^^^^^^^^^^ .. autosummary:: + :toctree: generated + LogNormal BoxCox BoxCoxShift YeoJohnson Modulus Manly - ----- """ from gstools.normalize.base import Normalizer diff --git a/gstools/normalize/normalizer.py b/gstools/normalize/normalizer.py index bbaf8e4ac..e0ec8f8b3 100644 --- a/gstools/normalize/normalizer.py +++ b/gstools/normalize/normalizer.py @@ -43,7 +43,7 @@ def derivative(self, values): class BoxCox(Normalizer): - r"""Box-Cox (1964) transformed fields [1]_. + r"""Box-Cox (1964) transformed fields. Parameters ---------- @@ -55,7 +55,7 @@ class BoxCox(Normalizer): Notes ----- - This transformation is given by: + This transformation is given by [1]_: .. math:: y=\begin{cases} @@ -91,7 +91,7 @@ def derivative(self, values): class BoxCoxShift(Normalizer): - r"""Box-Cox (1964) transformed fields including shifting [1]_. + r"""Box-Cox (1964) transformed fields including shifting. Parameters ---------- @@ -105,7 +105,7 @@ class BoxCoxShift(Normalizer): Notes ----- - This transformation is given by: + This transformation is given by [1]_: .. math:: y=\begin{cases} @@ -151,7 +151,7 @@ def derivative(self, values): class YeoJohnson(Normalizer): - r"""Yeo-Johnson (2000) transformed fields [1]_. + r"""Yeo-Johnson (2000) transformed fields. Parameters ---------- @@ -163,7 +163,7 @@ class YeoJohnson(Normalizer): Notes ----- - This transformation is given by: + This transformation is given by [1]_: .. math:: y=\begin{cases} @@ -235,7 +235,7 @@ def derivative(self, values): class Modulus(Normalizer): - r"""Modulus or John-Draper (1980) transformed fields [1]_. + r"""Modulus or John-Draper (1980) transformed fields. Parameters ---------- @@ -247,7 +247,7 @@ class Modulus(Normalizer): Notes ----- - This transformation is given by: + This transformation is given by [1]_: .. math:: y=\begin{cases} @@ -290,7 +290,7 @@ def derivative(self, values): class Manly(Normalizer): - r"""Manly (1971) transformed fields [1]_. + r"""Manly (1971) transformed fields. Parameters ---------- @@ -302,7 +302,7 @@ class Manly(Normalizer): Notes ----- - This transformation is given by: + This transformation is given by [1]_: .. math:: y=\begin{cases} From 1ebfcf69d7a58de764417414085e7b4b7a614240 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 6 Jan 2021 14:29:43 +0100 Subject: [PATCH 05/97] Covmodel: use 'sorted' to gain soterd opt-args --- gstools/covmodel/tools.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/gstools/covmodel/tools.py b/gstools/covmodel/tools.py index 96a5b0ae5..cf5233334 100644 --- a/gstools/covmodel/tools.py +++ b/gstools/covmodel/tools.py @@ -203,8 +203,7 @@ def set_opt_args(model, opt_arg): if def_arg not in opt_arg: opt_arg[def_arg] = default[def_arg] # save names of the optional arguments (sort them by name) - model._opt_arg = list(opt_arg.keys()) - model._opt_arg.sort() + model._opt_arg = sorted(opt_arg) # add the optional arguments as attributes to the class for opt_name in opt_arg: if opt_name in dir(model): # "dir" also respects properties From c0735e22698edcebfe71a450e86952da1b778f87 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 6 Jan 2021 14:30:15 +0100 Subject: [PATCH 06/97] Normalize: add more likelihood routines --- gstools/normalize/base.py | 72 ++++++++++++++++++++++++++------------- 1 file changed, 48 insertions(+), 24 deletions(-) diff --git a/gstools/normalize/base.py b/gstools/normalize/base.py index 5cd547076..7663f9e50 100644 --- a/gstools/normalize/base.py +++ b/gstools/normalize/base.py @@ -21,7 +21,7 @@ class Normalizer: Parameters ---------- data : array_like, optional - Input data to fit the transformation in order to gain normality. + Input data to fit the transformation to in order to gain normality. The default is None. **parameter Specified parameters given by name. If not given, default values @@ -92,27 +92,57 @@ def derivative(self, values): """ return spm.derivative(self.normalize, np.asanyarray(values), dx=1e-6) + def likelihood(self, data): + """Likelihood for given data with current parameters. + + Parameters + ---------- + data : array_like + Input data to fit the transformation to in order to gain normality. + + Returns + ------- + :class:`float` + Likelihood of the given data. + """ + return np.exp(self.loglikelihood(data)) + def loglikelihood(self, data): """Log-Likelihood for given data with current parameters. Parameters ---------- data : array_like - Input data to fit the transformation in order to gain normality. + Input data to fit the transformation to in order to gain normality. Returns ------- :class:`float` Log-Likelihood of the given data. + """ + add = -0.5 * np.size(data) * np.log(2 * np.pi) - 0.5 + return self.kernel_loglikelihood(data) + add + + def kernel_loglikelihood(self, data): + """Kernel Log-Likelihood for given data with current parameters. + + Parameters + ---------- + data : array_like + Input data to fit the transformation to in order to gain normality. + + Returns + ------- + :class:`float` + Kernel Log-Likelihood of the given data. Notes ----- - This loglikelihood function could be missing additive constants, + This loglikelihood function is neglecting additive constants, that are not needed for optimization. """ res = -0.5 * np.size(data) * np.log(np.var(self.normalize(data))) - res += np.sum(np.log(np.maximum(1e-16, self.derivative(data)))) - return res + return res + np.sum(np.log(np.maximum(1e-16, self.derivative(data)))) def fit(self, data, skip=None, **kwargs): """Fitting the transformation to data by maximizing Log-Likelihood. @@ -120,7 +150,7 @@ def fit(self, data, skip=None, **kwargs): Parameters ---------- data : array_like - Input data to fit the transformation in order to gain normality. + Input data to fit the transformation to in order to gain normality. skip : :class:`list` of :class:`str` or :any:`None`, optional Names of parameters to be skiped in fitting. The default is None. @@ -134,41 +164,35 @@ def fit(self, data, skip=None, **kwargs): Optimal paramters given by names. """ skip = [] if skip is None else skip - all_names = list(self.default_parameter().keys()) - para_names = [p for p in all_names if p not in skip] + all_names = sorted(self.default_parameter()) + para_names = [name for name in all_names if name not in skip] - def _neg_llf(par, dat): + def _neg_kllf(par, dat): for name, val in zip(para_names, np.atleast_1d(par)): setattr(self, name, val) - return -self.loglikelihood(dat) + return -self.kernel_loglikelihood(dat) if len(para_names) == 0: # transformations without para. (no opti.) - warnings.warn("Normalizer.fit: no parameters for fitting.") + warnings.warn(self.__class__.__name__ + ".fit: no parameters!") return {} if len(para_names) == 1: # one-para. transformations (simple opti.) # default bracket like in scipy's boxcox (if not given) kwargs.setdefault("bracket", (-2, 2)) - out = spo.minimize_scalar(_neg_llf, args=(data,), **kwargs) + out = spo.minimize_scalar(_neg_kllf, args=(data,), **kwargs) else: # general case # init guess from current values (if x0 not given) kwargs.setdefault("x0", [getattr(self, p) for p in para_names]) - out = spo.minimize(_neg_llf, args=(data,), **kwargs) + out = spo.minimize(_neg_kllf, args=(data,), **kwargs) # save optimization results self._opti = out for name, val in zip(para_names, np.atleast_1d(out.x)): setattr(self, name, val) return {name: getattr(self, name) for name in all_names} - def __str__(self): - """Return String representation.""" - return self.__repr__() - def __repr__(self): """Return String representation.""" - out_str = self.__class__.__name__ - para_strs = [] - for p in self.default_parameter(): - para_strs.append( - "{0}={1:.{2}}".format(p, float(getattr(self, p)), self._prec) - ) - return out_str + "(" + ", ".join(para_strs) + ")" + para_strs = [ + "{0}={1:.{2}}".format(p, float(getattr(self, p)), self._prec) + for p in sorted(self.default_parameter()) + ] + return self.__class__.__name__ + "(" + ", ".join(para_strs) + ")" From f42195b366c05e2b636c488747a9ebb983bebdf4 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 11 Jan 2021 18:28:17 +0100 Subject: [PATCH 07/97] Normalizer: add init value for self._opti --- gstools/normalize/base.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/gstools/normalize/base.py b/gstools/normalize/base.py index 7663f9e50..f7f299c21 100644 --- a/gstools/normalize/base.py +++ b/gstools/normalize/base.py @@ -29,11 +29,14 @@ class Normalizer: """ def __init__(self, data=None, **parameter): - # only use values, that have a provided default value + # only use parameter, that have a provided default value for key, value in self.default_parameter().items(): setattr(self, key, parameter.get(key, value)) + # fit parameters if data is given if data is not None: self.fit(data) + # optimization results + self._opti = None # precision for printing self._prec = 3 From 9f345bbb3c0d0653dec3c53c1ff4a61150af0b6d Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 11 Jan 2021 18:28:46 +0100 Subject: [PATCH 08/97] Tests: add stats tests for normalizer --- tests/test_normalize.py | 61 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 61 insertions(+) create mode 100644 tests/test_normalize.py diff --git a/tests/test_normalize.py b/tests/test_normalize.py new file mode 100644 index 000000000..ae033d37b --- /dev/null +++ b/tests/test_normalize.py @@ -0,0 +1,61 @@ +# -*- coding: utf-8 -*- +""" +This is the unittest of the Normalizer class. +""" + +import unittest +import numpy as np +import gstools as gs + + +def _rel_err(a, b): + return np.abs(a / ((a + b) / 2) - 1) + + +class TestNormalizer(unittest.TestCase): + def setUp(self): + self.seed = 20210111 + self.rng = gs.random.RNG(self.seed) + self.samp = self.rng.random.normal(11.1, 2.25, 1000) + self.lmb = 2.5 + + def test_fitting(self): + # boxcox with given data to init + bc_samples = gs.normalize.BoxCox(lmbda=self.lmb).denormalize(self.samp) + bc_norm = gs.normalize.BoxCox(data=bc_samples) + self.assertLess(_rel_err(self.lmb, bc_norm.lmbda), 1e-2) + self.assertAlmostEqual( + bc_norm.likelihood(bc_samples), + np.exp(bc_norm.loglikelihood(bc_samples)), + ) + # yeo-johnson with calling fit + yj_norm = gs.normalize.YeoJohnson(lmbda=self.lmb) + yj_samples = yj_norm.denormalize(self.samp) + yj_norm.fit(yj_samples) + self.assertLess(_rel_err(self.lmb, yj_norm.lmbda), 1e-2) + self.assertAlmostEqual( + yj_norm.likelihood(yj_samples), + np.exp(yj_norm.loglikelihood(yj_samples)), + ) + # modulus with calling fit + mo_norm = gs.normalize.Modulus(lmbda=self.lmb) + mo_samples = mo_norm.denormalize(self.samp) + mo_norm.fit(mo_samples) + self.assertLess(_rel_err(self.lmb, mo_norm.lmbda), 1e-2) + self.assertAlmostEqual( + mo_norm.likelihood(mo_samples), + np.exp(mo_norm.loglikelihood(mo_samples)), + ) + # manly with calling fit + ma_norm = gs.normalize.Manly(lmbda=self.lmb) + ma_samples = ma_norm.denormalize(self.samp) + ma_norm.fit(ma_samples) + self.assertLess(_rel_err(self.lmb, ma_norm.lmbda), 1e-2) + # self.assertAlmostEqual( + # ma_norm.likelihood(ma_samples), + # np.exp(ma_norm.loglikelihood(ma_samples)), + # ) # this is comparing infs + + +if __name__ == "__main__": + unittest.main() From d4528cfd17cc2bc1d31455f2a324265943d5c385 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 11 Jan 2021 23:07:21 +0100 Subject: [PATCH 09/97] Normalizer: fix bug in likelihood --- gstools/normalize/base.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gstools/normalize/base.py b/gstools/normalize/base.py index f7f299c21..b2b718958 100644 --- a/gstools/normalize/base.py +++ b/gstools/normalize/base.py @@ -123,7 +123,7 @@ def loglikelihood(self, data): :class:`float` Log-Likelihood of the given data. """ - add = -0.5 * np.size(data) * np.log(2 * np.pi) - 0.5 + add = -0.5 * np.size(data) * (np.log(2 * np.pi) + 1) return self.kernel_loglikelihood(data) + add def kernel_loglikelihood(self, data): From b47fb2d33ec069cdd5aff8d115ecd5d143bff0c0 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 11 Jan 2021 23:07:42 +0100 Subject: [PATCH 10/97] Normalizer: fix doc typo --- gstools/normalize/normalizer.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gstools/normalize/normalizer.py b/gstools/normalize/normalizer.py index e0ec8f8b3..8e5613e1c 100644 --- a/gstools/normalize/normalizer.py +++ b/gstools/normalize/normalizer.py @@ -171,7 +171,7 @@ class YeoJohnson(Normalizer): & x\geq 0,\, \lambda\neq 0 \\ \log(x+1) & x\geq 0,\, \lambda = 0 \\ - \frac{(|x|+1)^{2-\lambda} - 1}{2-\lambda} + -\frac{(|x|+1)^{2-\lambda} - 1}{2-\lambda} & x<0,\, \lambda\neq 2 \\ -\log(|x|+1) & x<0,\, \lambda = 2 From b2542af8eb7fea633c7e732ab1c147f167c7e723 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 11 Jan 2021 23:08:00 +0100 Subject: [PATCH 11/97] Tests: add more normalizer tests --- tests/test_normalize.py | 106 +++++++++++++++++++++++++++++++++++++++- 1 file changed, 104 insertions(+), 2 deletions(-) diff --git a/tests/test_normalize.py b/tests/test_normalize.py index ae033d37b..f8abd3432 100644 --- a/tests/test_normalize.py +++ b/tests/test_normalize.py @@ -16,8 +16,10 @@ class TestNormalizer(unittest.TestCase): def setUp(self): self.seed = 20210111 self.rng = gs.random.RNG(self.seed) - self.samp = self.rng.random.normal(11.1, 2.25, 1000) - self.lmb = 2.5 + self.mean = 11.1 + self.std = 2.25 + self.samp = self.rng.random.normal(self.mean, self.std, 1000) + self.lmb = 1.5 def test_fitting(self): # boxcox with given data to init @@ -56,6 +58,106 @@ def test_fitting(self): # np.exp(ma_norm.loglikelihood(ma_samples)), # ) # this is comparing infs + def test_boxcox(self): + # without shift + bc = gs.normalize.BoxCox(lmbda=0) + self.assertTrue( + np.all( + np.isclose(self.samp, bc.normalize(bc.denormalize(self.samp))) + ) + ) + bc.lmbda = self.lmb + self.assertTrue( + np.all( + np.isclose(self.samp, bc.normalize(bc.denormalize(self.samp))) + ) + ) + # with shift + bc = gs.normalize.BoxCoxShift(lmbda=0, shift=1.1) + self.assertTrue( + np.all( + np.isclose(self.samp, bc.normalize(bc.denormalize(self.samp))) + ) + ) + bc.lmbda = self.lmb + self.assertTrue( + np.all( + np.isclose(self.samp, bc.normalize(bc.denormalize(self.samp))) + ) + ) + + def test_yeojohnson(self): + yj = gs.normalize.YeoJohnson(lmbda=0) + self.assertTrue( + np.all( + np.isclose( + self.samp - self.mean, + yj.normalize(yj.denormalize(self.samp - self.mean)), + ) + ) + ) + yj.lmbda = 2 + self.assertTrue( + np.all( + np.isclose( + self.samp - self.mean, + yj.normalize(yj.denormalize(self.samp - self.mean)), + ) + ) + ) + # with shift + yj.lmbda = self.lmb + self.assertTrue( + np.all( + np.isclose( + self.samp - self.mean, + yj.normalize(yj.denormalize(self.samp - self.mean)), + ) + ) + ) + + def test_modulus(self): + mo = gs.normalize.Modulus(lmbda=0) + self.assertTrue( + np.all( + np.isclose(self.samp, mo.normalize(mo.denormalize(self.samp))) + ) + ) + mo.lmbda = self.lmb + self.assertTrue( + np.all( + np.isclose(self.samp, mo.normalize(mo.denormalize(self.samp))) + ) + ) + + def test_manly(self): + ma = gs.normalize.Manly(lmbda=0) + self.assertTrue( + np.all( + np.isclose(self.samp, ma.normalize(ma.denormalize(self.samp))) + ) + ) + ma.lmbda = self.lmb + self.assertTrue( + np.all( + np.isclose(self.samp, ma.normalize(ma.denormalize(self.samp))) + ) + ) + + def test_parameterless(self): + no = gs.normalize.LogNormal() + self.assertTrue( + np.all( + np.isclose(self.samp, no.normalize(no.denormalize(self.samp))) + ) + ) + no = gs.normalize.Normalizer() + self.assertTrue( + np.all( + np.isclose(self.samp, no.normalize(no.denormalize(self.samp))) + ) + ) + if __name__ == "__main__": unittest.main() From 3ccf7032822458c464bac01f7c42be6df6aab4c8 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 27 Jan 2021 14:02:59 +0100 Subject: [PATCH 12/97] Normalize: typo --- gstools/normalize/normalizer.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gstools/normalize/normalizer.py b/gstools/normalize/normalizer.py index 8e5613e1c..0fd72f2fc 100644 --- a/gstools/normalize/normalizer.py +++ b/gstools/normalize/normalizer.py @@ -55,7 +55,7 @@ class BoxCox(Normalizer): Notes ----- - This transformation is given by [1]_: + This transformation is given by [1]_: .. math:: y=\begin{cases} From 086780aa719bc66fbb300d81ba4c6c46e82f7b6b Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 27 Jan 2021 17:05:56 +0100 Subject: [PATCH 13/97] Examples: blackenend --- examples/09_spatio_temporal/02_precip_2d.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/examples/09_spatio_temporal/02_precip_2d.py b/examples/09_spatio_temporal/02_precip_2d.py index 7171260e6..a026d5504 100644 --- a/examples/09_spatio_temporal/02_precip_2d.py +++ b/examples/09_spatio_temporal/02_precip_2d.py @@ -52,13 +52,15 @@ ############################################################################### # plot the 2d precipitation field over time as an animation. + def _update_ani(time_step): im.set_array(srf.field[:, :, time_step].T) - return im, + return (im,) + fig, ax = plt.subplots() im = ax.imshow( - srf.field[:,:,0].T, + srf.field[:, :, 0].T, cmap="Blues", interpolation="bicubic", origin="lower", @@ -68,4 +70,6 @@ def _update_ani(time_step): ax.set_xlabel(r"$x$ / km") ax.set_ylabel(r"$y$ / km") -ani = animation.FuncAnimation(fig, _update_ani, len(t), interval=100, blit=True) +ani = animation.FuncAnimation( + fig, _update_ani, len(t), interval=100, blit=True +) From e632ef5042726b6eb089e188e50dad4972f5bed8 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Thu, 28 Jan 2021 12:19:19 +0100 Subject: [PATCH 14/97] Field: outsourcing helper functions --- gstools/field/tools.py | 235 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 235 insertions(+) create mode 100644 gstools/field/tools.py diff --git a/gstools/field/tools.py b/gstools/field/tools.py new file mode 100644 index 000000000..d0ba4ff3e --- /dev/null +++ b/gstools/field/tools.py @@ -0,0 +1,235 @@ +# -*- coding: utf-8 -*- +""" +GStools subpackage providing tools for Fields. + +.. currentmodule:: gstools.field.tools + +The following classes and functions are provided + +.. autosummary:: + to_vtk_helper + mesh_call +""" +import numpy as np +import meshio + +from gstools.tools.export import to_vtk, vtk_export + + +__all__ = ["to_vtk_helper", "mesh_call"] + + +def to_vtk_helper( + f_cls, filename=None, field_select="field", fieldname="field" +): # pragma: no cover + """Create a VTK/PyVista grid of the field or save it as a VTK file. + + This is an internal helper that will handle saving or creating objects + + Parameters + ---------- + f_cls : :any:`Field` + Field class in use. + filename : :class:`str` + Filename of the file to be saved, including the path. Note that an + ending (.vtr or .vtu) will be added to the name. If ``None`` is + passed, a PyVista dataset of the appropriate type will be returned. + field_select : :class:`str`, optional + Field that should be stored. Can be: + "field", "raw_field", "krige_field", "err_field" or "krige_var". + Default: "field" + fieldname : :class:`str`, optional + Name of the field in the VTK file. Default: "field" + """ + if f_cls.value_type == "vector": + if hasattr(f_cls, field_select): + field = getattr(f_cls, field_select) + else: + field = None + if not (f_cls.pos is None or field is None or f_cls.mesh_type is None): + suf = ["_X", "_Y", "_Z"] + fields = {} + for i in range(f_cls.model.dim): + fields[fieldname + suf[i]] = field[i] + if filename is None: + return to_vtk(f_cls.pos, fields, f_cls.mesh_type) + else: + return vtk_export(filename, f_cls.pos, fields, f_cls.mesh_type) + elif f_cls.value_type == "scalar": + if hasattr(f_cls, field_select): + field = getattr(f_cls, field_select) + else: + field = None + if not (f_cls.pos is None or field is None or f_cls.mesh_type is None): + if filename is None: + return to_vtk(f_cls.pos, {fieldname: field}, f_cls.mesh_type) + else: + return vtk_export( + filename, f_cls.pos, {fieldname: field}, f_cls.mesh_type + ) + else: + print("Field.to_vtk: '{}' not available.".format(field_select)) + else: + raise ValueError( + "Unknown field value type: {}".format(f_cls.value_type) + ) + + +def mesh_call( + f_cls, mesh, points="centroids", direction="all", name="field", **kwargs +): + """Generate a field on a given meshio or ogs5py mesh. + + Parameters + ---------- + mesh : meshio.Mesh or ogs5py.MSH or PyVista mesh + The given meshio, ogs5py, or PyVista mesh + points : :class:`str`, optional + The points to evaluate the field at. + Either the "centroids" of the mesh cells + (calculated as mean of the cell vertices) or the "points" + of the given mesh. + Default: "centroids" + direction : :class:`str` or :class:`list`, optional + Here you can state which direction should be choosen for + lower dimension. For example, if you got a 2D mesh in xz direction, + you have to pass "xz". By default, all directions are used. + One can also pass a list of indices. + Default: "all" + name : :class:`str` or :class:`list` of :class:`str`, optional + Name(s) to store the field(s) in the given mesh as point_data or + cell_data. If to few names are given, digits will be appended. + Default: "field" + **kwargs + Keyword arguments forwareded to `Field.__call__`. + + Notes + ----- + This will store the field in the given mesh under the given name, + if a meshio or PyVista mesh was given. + + See: https://github.com/nschloe/meshio + See: https://github.com/pyvista/pyvista + + See: :any:`Field.__call__` + """ + has_pyvista = False + has_ogs5py = False + + try: + import pyvista as pv + + has_pyvista = True + except ImportError: + pass + try: + import ogs5py as ogs + + has_ogs5py = True + except ImportError: + pass + + if isinstance(direction, str) and direction == "all": + select = list(range(f_cls.model.field_dim)) + elif isinstance(direction, str): + select = _get_select(direction)[: f_cls.model.field_dim] + else: + select = direction[: f_cls.model.field_dim] + if len(select) < f_cls.model.field_dim: + raise ValueError( + "Field.mesh: need at least {} direction(s), got '{}'".format( + f_cls.model.field_dim, direction + ) + ) + # convert pyvista mesh + if has_pyvista and pv.is_pyvista_dataset(mesh): + if points == "centroids": + pnts = mesh.cell_centers().points.T[select] + else: + pnts = mesh.points.T[select] + out = f_cls.unstructured(pos=pnts, **kwargs) + # Deal with the output + fields = [out] if isinstance(out, np.ndarray) else out + for f_name, field in zip(_names(name, len(fields)), fields): + mesh[f_name] = field + # convert ogs5py mesh + elif has_ogs5py and isinstance(mesh, ogs.MSH): + if points == "centroids": + pnts = mesh.centroids_flat.T[select] + else: + pnts = mesh.NODES.T[select] + out = f_cls.unstructured(pos=pnts, **kwargs) + # convert meshio mesh + elif isinstance(mesh, meshio.Mesh): + if points == "centroids": + # define unique order of cells + offset = [] + length = [] + mesh_dim = mesh.points.shape[1] + if mesh_dim < f_cls.model.field_dim: + raise ValueError("Field.mesh: mesh dimension too low!") + pnts = np.empty((0, mesh_dim), dtype=np.double) + for cell in mesh.cells: + pnt = np.mean(mesh.points[cell[1]], axis=1) + offset.append(pnts.shape[0]) + length.append(pnt.shape[0]) + pnts = np.vstack((pnts, pnt)) + # generate pos for __call__ + pnts = pnts.T[select] + out = f_cls.unstructured(pos=pnts, **kwargs) + fields = [out] if isinstance(out, np.ndarray) else out + f_lists = [] + for field in fields: + f_list = [] + for of, le in zip(offset, length): + f_list.append(field[of : of + le]) + f_lists.append(f_list) + for f_name, f_list in zip(_names(name, len(f_lists)), f_lists): + mesh.cell_data[f_name] = f_list + else: + out = f_cls.unstructured(pos=mesh.points.T[select], **kwargs) + fields = [out] if isinstance(out, np.ndarray) else out + for f_name, field in zip(_names(name, len(fields)), fields): + mesh.point_data[f_name] = field + else: + raise ValueError("Field.mesh: Unknown mesh format!") + return out + + +def _names(name, cnt): + name = [name] if isinstance(name, str) else list(name)[:cnt] + if len(name) < cnt: + name += [name[-1] + str(i + 1) for i in range(cnt - len(name))] + return name + + +def _get_select(direction): + select = [] + if not (0 < len(direction) < 4): + raise ValueError( + "Field.mesh: need 1 to 3 direction(s), got '{}'".format(direction) + ) + for axis in direction: + if axis == "x": + if 0 in select: + raise ValueError( + "Field.mesh: got duplicate directions {}".format(direction) + ) + select.append(0) + elif axis == "y": + if 1 in select: + raise ValueError( + "Field.mesh: got duplicate directions {}".format(direction) + ) + select.append(1) + elif axis == "z": + if 2 in select: + raise ValueError( + "Field.mesh: got duplicate directions {}".format(direction) + ) + select.append(2) + else: + raise ValueError( + "Field.mesh: got unknown direction {}".format(axis) + ) + return select From a85ee97318d1dc60435aead29088a0526aa3fee9 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Thu, 28 Jan 2021 14:15:00 +0100 Subject: [PATCH 15/97] Normalize: add name attribute to class --- gstools/normalize/base.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/gstools/normalize/base.py b/gstools/normalize/base.py index b2b718958..e0645c7d2 100644 --- a/gstools/normalize/base.py +++ b/gstools/normalize/base.py @@ -192,10 +192,15 @@ def _neg_kllf(par, dat): setattr(self, name, val) return {name: getattr(self, name) for name in all_names} + @property + def name(self): + """:class:`str`: The name of the normalizer class.""" + return self.__class__.__name__ + def __repr__(self): """Return String representation.""" para_strs = [ "{0}={1:.{2}}".format(p, float(getattr(self, p)), self._prec) for p in sorted(self.default_parameter()) ] - return self.__class__.__name__ + "(" + ", ".join(para_strs) + ")" + return self.name + "(" + ", ".join(para_strs) + ")" From 98974b10f4abd4ce7a2befe3a61668a2f2409dc0 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Thu, 28 Jan 2021 14:15:52 +0100 Subject: [PATCH 16/97] Tools: add eval_func to misc tools --- gstools/tools/misc.py | 70 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 70 insertions(+) diff --git a/gstools/tools/misc.py b/gstools/tools/misc.py index 93ef4cc0b..d9f861e01 100755 --- a/gstools/tools/misc.py +++ b/gstools/tools/misc.py @@ -9,6 +9,8 @@ .. autosummary:: list_format """ +import numpy as np +from gstools.tools.geometric import format_struct_pos_dim, gen_mesh def list_format(lst, prec): @@ -16,3 +18,71 @@ def list_format(lst, prec): return "[{}]".format( ", ".join("{x:.{p}}".format(x=x, p=prec) for x in lst) ) + + +def eval_func( + func_val, + pos, + dim, + mesh_type="unstructured", + value_type="scalar", + broadcast=False, +): + """ + Evaluate a function on a mesh. + + Parameters + ---------- + func_val : :any:`callable` or :class:`float` or :any:`None` + Function to be called or single value to be filled. + Should have the signiture f(x, [y, z, ...]) in case of callable. + In case of a float, the field will be filled with a single value and + in case of None, this value will be set to 0. + pos : :class:`list` + The position tuple, containing main direction and transversal + directions (x, [y, z, ...]). + dim : :class:`int` + The spatial dimension. + mesh_type : :class:`str`, optional + 'structured' / 'unstructured' + Default: 'unstructured' + value_type : :class:`str`, optional + Value type of the field. Either "scalar" or "vector". + The default is "scalar". + broadcast : :class:`bool`, optional + Whether to return a single value, if a single value was given. + Default: False + + Returns + ------- + :class:`numpy.ndarray` + Function values at the given points. + """ + # care about scalar inputs + if broadcast and not callable(func_val): + return 0.0 if func_val is None else float(func_val) + elif not callable(func_val): + func_val = _func_from_scalar(func_val) + # care about mesh and function call + if mesh_type != "unstructured": + pos, shape = format_struct_pos_dim(pos, dim) + pos = gen_mesh(pos) + else: + pos = np.array(pos, dtype=np.double).reshape(dim, -1) + shape = np.shape(pos[0]) + # prepend dimension if we have a vector field + if value_type == "vector": + shape = (dim,) + shape + return np.reshape(func_val(*pos), shape) + + +def _func_from_scalar(value, value_type="scalar"): + value = 0.0 if value is None else float(value) + + def _f(*pos): + if value_type == "vector": + return np.full_like(pos, value, dtype=np.double) + # scalar field has same shape like a single axis + return np.full_like(pos[0], value, dtype=np.double) + + return _f From 00257da589e9e11eb6b7ee80c97a3e74446d3dbd Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Thu, 28 Jan 2021 14:17:36 +0100 Subject: [PATCH 17/97] Field: add mean, normalizer, trend to Field class and refactoring --- gstools/field/base.py | 323 ++++++++++++++++++------------------------ 1 file changed, 134 insertions(+), 189 deletions(-) diff --git a/gstools/field/base.py b/gstools/field/base.py index 0c1cc4248..86019442f 100755 --- a/gstools/field/base.py +++ b/gstools/field/base.py @@ -10,20 +10,18 @@ Field """ # pylint: disable=C0103 - from functools import partial - import numpy as np -import meshio - from gstools.covmodel.base import CovModel from gstools.tools.geometric import format_struct_pos_dim, gen_mesh -from gstools.tools.export import to_vtk, vtk_export - +from gstools.tools.misc import eval_func +from gstools.normalize import Normalizer +from gstools.field.tools import mesh_call, to_vtk_helper __all__ = ["Field"] VALUE_TYPES = ["scalar", "vector"] +""":class:`list` of :class:`str`: valid field value types.""" class Field: @@ -33,18 +31,44 @@ class Field: ---------- model : :any:`CovModel` Covariance Model related to the field. + value_type : :class:`str`, optional + Value type of the field. Either "scalar" or "vector". + The default is "scalar". + mean : :any:`None` or :class:`float` or :any:`callable`, optional + Mean of the field if wanted. Could also be a callable. + The default is None. + normalizer : :any:`None` or :any:`Normalizer`, optional + Normalizer to be applied to the field. + The default is None. + trend : :any:`None` or :class:`float` or :any:`callable`, optional + Trend of the denormalized fields. If no normalizer is applied, + this behaves equal to 'mean'. + The default is None. """ - def __init__(self, model, value_type="scalar"): + def __init__( + self, + model, + value_type="scalar", + mean=None, + normalizer=None, + trend=None, + ): # initialize attributes self.pos = None self.mesh_type = None - self.field = None # initialize private attributes self._model = None - self.model = model self._value_type = None + self._mean = None + self._normalizer = None + self._trend = None + # set properties + self.model = model self.value_type = value_type + self.mean = mean + self.normalizer = normalizer + self.trend = trend def __call__(*args, **kwargs): """Generate the field.""" @@ -104,87 +128,7 @@ def mesh( See: :any:`Field.__call__` """ - has_pyvista = False - has_ogs5py = False - - try: - import pyvista as pv - - has_pyvista = True - except ImportError: - pass - try: - import ogs5py as ogs - - has_ogs5py = True - except ImportError: - pass - - if isinstance(direction, str) and direction == "all": - select = list(range(self.model.field_dim)) - elif isinstance(direction, str): - select = _get_select(direction)[: self.model.field_dim] - else: - select = direction[: self.model.field_dim] - if len(select) < self.model.field_dim: - raise ValueError( - "Field.mesh: need at least {} direction(s), got '{}'".format( - self.model.field_dim, direction - ) - ) - # convert pyvista mesh - if has_pyvista and pv.is_pyvista_dataset(mesh): - if points == "centroids": - pnts = mesh.cell_centers().points.T[select] - else: - pnts = mesh.points.T[select] - out = self.unstructured(pos=pnts, **kwargs) - # Deal with the output - fields = [out] if isinstance(out, np.ndarray) else out - for f_name, field in zip(_names(name, len(fields)), fields): - mesh[f_name] = field - # convert ogs5py mesh - elif has_ogs5py and isinstance(mesh, ogs.MSH): - if points == "centroids": - pnts = mesh.centroids_flat.T[select] - else: - pnts = mesh.NODES.T[select] - out = self.unstructured(pos=pnts, **kwargs) - # convert meshio mesh - elif isinstance(mesh, meshio.Mesh): - if points == "centroids": - # define unique order of cells - offset = [] - length = [] - mesh_dim = mesh.points.shape[1] - if mesh_dim < self.model.field_dim: - raise ValueError("Field.mesh: mesh dimension too low!") - pnts = np.empty((0, mesh_dim), dtype=np.double) - for cell in mesh.cells: - pnt = np.mean(mesh.points[cell[1]], axis=1) - offset.append(pnts.shape[0]) - length.append(pnt.shape[0]) - pnts = np.vstack((pnts, pnt)) - # generate pos for __call__ - pnts = pnts.T[select] - out = self.unstructured(pos=pnts, **kwargs) - fields = [out] if isinstance(out, np.ndarray) else out - f_lists = [] - for field in fields: - f_list = [] - for of, le in zip(offset, length): - f_list.append(field[of : of + le]) - f_lists.append(f_list) - for f_name, f_list in zip(_names(name, len(f_lists)), f_lists): - mesh.cell_data[f_name] = f_list - else: - out = self.unstructured(pos=mesh.points.T[select], **kwargs) - fields = [out] if isinstance(out, np.ndarray) else out - for f_name, field in zip(_names(name, len(fields)), fields): - mesh.point_data[f_name] = field - else: - raise ValueError("Field.mesh: Unknown mesh format!") - return out + return mesh_call(self, mesh, points, direction, name, **kwargs) def pre_pos(self, pos, mesh_type="unstructured"): """ @@ -226,64 +170,46 @@ def pre_pos(self, pos, mesh_type="unstructured"): # return isometrized pos tuple and resulting field shape return self.model.isometrize(pos), shape - def _to_vtk_helper( - self, filename=None, field_select="field", fieldname="field" - ): # pragma: no cover - """Create a VTK/PyVista grid of the field or save it as a VTK file. - - This is an internal helper that will handle saving or creating objects + def post_field(self, field, name="field", process=True, save=True): + """ + Postprocessing field values. Parameters ---------- - filename : :class:`str` - Filename of the file to be saved, including the path. Note that an - ending (.vtr or .vtu) will be added to the name. If ``None`` is - passed, a PyVista dataset of the appropriate type will be returned. - field_select : :class:`str`, optional - Field that should be stored. Can be: - "field", "raw_field", "krige_field", "err_field" or "krige_var". - Default: "field" - fieldname : :class:`str`, optional - Name of the field in the VTK file. Default: "field" + field : :class:`numpy.ndarray` + Field values. + name : :class:`str`, optional + Name. to store the field. + The default is "field". + process : :class:`bool`, optional + Whether to process field to apply mean, normalizer and trend. + The default is True. + save : :class:`bool`, optional + Whether to store the field under the given name. + The default is True. + + Returns + ------- + field : :class:`numpy.ndarray` + Processed field values. """ - if self.value_type == "vector": - if hasattr(self, field_select): - field = getattr(self, field_select) - else: - field = None - if not ( - self.pos is None or field is None or self.mesh_type is None - ): - suf = ["_X", "_Y", "_Z"] - fields = {} - for i in range(self.model.dim): - fields[fieldname + suf[i]] = field[i] - if filename is None: - return to_vtk(self.pos, fields, self.mesh_type) - else: - return vtk_export( - filename, self.pos, fields, self.mesh_type - ) - elif self.value_type == "scalar": - if hasattr(self, field_select): - field = getattr(self, field_select) - else: - field = None - if not ( - self.pos is None or field is None or self.mesh_type is None - ): - if filename is None: - return to_vtk(self.pos, {fieldname: field}, self.mesh_type) - else: - return vtk_export( - filename, self.pos, {fieldname: field}, self.mesh_type - ) - else: - print("Field.to_vtk: '{}' not available.".format(field_select)) - else: - raise ValueError( - "Unknown field value type: {}".format(self.value_type) + if process: + if self.pos is None: + raise ValueError("post_field: no 'pos' tuple set for field.") + kwargs = dict( + pos=self.pos, + dim=self.model.dim, + mesh_type=self.mesh_type, + value_type=self.value_type, + broadcast=True, ) + # apply mean - normalizer - trend + field += eval_func(func_val=self.mean, **kwargs) + field = self.normalizer.denormalize(field) + field += eval_func(self.trend, **kwargs) + if save: + setattr(self, name, field) + return field def to_pyvista( self, field_select="field", fieldname="field" @@ -299,8 +225,8 @@ def to_pyvista( fieldname : :class:`str`, optional Name of the field in the VTK file. Default: "field" """ - grid = self._to_vtk_helper( - filename=None, field_select=field_select, fieldname=fieldname + grid = to_vtk_helper( + self, filename=None, field_select=field_select, fieldname=fieldname ) return grid @@ -323,8 +249,11 @@ def vtk_export( """ if not isinstance(filename, str): raise TypeError("Please use a string filename.") - return self._to_vtk_helper( - filename=filename, field_select=field_select, fieldname=fieldname + return to_vtk_helper( + self, + filename=filename, + field_select=field_select, + fieldname=fieldname, ) def plot(self, field="field", fig=None, ax=None): # pragma: no cover @@ -334,8 +263,7 @@ def plot(self, field="field", fig=None, ax=None): # pragma: no cover Parameters ---------- field : :class:`str`, optional - Field that should be plotted. Can be: - "field", "raw_field", "krige_field", "err_field" or "krige_var". + Field that should be plotted. Default: "field" fig : :class:`Figure` or :any:`None` Figure to plot the axes on. If `None`, a new one will be created. @@ -384,6 +312,38 @@ def model(self, model): "Field: 'model' is not an instance of 'gstools.CovModel'" ) + @property + def mean(self): + """:class:`float` or :any:`callable`: The mean of the field.""" + return self._mean + + @mean.setter + def mean(self, mean): + self._mean = mean if (callable(mean) or mean is None) else float(mean) + + @property + def normalizer(self): + """:any:`Normalizer`: Normalizer of the field.""" + return self._normalizer + + @normalizer.setter + def normalizer(self, normalizer): + if isinstance(normalizer, Normalizer): + self._normalizer = normalizer + elif normalizer is None: + self._normalizer = Normalizer() + else: + raise ValueError("Field: 'normalizer' not of type 'Normalizer'.") + + @property + def trend(self): + """:class:`float` or :any:`callable`: The trend of the field.""" + return self._trend + + @trend.setter + def trend(self, tren): + self._trend = tren if (callable(tren) or tren is None) else float(tren) + @property def value_type(self): """:class:`str`: Type of the field values (scalar, vector).""" @@ -396,47 +356,32 @@ def value_type(self, value_type): raise ValueError("Field: value type not in {}".format(VALUE_TYPES)) self._value_type = value_type - def __repr__(self): - """Return String representation.""" - return "Field(model={0}, value_type={1})".format( - self.model, self.value_type - ) - + def _fmt_func_val(self, func_val): + if func_val is None: + return str(None) + if callable(func_val): + return "" + return "{0:.{p}}".format(float(func_val), p=self.model._prec) -def _names(name, cnt): - name = [name] if isinstance(name, str) else list(name)[:cnt] - if len(name) < cnt: - name += [name[-1] + str(i + 1) for i in range(cnt - len(name))] - return name + def _fmt_normalizer(self): + norm = self.normalizer + return str(None) if norm.__class__ is Normalizer else norm.name + @property + def name(self): + """:class:`str`: The name of the class.""" + return self.__class__.__name__ -def _get_select(direction): - select = [] - if not (0 < len(direction) < 4): - raise ValueError( - "Field.mesh: need 1 to 3 direction(s), got '{}'".format(direction) - ) - for axis in direction: - if axis == "x": - if 0 in select: - raise ValueError( - "Field.mesh: got duplicate directions {}".format(direction) - ) - select.append(0) - elif axis == "y": - if 1 in select: - raise ValueError( - "Field.mesh: got duplicate directions {}".format(direction) - ) - select.append(1) - elif axis == "z": - if 2 in select: - raise ValueError( - "Field.mesh: got duplicate directions {}".format(direction) - ) - select.append(2) - else: - raise ValueError( - "Field.mesh: got unknown direction {}".format(axis) + def __repr__(self): + """Return String representation.""" + return ( + "{0}(model={1}, value_type='{2}', " + "mean={3}, normalizer={4}, trend={5})".format( + self.name, + self.model.name, + self.value_type, + self._fmt_func_val(self.mean), + self._fmt_normalizer(), + self._fmt_func_val(self.trend), ) - return select + ) From 999e0e1e205e9deb22710c9f9413a4ec7b0c712f Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Thu, 28 Jan 2021 17:02:03 +0100 Subject: [PATCH 18/97] Field: use field_dim in post_field --- gstools/field/base.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/gstools/field/base.py b/gstools/field/base.py index 86019442f..ddea7c08b 100755 --- a/gstools/field/base.py +++ b/gstools/field/base.py @@ -57,6 +57,7 @@ def __init__( # initialize attributes self.pos = None self.mesh_type = None + self.field = None # initialize private attributes self._model = None self._value_type = None @@ -198,7 +199,7 @@ def post_field(self, field, name="field", process=True, save=True): raise ValueError("post_field: no 'pos' tuple set for field.") kwargs = dict( pos=self.pos, - dim=self.model.dim, + dim=self.model.field_dim, mesh_type=self.mesh_type, value_type=self.value_type, broadcast=True, From c06294e280437a98515fe01c33ff402bd7e2d80e Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Thu, 28 Jan 2021 17:13:34 +0100 Subject: [PATCH 19/97] SRF: add normalizer and trend --- gstools/field/srf.py | 46 +++++++++++++++++++++++++------------------- 1 file changed, 26 insertions(+), 20 deletions(-) diff --git a/gstools/field/srf.py b/gstools/field/srf.py index 2f8117800..072ae7c02 100644 --- a/gstools/field/srf.py +++ b/gstools/field/srf.py @@ -40,8 +40,16 @@ class SRF(Field): ---------- model : :any:`CovModel` Covariance Model of the spatial random field. - mean : :class:`float`, optional - mean value of the SRF + mean : :class:`float` or :any:`callable`, optional + Mean of the SRF (in normal form). Could also be a callable. + The default is 0.0. + normalizer : :any:`None` or :any:`Normalizer`, optional + Normalizer to be applied to the SRF to transform the field values. + The default is None. + trend : :any:`None` or :class:`float` or :any:`callable`, optional + Trend of the SRF (in transformed form). + If no normalizer is applied, this behaves equal to 'mean'. + The default is None. upscaling : :class:`str`, optional Method to be used for upscaling the variance at each point depending on the related element volume. @@ -77,18 +85,18 @@ def __init__( self, model, mean=0.0, + normalizer=None, + trend=None, upscaling="no_scaling", generator="RandMeth", **generator_kwargs ): - super().__init__(model) + super().__init__(model, mean=mean, normalizer=normalizer, trend=trend) # initialize private attributes self._generator = None self._upscaling = None self._upscaling_func = None - self._mean = None # initialize attributes - self.mean = mean self.upscaling = upscaling self.set_generator(generator, **generator_kwargs) @@ -125,15 +133,14 @@ def __call__( # get isometrized positions and the resulting field-shape iso_pos, shape = self.pre_pos(pos, mesh_type) # generate the field - self.field = np.reshape(self.generator(iso_pos), shape) + field = np.reshape(self.generator(iso_pos), shape) # upscaled variance if not np.isscalar(point_volumes) or not np.isclose(point_volumes, 0): scaled_var = self.upscaling_func(self.model, point_volumes) if np.size(scaled_var) > 1: scaled_var = np.reshape(scaled_var, shape) - self.field *= np.sqrt(scaled_var / self.model.sill) - self.field += self.mean - return self.field + field *= np.sqrt(scaled_var / self.model.sill) + return self.post_field(field) def upscaling_func(self, *args, **kwargs): """Upscaling method applied to the field variance.""" @@ -184,17 +191,16 @@ def upscaling(self, upscaling): "gstools.SRF: Unknown upscaling method: " + upscaling ) - @property - def mean(self): - """:class:`float`: The mean of the field.""" - return self._mean - - @mean.setter - def mean(self, mean): - self._mean = float(mean) - def __repr__(self): """Return String representation.""" - return "SRF(model={0}, mean={1:.{p}}, generator={2}".format( - self.model, self.mean, self.generator, p=self.model._prec + return ( + "{0}(model={1}, " + "mean={2}, normalizer={3}, trend={4}, generator={5})".format( + self.name, + self.model.name, + self._fmt_func_val(self.mean), + self._fmt_normalizer(), + self._fmt_func_val(self.trend), + self.generator.name, + ) ) From 8446d9aac90441f7dcc281aed28f519e460ea2e5 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Thu, 28 Jan 2021 17:23:34 +0100 Subject: [PATCH 20/97] Krige: add normalizer and trend to base-class; refactoring; add fit_normalizer; use post_field --- gstools/krige/base.py | 210 ++++++++++++++++++------------------------ 1 file changed, 90 insertions(+), 120 deletions(-) diff --git a/gstools/krige/base.py b/gstools/krige/base.py index d242b2611..e3ec485f6 100755 --- a/gstools/krige/base.py +++ b/gstools/krige/base.py @@ -20,12 +20,8 @@ calc_field_krige_and_variance, calc_field_krige, ) -from gstools.krige.tools import ( - set_condition, - get_drift_functions, - no_trend, - eval_func, -) +from gstools.krige.tools import set_condition, get_drift_functions +from gstools.tools.misc import eval_func __all__ = ["Krige"] @@ -55,15 +51,20 @@ class Krige(Field): ext_drift : :class:`numpy.ndarray` or :any:`None`, optional the external drift values at the given cond. positions. - trend_function : :any:`callable`, optional - A callable trend function. Should have the signiture: f(x, [y, z]) + mean : :class:`float`, optional + mean value used to shift normalized conditioning data. + Could also be a callable. The default is None. + normalizer : :any:`None` or :any:`Normalizer`, optional + Normalizer to be applied to the SRF to transform the field values. + The default is None. + trend : :any:`None` or :class:`float` or :any:`callable`, optional + A callable trend function. Should have the signiture: f(x, [y, z, ...]) This is used for detrended kriging, where the trended is subtracted from the conditions before kriging is applied. This can be used for regression kriging, where the trend function is determined by an external regression algorithm. - mean : :class:`float`, optional - mean value used for kriging, if the system is not unbiased. - If unbiased is `True`, the mean will be estimated. + If no normalizer is applied, this behaves equal to 'mean'. + The default is None. unbiased : :class:`bool`, optional Whether the kriging weights should sum up to 1, so the estimator is unbiased. If unbiased is `False` and no drifts are given, @@ -95,6 +96,9 @@ class Krige(Field): If you want to use another routine to invert the kriging matrix, you can pass a callable which takes a matrix and returns the inverse. Default: `1` + fit_normalizer : :class:`bool`, optional + Wheater to fit the data-normalizer to the given conditioning data. + Default: False """ def __init__( @@ -104,15 +108,17 @@ def __init__( cond_val, drift_functions=None, ext_drift=None, - trend_function=None, - mean=0.0, + mean=None, + normalizer=None, + trend=None, unbiased=True, exact=False, cond_err="nugget", pseudo_inv=True, pseudo_inv_type=1, + fit_normalizer=False, ): - super().__init__(model) + super().__init__(model, mean=mean, normalizer=normalizer, trend=trend) self.mean_field = None self.krige_var = None self._unbiased = bool(unbiased) @@ -120,23 +126,20 @@ def __init__( self._pseudo_inv = bool(pseudo_inv) self._pseudo_inv_type = None self.pseudo_inv_type = pseudo_inv_type - self._mean = None - self.mean = mean # initialize private attributes - self._value_type = "scalar" self._cond_pos = None self._cond_val = None self._cond_err = None self._krige_mat = None self._krige_pos = None self._cond_trend = None - self._trend_function = None - self.trend_function = trend_function self._cond_ext_drift = np.array([]) self._drift_functions = [] if drift_functions is not None: self.set_drift_functions(drift_functions) - self.set_condition(cond_pos, cond_val, ext_drift, cond_err) + self.set_condition( + cond_pos, cond_val, ext_drift, cond_err, fit_normalizer + ) def __call__( self, @@ -146,6 +149,7 @@ def __call__( chunk_size=None, only_mean=False, return_var=True, + post_process=True, ): """ Generate the kriging field. @@ -172,6 +176,9 @@ def __call__( return_var : :class:`bool`, optional Whether to return the variance along with the field. Default: `True` + post_process : :class:`bool`, optional + Whether to apply mean, normalizer and trend to the field. + Default: `True` Returns ------- @@ -189,7 +196,8 @@ def __call__( krige_var = np.empty(pnt_cnt, dtype=np.double) if return_var else None # set constant mean if present and wanted if only_mean and self.has_const_mean: - field[...] = self.get_mean() - self.mean # mean is added later + mean = 0.0 if self.mean is None else self.mean + field[...] = self.get_mean() - mean # mean is added later # execute the kriging routine else: # set chunk size @@ -212,23 +220,26 @@ def __call__( self._summate(field, krige_var, c_slice, k_vec, return_var) # reshape field if we got a structured mesh field = np.reshape(field, shape) - if only_mean: - self._post_field(field, "mean_field") - return self.mean_field - self._post_field(field) - if return_var: - krige_var = np.reshape(krige_var, shape) - self.krige_var = np.maximum(self.model.sill - krige_var, 0) - return self.field, self.krige_var + if only_mean: # care about 'kriging the mean' + return self.post_field(field, "mean_field", process=post_process) + # save field to class + field = self.post_field(field, "field", process=post_process) + if return_var: # care about the estimated error variance + krige_var = np.reshape( + np.maximum(self.model.sill - krige_var, 0), shape + ) + krige_var = self.post_field(krige_var, "krige_var", process=False) + return field, krige_var + # if we only calculate the field, overwrite the error variance self.krige_var = None - return self.field + return field def _summate(self, field, krige_var, c_slice, k_vec, return_var): - if return_var: + if return_var: # estimate error variance field[c_slice], krige_var[c_slice] = calc_field_krige_and_variance( self._krige_mat, k_vec, self._krige_cond ) - else: + else: # solely calculate the interpolated field field[c_slice] = calc_field_krige( self._krige_mat, k_vec, self._krige_cond ) @@ -252,8 +263,7 @@ def _get_krige_mat(self): self._get_dists(self._krige_pos) ) # apply the measurement error (nugget by default) - di = np.diag_indices(self.cond_no) - res[di] += self.cond_err + res[np.diag_indices(self.cond_no)] += self.cond_err # set unbias condition (weights have to sum up to 1) if self.unbiased: res[self.cond_no, : self.cond_no] = 1 @@ -343,27 +353,6 @@ def _pre_ext_drift(self, pnt_cnt, ext_drift=None, set_cond=False): raise ValueError("Krige: wrong number of ext. drift values.") return np.array([]) - def _post_field(self, field, name="field"): - """ - Postprocessing and saving of kriging field and error variance. - - Parameters - ---------- - field : :class:`numpy.ndarray` - Raw kriging field. - krige_var : :class:`numpy.ndarray` - Raw kriging error variance. - """ - field += self.mean - if self.trend_function is not no_trend: - field += eval_func( - self.trend_function, - self.pos, - self.model.field_dim, - self.mesh_type, - ) - setattr(self, name, field) - def _get_dists(self, pos1, pos2=None, pos2_slice=(0, None)): """ Calculate pairwise distances. @@ -401,22 +390,26 @@ def get_mean(self): # if there are drift-terms, no constant mean can be calculated -> None if not self.has_const_mean: return None + res = 0.0 # for simple kriging return the given mean + # correctly setting given mean + mean = 0.0 if self.mean is None else self.mean # for ordinary kriging return the estimated mean if self.unbiased: # set the right side of the kriging system to the limit of cov. mean_est = np.concatenate((np.full_like(self.cond_val, 0.0), [1])) # execute the kriging routine with einsum - return ( - np.einsum( - "i,ij,j", self._krige_cond, self._krige_mat, mean_est - ) - + self.mean + res = np.einsum( + "i,ij,j", self._krige_cond, self._krige_mat, mean_est ) - # for simple kriging return the given mean - return self.mean + return res + mean def set_condition( - self, cond_pos, cond_val, ext_drift=None, cond_err="nugget" + self, + cond_pos, + cond_val, + ext_drift=None, + cond_err="nugget", + fit_normalizer=False, ): """Set the conditions for kriging. @@ -437,11 +430,16 @@ def set_condition( The measurement error has to be <= nugget. The "exact=True" variant only works with "cond_err='nugget'". Default: "nugget" + fit_normalizer : :class:`bool`, optional + Wheater to fit the data-normalizer to the given conditioning data. + Default: False """ # correctly format cond_pos and cond_val self._cond_pos, self._cond_val = set_condition( cond_pos, cond_val, self.model.field_dim ) + if fit_normalizer: # fit normalizer to detrended data + self.normalizer.fit(self.cond_val - self.cond_trend) # set the measurement errors self.cond_err = cond_err # set the external drift values and the conditioning points @@ -449,7 +447,9 @@ def set_condition( self.cond_no, ext_drift, set_cond=True ) # upate the internal kriging settings - self.update() + self._krige_pos = self.model.isometrize(self.cond_pos) + # krige pos are the unrotated and isotropic condition positions + self._krige_mat = self._get_krige_mat() def set_drift_functions(self, drift_functions=None): """ @@ -488,26 +488,15 @@ def set_drift_functions(self, drift_functions=None): raise ValueError("Krige: Drift functions not callable") self._drift_functions = drift_functions - def update(self): - """Update the kriging settings.""" - self._krige_pos = self.model.isometrize(self.cond_pos) - # krige pos are the unrotated and isotropic condition positions - self._krige_mat = self._get_krige_mat() - if self.trend_function is no_trend: - self._cond_trend = 0.0 - else: - self._cond_trend = self.trend_function(*self.cond_pos) - @property def _krige_cond(self): """:class:`numpy.ndarray`: The prepared kriging conditions.""" pad_size = self.drift_no + int(self.unbiased) - return np.pad( - self.cond_val - self.mean - self.cond_trend, - (0, pad_size), - mode="constant", - constant_values=0, - ) + # detrend data and normalize + val = self.normalizer.normalize(self.cond_val - self.cond_trend) + # set to zero mean + val -= self.cond_mean + return np.pad(val, (0, pad_size), mode="constant", constant_values=0) @property def cond_pos(self): @@ -534,7 +523,7 @@ def cond_err(self, value): if self.exact: raise ValueError( "krige.cond_err: measurement errors can't be given, " - + "when interpolator should be exact." + "when interpolator should be exact." ) value = np.array(value, dtype=np.double).reshape(-1) if value.size == 1: @@ -551,6 +540,25 @@ def cond_no(self): """:class:`int`: The number of the conditions.""" return len(self._cond_val) + @property + def cond_ext_drift(self): + """:class:`numpy.ndarray`: The ext. drift at the conditions.""" + return self._cond_ext_drift + + @property + def cond_mean(self): + """:class:`numpy.ndarray`: Trend at the conditions.""" + return eval_func( + self.mean, self.cond_pos, self.model.field_dim, broadcast=True + ) + + @property + def cond_trend(self): + """:class:`numpy.ndarray`: Trend at the conditions.""" + return eval_func( + self.trend, self.cond_pos, self.model.field_dim, broadcast=True + ) + @property def unbiased(self): """:class:`bool`: Whether the kriging is unbiased or not.""" @@ -577,11 +585,6 @@ def pseudo_inv_type(self, val): raise ValueError("Krige: pseudo_inv_type needs to be in [1,2,3]") self._pseudo_inv_type = val - @property - def cond_ext_drift(self): - """:class:`numpy.ndarray`: The ext. drift at the conditions.""" - return self._cond_ext_drift - @property def drift_functions(self): """:class:`list` of :any:`callable`: The drift functions.""" @@ -590,7 +593,7 @@ def drift_functions(self): @property def has_const_mean(self): """:class:`bool`: Whether the field has a constant mean or not.""" - return self.drift_no == 0 + return self.drift_no == 0 and not callable(self.mean) @property def krige_size(self): @@ -612,39 +615,6 @@ def ext_drift_no(self): """:class:`int`: Number of external drift values per point.""" return self.cond_ext_drift.shape[0] - @property - def cond_trend(self): - """:class:`numpy.ndarray`: Trend at the conditions.""" - return self._cond_trend - - @property - def trend_function(self): - """:any:`callable`: The trend function.""" - return self._trend_function - - @trend_function.setter - def trend_function(self, trend_function): - if trend_function is None: - trend_function = no_trend - if not callable(trend_function): - raise ValueError("Detrended kriging: trend function not callable.") - self._trend_function = trend_function - # apply trend_function instantly (if cond_pos already set) - if self._cond_pos is not None: - if self._trend_function is no_trend: - self._cond_trend = 0.0 - else: - self._cond_trend = self._trend_function(*self.cond_pos) - - @property - def mean(self): - """:class:`float`: The mean of the field.""" - return self._mean - - @mean.setter - def mean(self, mean): - self._mean = float(mean) - @property def name(self): """:class:`str`: The name of the kriging class.""" From c39687306264ce076e4fc985ce975e44759bc655 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Thu, 28 Jan 2021 17:24:15 +0100 Subject: [PATCH 21/97] Krige: add normalizer and trend to methods (Detrended stays the same) --- gstools/krige/methods.py | 93 +++++++++++++++++++++++++++++++--------- 1 file changed, 72 insertions(+), 21 deletions(-) diff --git a/gstools/krige/methods.py b/gstools/krige/methods.py index 616f3d90d..26f3a3569 100644 --- a/gstools/krige/methods.py +++ b/gstools/krige/methods.py @@ -34,13 +34,19 @@ class Simple(Krige): cond_val : :class:`numpy.ndarray` the values of the conditions mean : :class:`float`, optional - mean value of the kriging field - trend_function : :any:`callable`, optional - A callable trend function. Should have the signiture: f(x, [y, z]) + mean value used to shift normalized conditioning data. + Could also be a callable. The default is None. + normalizer : :any:`None` or :any:`Normalizer`, optional + Normalizer to be applied to the SRF to transform the field values. + The default is None. + trend : :any:`None` or :class:`float` or :any:`callable`, optional + A callable trend function. Should have the signiture: f(x, [y, z, ...]) This is used for detrended kriging, where the trended is subtracted from the conditions before kriging is applied. This can be used for regression kriging, where the trend function is determined by an external regression algorithm. + If no normalizer is applied, this behaves equal to 'mean'. + The default is None. exact : :class:`bool`, optional Whether the interpolator should reproduce the exact input values. If `False`, `cond_err` is interpreted as measurement error @@ -68,6 +74,9 @@ class Simple(Krige): If you want to use another routine to invert the kriging matrix, you can pass a callable which takes a matrix and returns the inverse. Default: `1` + fit_normalizer : :class:`bool`, optional + Wheater to fit the data-normalizer to the given conditioning data. + Default: False """ def __init__( @@ -76,23 +85,27 @@ def __init__( cond_pos, cond_val, mean=0.0, - trend_function=None, + normalizer=None, + trend=None, exact=False, cond_err="nugget", pseudo_inv=True, pseudo_inv_type=1, + fit_normalizer=False, ): super().__init__( model, cond_pos, cond_val, mean=mean, - trend_function=trend_function, + normalizer=normalizer, + trend=trend, unbiased=False, exact=exact, cond_err=cond_err, pseudo_inv=pseudo_inv, pseudo_inv_type=pseudo_inv_type, + fit_normalizer=fit_normalizer, ) @@ -110,12 +123,17 @@ class Ordinary(Krige): tuple, containing the given condition positions (x, [y, z]) cond_val : :class:`numpy.ndarray` the values of the conditions - trend_function : :any:`callable`, optional - A callable trend function. Should have the signiture: f(x, [y, z]) + normalizer : :any:`None` or :any:`Normalizer`, optional + Normalizer to be applied to the SRF to transform the field values. + The default is None. + trend : :any:`None` or :class:`float` or :any:`callable`, optional + A callable trend function. Should have the signiture: f(x, [y, z, ...]) This is used for detrended kriging, where the trended is subtracted from the conditions before kriging is applied. This can be used for regression kriging, where the trend function is determined by an external regression algorithm. + If no normalizer is applied, this behaves equal to 'mean'. + The default is None. exact : :class:`bool`, optional Whether the interpolator should reproduce the exact input values. If `False`, `cond_err` is interpreted as measurement error @@ -143,6 +161,9 @@ class Ordinary(Krige): If you want to use another routine to invert the kriging matrix, you can pass a callable which takes a matrix and returns the inverse. Default: `1` + fit_normalizer : :class:`bool`, optional + Wheater to fit the data-normalizer to the given conditioning data. + Default: False """ def __init__( @@ -150,21 +171,25 @@ def __init__( model, cond_pos, cond_val, - trend_function=None, + normalizer=None, + trend=None, exact=False, cond_err="nugget", pseudo_inv=True, pseudo_inv_type=1, + fit_normalizer=False, ): super().__init__( model, cond_pos, cond_val, - trend_function=trend_function, + trend=trend, + normalizer=normalizer, exact=exact, cond_err=cond_err, pseudo_inv=pseudo_inv, pseudo_inv_type=pseudo_inv_type, + fit_normalizer=fit_normalizer, ) @@ -195,12 +220,17 @@ class Universal(Krige): * "linear" : regional linear drift (equals order=1) * "quadratic" : regional quadratic drift (equals order=2) - trend_function : :any:`callable`, optional - A callable trend function. Should have the signiture: f(x, [y, z]) + normalizer : :any:`None` or :any:`Normalizer`, optional + Normalizer to be applied to the SRF to transform the field values. + The default is None. + trend : :any:`None` or :class:`float` or :any:`callable`, optional + A callable trend function. Should have the signiture: f(x, [y, z, ...]) This is used for detrended kriging, where the trended is subtracted from the conditions before kriging is applied. This can be used for regression kriging, where the trend function is determined by an external regression algorithm. + If no normalizer is applied, this behaves equal to 'mean'. + The default is None. exact : :class:`bool`, optional Whether the interpolator should reproduce the exact input values. If `False`, `cond_err` is interpreted as measurement error @@ -228,6 +258,9 @@ class Universal(Krige): If you want to use another routine to invert the kriging matrix, you can pass a callable which takes a matrix and returns the inverse. Default: `1` + fit_normalizer : :class:`bool`, optional + Wheater to fit the data-normalizer to the given conditioning data. + Default: False """ def __init__( @@ -236,22 +269,26 @@ def __init__( cond_pos, cond_val, drift_functions, - trend_function=None, + normalizer=None, + trend=None, exact=False, cond_err="nugget", pseudo_inv=True, pseudo_inv_type=1, + fit_normalizer=False, ): super().__init__( model, cond_pos, cond_val, drift_functions=drift_functions, - trend_function=trend_function, + normalizer=normalizer, + trend=trend, exact=exact, cond_err=cond_err, pseudo_inv=pseudo_inv, pseudo_inv_type=pseudo_inv_type, + fit_normalizer=fit_normalizer, ) @@ -277,12 +314,17 @@ class ExtDrift(Krige): the values of the conditions ext_drift : :class:`numpy.ndarray` the external drift values at the given condition positions. - trend_function : :any:`callable`, optional - A callable trend function. Should have the signiture: f(x, [y, z]) + normalizer : :any:`None` or :any:`Normalizer`, optional + Normalizer to be applied to the SRF to transform the field values. + The default is None. + trend : :any:`None` or :class:`float` or :any:`callable`, optional + A callable trend function. Should have the signiture: f(x, [y, z, ...]) This is used for detrended kriging, where the trended is subtracted from the conditions before kriging is applied. This can be used for regression kriging, where the trend function is determined by an external regression algorithm. + If no normalizer is applied, this behaves equal to 'mean'. + The default is None. exact : :class:`bool`, optional Whether the interpolator should reproduce the exact input values. If `False`, `cond_err` is interpreted as measurement error @@ -310,6 +352,9 @@ class ExtDrift(Krige): If you want to use another routine to invert the kriging matrix, you can pass a callable which takes a matrix and returns the inverse. Default: `1` + fit_normalizer : :class:`bool`, optional + Wheater to fit the data-normalizer to the given conditioning data. + Default: False """ def __init__( @@ -318,22 +363,26 @@ def __init__( cond_pos, cond_val, ext_drift, - trend_function=None, + normalizer=None, + trend=None, exact=False, cond_err="nugget", pseudo_inv=True, pseudo_inv_type=1, + fit_normalizer=False, ): super().__init__( model, cond_pos, cond_val, ext_drift=ext_drift, - trend_function=trend_function, + normalizer=normalizer, + trend=trend, exact=exact, cond_err=cond_err, pseudo_inv=pseudo_inv, pseudo_inv_type=pseudo_inv_type, + fit_normalizer=fit_normalizer, ) @@ -348,8 +397,10 @@ class Detrended(Krige): This can be used for regression kriging, where the trend function is determined by an external regression algorithm. - This is just a shortcut for simple kriging with a given trend function - and zero mean. A trend can be given with EVERY provided kriging routine. + This is just a shortcut for simple kriging with a given trend function, + zero mean and no normalizer. + + A trend can be given with EVERY provided kriging routine. Parameters ---------- @@ -395,7 +446,7 @@ def __init__( model, cond_pos, cond_val, - trend_function, + trend, exact=False, cond_err="nugget", pseudo_inv=True, @@ -405,7 +456,7 @@ def __init__( model, cond_pos, cond_val, - trend_function=trend_function, + trend=trend, unbiased=False, exact=exact, cond_err=cond_err, From 185bc0c53083effe6ce336ed69a8cf5b3dd0edef Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Thu, 28 Jan 2021 17:24:42 +0100 Subject: [PATCH 22/97] Krige: remove moved routines --- gstools/krige/tools.py | 57 ++---------------------------------------- 1 file changed, 2 insertions(+), 55 deletions(-) diff --git a/gstools/krige/tools.py b/gstools/krige/tools.py index 64448f547..5d220f034 100644 --- a/gstools/krige/tools.py +++ b/gstools/krige/tools.py @@ -9,66 +9,13 @@ .. autosummary:: set_condition get_drift_functions - no_trend - eval_func """ # pylint: disable=C0103 from itertools import combinations_with_replacement import numpy as np -from gstools.tools.geometric import format_struct_pos_dim, gen_mesh -__all__ = ["no_trend", "eval_func", "set_condition", "get_drift_functions"] - - -def no_trend(*args, **kwargs): - """ - Zero trend dummy function. - - Parameters - ---------- - *args : any - Ignored arguments. - **kwargs : any - Ignored keyword arguments. - - Returns - ------- - float - A zero trend given as single float. - - """ - return 0.0 - - -def eval_func(func, pos, dim, mesh_type="structured"): - """ - Evaluate a function on a mesh. - - Parameters - ---------- - func : :any:`callable` - The function to be called. Should have the signiture f(x, [y, z, ...]). - pos : :class:`list` - The position tuple, containing main direction and transversal - directions (x, [y, z, ...]). - dim : :class:`int` - The spatial dimension. - mesh_type : :class:`str`, optional - 'structured' / 'unstructured' - - Returns - ------- - :class:`numpy.ndarray` - Function values at the given points. - """ - if mesh_type != "unstructured": - pos, shape = format_struct_pos_dim(pos, dim) - pos = gen_mesh(pos) - else: - pos = np.array(pos, dtype=np.double).reshape(dim, -1) - shape = np.shape(pos[0]) - return np.reshape(func(*pos), shape) +__all__ = ["set_condition", "get_drift_functions"] def set_condition(cond_pos, cond_val, dim): @@ -102,7 +49,7 @@ def set_condition(cond_pos, cond_val, dim): if len(cond_pos[0]) != len(cond_val): raise ValueError( "Please check your 'cond_pos' and 'cond_val' parameters. " - + "The shapes do not match." + "The shapes do not match." ) return cond_pos, cond_val From 9b8e78266b2163c87f1b19b55d7bd7439df1299d Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Thu, 28 Jan 2021 17:25:26 +0100 Subject: [PATCH 23/97] Tests: change calling sequence in Kriging routines (trend was moved) --- tests/test_krige.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/tests/test_krige.py b/tests/test_krige.py index e0e8a762a..6d4d935d5 100644 --- a/tests/test_krige.py +++ b/tests/test_krige.py @@ -84,7 +84,10 @@ def test_ordinary(self): angles=[2, 1, 0.5], ) ordinary = krige.Ordinary( - model, self.cond_pos[:dim], self.cond_val, trend_func + model, + self.cond_pos[:dim], + self.cond_val, + trend=trend_func, ) field_1, __ = ordinary.unstructured(self.grids[dim - 1]) field_1 = field_1.reshape(self.grid_shape[:dim]) From 4bacbbb77ad5c5b9003a5f0a8191cb521a14a98b Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Thu, 28 Jan 2021 18:01:41 +0100 Subject: [PATCH 24/97] CondSRF: add support for normalizer and trend --- gstools/field/cond_srf.py | 38 +++++++++++++++++++++++++++++++++++--- 1 file changed, 35 insertions(+), 3 deletions(-) diff --git a/gstools/field/cond_srf.py b/gstools/field/cond_srf.py index 6bff44972..46493ae75 100644 --- a/gstools/field/cond_srf.py +++ b/gstools/field/cond_srf.py @@ -56,7 +56,11 @@ def __init__(self, krige, generator="RandMeth", **generator_kwargs): # initialize private attributes self._value_type = None self._generator = None + self._mean = None + self._normalizer = None + self._trend = None # initialize attributes + self.normalizer = None self.set_generator(generator, **generator_kwargs) def __call__(self, pos, seed=np.nan, mesh_type="unstructured", **kwargs): @@ -84,6 +88,7 @@ def __call__(self, pos, seed=np.nan, mesh_type="unstructured", **kwargs): kwargs["mesh_type"] = mesh_type kwargs["only_mean"] = False # overwrite if given kwargs["return_var"] = True # overwrite if given + kwargs["post_process"] = False # overwrite if given # update the model/seed in the generator if any changes were made self.generator.update(self.model, seed) # get isometrized positions and the resulting field-shape @@ -94,8 +99,12 @@ def __call__(self, pos, seed=np.nan, mesh_type="unstructured", **kwargs): ) field, krige_var = self.krige(pos, **kwargs) var_scale, nugget = self.get_scaling(krige_var, shape) - self.field = field + var_scale * self.raw_field + nugget - return self.field + field = self.krige.post_field( + field=field + var_scale * self.raw_field + nugget, save=False + ) + self.krige.post_field(self.krige.field) + # processing is done in the kriging class + return self.post_field(field, process=False) def get_scaling(self, krige_var, shape): """ @@ -159,8 +168,31 @@ def model(self): def model(self, model): pass + @property + def mean(self): + """:class:`float` or :any:`callable`: The mean of the field.""" + return self.krige.mean + + @mean.setter + def mean(self, mean): + pass + + @property + def normalizer(self): + """:any:`Normalizer`: Normalizer of the field.""" + return self.krige.normalizer + + @normalizer.setter + def normalizer(self, normalizer): + pass + + @property + def trend(self): + """:class:`float` or :any:`callable`: The trend of the field.""" + return self._trend + def __repr__(self): """Return String representation.""" return "CondSRF(krige={0}, generator={1})".format( - self.krige, self.generator + self.krige, self.generator.name ) From 78357692861d68a236c06ee943331889ed291d87 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Thu, 28 Jan 2021 18:02:00 +0100 Subject: [PATCH 25/97] Docs: update docs accordingly --- README.md | 8 ++++---- docs/source/index.rst | 8 ++++---- examples/05_kriging/README.rst | 19 ++++++++++++++----- 3 files changed, 22 insertions(+), 13 deletions(-) diff --git a/README.md b/README.md index 54b9412bc..94e253331 100644 --- a/README.md +++ b/README.md @@ -219,15 +219,15 @@ cond_val = [0.47, 0.56, 0.74, 1.47, 1.74] gridx = np.linspace(0.0, 15.0, 151) -# spatial random field class +# conditioned spatial random field class model = gs.Gaussian(dim=1, var=0.5, len_scale=2) -srf = gs.SRF(model) -srf.set_condition(cond_pos, cond_val, "ordinary") +krige = gs.krige.Ordinary(model, cond_pos, cond_val) +cond_srf = gs.CondSRF(krige) # generate the ensemble of field realizations fields = [] for i in range(100): - fields.append(srf(gridx, seed=i)) + fields.append(cond_srf(gridx, seed=i)) plt.plot(gridx, fields[i], color="k", alpha=0.1) plt.scatter(cond_pos, cond_val, color="k") plt.show() diff --git a/docs/source/index.rst b/docs/source/index.rst index 073df0ff8..585935efd 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -253,15 +253,15 @@ generate 100 realizations and plot them: gridx = np.linspace(0.0, 15.0, 151) - # spatial random field class + # conditioned spatial random field class model = gs.Gaussian(dim=1, var=0.5, len_scale=2) - srf = gs.SRF(model) - srf.set_condition(cond_pos, cond_val, "ordinary") + krige = gs.krige.Ordinary(model, cond_pos, cond_val) + cond_srf = gs.CondSRF(krige) # generate the ensemble of field realizations fields = [] for i in range(100): - fields.append(srf(gridx, seed=i)) + fields.append(cond_srf(gridx, seed=i)) plt.plot(gridx, fields[i], color="k", alpha=0.1) plt.scatter(cond_pos, cond_val, color="k") plt.show() diff --git a/examples/05_kriging/README.rst b/examples/05_kriging/README.rst index 07ff852d2..bfbc7e461 100644 --- a/examples/05_kriging/README.rst +++ b/examples/05_kriging/README.rst @@ -33,11 +33,20 @@ the features you want: In contrast to the internal drift, that is evaluated at the desired points with the given functions, the external drift has to given for each point form an "external" source. This results in :any:`ExtDrift` kriging. -* `trend_function`: If you already have fitted a trend model, that is provided as a - callable, you can give it to the kriging routine. This trend is subtracted from the - conditional values before the kriging is done, meaning, that only the residuals are - used for kriging. This can be used with separate regression of your data. - This results in :any:`Detrended` kriging. +* `trend`, `mean`, `normalizer`: These are used to pre- and post-process data. + If you already have fitted a trend model, that is provided as a callable function, + you can give it to the kriging routine. Normalizer are power-transformations + to gain normality. + `mean` behaves similar to `trend` but is applied at another position: + + 1. conditioning data is de-trended (substracting trend) + 2. detrended conditioning data is then normalized (in order to follow a normal distribution) + 3. normalized conditioning data is set to zero mean (subsracting mean) + + Cosequently, when there is no normalizer given, trend and mean are the same thing + and only one should be used. + :any:`Detrended` kriging is a shortcut to provide only a trend and simple kriging + with normal data. * `exact` and `cond_err`: To incorporate the nugget effect and/or measurement errors, one can set `exact` to `False` and provide either individual measurement errors for each point or set the nugget as a constant measurement error everywhere. From 6cda28fa45d1d157e5e9b45977d2181dc50e68ae Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Thu, 28 Jan 2021 18:15:50 +0100 Subject: [PATCH 26/97] Examples: fix call signature in ordinary kriging --- examples/05_kriging/07_detrended_ordinary_kriging.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/05_kriging/07_detrended_ordinary_kriging.py b/examples/05_kriging/07_detrended_ordinary_kriging.py index 585fdcd10..9e468e634 100755 --- a/examples/05_kriging/07_detrended_ordinary_kriging.py +++ b/examples/05_kriging/07_detrended_ordinary_kriging.py @@ -21,7 +21,7 @@ def trend(x): drift_field = drift(gridx) + trend(gridx) # kriging model = Gaussian(dim=1, var=0.1, len_scale=2) -krig_trend = krige.Ordinary(model, cond_pos, cond_val, trend) +krig_trend = krige.Ordinary(model, cond_pos, cond_val, trend=trend) krig_trend(gridx) ax = krig_trend.plot() ax.scatter(cond_pos, cond_val, color="k", zorder=10, label="Conditions") From f626c02d25754bce1edbc189d429ab9f642feebd Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Fri, 29 Jan 2021 11:34:02 +0100 Subject: [PATCH 27/97] Transform: better description to distinguish from normalizer --- gstools/transform/__init__.py | 2 +- gstools/transform/field.py | 7 +++---- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/gstools/transform/__init__.py b/gstools/transform/__init__.py index 1937e6a73..0c52fd6b2 100644 --- a/gstools/transform/__init__.py +++ b/gstools/transform/__init__.py @@ -1,6 +1,6 @@ # -*- coding: utf-8 -*- """ -GStools subpackage providing transformations. +GStools subpackage providing transformations to post-process normal fields. .. currentmodule:: gstools.transform diff --git a/gstools/transform/field.py b/gstools/transform/field.py index 89e88704c..f912d4a16 100644 --- a/gstools/transform/field.py +++ b/gstools/transform/field.py @@ -108,15 +108,14 @@ def discrete(fld, values, thresholds="arithmetic"): if len(values) != len(thresholds) + 1: raise ValueError( "discrete transformation: " - + "len(values) != len(thresholds) + 1" + "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." + "discrete transformation: thresholds need to be ascending." ) # use a separate result so the intermediate results are not affected result = np.empty_like(fld.field) @@ -135,7 +134,7 @@ def discrete(fld, values, thresholds="arithmetic"): def boxcox(fld, lmbda=1, shift=0): """ - Box-Cox transformation. + (Inverse) Box-Cox transformation to denormalize data. After this transformation, the again Box-Cox transformed field is normal distributed. From 19519b20649d52bd480addf15611823667fbbbc6 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Fri, 29 Jan 2021 11:34:45 +0100 Subject: [PATCH 28/97] Normalizer: rename subpackage to 'normalizer' --- docs/source/normalize.rst | 8 -------- docs/source/normalizer.rst | 8 ++++++++ docs/source/package.rst | 2 +- gstools/__init__.py | 6 +++--- gstools/field/base.py | 2 +- gstools/{normalize => normalizer}/__init__.py | 6 +++--- gstools/{normalize => normalizer}/base.py | 2 +- .../{normalize/normalizer.py => normalizer/methods.py} | 4 ++-- 8 files changed, 19 insertions(+), 19 deletions(-) delete mode 100644 docs/source/normalize.rst create mode 100644 docs/source/normalizer.rst rename gstools/{normalize => normalizer}/__init__.py (81%) rename gstools/{normalize => normalizer}/base.py (99%) rename gstools/{normalize/normalizer.py => normalizer/methods.py} (99%) diff --git a/docs/source/normalize.rst b/docs/source/normalize.rst deleted file mode 100644 index 160356404..000000000 --- a/docs/source/normalize.rst +++ /dev/null @@ -1,8 +0,0 @@ -gstools.normalize -================= - -.. automodule:: gstools.normalize - -.. raw:: latex - - \clearpage diff --git a/docs/source/normalizer.rst b/docs/source/normalizer.rst new file mode 100644 index 000000000..396c7ceaf --- /dev/null +++ b/docs/source/normalizer.rst @@ -0,0 +1,8 @@ +gstools.normalizer +================== + +.. automodule:: gstools.normalizer + +.. raw:: latex + + \clearpage diff --git a/docs/source/package.rst b/docs/source/package.rst index 487351a91..792fa9731 100644 --- a/docs/source/package.rst +++ b/docs/source/package.rst @@ -18,4 +18,4 @@ GSTools API random.rst tools.rst transform.rst - normalize.rst + normalizer.rst diff --git a/gstools/__init__.py b/gstools/__init__.py index 200233081..ba6ef832a 100644 --- a/gstools/__init__.py +++ b/gstools/__init__.py @@ -21,7 +21,7 @@ random tools transform - normalize + normalizer Classes ======= @@ -123,7 +123,7 @@ tools, krige, transform, - normalize, + normalizer, ) from gstools.field import SRF, CondSRF from gstools.tools import ( @@ -171,7 +171,7 @@ __all__ = ["__version__"] __all__ += ["covmodel", "field", "variogram", "krige", "random", "tools"] -__all__ += ["transform", "normalize"] +__all__ += ["transform", "normalizer"] __all__ += [ "CovModel", "Gaussian", diff --git a/gstools/field/base.py b/gstools/field/base.py index ddea7c08b..00c20faef 100755 --- a/gstools/field/base.py +++ b/gstools/field/base.py @@ -15,7 +15,7 @@ from gstools.covmodel.base import CovModel from gstools.tools.geometric import format_struct_pos_dim, gen_mesh from gstools.tools.misc import eval_func -from gstools.normalize import Normalizer +from gstools.normalizer import Normalizer from gstools.field.tools import mesh_call, to_vtk_helper __all__ = ["Field"] diff --git a/gstools/normalize/__init__.py b/gstools/normalizer/__init__.py similarity index 81% rename from gstools/normalize/__init__.py rename to gstools/normalizer/__init__.py index 2dab088b2..d25487b2b 100644 --- a/gstools/normalize/__init__.py +++ b/gstools/normalizer/__init__.py @@ -2,7 +2,7 @@ """ GStools subpackage providing normalization routines. -.. currentmodule:: gstools.normalize +.. currentmodule:: gstools.normalizer Base-Normalizer ^^^^^^^^^^^^^^^ @@ -26,8 +26,8 @@ Manly """ -from gstools.normalize.base import Normalizer -from gstools.normalize.normalizer import ( +from gstools.normalizer.base import Normalizer +from gstools.normalizer.methods import ( LogNormal, BoxCox, BoxCoxShift, diff --git a/gstools/normalize/base.py b/gstools/normalizer/base.py similarity index 99% rename from gstools/normalize/base.py rename to gstools/normalizer/base.py index e0645c7d2..12814162c 100644 --- a/gstools/normalize/base.py +++ b/gstools/normalizer/base.py @@ -2,7 +2,7 @@ """ GStools subpackage providing the base class for normalizers. -.. currentmodule:: gstools.normalize.base +.. currentmodule:: gstools.normalizer.base The following classes are provided diff --git a/gstools/normalize/normalizer.py b/gstools/normalizer/methods.py similarity index 99% rename from gstools/normalize/normalizer.py rename to gstools/normalizer/methods.py index 0fd72f2fc..fab09cd1d 100644 --- a/gstools/normalize/normalizer.py +++ b/gstools/normalizer/methods.py @@ -2,7 +2,7 @@ """ GStools subpackage providing different normalizer transformations. -.. currentmodule:: gstools.normalize.normalizer +.. currentmodule:: gstools.normalizer.methods The following classes are provided @@ -15,7 +15,7 @@ Manly """ import numpy as np -from gstools.normalize.base import Normalizer +from gstools.normalizer.base import Normalizer class LogNormal(Normalizer): From 5d541b619a7f2ab584df7a9bd8cdf8c1a80e4cb5 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Fri, 29 Jan 2021 11:38:20 +0100 Subject: [PATCH 29/97] Tests: import fix for normalizer --- tests/test_normalize.py | 64 ++++++++++++++++++++--------------------- 1 file changed, 32 insertions(+), 32 deletions(-) diff --git a/tests/test_normalize.py b/tests/test_normalize.py index f8abd3432..a722f64be 100644 --- a/tests/test_normalize.py +++ b/tests/test_normalize.py @@ -18,21 +18,21 @@ def setUp(self): self.rng = gs.random.RNG(self.seed) self.mean = 11.1 self.std = 2.25 - self.samp = self.rng.random.normal(self.mean, self.std, 1000) + self.smp = self.rng.random.normal(self.mean, self.std, 1000) self.lmb = 1.5 def test_fitting(self): # boxcox with given data to init - bc_samples = gs.normalize.BoxCox(lmbda=self.lmb).denormalize(self.samp) - bc_norm = gs.normalize.BoxCox(data=bc_samples) + bc_samples = gs.normalizer.BoxCox(lmbda=self.lmb).denormalize(self.smp) + bc_norm = gs.normalizer.BoxCox(data=bc_samples) self.assertLess(_rel_err(self.lmb, bc_norm.lmbda), 1e-2) self.assertAlmostEqual( bc_norm.likelihood(bc_samples), np.exp(bc_norm.loglikelihood(bc_samples)), ) # yeo-johnson with calling fit - yj_norm = gs.normalize.YeoJohnson(lmbda=self.lmb) - yj_samples = yj_norm.denormalize(self.samp) + yj_norm = gs.normalizer.YeoJohnson(lmbda=self.lmb) + yj_samples = yj_norm.denormalize(self.smp) yj_norm.fit(yj_samples) self.assertLess(_rel_err(self.lmb, yj_norm.lmbda), 1e-2) self.assertAlmostEqual( @@ -40,8 +40,8 @@ def test_fitting(self): np.exp(yj_norm.loglikelihood(yj_samples)), ) # modulus with calling fit - mo_norm = gs.normalize.Modulus(lmbda=self.lmb) - mo_samples = mo_norm.denormalize(self.samp) + mo_norm = gs.normalizer.Modulus(lmbda=self.lmb) + mo_samples = mo_norm.denormalize(self.smp) mo_norm.fit(mo_samples) self.assertLess(_rel_err(self.lmb, mo_norm.lmbda), 1e-2) self.assertAlmostEqual( @@ -49,8 +49,8 @@ def test_fitting(self): np.exp(mo_norm.loglikelihood(mo_samples)), ) # manly with calling fit - ma_norm = gs.normalize.Manly(lmbda=self.lmb) - ma_samples = ma_norm.denormalize(self.samp) + ma_norm = gs.normalizer.Manly(lmbda=self.lmb) + ma_samples = ma_norm.denormalize(self.smp) ma_norm.fit(ma_samples) self.assertLess(_rel_err(self.lmb, ma_norm.lmbda), 1e-2) # self.assertAlmostEqual( @@ -60,39 +60,39 @@ def test_fitting(self): def test_boxcox(self): # without shift - bc = gs.normalize.BoxCox(lmbda=0) + bc = gs.normalizer.BoxCox(lmbda=0) self.assertTrue( np.all( - np.isclose(self.samp, bc.normalize(bc.denormalize(self.samp))) + np.isclose(self.smp, bc.normalize(bc.denormalize(self.smp))) ) ) bc.lmbda = self.lmb self.assertTrue( np.all( - np.isclose(self.samp, bc.normalize(bc.denormalize(self.samp))) + np.isclose(self.smp, bc.normalize(bc.denormalize(self.smp))) ) ) # with shift - bc = gs.normalize.BoxCoxShift(lmbda=0, shift=1.1) + bc = gs.normalizer.BoxCoxShift(lmbda=0, shift=1.1) self.assertTrue( np.all( - np.isclose(self.samp, bc.normalize(bc.denormalize(self.samp))) + np.isclose(self.smp, bc.normalize(bc.denormalize(self.smp))) ) ) bc.lmbda = self.lmb self.assertTrue( np.all( - np.isclose(self.samp, bc.normalize(bc.denormalize(self.samp))) + np.isclose(self.smp, bc.normalize(bc.denormalize(self.smp))) ) ) def test_yeojohnson(self): - yj = gs.normalize.YeoJohnson(lmbda=0) + yj = gs.normalizer.YeoJohnson(lmbda=0) self.assertTrue( np.all( np.isclose( - self.samp - self.mean, - yj.normalize(yj.denormalize(self.samp - self.mean)), + self.smp - self.mean, + yj.normalize(yj.denormalize(self.smp - self.mean)), ) ) ) @@ -100,8 +100,8 @@ def test_yeojohnson(self): self.assertTrue( np.all( np.isclose( - self.samp - self.mean, - yj.normalize(yj.denormalize(self.samp - self.mean)), + self.smp - self.mean, + yj.normalize(yj.denormalize(self.smp - self.mean)), ) ) ) @@ -110,51 +110,51 @@ def test_yeojohnson(self): self.assertTrue( np.all( np.isclose( - self.samp - self.mean, - yj.normalize(yj.denormalize(self.samp - self.mean)), + self.smp - self.mean, + yj.normalize(yj.denormalize(self.smp - self.mean)), ) ) ) def test_modulus(self): - mo = gs.normalize.Modulus(lmbda=0) + mo = gs.normalizer.Modulus(lmbda=0) self.assertTrue( np.all( - np.isclose(self.samp, mo.normalize(mo.denormalize(self.samp))) + np.isclose(self.smp, mo.normalize(mo.denormalize(self.smp))) ) ) mo.lmbda = self.lmb self.assertTrue( np.all( - np.isclose(self.samp, mo.normalize(mo.denormalize(self.samp))) + np.isclose(self.smp, mo.normalize(mo.denormalize(self.smp))) ) ) def test_manly(self): - ma = gs.normalize.Manly(lmbda=0) + ma = gs.normalizer.Manly(lmbda=0) self.assertTrue( np.all( - np.isclose(self.samp, ma.normalize(ma.denormalize(self.samp))) + np.isclose(self.smp, ma.normalize(ma.denormalize(self.smp))) ) ) ma.lmbda = self.lmb self.assertTrue( np.all( - np.isclose(self.samp, ma.normalize(ma.denormalize(self.samp))) + np.isclose(self.smp, ma.normalize(ma.denormalize(self.smp))) ) ) def test_parameterless(self): - no = gs.normalize.LogNormal() + no = gs.normalizer.LogNormal() self.assertTrue( np.all( - np.isclose(self.samp, no.normalize(no.denormalize(self.samp))) + np.isclose(self.smp, no.normalize(no.denormalize(self.smp))) ) ) - no = gs.normalize.Normalizer() + no = gs.normalizer.Normalizer() self.assertTrue( np.all( - np.isclose(self.samp, no.normalize(no.denormalize(self.samp))) + np.isclose(self.smp, no.normalize(no.denormalize(self.smp))) ) ) From d5d2454a4295ac4c6291e04a49105d5adecb58dd Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Fri, 29 Jan 2021 12:27:10 +0100 Subject: [PATCH 30/97] Examples: add first normalizer example --- .../10_normalizer/00_lognormal_kriging.py | 53 +++++++++++++++++++ examples/10_normalizer/README.rst | 36 +++++++++++++ 2 files changed, 89 insertions(+) create mode 100644 examples/10_normalizer/00_lognormal_kriging.py create mode 100644 examples/10_normalizer/README.rst diff --git a/examples/10_normalizer/00_lognormal_kriging.py b/examples/10_normalizer/00_lognormal_kriging.py new file mode 100644 index 000000000..213a6912c --- /dev/null +++ b/examples/10_normalizer/00_lognormal_kriging.py @@ -0,0 +1,53 @@ +r""" +Log-Normal Kriging +------------------ + +Log Normal kriging is a term to describe a special workflow for kriging to +deal with log-normal data, like conductivity or transmissivity in hydrogeology. + +It simply means to first convert the input data to a normal distribution, i.e. +applying a logarithic function, then interpolating these values with kriging +and transforming the result back with the exponential function. + +The resulting kriging variance describes the error variance of the log-values +of the target variable. + +In this example we will use ordinary kriging. +""" +import numpy as np +import gstools as gs + +# condtions +cond_pos = [0.3, 1.9, 1.1, 3.3, 4.7] +cond_val = [0.47, 0.56, 0.74, 1.47, 1.74] +# resulting grid +gridx = np.linspace(0.0, 15.0, 151) +# stable covariance model +model = gs.Stable(dim=1, var=0.5, len_scale=2.56, alpha=1.9) + +############################################################################### +# In order to result in log-normal kriging, we will use the :any:`LogNormal` +# Normalizer. This is a parameter-less normalizer, so we don't have to fit it. +normalizer = gs.normalizer.LogNormal() + +############################################################################### +# Now we generate the interpolated field as well as the mean field. +# This can be done by setting `only_mean=True` in :any:`Krige.__call__`. +# The result is then stored as `mean_field`. +# +# In terms of log-normal kriging, this mean represents the geometric mean of +# the field. +krige = gs.krige.Ordinary(model, cond_pos, cond_val, normalizer=normalizer) +# interpolate the field +krige(gridx) +# also generate the mean field +krige(gridx, only_mean=True) + +############################################################################### +# And that's it. Let's have a look at the results. +ax = krige.plot() +# plotting the geometric mean +krige.plot("mean_field", ax=ax) +# plotting the conditioning data +ax.scatter(cond_pos, cond_val, color="k", zorder=10, label="Conditions") +ax.legend() diff --git a/examples/10_normalizer/README.rst b/examples/10_normalizer/README.rst new file mode 100644 index 000000000..a3eb89bab --- /dev/null +++ b/examples/10_normalizer/README.rst @@ -0,0 +1,36 @@ +Normalizing Data +================ + +When dealing with real-world data, one can't assume it to be normal distributed. +In fact, many properties are modeled by applying different transformations, +for example conductivity is often assumed to be log-normal or precipitation +is transformed using the famous box-cox power transformation. + +These "normalizers" are often represented as parameteric power transforms and +one is interested in finding the best parameter to gain normality in the input +data. + +This is of special interest when kriging should be applied, since the target +variable of the kriging interpolation is assumed to be normal distributed. + +GSTools provides a set of Normalizers and routines to automatically fit these +to input data by minimizing the likelihood function. + +Provided Normalizers +-------------------- + +The following normalizers can be passed to all Field-classes +(:any:`SRF` or :any:`Krige`) or can be used as standalone tools to analyse data. + +.. currentmodule:: gstools.normalizer + +.. autosummary:: + LogNormal + BoxCox + BoxCoxShift + YeoJohnson + Modulus + Manly + +Examples +-------- From ea759ac388c4cef359d19defcc431d6183802cc4 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Fri, 29 Jan 2021 12:46:12 +0100 Subject: [PATCH 31/97] Examples: update Normalizers --- docs/source/conf.py | 2 ++ examples/10_normalizer/README.rst | 19 +++++++++++++++++++ 2 files changed, 21 insertions(+) diff --git a/docs/source/conf.py b/docs/source/conf.py index fbdfd2fc0..b09889486 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -292,6 +292,7 @@ def setup(app): "../../examples/07_transformations/", "../../examples/08_geo_coordinates/", "../../examples/09_spatio_temporal/", + "../../examples/10_normalizer/", ], # path where to save gallery generated examples "gallery_dirs": [ @@ -305,6 +306,7 @@ def setup(app): "examples/07_transformations/", "examples/08_geo_coordinates/", "examples/09_spatio_temporal/", + "examples/10_normalizer/", ], # Pattern to search for example files "filename_pattern": r"\.py", diff --git a/examples/10_normalizer/README.rst b/examples/10_normalizer/README.rst index a3eb89bab..9f86259e5 100644 --- a/examples/10_normalizer/README.rst +++ b/examples/10_normalizer/README.rst @@ -16,6 +16,25 @@ variable of the kriging interpolation is assumed to be normal distributed. GSTools provides a set of Normalizers and routines to automatically fit these to input data by minimizing the likelihood function. +Mean, Trend and Normalizers +--------------------------- + +All Field classes (:any:`SRF` or :any:`Krige`) provide the input of `mean`, +`normalizer` and `trend`. + +A `trend` can be a callable function, that represents a trend in input data. +For example a linear decrease of temperature with height. + +The `normalizer` will be applied after the data was detrended, i.e. the trend +was substracted from the data. + +The `mean` is now interpreted as the mean of the normalized data. Users +could also provide a callable mean, but it is mostly meant to be constant. + +When no normalizer is given, `trend` and `mean` basically behave the same. +We just decided that a trend is associated with raw data and a mean is used +in the context of normally distributed data. + Provided Normalizers -------------------- From 277489968b404b5679d84c4989e746969c6c8f08 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Fri, 29 Jan 2021 12:52:18 +0100 Subject: [PATCH 32/97] Examples: use bullet points to describe mean, trend, normalizer --- examples/10_normalizer/README.rst | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/10_normalizer/README.rst b/examples/10_normalizer/README.rst index 9f86259e5..f13432992 100644 --- a/examples/10_normalizer/README.rst +++ b/examples/10_normalizer/README.rst @@ -20,15 +20,15 @@ Mean, Trend and Normalizers --------------------------- All Field classes (:any:`SRF` or :any:`Krige`) provide the input of `mean`, -`normalizer` and `trend`. +`normalizer` and `trend`: -A `trend` can be a callable function, that represents a trend in input data. +* A `trend` can be a callable function, that represents a trend in input data. For example a linear decrease of temperature with height. -The `normalizer` will be applied after the data was detrended, i.e. the trend +* The `normalizer` will be applied after the data was detrended, i.e. the trend was substracted from the data. -The `mean` is now interpreted as the mean of the normalized data. Users +* The `mean` is now interpreted as the mean of the normalized data. Users could also provide a callable mean, but it is mostly meant to be constant. When no normalizer is given, `trend` and `mean` basically behave the same. From ab3969a255b69e21d102511a04309981ddec3e04 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Fri, 29 Jan 2021 12:52:49 +0100 Subject: [PATCH 33/97] CondSRF: bug fix for trend property --- gstools/field/cond_srf.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/gstools/field/cond_srf.py b/gstools/field/cond_srf.py index 46493ae75..917a5a680 100644 --- a/gstools/field/cond_srf.py +++ b/gstools/field/cond_srf.py @@ -189,7 +189,11 @@ def normalizer(self, normalizer): @property def trend(self): """:class:`float` or :any:`callable`: The trend of the field.""" - return self._trend + return self.krige.trend + + @trend.setter + def trend(self, tren): + pass def __repr__(self): """Return String representation.""" From db681d8e61903b8894bf0e7386671c605b31a49d Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Fri, 29 Jan 2021 13:53:16 +0100 Subject: [PATCH 34/97] CondSRF: consequently forward properties to kriging class --- gstools/field/base.py | 1 - gstools/field/cond_srf.py | 17 +++++++++++++---- 2 files changed, 13 insertions(+), 5 deletions(-) diff --git a/gstools/field/base.py b/gstools/field/base.py index 00c20faef..9c22d12cb 100755 --- a/gstools/field/base.py +++ b/gstools/field/base.py @@ -352,7 +352,6 @@ def value_type(self): @value_type.setter def value_type(self, value_type): - """:class:`str`: Type of the field values (scalar, vector).""" if value_type not in VALUE_TYPES: raise ValueError("Field: value type not in {}".format(VALUE_TYPES)) self._value_type = value_type diff --git a/gstools/field/cond_srf.py b/gstools/field/cond_srf.py index 917a5a680..74db5afa9 100644 --- a/gstools/field/cond_srf.py +++ b/gstools/field/cond_srf.py @@ -166,7 +166,7 @@ def model(self): @model.setter def model(self, model): - pass + self.krige.model = model @property def mean(self): @@ -175,7 +175,7 @@ def mean(self): @mean.setter def mean(self, mean): - pass + self.krige.mean = mean @property def normalizer(self): @@ -184,7 +184,7 @@ def normalizer(self): @normalizer.setter def normalizer(self, normalizer): - pass + self.krige.normalizer = normalizer @property def trend(self): @@ -193,7 +193,16 @@ def trend(self): @trend.setter def trend(self, tren): - pass + self.krige.trend = tren + + @property + def value_type(self): + """:class:`str`: Type of the field values (scalar, vector).""" + return self.krige.value_type + + @value_type.setter + def value_type(self, value_type): + self.krige.value_type = value_type def __repr__(self): """Return String representation.""" From bbf19e7e7be73a65ced18dd592a7f41bf608ee5f Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Fri, 29 Jan 2021 13:55:10 +0100 Subject: [PATCH 35/97] Krige: fix external drift bug (when using ext+int drift); adding default values for set_condition; add deleter for cond_ext_drift --- gstools/krige/base.py | 42 ++++++++++++++++++++++++++++++++---------- 1 file changed, 32 insertions(+), 10 deletions(-) diff --git a/gstools/krige/base.py b/gstools/krige/base.py index e3ec485f6..2670116bb 100755 --- a/gstools/krige/base.py +++ b/gstools/krige/base.py @@ -343,8 +343,8 @@ def _pre_ext_drift(self, pnt_cnt, ext_drift=None, set_cond=False): return ext_drift ext_drift = np.array(ext_drift, dtype=np.double, ndmin=2) ext_shape = np.shape(ext_drift) - shape = (self.drift_no, pnt_cnt) - if self.drift_no > 1 and ext_shape[0] != self.drift_no: + shape = (self.ext_drift_no, pnt_cnt) + if self.drift_no > 1 and ext_shape[0] != self.ext_drift_no: raise ValueError("Krige: wrong number of external drifts.") if np.prod(ext_shape) != np.prod(shape): raise ValueError("Krige: wrong number of ext. drift values.") @@ -405,24 +405,28 @@ def get_mean(self): def set_condition( self, - cond_pos, - cond_val, + cond_pos=None, + cond_val=None, ext_drift=None, - cond_err="nugget", + cond_err=None, fit_normalizer=False, ): """Set the conditions for kriging. + This method could also be used to update the kriging setup, when + properties were changed. Then you can call it without arguments. + Parameters ---------- - cond_pos : :class:`list` - the position tuple of the conditions (x, [y, z]) - cond_val : :class:`numpy.ndarray` - the values of the conditions + cond_pos : :class:`list`, optional + the position tuple of the conditions (x, [y, z]). Default: current. + cond_val : :class:`numpy.ndarray`, optional + the values of the conditions. Default: current. ext_drift : :class:`numpy.ndarray` or :any:`None`, optional the external drift values at the given conditions (only for EDK) For multiple external drifts, the first dimension - should be the index of the drift term. + should be the index of the drift term. When passing `None`, the + extisting external drift will be used. cond_err : :class:`str`, :class :class:`float`, :class:`list`, optional The measurement error at the conditioning points. Either "nugget" to apply the model-nugget, a single value applied @@ -433,7 +437,21 @@ def set_condition( fit_normalizer : :class:`bool`, optional Wheater to fit the data-normalizer to the given conditioning data. Default: False + + Notes + ----- + In order to remove an existing external drift use: + + >>> del Krige.cond_ext_drift """ + # Default values + cond_pos = self._cond_pos if cond_pos is None else cond_pos + cond_val = self._cond_val if cond_val is None else cond_val + ext_drift = self._cond_ext_drift if ext_drift is None else ext_drift + cond_err = self._cond_err if cond_err is None else cond_err + cond_err = "nugget" if cond_err is None else cond_err + if cond_pos is None or cond_val is None: + raise ValueError("Krige.set_condition: missing cond_pos/cond_val.") # correctly format cond_pos and cond_val self._cond_pos, self._cond_val = set_condition( cond_pos, cond_val, self.model.field_dim @@ -545,6 +563,10 @@ def cond_ext_drift(self): """:class:`numpy.ndarray`: The ext. drift at the conditions.""" return self._cond_ext_drift + @cond_ext_drift.deleter + def cond_ext_drift(self): + self._cond_ext_drift = np.array([]) + @property def cond_mean(self): """:class:`numpy.ndarray`: Trend at the conditions.""" From 5e169a5af4e421a7fd8244fe86049f4406207196 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Fri, 29 Jan 2021 14:18:09 +0100 Subject: [PATCH 36/97] Krige: add fit_variogram argument to automatically fit covmodel to given data --- gstools/krige/base.py | 43 +++++++++++++++++++++++++++++++++++++++- gstools/krige/methods.py | 40 +++++++++++++++++++++++++++++++++++++ 2 files changed, 82 insertions(+), 1 deletion(-) diff --git a/gstools/krige/base.py b/gstools/krige/base.py index 2670116bb..caae4473c 100755 --- a/gstools/krige/base.py +++ b/gstools/krige/base.py @@ -22,6 +22,8 @@ ) from gstools.krige.tools import set_condition, get_drift_functions from gstools.tools.misc import eval_func +from gstools.tools.geometric import rotated_main_axes +from gstools.variogram import vario_estimate __all__ = ["Krige"] @@ -99,6 +101,17 @@ class Krige(Field): fit_normalizer : :class:`bool`, optional Wheater to fit the data-normalizer to the given conditioning data. Default: False + fit_variogram : :class:`bool`, optional + Wheater to fit the given variogram model to the data. + This is done by using isotropy settings of the given model, + assuming the sill to be the data variance and with the + standard bins provided by the :any:`variogram` submodule. + Default: False + + Notes + ----- + If you have changed any properties in the class, you can update the kriging + setup by calling :any:`Krige.set_condition` without any arguments. """ def __init__( @@ -117,6 +130,7 @@ def __init__( pseudo_inv=True, pseudo_inv_type=1, fit_normalizer=False, + fit_variogram=False, ): super().__init__(model, mean=mean, normalizer=normalizer, trend=trend) self.mean_field = None @@ -138,7 +152,12 @@ def __init__( if drift_functions is not None: self.set_drift_functions(drift_functions) self.set_condition( - cond_pos, cond_val, ext_drift, cond_err, fit_normalizer + cond_pos, + cond_val, + ext_drift, + cond_err, + fit_normalizer, + fit_variogram, ) def __call__( @@ -410,6 +429,7 @@ def set_condition( ext_drift=None, cond_err=None, fit_normalizer=False, + fit_variogram=False, ): """Set the conditions for kriging. @@ -437,6 +457,12 @@ def set_condition( fit_normalizer : :class:`bool`, optional Wheater to fit the data-normalizer to the given conditioning data. Default: False + fit_variogram : :class:`bool`, optional + Wheater to fit the given variogram model to the data. + This is done by using isotropy settings of the given model, + assuming the sill to be the data variance and with the + standard bins provided by the :any:`variogram` submodule. + Default: False Notes ----- @@ -458,6 +484,21 @@ def set_condition( ) if fit_normalizer: # fit normalizer to detrended data self.normalizer.fit(self.cond_val - self.cond_trend) + if fit_variogram: # fitting model to empirical variogram of data + # normalize field + field = self.normalizer.normalize(self.cond_val - self.cond_trend) + field -= self.cond_mean + sill = np.var(field) + if self.model.is_isotropic: + emp_vario = vario_estimate( + self.cond_pos, field, latlon=self.model.latlon + ) + else: + axes = rotated_main_axes(self.model.dim, self.model.angles) + emp_vario = vario_estimate( + self.cond_pos, field, direction=axes + ) + self.model.fit_variogram(*emp_vario, sill=sill) # set the measurement errors self.cond_err = cond_err # set the external drift values and the conditioning points diff --git a/gstools/krige/methods.py b/gstools/krige/methods.py index 26f3a3569..def09a7d5 100644 --- a/gstools/krige/methods.py +++ b/gstools/krige/methods.py @@ -77,6 +77,12 @@ class Simple(Krige): fit_normalizer : :class:`bool`, optional Wheater to fit the data-normalizer to the given conditioning data. Default: False + fit_variogram : :class:`bool`, optional + Wheater to fit the given variogram model to the data. + This is done by using isotropy settings of the given model, + assuming the sill to be the data variance and with the + standard bins provided by the :any:`variogram` submodule. + Default: False """ def __init__( @@ -92,6 +98,7 @@ def __init__( pseudo_inv=True, pseudo_inv_type=1, fit_normalizer=False, + fit_variogram=False, ): super().__init__( model, @@ -106,6 +113,7 @@ def __init__( pseudo_inv=pseudo_inv, pseudo_inv_type=pseudo_inv_type, fit_normalizer=fit_normalizer, + fit_variogram=fit_variogram, ) @@ -164,6 +172,12 @@ class Ordinary(Krige): fit_normalizer : :class:`bool`, optional Wheater to fit the data-normalizer to the given conditioning data. Default: False + fit_variogram : :class:`bool`, optional + Wheater to fit the given variogram model to the data. + This is done by using isotropy settings of the given model, + assuming the sill to be the data variance and with the + standard bins provided by the :any:`variogram` submodule. + Default: False """ def __init__( @@ -178,6 +192,7 @@ def __init__( pseudo_inv=True, pseudo_inv_type=1, fit_normalizer=False, + fit_variogram=False, ): super().__init__( model, @@ -190,6 +205,7 @@ def __init__( pseudo_inv=pseudo_inv, pseudo_inv_type=pseudo_inv_type, fit_normalizer=fit_normalizer, + fit_variogram=fit_variogram, ) @@ -261,6 +277,12 @@ class Universal(Krige): fit_normalizer : :class:`bool`, optional Wheater to fit the data-normalizer to the given conditioning data. Default: False + fit_variogram : :class:`bool`, optional + Wheater to fit the given variogram model to the data. + This is done by using isotropy settings of the given model, + assuming the sill to be the data variance and with the + standard bins provided by the :any:`variogram` submodule. + Default: False """ def __init__( @@ -276,6 +298,7 @@ def __init__( pseudo_inv=True, pseudo_inv_type=1, fit_normalizer=False, + fit_variogram=False, ): super().__init__( model, @@ -289,6 +312,7 @@ def __init__( pseudo_inv=pseudo_inv, pseudo_inv_type=pseudo_inv_type, fit_normalizer=fit_normalizer, + fit_variogram=fit_variogram, ) @@ -355,6 +379,12 @@ class ExtDrift(Krige): fit_normalizer : :class:`bool`, optional Wheater to fit the data-normalizer to the given conditioning data. Default: False + fit_variogram : :class:`bool`, optional + Wheater to fit the given variogram model to the data. + This is done by using isotropy settings of the given model, + assuming the sill to be the data variance and with the + standard bins provided by the :any:`variogram` submodule. + Default: False """ def __init__( @@ -370,6 +400,7 @@ def __init__( pseudo_inv=True, pseudo_inv_type=1, fit_normalizer=False, + fit_variogram=False, ): super().__init__( model, @@ -383,6 +414,7 @@ def __init__( pseudo_inv=pseudo_inv, pseudo_inv_type=pseudo_inv_type, fit_normalizer=fit_normalizer, + fit_variogram=fit_variogram, ) @@ -439,6 +471,12 @@ class Detrended(Krige): If you want to use another routine to invert the kriging matrix, you can pass a callable which takes a matrix and returns the inverse. Default: `1` + fit_variogram : :class:`bool`, optional + Wheater to fit the given variogram model to the data. + This is done by using isotropy settings of the given model, + assuming the sill to be the data variance and with the + standard bins provided by the :any:`variogram` submodule. + Default: False """ def __init__( @@ -451,6 +489,7 @@ def __init__( cond_err="nugget", pseudo_inv=True, pseudo_inv_type=1, + fit_variogram=False, ): super().__init__( model, @@ -462,4 +501,5 @@ def __init__( cond_err=cond_err, pseudo_inv=pseudo_inv, pseudo_inv_type=pseudo_inv_type, + fit_variogram=fit_variogram, ) From 542c0cb0290a3ffe9145724f3e0cf8daad347e03 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Fri, 29 Jan 2021 14:28:18 +0100 Subject: [PATCH 37/97] Krige: bugfix in ext_drift setting --- gstools/krige/base.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/gstools/krige/base.py b/gstools/krige/base.py index caae4473c..aefceb8b2 100755 --- a/gstools/krige/base.py +++ b/gstools/krige/base.py @@ -57,7 +57,7 @@ class Krige(Field): mean value used to shift normalized conditioning data. Could also be a callable. The default is None. normalizer : :any:`None` or :any:`Normalizer`, optional - Normalizer to be applied to the SRF to transform the field values. + Normalizer to be applied to the input data to gain normality. The default is None. trend : :any:`None` or :class:`float` or :any:`callable`, optional A callable trend function. Should have the signiture: f(x, [y, z, ...]) @@ -355,12 +355,13 @@ def _pre_ext_drift(self, pnt_cnt, ext_drift=None, set_cond=False): the drift values at the given positions """ if ext_drift is not None: + ext_drift = np.array(ext_drift, dtype=np.double, ndmin=2) + if ext_drift.size == 0: # treat empty array as no ext_drift + return np.array([]) if set_cond: - ext_drift = np.array(ext_drift, dtype=np.double, ndmin=2) if len(ext_drift.shape) > 2 or ext_drift.shape[1] != pnt_cnt: raise ValueError("Krige: wrong number of ext. drifts.") return ext_drift - ext_drift = np.array(ext_drift, dtype=np.double, ndmin=2) ext_shape = np.shape(ext_drift) shape = (self.ext_drift_no, pnt_cnt) if self.drift_no > 1 and ext_shape[0] != self.ext_drift_no: From 50dfe664d393e91fc15d7b1d4d555e635b9b6b60 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Fri, 29 Jan 2021 14:28:27 +0100 Subject: [PATCH 38/97] Krige: better doc for normalizer --- gstools/krige/methods.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/gstools/krige/methods.py b/gstools/krige/methods.py index def09a7d5..ee3534a9d 100644 --- a/gstools/krige/methods.py +++ b/gstools/krige/methods.py @@ -37,7 +37,7 @@ class Simple(Krige): mean value used to shift normalized conditioning data. Could also be a callable. The default is None. normalizer : :any:`None` or :any:`Normalizer`, optional - Normalizer to be applied to the SRF to transform the field values. + Normalizer to be applied to the input data to gain normality. The default is None. trend : :any:`None` or :class:`float` or :any:`callable`, optional A callable trend function. Should have the signiture: f(x, [y, z, ...]) @@ -132,7 +132,7 @@ class Ordinary(Krige): cond_val : :class:`numpy.ndarray` the values of the conditions normalizer : :any:`None` or :any:`Normalizer`, optional - Normalizer to be applied to the SRF to transform the field values. + Normalizer to be applied to the input data to gain normality. The default is None. trend : :any:`None` or :class:`float` or :any:`callable`, optional A callable trend function. Should have the signiture: f(x, [y, z, ...]) @@ -237,7 +237,7 @@ class Universal(Krige): * "quadratic" : regional quadratic drift (equals order=2) normalizer : :any:`None` or :any:`Normalizer`, optional - Normalizer to be applied to the SRF to transform the field values. + Normalizer to be applied to the input data to gain normality. The default is None. trend : :any:`None` or :class:`float` or :any:`callable`, optional A callable trend function. Should have the signiture: f(x, [y, z, ...]) @@ -339,7 +339,7 @@ class ExtDrift(Krige): ext_drift : :class:`numpy.ndarray` the external drift values at the given condition positions. normalizer : :any:`None` or :any:`Normalizer`, optional - Normalizer to be applied to the SRF to transform the field values. + Normalizer to be applied to the input data to gain normality. The default is None. trend : :any:`None` or :class:`float` or :any:`callable`, optional A callable trend function. Should have the signiture: f(x, [y, z, ...]) From 0fed22f133fe26696b1909f3306c94a72472237f Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Fri, 29 Jan 2021 14:35:29 +0100 Subject: [PATCH 39/97] Docs: add normalizer tutorial to toctree --- docs/source/tutorials.rst | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/source/tutorials.rst b/docs/source/tutorials.rst index bf89d159c..28d39e563 100644 --- a/docs/source/tutorials.rst +++ b/docs/source/tutorials.rst @@ -21,4 +21,5 @@ explore its whole beauty and power. examples/07_transformations/index examples/08_geo_coordinates/index examples/09_spatio_temporal/index + examples/10_normalizer/index examples/00_misc/index From 1f93a8e06442f3ea4261cac93d1dab415750020a Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Fri, 29 Jan 2021 15:43:13 +0100 Subject: [PATCH 40/97] Examples: update normalizer readme --- examples/10_normalizer/README.rst | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/10_normalizer/README.rst b/examples/10_normalizer/README.rst index f13432992..972a1c3d8 100644 --- a/examples/10_normalizer/README.rst +++ b/examples/10_normalizer/README.rst @@ -23,13 +23,13 @@ All Field classes (:any:`SRF` or :any:`Krige`) provide the input of `mean`, `normalizer` and `trend`: * A `trend` can be a callable function, that represents a trend in input data. -For example a linear decrease of temperature with height. + For example a linear decrease of temperature with height. * The `normalizer` will be applied after the data was detrended, i.e. the trend -was substracted from the data. + was substracted from the data, in order to gain normality. -* The `mean` is now interpreted as the mean of the normalized data. Users -could also provide a callable mean, but it is mostly meant to be constant. +* The `mean` is now interpreted as the mean of the normalized data. The user + could also provide a callable mean, but it is mostly meant to be constant. When no normalizer is given, `trend` and `mean` basically behave the same. We just decided that a trend is associated with raw data and a mean is used From 788e87b7547c0209cfeba4d6f23c5a313e9502ad Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Fri, 29 Jan 2021 15:43:56 +0100 Subject: [PATCH 41/97] Examples: add second normalizer example on auto fitting --- examples/10_normalizer/01_auto_fit.py | 83 +++++++++++++++++++++++++++ 1 file changed, 83 insertions(+) create mode 100644 examples/10_normalizer/01_auto_fit.py diff --git a/examples/10_normalizer/01_auto_fit.py b/examples/10_normalizer/01_auto_fit.py new file mode 100644 index 000000000..0461c0f5b --- /dev/null +++ b/examples/10_normalizer/01_auto_fit.py @@ -0,0 +1,83 @@ +""" +Automatically fitting +--------------------- + +Here we generate synthetic lognormal data, that should be interpolated +Here we transform a field to a log-normal distribution: +""" +import numpy as np +import gstools as gs +import matplotlib.pyplot as plt + +# structured field with a size of 60x60 +x = y = range(60) +pos = gs.tools.geometric.gen_mesh([x, y]) +model = gs.Gaussian(dim=2, var=1, len_scale=10) +srf = gs.SRF(model, seed=20170519, normalizer=gs.normalizer.LogNormal()) +srf(pos) + +############################################################################### +# Now we want to interpolate the given field to finer resolution (down-scaling) +# and we want to normalize the given data with the BoxCox transformation. +# +# First, we convert the given condtions to be unstructured +# and set the output resolution. +# +# To make the data more scarce, we will sample 60 points. + +ids = np.arange(srf.field.size) +sampled = np.random.RandomState(20210201).choice(ids, size=60, replace=False) + +# sample conditioning points from generated field +cond_pos = pos[:, sampled] +cond_val = srf.field[sampled] + +############################################################################### +# Now we set up the kriging routine. We use a :any:`Stable` model, that should +# be fitted to the given data (by setting ``fit_variogram=True``) and we +# use the :any:`BoxCox` normalizer in order to gain normality. +# +# The BoxCox normalizer will be fitted automatically to the data, +# by setting ``fit_normalizer=True``. + +krige = gs.krige.Ordinary( + model=gs.Stable(dim=2), + cond_pos=cond_pos, + cond_val=cond_val, + normalizer=gs.normalizer.BoxCox(), + fit_normalizer=True, + fit_variogram=True, +) + +############################################################################### +# First, let's have a look at the fitting results: + +print(krige.model) +print(krige.normalizer) + +############################################################################### +# As we see, it went quite well. Variance is a bit underestimated, but +# length scale and nugget are good. The shape parameter of the stable model +# is correctly estimated to be close to `2`, +# so we result in a gaussian like model. +# +# The BoxCox parameter `lmbda` was estimated to be almost 0, which means, +# the log-normal distribution was correctly fitted. +# +# Now let's run the kriging interpolation and let's have a look at the +# resulting field. We will also generate the original field for comparison. + +# interpolate +krige(pos) + +# plot +fig, ax = plt.subplots(1, 3, figsize=[8, 3]) +ax[0].imshow(srf.field.reshape(60, 60).T, origin="lower") +ax[1].scatter(*cond_pos, c=cond_val) +ax[2].imshow(krige.field.reshape(60, 60).T, origin="lower") +# titles +ax[0].set_title("original field") +ax[1].set_title("sampled field") +ax[2].set_title("interpolated field") +# set aspect ratio to equal in all plots +[ax[i].set_aspect("equal") for i in range(3)] From e4a1be5097c8386863ce0e55850ba418c867566b Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Fri, 29 Jan 2021 15:53:50 +0100 Subject: [PATCH 42/97] Examples: update auto-fit normalizer example --- examples/10_normalizer/01_auto_fit.py | 28 +++++++++++++++------------ 1 file changed, 16 insertions(+), 12 deletions(-) diff --git a/examples/10_normalizer/01_auto_fit.py b/examples/10_normalizer/01_auto_fit.py index 0461c0f5b..df40672a5 100644 --- a/examples/10_normalizer/01_auto_fit.py +++ b/examples/10_normalizer/01_auto_fit.py @@ -1,9 +1,14 @@ """ -Automatically fitting ---------------------- +Automatic fitting +----------------- -Here we generate synthetic lognormal data, that should be interpolated -Here we transform a field to a log-normal distribution: +In order to demonstrate how to automatically fit normalizer and variograms, +we generate synthetic lognormal data, that should be interpolated with +ordinary kriging. + +Normalizers are fitted by minimizing the likelihood function and variograms +are fitted by estimating the empirical variogram with automatic binning and +fitting the theoretical model to it. """ import numpy as np import gstools as gs @@ -14,23 +19,22 @@ pos = gs.tools.geometric.gen_mesh([x, y]) model = gs.Gaussian(dim=2, var=1, len_scale=10) srf = gs.SRF(model, seed=20170519, normalizer=gs.normalizer.LogNormal()) +# generate the original field srf(pos) ############################################################################### -# Now we want to interpolate the given field to finer resolution (down-scaling) +# Now we want to interpolate sampled data from the given field (in order to +# pretend, we got any measured real-world data) # and we want to normalize the given data with the BoxCox transformation. # -# First, we convert the given condtions to be unstructured -# and set the output resolution. -# -# To make the data more scarce, we will sample 60 points. +# Here, we will sample 60 points and set the conditioning points and values. ids = np.arange(srf.field.size) -sampled = np.random.RandomState(20210201).choice(ids, size=60, replace=False) +samples = np.random.RandomState(20210201).choice(ids, size=60, replace=False) # sample conditioning points from generated field -cond_pos = pos[:, sampled] -cond_val = srf.field[sampled] +cond_pos = pos[:, samples] +cond_val = srf.field[samples] ############################################################################### # Now we set up the kriging routine. We use a :any:`Stable` model, that should From 96d4a077bb5f76261fc1405ef4359ab5e227c312 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Fri, 29 Jan 2021 16:00:56 +0100 Subject: [PATCH 43/97] Examples: update auto-fit normalizer example again --- examples/10_normalizer/01_auto_fit.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/10_normalizer/01_auto_fit.py b/examples/10_normalizer/01_auto_fit.py index df40672a5..0681cf0cc 100644 --- a/examples/10_normalizer/01_auto_fit.py +++ b/examples/10_normalizer/01_auto_fit.py @@ -69,7 +69,7 @@ # the log-normal distribution was correctly fitted. # # Now let's run the kriging interpolation and let's have a look at the -# resulting field. We will also generate the original field for comparison. +# resulting field. We'll compare the original, sampled and interpolated field. # interpolate krige(pos) From 9520b35dde4227ec939a292f49dba2fbe945339f Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Fri, 29 Jan 2021 17:15:04 +0100 Subject: [PATCH 44/97] Examples: update auto-fit normalizer example again --- examples/10_normalizer/01_auto_fit.py | 39 ++++++++++++++++++--------- 1 file changed, 27 insertions(+), 12 deletions(-) diff --git a/examples/10_normalizer/01_auto_fit.py b/examples/10_normalizer/01_auto_fit.py index 0681cf0cc..bf3db7744 100644 --- a/examples/10_normalizer/01_auto_fit.py +++ b/examples/10_normalizer/01_auto_fit.py @@ -9,6 +9,13 @@ Normalizers are fitted by minimizing the likelihood function and variograms are fitted by estimating the empirical variogram with automatic binning and fitting the theoretical model to it. + +Artificial data +^^^^^^^^^^^^^^^ + +Here we generate log-normal data following a gaussian covariance model. +We will generate the "original" field on a 60x60 mesh, where we will take +samples from in order to pretend a data-scarcity situation. """ import numpy as np import gstools as gs @@ -23,11 +30,7 @@ srf(pos) ############################################################################### -# Now we want to interpolate sampled data from the given field (in order to -# pretend, we got any measured real-world data) -# and we want to normalize the given data with the BoxCox transformation. -# -# Here, we will sample 60 points and set the conditioning points and values. +# Here, we sample 60 points and set the conditioning points and values. ids = np.arange(srf.field.size) samples = np.random.RandomState(20210201).choice(ids, size=60, replace=False) @@ -37,11 +40,18 @@ cond_val = srf.field[samples] ############################################################################### +# Fitting and Interpolation +# ^^^^^^^^^^^^^^^^^^^^^^^^^ +# +# Now we want to interpolate sampled data from the given field (in order to +# pretend, we got any measured real-world data) +# and we want to normalize the given data with the BoxCox transformation. +# # Now we set up the kriging routine. We use a :any:`Stable` model, that should -# be fitted to the given data (by setting ``fit_variogram=True``) and we -# use the :any:`BoxCox` normalizer in order to gain normality. +# be fitted automatically to the given data (by setting ``fit_variogram=True``) +# and we use the :any:`BoxCox` normalizer in order to gain normality. # -# The BoxCox normalizer will be fitted automatically to the data, +# The normalizer will be fitted automatically to the data, # by setting ``fit_normalizer=True``. krige = gs.krige.Ordinary( @@ -68,13 +78,18 @@ # The BoxCox parameter `lmbda` was estimated to be almost 0, which means, # the log-normal distribution was correctly fitted. # -# Now let's run the kriging interpolation and let's have a look at the -# resulting field. We'll compare the original, sampled and interpolated field. - +# Now let's run the kriging interpolation and # interpolate krige(pos) -# plot +############################################################################### +# Plotting +# ^^^^^^^^ +# +# Now let's compare the original, sampled and interpolated fields. +# As we'll see, there is a lot of information in the covariance structure +# of the measurement samples and the field is reconstructed quite accurately. + fig, ax = plt.subplots(1, 3, figsize=[8, 3]) ax[0].imshow(srf.field.reshape(60, 60).T, origin="lower") ax[1].scatter(*cond_pos, c=cond_val) From 8b4e4c288c2782067fc43e6a583b0e3b386922bd Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Fri, 29 Jan 2021 17:48:29 +0100 Subject: [PATCH 45/97] Examples: update auto-fit normalizer example again --- examples/10_normalizer/01_auto_fit.py | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/examples/10_normalizer/01_auto_fit.py b/examples/10_normalizer/01_auto_fit.py index bf3db7744..a743064ae 100644 --- a/examples/10_normalizer/01_auto_fit.py +++ b/examples/10_normalizer/01_auto_fit.py @@ -14,8 +14,8 @@ ^^^^^^^^^^^^^^^ Here we generate log-normal data following a gaussian covariance model. -We will generate the "original" field on a 60x60 mesh, where we will take -samples from in order to pretend a data-scarcity situation. +We will generate the "original" field on a 60x60 mesh, from which we will take +samples in order to pretend a situation of data-scarcity. """ import numpy as np import gstools as gs @@ -43,16 +43,18 @@ # Fitting and Interpolation # ^^^^^^^^^^^^^^^^^^^^^^^^^ # -# Now we want to interpolate sampled data from the given field (in order to -# pretend, we got any measured real-world data) +# Now we want to interpolate the "meassured" samples # and we want to normalize the given data with the BoxCox transformation. # -# Now we set up the kriging routine. We use a :any:`Stable` model, that should -# be fitted automatically to the given data (by setting ``fit_variogram=True``) -# and we use the :any:`BoxCox` normalizer in order to gain normality. +# Here we set up the kriging routine and use a :any:`Stable` model, that should +# be fitted automatically to the given data +# and we pass the :any:`BoxCox` normalizer in order to gain normality. # # The normalizer will be fitted automatically to the data, # by setting ``fit_normalizer=True``. +# +# The covariance/variogram model will be fitted by an automatic workflow +# by setting ``fit_variogram=True``. krige = gs.krige.Ordinary( model=gs.Stable(dim=2), From 3ed157a853d1b2c8e038cedbc472b8478e68f21f Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Fri, 29 Jan 2021 17:49:45 +0100 Subject: [PATCH 46/97] Examples: update auto-fit normalizer example again --- examples/10_normalizer/01_auto_fit.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/10_normalizer/01_auto_fit.py b/examples/10_normalizer/01_auto_fit.py index a743064ae..0a809565d 100644 --- a/examples/10_normalizer/01_auto_fit.py +++ b/examples/10_normalizer/01_auto_fit.py @@ -80,15 +80,15 @@ # The BoxCox parameter `lmbda` was estimated to be almost 0, which means, # the log-normal distribution was correctly fitted. # -# Now let's run the kriging interpolation and -# interpolate +# Now let's run the kriging interpolation. + krige(pos) ############################################################################### # Plotting # ^^^^^^^^ # -# Now let's compare the original, sampled and interpolated fields. +# Finally let's compare the original, sampled and interpolated fields. # As we'll see, there is a lot of information in the covariance structure # of the measurement samples and the field is reconstructed quite accurately. From 4936a844f4bd88175ddfcc5ae36ce91dcb8171c8 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Fri, 29 Jan 2021 18:38:39 +0100 Subject: [PATCH 47/97] Examples: update auto-fit variogram example --- examples/03_variogram/06_auto_bin_latlon.py | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/examples/03_variogram/06_auto_bin_latlon.py b/examples/03_variogram/06_auto_bin_latlon.py index 4ca029259..603e090ae 100644 --- a/examples/03_variogram/06_auto_bin_latlon.py +++ b/examples/03_variogram/06_auto_bin_latlon.py @@ -72,5 +72,16 @@ ############################################################################### # Looks good, doesn't it? # +# This workflow is also implemented in the :any:`Krige` class, by setting +# ``fit_variogram=True``. Then the whole procedure shortens: + +krige = gs.krige.Ordinary(sph, pos, field, fit_variogram=True) +krige.structured((grid_lat, grid_lon)) + +# plot the result +krige.plot() +# show the fitting results +print(krige.model) + # This example shows, that setting up variogram estimation and kriging routines -# is straight forward with GSTools. ;-) +# is straight forward with GSTools! From ecbf9377f9b47724aa6b9bfd4180ac9b021b5dd6 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Fri, 29 Jan 2021 19:00:47 +0100 Subject: [PATCH 48/97] DOC: update --- README.md | 7 +++++-- examples/03_variogram/06_auto_bin_latlon.py | 1 + 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 94e253331..7c3277793 100644 --- a/README.md +++ b/README.md @@ -19,11 +19,14 @@ GeoStatTools provides geostatistical tools for various purposes: - random field generation +- simple, ordinary, universal and external drift kriging - conditioned field generation - incompressible random vector field generation -- simple and ordinary kriging -- variogram estimation and fitting +- (automatted) variogram estimation and fitting +- directional variogram estimation and modelling +- data normalization and transformation - many readily provided and even user-defined covariance models +- metric spatio-temporal modelling - plotting and exporting routines diff --git a/examples/03_variogram/06_auto_bin_latlon.py b/examples/03_variogram/06_auto_bin_latlon.py index 603e090ae..22ccc377b 100644 --- a/examples/03_variogram/06_auto_bin_latlon.py +++ b/examples/03_variogram/06_auto_bin_latlon.py @@ -83,5 +83,6 @@ # show the fitting results print(krige.model) +############################################################################### # This example shows, that setting up variogram estimation and kriging routines # is straight forward with GSTools! From e6b7b48f3aa6504efd84d5b9aed85b959fc80296 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Fri, 29 Jan 2021 19:07:27 +0100 Subject: [PATCH 49/97] Docs: update index --- docs/source/index.rst | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/docs/source/index.rst b/docs/source/index.rst index 585935efd..7df0543e7 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -6,9 +6,17 @@ GSTools Quickstart :width: 150px :align: center -GeoStatTools provides geostatistical tools for random field generation and -variogram estimation based on many readily provided and even user-defined -covariance models. +GeoStatTools provides geostatistical tools for various purposes: +- random field generation +- simple, ordinary, universal and external drift kriging +- conditioned field generation +- incompressible random vector field generation +- (automatted) variogram estimation and fitting +- directional variogram estimation and modelling +- data normalization and transformation +- many readily provided and even user-defined covariance models +- metric spatio-temporal modelling +- plotting and exporting routines Installation @@ -99,6 +107,7 @@ showing the most important use cases of GSTools, which are - `Field transformations `__ - `Geographic Coordinates `__ - `Spatio-Temporal Modelling `__ +- `Normalizing Data `__ - `Miscellaneous examples `__ Some more examples are provided in the examples folder. From 6b9b535b16a4095ab5142ace1b1dde1cbe5992b5 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Fri, 29 Jan 2021 19:34:32 +0100 Subject: [PATCH 50/97] README: link normalizer models --- README.md | 2 ++ docs/source/index.rst | 3 +-- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 7c3277793..3fb1dfb18 100644 --- a/README.md +++ b/README.md @@ -86,6 +86,7 @@ The documentation also includes some [tutorials][tut_link], showing the most imp - [Field transformations][tut7_link] - [Geographic Coordinates][tut8_link] - [Spatio-Temporal Modelling][tut9_link] +- [Normalizing Data][tut10_link] - [Miscellaneous examples][tut0_link] The associated python scripts are provided in the `examples` folder. @@ -366,6 +367,7 @@ You can contact us via . [tut7_link]: https://geostat-framework.readthedocs.io/projects/gstools/en/stable/examples/07_transformations/index.html [tut8_link]: https://geostat-framework.readthedocs.io/projects/gstools/en/stable/examples/08_geo_coordinates/index.html [tut9_link]: https://geostat-framework.readthedocs.io/projects/gstools/en/stable/examples/09_spatio_temporal/index.html +[tut9_link]: https://geostat-framework.readthedocs.io/projects/gstools/en/stable/examples/10_normalizer/index.html [tut0_link]: https://geostat-framework.readthedocs.io/projects/gstools/en/stable/examples/00_misc/index.html [cor_link]: https://en.wikipedia.org/wiki/Autocovariance#Normalization [vtk_link]: https://www.vtk.org/ diff --git a/docs/source/index.rst b/docs/source/index.rst index 7df0543e7..372f7a060 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -7,6 +7,7 @@ GSTools Quickstart :align: center GeoStatTools provides geostatistical tools for various purposes: + - random field generation - simple, ordinary, universal and external drift kriging - conditioned field generation @@ -110,8 +111,6 @@ showing the most important use cases of GSTools, which are - `Normalizing Data `__ - `Miscellaneous examples `__ -Some more examples are provided in the examples folder. - Spatial Random Field Generation =============================== From f49b4da3069ffc5f2845d7c81ee78724db5d1800 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sat, 30 Jan 2021 16:19:01 +0100 Subject: [PATCH 51/97] README: fix link --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 3fb1dfb18..2df5aef2e 100644 --- a/README.md +++ b/README.md @@ -367,7 +367,7 @@ You can contact us via . [tut7_link]: https://geostat-framework.readthedocs.io/projects/gstools/en/stable/examples/07_transformations/index.html [tut8_link]: https://geostat-framework.readthedocs.io/projects/gstools/en/stable/examples/08_geo_coordinates/index.html [tut9_link]: https://geostat-framework.readthedocs.io/projects/gstools/en/stable/examples/09_spatio_temporal/index.html -[tut9_link]: https://geostat-framework.readthedocs.io/projects/gstools/en/stable/examples/10_normalizer/index.html +[tut10_link]: https://geostat-framework.readthedocs.io/projects/gstools/en/stable/examples/10_normalizer/index.html [tut0_link]: https://geostat-framework.readthedocs.io/projects/gstools/en/stable/examples/00_misc/index.html [cor_link]: https://en.wikipedia.org/wiki/Autocovariance#Normalization [vtk_link]: https://www.vtk.org/ From 058f27b19130c8e7f25b4e8363ea9b4b16f9b6b3 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sat, 30 Jan 2021 16:20:08 +0100 Subject: [PATCH 52/97] Field: allow passing uninitialized Normalizer classes --- gstools/field/base.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/gstools/field/base.py b/gstools/field/base.py index 9c22d12cb..94623e087 100755 --- a/gstools/field/base.py +++ b/gstools/field/base.py @@ -329,6 +329,8 @@ def normalizer(self): @normalizer.setter def normalizer(self, normalizer): + if isinstance(normalizer, type) and issubclass(normalizer, Normalizer): + self._normalizer = normalizer() if isinstance(normalizer, Normalizer): self._normalizer = normalizer elif normalizer is None: From cbc305abcda399b276a8f9e20a63a93eadf6ea0d Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sat, 30 Jan 2021 16:20:47 +0100 Subject: [PATCH 53/97] Tools: add eval_func to doc; minor fix --- gstools/tools/misc.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/gstools/tools/misc.py b/gstools/tools/misc.py index d9f861e01..7fc549b19 100755 --- a/gstools/tools/misc.py +++ b/gstools/tools/misc.py @@ -8,11 +8,15 @@ .. autosummary:: list_format + eval_func """ import numpy as np from gstools.tools.geometric import format_struct_pos_dim, gen_mesh +__all__ = ["list_format", "eval_func"] + + def list_format(lst, prec): """Format a list of floats.""" return "[{}]".format( @@ -61,7 +65,7 @@ def eval_func( # care about scalar inputs if broadcast and not callable(func_val): return 0.0 if func_val is None else float(func_val) - elif not callable(func_val): + if not callable(func_val): func_val = _func_from_scalar(func_val) # care about mesh and function call if mesh_type != "unstructured": From b8572cd0ab491dfd9f1dbade34c29333c0571859 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sat, 30 Jan 2021 16:21:50 +0100 Subject: [PATCH 54/97] Krige: get_mean respects normalizer now; set_condition better handles existing external drift --- gstools/krige/base.py | 20 ++++++++++++++------ 1 file changed, 14 insertions(+), 6 deletions(-) diff --git a/gstools/krige/base.py b/gstools/krige/base.py index aefceb8b2..91137802c 100755 --- a/gstools/krige/base.py +++ b/gstools/krige/base.py @@ -396,7 +396,7 @@ def _get_dists(self, pos1, pos2=None, pos2_slice=(0, None)): return cdist(pos1.T, pos2.T[slice(*pos2_slice), ...]) def get_mean(self): - """Calculate the estimated mean. + """Calculate the estimated mean of the detrended field. Returns ------- @@ -405,7 +405,9 @@ def get_mean(self): Notes ----- - Only if the Kriging System has a constant mean. + Only not ``None`` if the Kriging System has a constant mean. + This means, no drift is given and the given field-mean is constant. + The result is neglecting a potential given trend. """ # if there are drift-terms, no constant mean can be calculated -> None if not self.has_const_mean: @@ -421,7 +423,7 @@ def get_mean(self): res = np.einsum( "i,ij,j", self._krige_cond, self._krige_mat, mean_est ) - return res + mean + return self.normalizer.denormalize(res + mean) def set_condition( self, @@ -471,12 +473,17 @@ def set_condition( >>> del Krige.cond_ext_drift """ - # Default values + # only use existing external drift, if no new positions are given + ext_drift = ( + self._cond_ext_drift + if (ext_drift is None and cond_pos is None) + else ext_drift + ) + # use existing values or set default cond_pos = self._cond_pos if cond_pos is None else cond_pos cond_val = self._cond_val if cond_val is None else cond_val - ext_drift = self._cond_ext_drift if ext_drift is None else ext_drift cond_err = self._cond_err if cond_err is None else cond_err - cond_err = "nugget" if cond_err is None else cond_err + cond_err = "nugget" if cond_err is None else cond_err # default if cond_pos is None or cond_val is None: raise ValueError("Krige.set_condition: missing cond_pos/cond_val.") # correctly format cond_pos and cond_val @@ -499,6 +506,7 @@ def set_condition( emp_vario = vario_estimate( self.cond_pos, field, direction=axes ) + # set the sill to the field variance self.model.fit_variogram(*emp_vario, sill=sill) # set the measurement errors self.cond_err = cond_err From 44aaf728daf9bdeb9e8cccbd0c5a16549676b249 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sat, 30 Jan 2021 17:13:48 +0100 Subject: [PATCH 55/97] Filed: fix normalizer setting with class --- gstools/field/base.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gstools/field/base.py b/gstools/field/base.py index 94623e087..5388fa7f4 100755 --- a/gstools/field/base.py +++ b/gstools/field/base.py @@ -331,7 +331,7 @@ def normalizer(self): def normalizer(self, normalizer): if isinstance(normalizer, type) and issubclass(normalizer, Normalizer): self._normalizer = normalizer() - if isinstance(normalizer, Normalizer): + elif isinstance(normalizer, Normalizer): self._normalizer = normalizer elif normalizer is None: self._normalizer = Normalizer() From bc3e893ec1cc0ae55bd7406bb2116ad88275468f Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sat, 30 Jan 2021 17:14:40 +0100 Subject: [PATCH 56/97] Krige: 'get_mean' add switch to enable processing; fix kriging the mean --- gstools/krige/base.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/gstools/krige/base.py b/gstools/krige/base.py index 91137802c..94616f9b7 100755 --- a/gstools/krige/base.py +++ b/gstools/krige/base.py @@ -215,8 +215,7 @@ def __call__( krige_var = np.empty(pnt_cnt, dtype=np.double) if return_var else None # set constant mean if present and wanted if only_mean and self.has_const_mean: - mean = 0.0 if self.mean is None else self.mean - field[...] = self.get_mean() - mean # mean is added later + field[...] = self.get_mean(post_process=False) # execute the kriging routine else: # set chunk size @@ -395,13 +394,16 @@ def _get_dists(self, pos1, pos2=None, pos2_slice=(0, None)): return cdist(pos1.T, pos1.T) return cdist(pos1.T, pos2.T[slice(*pos2_slice), ...]) - def get_mean(self): + def get_mean(self, post_process=True): """Calculate the estimated mean of the detrended field. Returns ------- mean : :class:`float` or :any:`None` Mean of the Kriging System. + post_process : :class:`bool`, optional + Whether to apply field-mean and normalizer. + Default: `True` Notes ----- @@ -423,7 +425,7 @@ def get_mean(self): res = np.einsum( "i,ij,j", self._krige_cond, self._krige_mat, mean_est ) - return self.normalizer.denormalize(res + mean) + return self.normalizer.denormalize(res + mean) if post_process else res def set_condition( self, From 83126c9a879dce41c33a6857934542dc25d3fb60 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sat, 30 Jan 2021 17:15:15 +0100 Subject: [PATCH 57/97] CondSRF: cleanup init; make krige a property --- gstools/field/cond_srf.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/gstools/field/cond_srf.py b/gstools/field/cond_srf.py index 74db5afa9..adb64b7bc 100644 --- a/gstools/field/cond_srf.py +++ b/gstools/field/cond_srf.py @@ -47,20 +47,15 @@ class CondSRF(Field): def __init__(self, krige, generator="RandMeth", **generator_kwargs): if not isinstance(krige, Krige): raise ValueError("CondSRF: krige should be an instance of Krige.") - self.krige = krige + self._krige = krige # initialize attributes self.pos = None self.mesh_type = None self.field = None self.raw_field = None # initialize private attributes - self._value_type = None self._generator = None - self._mean = None - self._normalizer = None - self._trend = None # initialize attributes - self.normalizer = None self.set_generator(generator, **generator_kwargs) def __call__(self, pos, seed=np.nan, mesh_type="unstructured", **kwargs): @@ -148,12 +143,17 @@ def set_generator(self, generator, **generator_kwargs): if generator in GENERATOR: gen = GENERATOR[generator] self._generator = gen(self.model, **generator_kwargs) - self.value_type = self._generator.value_type + self.value_type = self.generator.value_type else: raise ValueError("gstools.CondSRF: Unknown generator " + generator) if self.value_type != "scalar": raise ValueError("CondSRF: only scalar field value type allowed.") + @property + def krige(self): + """:any:`Krige`: The underlying kriging class.""" + return self._krige + @property def generator(self): """:any:`callable`: The generator of the field.""" From 0fa29dfe9666cba096aca7d386eae609cb99b7b8 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sat, 30 Jan 2021 17:17:54 +0100 Subject: [PATCH 58/97] Examples: minor updates in normalizer examples --- examples/10_normalizer/00_lognormal_kriging.py | 2 +- examples/10_normalizer/01_auto_fit.py | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/examples/10_normalizer/00_lognormal_kriging.py b/examples/10_normalizer/00_lognormal_kriging.py index 213a6912c..2ce268f00 100644 --- a/examples/10_normalizer/00_lognormal_kriging.py +++ b/examples/10_normalizer/00_lognormal_kriging.py @@ -28,7 +28,7 @@ ############################################################################### # In order to result in log-normal kriging, we will use the :any:`LogNormal` # Normalizer. This is a parameter-less normalizer, so we don't have to fit it. -normalizer = gs.normalizer.LogNormal() +normalizer = gs.normalizer.LogNormal ############################################################################### # Now we generate the interpolated field as well as the mean field. diff --git a/examples/10_normalizer/01_auto_fit.py b/examples/10_normalizer/01_auto_fit.py index 0a809565d..78138565d 100644 --- a/examples/10_normalizer/01_auto_fit.py +++ b/examples/10_normalizer/01_auto_fit.py @@ -8,7 +8,8 @@ Normalizers are fitted by minimizing the likelihood function and variograms are fitted by estimating the empirical variogram with automatic binning and -fitting the theoretical model to it. +fitting the theoretical model to it. Thereby the sill is constrained to match +the field variance. Artificial data ^^^^^^^^^^^^^^^ From f3fc351103fd13c7c0576d8339f6ef07bab66546 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sun, 31 Jan 2021 01:06:57 +0100 Subject: [PATCH 59/97] Krige: better handling of field-mean in kriging the mean --- gstools/krige/base.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/gstools/krige/base.py b/gstools/krige/base.py index 94616f9b7..d82643536 100755 --- a/gstools/krige/base.py +++ b/gstools/krige/base.py @@ -214,7 +214,7 @@ def __call__( field = np.empty(pnt_cnt, dtype=np.double) krige_var = np.empty(pnt_cnt, dtype=np.double) if return_var else None # set constant mean if present and wanted - if only_mean and self.has_const_mean: + if only_mean and self.drift_no == 0: field[...] = self.get_mean(post_process=False) # execute the kriging routine else: @@ -412,7 +412,8 @@ def get_mean(self, post_process=True): The result is neglecting a potential given trend. """ # if there are drift-terms, no constant mean can be calculated -> None - if not self.has_const_mean: + # if mean should not be post-processed, it exists when no drift given + if not self.has_const_mean and (post_process or self.drift_no > 0): return None res = 0.0 # for simple kriging return the given mean # correctly setting given mean From f75ee1c6bc973193143f042c867cb6521a11c721 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sun, 31 Jan 2021 18:49:16 +0100 Subject: [PATCH 60/97] Field: better formatting of mean, normalizer and trend in repr --- gstools/field/base.py | 33 ++++++++++----------------------- gstools/field/srf.py | 15 +++++---------- gstools/field/tools.py | 26 +++++++++++++++++++++++++- gstools/krige/base.py | 3 ++- 4 files changed, 42 insertions(+), 35 deletions(-) diff --git a/gstools/field/base.py b/gstools/field/base.py index 5388fa7f4..2f50a7587 100755 --- a/gstools/field/base.py +++ b/gstools/field/base.py @@ -16,7 +16,7 @@ from gstools.tools.geometric import format_struct_pos_dim, gen_mesh from gstools.tools.misc import eval_func from gstools.normalizer import Normalizer -from gstools.field.tools import mesh_call, to_vtk_helper +from gstools.field.tools import mesh_call, to_vtk_helper, fmt_mean_norm_trend __all__ = ["Field"] @@ -207,7 +207,7 @@ def post_field(self, field, name="field", process=True, save=True): # apply mean - normalizer - trend field += eval_func(func_val=self.mean, **kwargs) field = self.normalizer.denormalize(field) - field += eval_func(self.trend, **kwargs) + field += eval_func(func_val=self.trend, **kwargs) if save: setattr(self, name, field) return field @@ -358,32 +358,19 @@ def value_type(self, value_type): raise ValueError("Field: value type not in {}".format(VALUE_TYPES)) self._value_type = value_type - def _fmt_func_val(self, func_val): - if func_val is None: - return str(None) - if callable(func_val): - return "" - return "{0:.{p}}".format(float(func_val), p=self.model._prec) - - def _fmt_normalizer(self): - norm = self.normalizer - return str(None) if norm.__class__ is Normalizer else norm.name - @property def name(self): """:class:`str`: The name of the class.""" return self.__class__.__name__ + def _fmt_mean_norm_trend(self): + return fmt_mean_norm_trend(self) + def __repr__(self): """Return String representation.""" - return ( - "{0}(model={1}, value_type='{2}', " - "mean={3}, normalizer={4}, trend={5})".format( - self.name, - self.model.name, - self.value_type, - self._fmt_func_val(self.mean), - self._fmt_normalizer(), - self._fmt_func_val(self.trend), - ) + return "{0}(model={1}, value_type='{2}'{3})".format( + self.name, + self.model.name, + self.value_type, + self._fmt_mean_norm_trend(), ) diff --git a/gstools/field/srf.py b/gstools/field/srf.py index 072ae7c02..d92a46c82 100644 --- a/gstools/field/srf.py +++ b/gstools/field/srf.py @@ -193,14 +193,9 @@ def upscaling(self, upscaling): def __repr__(self): """Return String representation.""" - return ( - "{0}(model={1}, " - "mean={2}, normalizer={3}, trend={4}, generator={5})".format( - self.name, - self.model.name, - self._fmt_func_val(self.mean), - self._fmt_normalizer(), - self._fmt_func_val(self.trend), - self.generator.name, - ) + return "{0}(model={1}{2}, generator={3})".format( + self.name, + self.model.name, + self._fmt_mean_norm_trend(), + self.generator.name, ) diff --git a/gstools/field/tools.py b/gstools/field/tools.py index d0ba4ff3e..ef57285db 100644 --- a/gstools/field/tools.py +++ b/gstools/field/tools.py @@ -13,10 +13,34 @@ import numpy as np import meshio +from gstools.normalizer import Normalizer from gstools.tools.export import to_vtk, vtk_export -__all__ = ["to_vtk_helper", "mesh_call"] +__all__ = ["fmt_mean_norm_trend", "to_vtk_helper", "mesh_call"] + + +def _fmt_func_val(f_cls, func_val): + if func_val is None: + return str(None) + if callable(func_val): + return "" # or format(func_val.__name__) + return "{0:.{p}}".format(float(func_val), p=f_cls.model._prec) + + +def _fmt_normalizer(f_cls): + norm = f_cls.normalizer + return str(None) if norm.__class__ is Normalizer else norm.name + + +def fmt_mean_norm_trend(f_cls): + """Format string repr. for mean, normalizer and trend of a field.""" + args = [ + "mean=" + _fmt_func_val(f_cls, f_cls.mean), + "normalizer=" + _fmt_normalizer(f_cls), + "trend=" + _fmt_func_val(f_cls, f_cls.trend), + ] + return "".join([", " + arg for arg in args if not arg.endswith("None")]) def to_vtk_helper( diff --git a/gstools/krige/base.py b/gstools/krige/base.py index d82643536..78df94bc4 100755 --- a/gstools/krige/base.py +++ b/gstools/krige/base.py @@ -697,8 +697,9 @@ def name(self): def __repr__(self): """Return String representation.""" - return "{0}(model={1}, cond_no={2})".format( + return "{0}(model={1}, cond_no={2}{3})".format( self.name, self.model.name, self.cond_no, + self._fmt_mean_norm_trend(), ) From 46db911faac96330953adce7e66b96dd7193dc73 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 1 Feb 2021 21:45:14 +0100 Subject: [PATCH 61/97] Normalizer: ignore nan-data; add range check for data; better subclassing; use 'data' as arg name --- gstools/normalizer/base.py | 106 +++++++++++------ gstools/normalizer/methods.py | 209 ++++++++++++++++++---------------- 2 files changed, 184 insertions(+), 131 deletions(-) diff --git a/gstools/normalizer/base.py b/gstools/normalizer/base.py index 12814162c..d4ce96613 100644 --- a/gstools/normalizer/base.py +++ b/gstools/normalizer/base.py @@ -24,13 +24,20 @@ class Normalizer: Input data to fit the transformation to in order to gain normality. The default is None. **parameter - Specified parameters given by name. If not given, default values - will be applied. + Specified parameters given by name. If not given, default parameters + will be used. """ + default_parameter = {} + """:class:`dict`: Default parameters of the Normalizer.""" + normalize_range = (-np.inf, np.inf) + """:class:`tuple`: Valid range for input data.""" + denormalize_range = (-np.inf, np.inf) + """:class:`tuple`: Valid range for output/normal data.""" + def __init__(self, data=None, **parameter): # only use parameter, that have a provided default value - for key, value in self.default_parameter().items(): + for key, value in self.default_parameter.items(): setattr(self, key, parameter.get(key, value)) # fit parameters if data is given if data is not None: @@ -40,60 +47,93 @@ def __init__(self, data=None, **parameter): # precision for printing self._prec = 3 - def default_parameter(self): - """Get default parameters for the transformation. + def _denormalize(self, data): + return data - Returns - ------- - :class:`dict` - Default parameters. - """ - return {} + def _normalize(self, data): + return data + + def _derivative(self, data): + return spm.derivative(self._normalize, data, dx=1e-6) - def denormalize(self, values): + def _loglikelihood(self, data): + add = -0.5 * np.size(data) * (np.log(2 * np.pi) + 1) + return self._kernel_loglikelihood(data) + add + + def _kernel_loglikelihood(self, data): + res = -0.5 * np.size(data) * np.log(np.var(self._normalize(data))) + return res + np.sum(np.log(np.maximum(1e-16, self._derivative(data)))) + + def _check_input(self, data, data_range=None, return_output_template=True): + is_data = np.array(np.logical_not(np.isnan(data))) + if return_output_template: + out = np.full_like(data, np.nan, dtype=np.double) + data = np.array(data, dtype=np.double)[is_data] + if data_range is not None and np.min(np.abs(data_range)) < np.inf: + dat_in = np.logical_and(data > data_range[0], data < data_range[1]) + if not np.all(dat_in): + warnings.warn( + "{0}: data (min: {1}, max: {2}) out of range: {3}. " + "Affected values will be treated as NaN.".format( + self.name, np.min(data), np.max(data), data_range + ) + ) + is_data[is_data] &= dat_in + data = data[dat_in] + if return_output_template: + return data, is_data, out + return data + + def denormalize(self, data): """Transform to input distribution. Parameters ---------- - values : array_like - Input values (normal distributed). + data : array_like + Input data (normal distributed). Returns ------- :class:`numpy.ndarray` - Denormalized values. + Denormalized data. """ - return values + data, is_data, out = self._check_input(data, self.denormalize_range) + out[is_data] = self._denormalize(data) + return out - def normalize(self, values): + def normalize(self, data): """Transform to normal distribution. Parameters ---------- - values : array_like - Input values (not normal distributed). + data : array_like + Input data (not normal distributed). Returns ------- :class:`numpy.ndarray` - Normalized values. + Normalized data. """ - return values + data, is_data, out = self._check_input(data, self.normalize_range) + out[is_data] = self._normalize(data) + return out - def derivative(self, values): + def derivative(self, data): """Factor for normal PDF to gain target PDF. Parameters ---------- - values : array_like - Input values. + data : array_like + Input data (not normal distributed). Returns ------- :class:`numpy.ndarray` Derivative of the normalization transformation function. """ - return spm.derivative(self.normalize, np.asanyarray(values), dx=1e-6) + data, is_data, out = self._check_input(data, self.normalize_range) + out[is_data] = self._derivative(data) + return out def likelihood(self, data): """Likelihood for given data with current parameters. @@ -123,8 +163,8 @@ def loglikelihood(self, data): :class:`float` Log-Likelihood of the given data. """ - add = -0.5 * np.size(data) * (np.log(2 * np.pi) + 1) - return self.kernel_loglikelihood(data) + add + data = self._check_input(data, self.normalize_range, False) + return self._loglikelihood(data) def kernel_loglikelihood(self, data): """Kernel Log-Likelihood for given data with current parameters. @@ -144,8 +184,8 @@ def kernel_loglikelihood(self, data): This loglikelihood function is neglecting additive constants, that are not needed for optimization. """ - res = -0.5 * np.size(data) * np.log(np.var(self.normalize(data))) - return res + np.sum(np.log(np.maximum(1e-16, self.derivative(data)))) + data = self._check_input(data, self.normalize_range, False) + return self._kernel_loglikelihood(data) def fit(self, data, skip=None, **kwargs): """Fitting the transformation to data by maximizing Log-Likelihood. @@ -167,7 +207,7 @@ def fit(self, data, skip=None, **kwargs): Optimal paramters given by names. """ skip = [] if skip is None else skip - all_names = sorted(self.default_parameter()) + all_names = sorted(self.default_parameter) para_names = [name for name in all_names if name not in skip] def _neg_kllf(par, dat): @@ -176,14 +216,14 @@ def _neg_kllf(par, dat): return -self.kernel_loglikelihood(dat) if len(para_names) == 0: # transformations without para. (no opti.) - warnings.warn(self.__class__.__name__ + ".fit: no parameters!") + warnings.warn(self.name + ".fit: no parameters!") return {} if len(para_names) == 1: # one-para. transformations (simple opti.) # default bracket like in scipy's boxcox (if not given) kwargs.setdefault("bracket", (-2, 2)) out = spo.minimize_scalar(_neg_kllf, args=(data,), **kwargs) else: # general case - # init guess from current values (if x0 not given) + # init guess from current parameters (if x0 not given) kwargs.setdefault("x0", [getattr(self, p) for p in para_names]) out = spo.minimize(_neg_kllf, args=(data,), **kwargs) # save optimization results @@ -201,6 +241,6 @@ def __repr__(self): """Return String representation.""" para_strs = [ "{0}={1:.{2}}".format(p, float(getattr(self, p)), self._prec) - for p in sorted(self.default_parameter()) + for p in sorted(self.default_parameter) ] return self.name + "(" + ", ".join(para_strs) + ")" diff --git a/gstools/normalizer/methods.py b/gstools/normalizer/methods.py index fab09cd1d..6a9d9a12e 100644 --- a/gstools/normalizer/methods.py +++ b/gstools/normalizer/methods.py @@ -29,17 +29,17 @@ class LogNormal(Normalizer): y=\log(x) """ - def denormalize(self, values): - """Transform to log-normal distribution.""" - return np.exp(values) + normalize_range = (0.0, np.inf) + """Valid range for input data.""" - def normalize(self, values): - """Transform to normal distribution.""" - return np.log(values) + def _denormalize(self, data): + return np.exp(data) - def derivative(self, values): - """Factor for normal PDF to gain target PDF.""" - return np.power(values, -1) + def _normalize(self, data): + return np.log(data) + + def _derivative(self, data): + return np.power(data, -1) class BoxCox(Normalizer): @@ -69,25 +69,32 @@ class BoxCox(Normalizer): of the Royal Statistical Society B, 26, 211-252 (1964). """ - def default_parameter(self): - """Get default parameters.""" - return {"lmbda": 1} + default_parameter = {"lmbda": 1} + """:class:`dict`: Default parameter of the BoxCox-Normalizer.""" + normalize_range = (0.0, np.inf) + """:class:`tuple`: Valid range for input data.""" + + @property + def denormalize_range(self): + """:class:`tuple`: Valid range for output data depending on lmbda. + + (-1/lmbda, inf) + """ + with np.errstate(divide="ignore"): + return (-np.divide(1, self.lmbda), np.inf) - def denormalize(self, values): - """Transform to target distribution.""" + def _denormalize(self, data): if np.isclose(self.lmbda, 0): - return np.exp(values) - return (1 + np.multiply(values, self.lmbda)) ** (1 / self.lmbda) + return np.exp(data) + return (1 + np.multiply(data, self.lmbda)) ** (1 / self.lmbda) - def normalize(self, values): - """Transform to normal distribution.""" + def _normalize(self, data): if np.isclose(self.lmbda, 0): - return np.log(values) - return (np.power(values, self.lmbda) - 1) / self.lmbda + return np.log(data) + return (np.power(data, self.lmbda) - 1) / self.lmbda - def derivative(self, values): - """Factor for normal PDF to gain target PDF.""" - return np.power(values, self.lmbda - 1) + def _derivative(self, data): + return np.power(data, self.lmbda - 1) class BoxCoxShift(Normalizer): @@ -127,27 +134,40 @@ class BoxCoxShift(Normalizer): of the Royal Statistical Society B, 26, 211-252 (1964). """ - def default_parameter(self): - """Get default parameters.""" - return {"shift": 0, "lmbda": 1} + default_parameter = {"shift": 0, "lmbda": 1} + """:class:`dict`: Default parameters of the BoxCoxShift-Normalizer.""" + + @property + def normalize_range(self): + """:class:`tuple`: Valid range for input data depending on shift. + + (-shift, inf) + """ + return (-self.shift, np.inf) - def denormalize(self, values): - """Transform to target distribution.""" + @property + def denormalize_range(self): + """:class:`tuple`: Valid range for output data depending on lmbda. + + (-1/lmbda, inf) + """ + with np.errstate(divide="ignore"): + return (-np.divide(1, self.lmbda), np.inf) + + def _denormalize(self, data): if np.isclose(self.lmbda, 0): - return np.exp(values) - self.shift - return (1 + np.multiply(values, self.lmbda)) ** ( + return np.exp(data) - self.shift + return (1 + np.multiply(data, self.lmbda)) ** ( 1 / self.lmbda ) - self.shift - def normalize(self, values): - """Transform to normal distribution.""" + def _normalize(self, data): if np.isclose(self.lmbda, 0): - return np.log(np.add(values, self.shift)) - return (np.add(values, self.shift) ** self.lmbda - 1) / self.lmbda + return np.log(np.add(data, self.shift)) + return (np.add(data, self.shift) ** self.lmbda - 1) / self.lmbda - def derivative(self, values): - """Factor for normal PDF to gain target PDF.""" - return np.power(np.add(values, self.shift), self.lmbda - 1) + def _derivative(self, data): + return np.power(np.add(data, self.shift), self.lmbda - 1) class YeoJohnson(Normalizer): @@ -185,53 +205,47 @@ class YeoJohnson(Normalizer): (2000). """ - def default_parameter(self): - """Get default parameters.""" - return {"lmbda": 1} + default_parameter = {"lmbda": 1} + """:class:`dict`: Default parameter of the YeoJohnson-Normalizer.""" - def denormalize(self, values): - """Transform to target distribution.""" - values = np.asanyarray(values) - res = np.zeros_like(values, dtype=np.double) - pos = values >= 0 - # when values >= 0 + def _denormalize(self, data): + data = np.asanyarray(data) + res = np.zeros_like(data, dtype=np.double) + pos = data >= 0 + # when data >= 0 if np.isclose(self.lmbda, 0): - res[pos] = np.expm1(values[pos]) + res[pos] = np.expm1(data[pos]) else: # self.lmbda != 0 - res[pos] = ( - np.power(values[pos] * self.lmbda + 1, 1 / self.lmbda) - 1 - ) - # when values < 0 + res[pos] = np.power(data[pos] * self.lmbda + 1, 1 / self.lmbda) - 1 + # when data < 0 if np.isclose(self.lmbda, 2): - res[~pos] = -np.expm1(-values[~pos]) + res[~pos] = -np.expm1(-data[~pos]) else: # self.lmbda != 2 res[~pos] = 1 - np.power( - -(2 - self.lmbda) * values[~pos] + 1, 1 / (2 - self.lmbda) + -(2 - self.lmbda) * data[~pos] + 1, 1 / (2 - self.lmbda) ) return res - def normalize(self, values): - """Transform to normal distribution.""" - values = np.asanyarray(values) - res = np.zeros_like(values, dtype=np.double) - pos = values >= 0 - # when values >= 0 + def _normalize(self, data): + data = np.asanyarray(data) + res = np.zeros_like(data, dtype=np.double) + pos = data >= 0 + # when data >= 0 if np.isclose(self.lmbda, 0): - res[pos] = np.log1p(values[pos]) + res[pos] = np.log1p(data[pos]) else: # self.lmbda != 0 - res[pos] = (np.power(values[pos] + 1, self.lmbda) - 1) / self.lmbda - # when values < 0 + res[pos] = (np.power(data[pos] + 1, self.lmbda) - 1) / self.lmbda + # when data < 0 if np.isclose(self.lmbda, 2): - res[~pos] = -np.log1p(-values[~pos]) + res[~pos] = -np.log1p(-data[~pos]) else: # self.lmbda != 2 - res[~pos] = -(np.power(-values[~pos] + 1, 2 - self.lmbda) - 1) / ( + res[~pos] = -(np.power(-data[~pos] + 1, 2 - self.lmbda) - 1) / ( 2 - self.lmbda ) return res - def derivative(self, values): - """Factor for normal PDF to gain target PDF.""" - return (np.abs(values) + 1) ** (np.sign(values) * (self.lmbda - 1)) + def _derivative(self, data): + return (np.abs(data) + 1) ** (np.sign(data) * (self.lmbda - 1)) class Modulus(Normalizer): @@ -262,31 +276,25 @@ class Modulus(Normalizer): of the Royal Statistical Society C, 29.2, 190-197, (1980) """ - def default_parameter(self): - """Get default parameters.""" - return {"lmbda": 1} + default_parameter = {"lmbda": 1} + """:class:`dict`: Default parameter of the Modulus-Normalizer.""" - def denormalize(self, values): - """Transform to target distribution.""" + def _denormalize(self, data): if np.isclose(self.lmbda, 0): - return np.sign(values) * np.expm1(np.abs(values)) - return np.sign(values) * ( - (1 + self.lmbda * np.abs(values)) ** (1 / self.lmbda) - 1 + return np.sign(data) * np.expm1(np.abs(data)) + return np.sign(data) * ( + (1 + self.lmbda * np.abs(data)) ** (1 / self.lmbda) - 1 ) - def normalize(self, values): - """Transform to normal distribution.""" + def _normalize(self, data): if np.isclose(self.lmbda, 0): - return np.sign(values) * np.log1p(np.abs(values)) + return np.sign(data) * np.log1p(np.abs(data)) return ( - np.sign(values) - * ((np.abs(values) + 1) ** self.lmbda - 1) - / self.lmbda + np.sign(data) * ((np.abs(data) + 1) ** self.lmbda - 1) / self.lmbda ) - def derivative(self, values): - """Factor for normal PDF to gain target PDF.""" - return np.power(np.abs(values) + 1, self.lmbda - 1) + def _derivative(self, data): + return np.power(np.abs(data) + 1, self.lmbda - 1) class Manly(Normalizer): @@ -316,22 +324,27 @@ class Manly(Normalizer): Journal of the Royal Statistical Society D, 25.1, 37-42 (1976). """ - def default_parameter(self): - """Get default parameters.""" - return {"lmbda": 1} + default_parameter = {"lmbda": 1} + """:class:`dict`: Default parameter of the Manly-Normalizer.""" + + @property + def denormalize_range(self): + """:class:`tuple`: Valid range for output data depending on lmbda. + + (-1/lmbda, inf) + """ + with np.errstate(divide="ignore"): + return (-np.divide(1, self.lmbda), np.inf) - def denormalize(self, values): - """Transform to target distribution.""" + def _denormalize(self, data): if np.isclose(self.lmbda, 0): - return values - return np.log1p(np.multiply(values, self.lmbda)) / self.lmbda + return data + return np.log1p(np.multiply(data, self.lmbda)) / self.lmbda - def normalize(self, values): - """Transform to normal distribution.""" + def _normalize(self, data): if np.isclose(self.lmbda, 0): - return values - return np.expm1(np.multiply(values, self.lmbda)) / self.lmbda + return data + return np.expm1(np.multiply(data, self.lmbda)) / self.lmbda - def derivative(self, values): - """Factor for normal PDF to gain target PDF.""" - return np.exp(np.multiply(values, self.lmbda)) + def _derivative(self, data): + return np.exp(np.multiply(data, self.lmbda)) From 9092ef9834e74773d43f2b4a8252e4d65453849c Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 1 Feb 2021 22:03:08 +0100 Subject: [PATCH 62/97] Docs: cleanup doc; Sphinx bugfixes (wrong refs) --- docs/source/field.base.rst | 10 ---------- docs/source/field.rst | 5 ----- docs/source/variogram.rst | 2 -- gstools/field/__init__.py | 5 ++++- gstools/field/cond_srf.py | 6 +++--- gstools/field/generator.py | 4 ++-- gstools/krige/base.py | 4 ++-- gstools/krige/methods.py | 10 +++++----- gstools/variogram/__init__.py | 4 ++++ 9 files changed, 20 insertions(+), 30 deletions(-) delete mode 100755 docs/source/field.base.rst diff --git a/docs/source/field.base.rst b/docs/source/field.base.rst deleted file mode 100755 index 5e89c99ec..000000000 --- a/docs/source/field.base.rst +++ /dev/null @@ -1,10 +0,0 @@ -gstools.field.base ------------------- - -.. automodule:: gstools.field.base - :members: - :undoc-members: - -.. raw:: latex - - \clearpage diff --git a/docs/source/field.rst b/docs/source/field.rst index 24ef6357f..c37b49660 100644 --- a/docs/source/field.rst +++ b/docs/source/field.rst @@ -2,10 +2,6 @@ gstools.field ============= .. automodule:: gstools.field - :members: - :undoc-members: - :inherited-members: - :show-inheritance: .. raw:: latex @@ -16,4 +12,3 @@ gstools.field field.generator.rst field.upscaling.rst - field.base.rst diff --git a/docs/source/variogram.rst b/docs/source/variogram.rst index cb6bd288d..2f50d2669 100644 --- a/docs/source/variogram.rst +++ b/docs/source/variogram.rst @@ -2,8 +2,6 @@ gstools.variogram ================= .. automodule:: gstools.variogram - :members: - :undoc-members: .. raw:: latex diff --git a/gstools/field/__init__.py b/gstools/field/__init__.py index d92cdbc4b..01a7fc681 100644 --- a/gstools/field/__init__.py +++ b/gstools/field/__init__.py @@ -10,12 +10,13 @@ .. autosummary:: generator upscaling - base Spatial Random Field ^^^^^^^^^^^^^^^^^^^^ .. autosummary:: + :toctree: generated + SRF CondSRF @@ -23,6 +24,8 @@ ^^^^^^^^^^^^^^^^ .. autosummary:: + :toctree: generated + Field ---- diff --git a/gstools/field/cond_srf.py b/gstools/field/cond_srf.py index adb64b7bc..33b5b73a6 100644 --- a/gstools/field/cond_srf.py +++ b/gstools/field/cond_srf.py @@ -29,8 +29,8 @@ class CondSRF(Field): Parameters ---------- - model : :any:`CovModel` - Covariance Model of the spatial random field. + krige : :any:`Krige` + Kriging setup to condition the spatial random field. generator : :class:`str`, optional Name of the field generator to be used. At the moment, only the following generator is provided: @@ -116,7 +116,7 @@ def get_scaling(self, krige_var, shape): ------- var_scale : :class:`numpy.ndarray` Variance scaling factor for the random field. - nugget : :class:`numpy.ndarray` or `class:`int + nugget : :class:`numpy.ndarray` or :class:`int` Nugget to be added to the field. """ if self.model.nugget > 0: diff --git a/gstools/field/generator.py b/gstools/field/generator.py index 8b5a1f36e..5bcdd8d36 100644 --- a/gstools/field/generator.py +++ b/gstools/field/generator.py @@ -41,7 +41,7 @@ class RandMeth: verbose : :class:`bool`, optional Be chatty during the generation. Default: :any:`False` - sampling :class:`str`, optional + sampling : :class:`str`, optional Sampling strategy. Either * "auto": select best strategy depending on given model @@ -331,7 +331,7 @@ class IncomprRandMeth(RandMeth): verbose : :class:`bool`, optional State if there should be output during the generation. Default: :any:`False` - sampling :class:`str`, optional + sampling : :class:`str`, optional Sampling strategy. Either * "auto": select best strategy depending on given model diff --git a/gstools/krige/base.py b/gstools/krige/base.py index 78df94bc4..4b60c24dc 100755 --- a/gstools/krige/base.py +++ b/gstools/krige/base.py @@ -105,7 +105,7 @@ class Krige(Field): Wheater to fit the given variogram model to the data. This is done by using isotropy settings of the given model, assuming the sill to be the data variance and with the - standard bins provided by the :any:`variogram` submodule. + standard bins provided by the :any:`standard_bins` routine. Default: False Notes @@ -467,7 +467,7 @@ def set_condition( Wheater to fit the given variogram model to the data. This is done by using isotropy settings of the given model, assuming the sill to be the data variance and with the - standard bins provided by the :any:`variogram` submodule. + standard bins provided by the :any:`standard_bins` routine. Default: False Notes diff --git a/gstools/krige/methods.py b/gstools/krige/methods.py index ee3534a9d..48aea0a20 100644 --- a/gstools/krige/methods.py +++ b/gstools/krige/methods.py @@ -81,7 +81,7 @@ class Simple(Krige): Wheater to fit the given variogram model to the data. This is done by using isotropy settings of the given model, assuming the sill to be the data variance and with the - standard bins provided by the :any:`variogram` submodule. + standard bins provided by the :any:`standard_bins` routine. Default: False """ @@ -176,7 +176,7 @@ class Ordinary(Krige): Wheater to fit the given variogram model to the data. This is done by using isotropy settings of the given model, assuming the sill to be the data variance and with the - standard bins provided by the :any:`variogram` submodule. + standard bins provided by the :any:`standard_bins` routine. Default: False """ @@ -281,7 +281,7 @@ class Universal(Krige): Wheater to fit the given variogram model to the data. This is done by using isotropy settings of the given model, assuming the sill to be the data variance and with the - standard bins provided by the :any:`variogram` submodule. + standard bins provided by the :any:`standard_bins` routine. Default: False """ @@ -383,7 +383,7 @@ class ExtDrift(Krige): Wheater to fit the given variogram model to the data. This is done by using isotropy settings of the given model, assuming the sill to be the data variance and with the - standard bins provided by the :any:`variogram` submodule. + standard bins provided by the :any:`standard_bins` routine. Default: False """ @@ -475,7 +475,7 @@ class Detrended(Krige): Wheater to fit the given variogram model to the data. This is done by using isotropy settings of the given model, assuming the sill to be the data variance and with the - standard bins provided by the :any:`variogram` submodule. + standard bins provided by the :any:`standard_bins` routine. Default: False """ diff --git a/gstools/variogram/__init__.py b/gstools/variogram/__init__.py index 22a362b1a..fe2c148fc 100644 --- a/gstools/variogram/__init__.py +++ b/gstools/variogram/__init__.py @@ -8,6 +8,8 @@ ^^^^^^^^^^^^^^^^^^^^ .. autosummary:: + :toctree: generated + vario_estimate vario_estimate_axis @@ -15,6 +17,8 @@ ^^^^^^^ .. autosummary:: + :toctree: generated + standard_bins ---- From 27a691b0081975df7e8ee4fb3c3aa539c652b3fd Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Tue, 2 Feb 2021 21:54:53 +0100 Subject: [PATCH 63/97] Normalizer: add routines to add/remove mean-normalizer-trend --- gstools/normalizer/__init__.py | 15 +++ gstools/normalizer/tools.py | 182 +++++++++++++++++++++++++++++++++ 2 files changed, 197 insertions(+) create mode 100644 gstools/normalizer/tools.py diff --git a/gstools/normalizer/__init__.py b/gstools/normalizer/__init__.py index d25487b2b..af3a5f243 100644 --- a/gstools/normalizer/__init__.py +++ b/gstools/normalizer/__init__.py @@ -24,6 +24,15 @@ YeoJohnson Modulus Manly + +Convenience Routines +^^^^^^^^^^^^^^^^^^^^ + +.. autosummary:: + :toctree: generated + + apply_mean_norm_trend + remove_trend_norm_mean """ from gstools.normalizer.base import Normalizer @@ -35,6 +44,10 @@ Modulus, Manly, ) +from gstools.normalizer.tools import ( + apply_mean_norm_trend, + remove_trend_norm_mean, +) __all__ = [ "Normalizer", @@ -44,4 +57,6 @@ "YeoJohnson", "Modulus", "Manly", + "apply_mean_norm_trend", + "remove_trend_norm_mean", ] diff --git a/gstools/normalizer/tools.py b/gstools/normalizer/tools.py new file mode 100644 index 000000000..094b5a748 --- /dev/null +++ b/gstools/normalizer/tools.py @@ -0,0 +1,182 @@ +# -*- coding: utf-8 -*- +""" +GStools subpackage providing tools for Normalizers. + +.. currentmodule:: gstools.normalizer.tools + +The following classes and functions are provided + +.. autosummary:: + apply_mean_norm_trend + remove_trend_norm_mean +""" +import numpy as np + +from gstools.normalizer import Normalizer +from gstools.tools.misc import eval_func +from gstools.tools.geometric import ( + format_struct_pos_shape, + format_unstruct_pos_shape, +) + +__all__ = ["apply_mean_norm_trend", "remove_trend_norm_mean"] + + +def _check_normalizer(normalizer): + if isinstance(normalizer, type) and issubclass(normalizer, Normalizer): + normalizer = normalizer() + elif normalizer is None: + normalizer = Normalizer() + elif not isinstance(normalizer, Normalizer): + raise ValueError("Check: 'normalizer' not of type 'Normalizer'.") + return normalizer + + +def apply_mean_norm_trend( + pos, + field, + mean=None, + normalizer=None, + trend=None, + mesh_type="unstructured", + value_type="scalar", + check_shape=True, + stacked=False, +): + """ + Apply mean, de-normalization and trend to given field. + + Parameters + ---------- + pos : :any:`iterable` + Position tuple, containing main direction and transversal directions. + field : :class:`numpy.ndarray` or :class:`list` of :class:`numpy.ndarray` + The spatially distributed data. + You can pass a list of fields, that will be used simultaneously. + Then you need to set ``stacked=True``. + mean : :any:`None` or :class:`float` or :any:`callable`, optional + Mean of the field if wanted. Could also be a callable. + The default is None. + normalizer : :any:`None` or :any:`Normalizer`, optional + Normalizer to be applied to the field. + The default is None. + trend : :any:`None` or :class:`float` or :any:`callable`, optional + Trend of the denormalized fields. If no normalizer is applied, + this behaves equal to 'mean'. + The default is None. + mesh_type : :class:`str`, optional + 'structured' / 'unstructured' + Default: 'unstructured' + value_type : :class:`str`, optional + Value type of the field. Either "scalar" or "vector". + The default is "scalar". + check_shape : :class:`bool`, optional + Wheather to check pos and field shapes. The default is True. + stacked : :class:`bool`, optional + Wheather the field is stacked or not. The default is False. + + Returns + ------- + :class:`numpy.ndarray` + The transformed field. + """ + normalizer = _check_normalizer(normalizer) + if check_shape: + if mesh_type != "unstructured": + pos, shape, dim = format_struct_pos_shape( + pos, field.shape, check_stacked_shape=stacked + ) + else: + pos, shape, dim = format_unstruct_pos_shape( + pos, field.shape, check_stacked_shape=stacked + ) + field = np.array(field, dtype=np.double).reshape(shape) + else: + dim = len(pos) + if not stacked: + field = [field] + field_cnt = len(field) + for i in range(field_cnt): + field[i] += eval_func(mean, pos, dim, mesh_type, value_type, True) + field = normalizer.denormalize(field) + for i in range(field_cnt): + field[i] += eval_func(trend, pos, dim, mesh_type, value_type, True) + return field if stacked else field[0] + + +def remove_trend_norm_mean( + pos, + field, + mean=None, + normalizer=None, + trend=None, + mesh_type="unstructured", + value_type="scalar", + check_shape=True, + stacked=False, + fit_normalizer=False, +): + """ + Remove trend, de-normalization and mean from given field. + + Parameters + ---------- + pos : :any:`iterable` + Position tuple, containing main direction and transversal directions. + field : :class:`numpy.ndarray` or :class:`list` of :class:`numpy.ndarray` + The spatially distributed data. + You can pass a list of fields, that will be used simultaneously. + Then you need to set ``stacked=True``. + mean : :any:`None` or :class:`float` or :any:`callable`, optional + Mean of the field if wanted. Could also be a callable. + The default is None. + normalizer : :any:`None` or :any:`Normalizer`, optional + Normalizer to be applied to the field. + The default is None. + trend : :any:`None` or :class:`float` or :any:`callable`, optional + Trend of the denormalized fields. If no normalizer is applied, + this behaves equal to 'mean'. + The default is None. + mesh_type : :class:`str`, optional + 'structured' / 'unstructured' + Default: 'unstructured' + value_type : :class:`str`, optional + Value type of the field. Either "scalar" or "vector". + The default is "scalar". + check_shape : :class:`bool`, optional + Wheather to check pos and field shapes. The default is True. + stacked : :class:`bool`, optional + Wheather the field is stacked or not. The default is False. + fit_normalizer : :class:`bool`, optional + Wheater to fit the data-normalizer to the given (detrended) field. + Default: False + + Returns + ------- + :class:`numpy.ndarray` + The cleaned field. + """ + normalizer = _check_normalizer(normalizer) + if check_shape: + if mesh_type != "unstructured": + pos, shape, dim = format_struct_pos_shape( + pos, field.shape, check_stacked_shape=stacked + ) + else: + pos, shape, dim = format_unstruct_pos_shape( + pos, field.shape, check_stacked_shape=stacked + ) + field = np.array(field, dtype=np.double).reshape(shape) + else: + dim = len(pos) + if not stacked: + field = [field] + field_cnt = len(field) + for i in range(field_cnt): + field[i] -= eval_func(trend, pos, dim, mesh_type, value_type, True) + if fit_normalizer: + normalizer.fit(field) + field = normalizer.normalize(field) + for i in range(field_cnt): + field[i] -= eval_func(mean, pos, dim, mesh_type, value_type, True) + return field if stacked else field[0] From 614a6f82739940bbef94341dd8ef29472752090b Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Tue, 2 Feb 2021 21:56:47 +0100 Subject: [PATCH 64/97] Variogram: add mean, normalizer, trend and fit_normalizer to variogram estimation --- gstools/variogram/variogram.py | 31 +++++++++++++++++++++++++++++-- 1 file changed, 29 insertions(+), 2 deletions(-) diff --git a/gstools/variogram/variogram.py b/gstools/variogram/variogram.py index ddf9835c7..a06b6be5b 100644 --- a/gstools/variogram/variogram.py +++ b/gstools/variogram/variogram.py @@ -27,6 +27,7 @@ directional, ) from gstools.variogram.binning import standard_bins +from gstools.normalizer.tools import remove_trend_norm_mean __all__ = [ "vario_estimate", @@ -81,6 +82,10 @@ def vario_estimate( mask=np.ma.nomask, mesh_type="unstructured", return_counts=False, + mean=None, + normalizer=None, + trend=None, + fit_normalizer=False, ): r""" Estimates the empirical variogram. @@ -126,8 +131,10 @@ def vario_estimate( You can pass a list of fields, that will be used simultaneously. This could be helpful, when there are multiple realizations at the same points, with the same statistical properties. - bin_edges : :class:`numpy.ndarray` - the bins on which the variogram will be calculated + bin_edges : :class:`numpy.ndarray`, optional + the bins on which the variogram will be calculated. + If None are given, standard bins provided by the :any:`standard_bins` + routine will be used. Default: :any:`None` sampling_size : :class:`int` or :any:`None`, optional for large input data, this method can take a long time to compute the variogram, therefore this argument specifies @@ -196,6 +203,19 @@ def vario_estimate( if set to true, this function will also return the number of data points found at each lag distance as a third return value Default: False + mean : :class:`float`, optional + mean value used to shift normalized input data. + Could also be a callable. The default is None. + normalizer : :any:`None` or :any:`Normalizer`, optional + Normalizer to be applied to the input data to gain normality. + The default is None. + trend : :any:`None` or :class:`float` or :any:`callable`, optional + A callable trend function. Should have the signiture: f(x, [y, z, ...]) + If no normalizer is applied, this behaves equal to 'mean'. + The default is None. + fit_normalizer : :class:`bool`, optional + Wheater to fit the data-normalizer to the given (detrended) field. + Default: False Returns ------- @@ -296,6 +316,13 @@ def vario_estimate( if bin_edges is None: bin_edges = standard_bins(pos, dim, latlon) bin_centres = (bin_edges[:-1] + bin_edges[1:]) / 2.0 + # normalize field + field = remove_trend_norm_mean( + *(pos, field, mean, normalizer, trend), + check_shape=False, + stacked=True, + fit_normalizer=fit_normalizer, + ) # select variogram estimator cython_estimator = _set_estimator(estimator) # run From 48ed439d36280024eeffc2be19386149c64b2d8f Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Tue, 2 Feb 2021 22:00:35 +0100 Subject: [PATCH 65/97] Field: use apply_mean_normalizer_trend routine --- gstools/field/base.py | 24 +++++++++--------------- 1 file changed, 9 insertions(+), 15 deletions(-) diff --git a/gstools/field/base.py b/gstools/field/base.py index 2f50a7587..8a1c352ea 100755 --- a/gstools/field/base.py +++ b/gstools/field/base.py @@ -17,6 +17,7 @@ from gstools.tools.misc import eval_func from gstools.normalizer import Normalizer from gstools.field.tools import mesh_call, to_vtk_helper, fmt_mean_norm_trend +from gstools.normalizer.tools import apply_mean_norm_trend, _check_normalizer __all__ = ["Field"] @@ -197,17 +198,17 @@ def post_field(self, field, name="field", process=True, save=True): if process: if self.pos is None: raise ValueError("post_field: no 'pos' tuple set for field.") - kwargs = dict( + field = apply_mean_norm_trend( pos=self.pos, - dim=self.model.field_dim, + field=field, mesh_type=self.mesh_type, value_type=self.value_type, - broadcast=True, + mean=self.mean, + normalizer=self.normalizer, + trend=self.trend, + check_shape=False, + stacked=False, ) - # apply mean - normalizer - trend - field += eval_func(func_val=self.mean, **kwargs) - field = self.normalizer.denormalize(field) - field += eval_func(func_val=self.trend, **kwargs) if save: setattr(self, name, field) return field @@ -329,14 +330,7 @@ def normalizer(self): @normalizer.setter def normalizer(self, normalizer): - if isinstance(normalizer, type) and issubclass(normalizer, Normalizer): - self._normalizer = normalizer() - elif isinstance(normalizer, Normalizer): - self._normalizer = normalizer - elif normalizer is None: - self._normalizer = Normalizer() - else: - raise ValueError("Field: 'normalizer' not of type 'Normalizer'.") + self._normalizer = _check_normalizer(normalizer) @property def trend(self): From 2c9dcd95194330364b2ec135245eb0f317781d9e Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Tue, 2 Feb 2021 22:00:48 +0100 Subject: [PATCH 66/97] Doc: minor fix --- gstools/field/__init__.py | 2 -- gstools/krige/__init__.py | 2 -- 2 files changed, 4 deletions(-) diff --git a/gstools/field/__init__.py b/gstools/field/__init__.py index 01a7fc681..f84d87fe4 100644 --- a/gstools/field/__init__.py +++ b/gstools/field/__init__.py @@ -27,8 +27,6 @@ :toctree: generated Field - ----- """ from gstools.field.base import Field diff --git a/gstools/krige/__init__.py b/gstools/krige/__init__.py index 420b3f77f..76670d77a 100644 --- a/gstools/krige/__init__.py +++ b/gstools/krige/__init__.py @@ -16,8 +16,6 @@ Universal ExtDrift Detrended - ----- """ from gstools.krige.base import Krige from gstools.krige.methods import ( From a9771ace24e811fe1a64b405b1dcf5562b454b3a Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Tue, 2 Feb 2021 23:24:37 +0100 Subject: [PATCH 67/97] Normalizer: fix some wrong data ranges for negative lambdas --- gstools/normalizer/methods.py | 29 +++++++++++++++++++---------- 1 file changed, 19 insertions(+), 10 deletions(-) diff --git a/gstools/normalizer/methods.py b/gstools/normalizer/methods.py index 6a9d9a12e..1a08bbf7c 100644 --- a/gstools/normalizer/methods.py +++ b/gstools/normalizer/methods.py @@ -78,10 +78,13 @@ class BoxCox(Normalizer): def denormalize_range(self): """:class:`tuple`: Valid range for output data depending on lmbda. - (-1/lmbda, inf) + `(-1/lmbda, inf)` or `(-inf, -1/lmbda)` """ - with np.errstate(divide="ignore"): - return (-np.divide(1, self.lmbda), np.inf) + if np.isclose(self.lmbda, 0): + return (-np.inf, np.inf) + if self.lmbda < 0: + return (-np.inf, np.divide(1, self.lmbda)) + return (-np.divide(1, self.lmbda), np.inf) def _denormalize(self, data): if np.isclose(self.lmbda, 0): @@ -141,7 +144,7 @@ class BoxCoxShift(Normalizer): def normalize_range(self): """:class:`tuple`: Valid range for input data depending on shift. - (-shift, inf) + `(-shift, inf)` """ return (-self.shift, np.inf) @@ -149,10 +152,13 @@ def normalize_range(self): def denormalize_range(self): """:class:`tuple`: Valid range for output data depending on lmbda. - (-1/lmbda, inf) + `(-1/lmbda, inf)` or `(-inf, -1/lmbda)` """ - with np.errstate(divide="ignore"): - return (-np.divide(1, self.lmbda), np.inf) + if np.isclose(self.lmbda, 0): + return (-np.inf, np.inf) + if self.lmbda < 0: + return (-np.inf, np.divide(1, self.lmbda)) + return (-np.divide(1, self.lmbda), np.inf) def _denormalize(self, data): if np.isclose(self.lmbda, 0): @@ -331,10 +337,13 @@ class Manly(Normalizer): def denormalize_range(self): """:class:`tuple`: Valid range for output data depending on lmbda. - (-1/lmbda, inf) + `(-1/lmbda, inf)` or `(-inf, -1/lmbda)` """ - with np.errstate(divide="ignore"): - return (-np.divide(1, self.lmbda), np.inf) + if np.isclose(self.lmbda, 0): + return (-np.inf, np.inf) + if self.lmbda < 0: + return (-np.inf, np.divide(1, self.lmbda)) + return (-np.divide(1, self.lmbda), np.inf) def _denormalize(self, data): if np.isclose(self.lmbda, 0): From 7b7278d372b3504b0816ae30c067244c78c7a6c1 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Tue, 2 Feb 2021 23:25:02 +0100 Subject: [PATCH 68/97] Examples: add normalizer comparison example --- examples/10_normalizer/02_compare.py | 65 ++++++++++++++++++++++++++++ 1 file changed, 65 insertions(+) create mode 100644 examples/10_normalizer/02_compare.py diff --git a/examples/10_normalizer/02_compare.py b/examples/10_normalizer/02_compare.py new file mode 100644 index 000000000..2dd2cf4c6 --- /dev/null +++ b/examples/10_normalizer/02_compare.py @@ -0,0 +1,65 @@ +""" +Normalizer Comparison +--------------------- + +Let's compare the transformation behavior of the provided normalizers. + +But first, we define a convenience routine and make some imports as always. +""" +import numpy as np +import gstools as gs +import matplotlib.pyplot as plt + + +def dashes(i=1, max_n=12, width=1): + """Return line dashes.""" + return i * [width, width] + [max_n * 2 * width - 2 * i * width, width] + + +############################################################################### +# We select 4 normalizers depending on a single parameter lambda and +# plot their transformation behavior within the interval [-5, 5]. +# +# For the shape parameter lambda, we create a list of 8 values ranging from +# -1 to 2.5. + +lmbdas = [i * 0.5 for i in range(-2, 6)] +normalizers = [ + gs.normalizer.BoxCox, + gs.normalizer.YeoJohnson, + gs.normalizer.Modulus, + gs.normalizer.Manly, +] + +############################################################################### +# Let's plot them! + +fig, ax = plt.subplots(2, 2, figsize=[8, 8]) +for i, norm in enumerate(normalizers): + # correctly setting the data range + x_rng = norm.normalize_range + x = np.linspace(max(-5, x_rng[0] + 0.01), min(5, x_rng[1] - 0.01)) + for j, lmbda in enumerate(lmbdas): + ax.flat[i].plot( + x, + norm(lmbda=lmbda).normalize(x), + label=r"$\lambda=" + str(lmbda) + "$", + color="k", + alpha=0.2 + j * 0.1, + dashes=dashes(j), + ) + # axis formatting + ax.flat[i].grid(which="both", color="grey", linestyle="-", alpha=0.2) + ax.flat[i].set_ylim((-5, 5)) + ax.flat[i].set_xlim((-5, 5)) + ax.flat[i].set_title(norm().name) +# figure formatting +handles, labels = ax.flat[-1].get_legend_handles_labels() +fig.legend(handles, labels, loc="lower center", ncol=4, handlelength=3.0) +fig.suptitle("Normalizer Comparison", fontsize=20) +fig.show() + +############################################################################### +# The missing :any:`LogNormal` transformation is covered by the :any:`BoxCox` +# transformation for lambda=0. The :any:`BoxCoxShift` transformation is +# simply the :any:`BoxCox` transformation shifted on the X-axis. From f3059959af800a0c6174b1be935b052d2757bdab Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Tue, 2 Feb 2021 23:27:45 +0100 Subject: [PATCH 69/97] Examples: saver call of normalize_range --- examples/10_normalizer/02_compare.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/10_normalizer/02_compare.py b/examples/10_normalizer/02_compare.py index 2dd2cf4c6..fe0bff510 100644 --- a/examples/10_normalizer/02_compare.py +++ b/examples/10_normalizer/02_compare.py @@ -37,7 +37,7 @@ def dashes(i=1, max_n=12, width=1): fig, ax = plt.subplots(2, 2, figsize=[8, 8]) for i, norm in enumerate(normalizers): # correctly setting the data range - x_rng = norm.normalize_range + x_rng = norm().normalize_range x = np.linspace(max(-5, x_rng[0] + 0.01), min(5, x_rng[1] - 0.01)) for j, lmbda in enumerate(lmbdas): ax.flat[i].plot( From 9b55002fae4c10064a0aeec58443bef444e73659 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Fri, 5 Feb 2021 14:55:38 +0100 Subject: [PATCH 70/97] Field: add 'dim' property that's identical to model.field_dim --- gstools/field/base.py | 12 ++++++++---- gstools/field/plot.py | 6 +++--- gstools/field/tools.py | 12 ++++++------ gstools/krige/base.py | 21 ++++++++++----------- 4 files changed, 27 insertions(+), 24 deletions(-) diff --git a/gstools/field/base.py b/gstools/field/base.py index 8a1c352ea..1d7fbfc58 100755 --- a/gstools/field/base.py +++ b/gstools/field/base.py @@ -154,19 +154,18 @@ def pre_pos(self, pos, mesh_type="unstructured"): """ # save mesh-type self.mesh_type = mesh_type - dim = self.model.field_dim # save pos tuple if mesh_type != "unstructured": - pos, shape = format_struct_pos_dim(pos, dim) + pos, shape = format_struct_pos_dim(pos, self.dim) self.pos = pos pos = gen_mesh(pos) else: - pos = np.array(pos, dtype=np.double).reshape(dim, -1) + pos = np.array(pos, dtype=np.double).reshape(self.dim, -1) self.pos = pos shape = np.shape(pos[0]) # prepend dimension if we have a vector field if self.value_type == "vector": - shape = (self.model.dim,) + shape + shape = (self.dim,) + shape if self.model.latlon: raise ValueError("Field: Vector fields not allowed for latlon") # return isometrized pos tuple and resulting field shape @@ -352,6 +351,11 @@ def value_type(self, value_type): raise ValueError("Field: value type not in {}".format(VALUE_TYPES)) self._value_type = value_type + @property + def dim(self): + """:class:`int`: Dimension of the field.""" + return self.model.field_dim + @property def name(self): """:class:`str`: The name of the class.""" diff --git a/gstools/field/plot.py b/gstools/field/plot.py index edfbd6ee3..572d88a80 100644 --- a/gstools/field/plot.py +++ b/gstools/field/plot.py @@ -42,13 +42,13 @@ def plot_field(fld, field="field", fig=None, ax=None): # pragma: no cover """ plot_field = getattr(fld, field) assert not (fld.pos is None or plot_field is None) - if fld.model.field_dim == 1: + if fld.dim == 1: ax = _plot_1d(fld.pos, plot_field, fig, ax) - elif fld.model.field_dim == 2: + elif fld.dim == 2: ax = _plot_2d( fld.pos, plot_field, fld.mesh_type, fig, ax, fld.model.latlon ) - elif fld.model.field_dim == 3: + elif fld.dim == 3: ax = _plot_3d(fld.pos, plot_field, fld.mesh_type, fig, ax) else: raise ValueError("Field.plot: only possible for dim=1,2,3!") diff --git a/gstools/field/tools.py b/gstools/field/tools.py index ef57285db..daaa4693f 100644 --- a/gstools/field/tools.py +++ b/gstools/field/tools.py @@ -154,15 +154,15 @@ def mesh_call( pass if isinstance(direction, str) and direction == "all": - select = list(range(f_cls.model.field_dim)) + select = list(range(f_cls.dim)) elif isinstance(direction, str): - select = _get_select(direction)[: f_cls.model.field_dim] + select = _get_select(direction)[: f_cls.dim] else: - select = direction[: f_cls.model.field_dim] - if len(select) < f_cls.model.field_dim: + select = direction[: f_cls.dim] + if len(select) < f_cls.dim: raise ValueError( "Field.mesh: need at least {} direction(s), got '{}'".format( - f_cls.model.field_dim, direction + f_cls.dim, direction ) ) # convert pyvista mesh @@ -190,7 +190,7 @@ def mesh_call( offset = [] length = [] mesh_dim = mesh.points.shape[1] - if mesh_dim < f_cls.model.field_dim: + if mesh_dim < f_cls.dim: raise ValueError("Field.mesh: mesh dimension too low!") pnts = np.empty((0, mesh_dim), dtype=np.double) for cell in mesh.cells: diff --git a/gstools/krige/base.py b/gstools/krige/base.py index 4b60c24dc..e60196fdb 100755 --- a/gstools/krige/base.py +++ b/gstools/krige/base.py @@ -397,13 +397,16 @@ def _get_dists(self, pos1, pos2=None, pos2_slice=(0, None)): def get_mean(self, post_process=True): """Calculate the estimated mean of the detrended field. + Parameters + ---------- + post_process : :class:`bool`, optional + Whether to apply field-mean and normalizer. + Default: `True` + Returns ------- mean : :class:`float` or :any:`None` Mean of the Kriging System. - post_process : :class:`bool`, optional - Whether to apply field-mean and normalizer. - Default: `True` Notes ----- @@ -491,7 +494,7 @@ def set_condition( raise ValueError("Krige.set_condition: missing cond_pos/cond_val.") # correctly format cond_pos and cond_val self._cond_pos, self._cond_val = set_condition( - cond_pos, cond_val, self.model.field_dim + cond_pos, cond_val, self.dim ) if fit_normalizer: # fit normalizer to detrended data self.normalizer.fit(self.cond_val - self.cond_trend) @@ -544,7 +547,7 @@ def set_drift_functions(self, drift_functions=None): self._drift_functions = [] elif isinstance(drift_functions, (str, int)): self._drift_functions = get_drift_functions( - self.model.field_dim, drift_functions + self.dim, drift_functions ) else: if isinstance(drift_functions, collections.abc.Iterator): @@ -623,16 +626,12 @@ def cond_ext_drift(self): @property def cond_mean(self): """:class:`numpy.ndarray`: Trend at the conditions.""" - return eval_func( - self.mean, self.cond_pos, self.model.field_dim, broadcast=True - ) + return eval_func(self.mean, self.cond_pos, self.dim, broadcast=True) @property def cond_trend(self): """:class:`numpy.ndarray`: Trend at the conditions.""" - return eval_func( - self.trend, self.cond_pos, self.model.field_dim, broadcast=True - ) + return eval_func(self.trend, self.cond_pos, self.dim, broadcast=True) @property def unbiased(self): From 9bb8daafc08880eac38a15c3b86a0e540bcf9785 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Fri, 12 Feb 2021 20:00:11 +0100 Subject: [PATCH 71/97] Field: allow vector valued mean/trend for vector valued fields --- gstools/field/base.py | 17 ++++++++++++----- gstools/field/srf.py | 4 ++++ gstools/tools/misc.py | 30 +++++++++++++++++++++--------- 3 files changed, 37 insertions(+), 14 deletions(-) diff --git a/gstools/field/base.py b/gstools/field/base.py index 1d7fbfc58..12a6e5915 100755 --- a/gstools/field/base.py +++ b/gstools/field/base.py @@ -14,8 +14,6 @@ import numpy as np from gstools.covmodel.base import CovModel from gstools.tools.geometric import format_struct_pos_dim, gen_mesh -from gstools.tools.misc import eval_func -from gstools.normalizer import Normalizer from gstools.field.tools import mesh_call, to_vtk_helper, fmt_mean_norm_trend from gstools.normalizer.tools import apply_mean_norm_trend, _check_normalizer @@ -25,6 +23,15 @@ """:class:`list` of :class:`str`: valid field value types.""" +def _set_mean_trend(value, dim): + if callable(value) or value is None: + return value + value = np.array(value, dtype=np.double).ravel() + if value.size > 1 and value.size != dim: # vector mean + raise ValueError("Mean/Trend: Wrong size ({})".format(value)) + return value if value.size > 1 else value.item() + + class Field: """A base class for random fields, kriging fields, etc. @@ -320,7 +327,7 @@ def mean(self): @mean.setter def mean(self, mean): - self._mean = mean if (callable(mean) or mean is None) else float(mean) + self._mean = _set_mean_trend(mean, self.dim) @property def normalizer(self): @@ -337,8 +344,8 @@ def trend(self): return self._trend @trend.setter - def trend(self, tren): - self._trend = tren if (callable(tren) or tren is None) else float(tren) + def trend(self, trend): + self._trend = _set_mean_trend(trend, self.dim) @property def value_type(self): diff --git a/gstools/field/srf.py b/gstools/field/srf.py index d92a46c82..50ceed654 100644 --- a/gstools/field/srf.py +++ b/gstools/field/srf.py @@ -163,6 +163,10 @@ def set_generator(self, generator, **generator_kwargs): self.value_type = self._generator.value_type else: raise ValueError("gstools.SRF: Unknown generator: " + generator) + for val in [self.mean, self.trend]: + if not callable(val) and val is not None: + if np.size(val) > 1 and self.value_type == "scalar": + raise ValueError("Mean/Trend: Wrong size ({})".format(val)) @property def generator(self): diff --git a/gstools/tools/misc.py b/gstools/tools/misc.py index 7fc549b19..6edacd688 100755 --- a/gstools/tools/misc.py +++ b/gstools/tools/misc.py @@ -63,10 +63,11 @@ def eval_func( Function values at the given points. """ # care about scalar inputs - if broadcast and not callable(func_val): - return 0.0 if func_val is None else float(func_val) + func_val = 0 if func_val is None else func_val + if broadcast and not callable(func_val) and np.size(func_val) == 1: + return np.array(func_val, dtype=np.double).item() if not callable(func_val): - func_val = _func_from_scalar(func_val) + func_val = _func_from_single_val(func_val, dim, value_type=value_type) # care about mesh and function call if mesh_type != "unstructured": pos, shape = format_struct_pos_dim(pos, dim) @@ -80,13 +81,24 @@ def eval_func( return np.reshape(func_val(*pos), shape) -def _func_from_scalar(value, value_type="scalar"): - value = 0.0 if value is None else float(value) +def _func_from_single_val(value, dim=None, value_type="scalar"): + # care about broadcasting vector values for each dim + v_d = dim if value_type == "vector" else 1 # value dim + if v_d is None: + raise ValueError("_func_from_single_val: dim needed for vector value.") + value = np.array(value, dtype=np.double).ravel()[:v_d] + # fill up vector valued output to dimension with last value + value = np.pad( + value, (0, v_d - len(value)), "constant", constant_values=value[-1] + ) def _f(*pos): - if value_type == "vector": - return np.full_like(pos, value, dtype=np.double) - # scalar field has same shape like a single axis - return np.full_like(pos[0], value, dtype=np.double) + # zip uses shortes len of iterables given (correct for scalar value) + return np.concatenate( + [ + np.full_like(p, val, dtype=np.double) + for p, val in zip(pos, value) + ] + ) return _f From 56ec497e2d52079c55f74ff8012d0fa9876e7661 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Fri, 12 Feb 2021 20:04:22 +0100 Subject: [PATCH 72/97] Field: bugfix in mesh method: transpose vector fields --- gstools/field/tools.py | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/gstools/field/tools.py b/gstools/field/tools.py index daaa4693f..074a44b2e 100644 --- a/gstools/field/tools.py +++ b/gstools/field/tools.py @@ -102,7 +102,7 @@ def to_vtk_helper( def mesh_call( f_cls, mesh, points="centroids", direction="all", name="field", **kwargs ): - """Generate a field on a given meshio or ogs5py mesh. + """Generate a field on a given meshio, ogs5py or pyvista mesh. Parameters ---------- @@ -133,9 +133,10 @@ def mesh_call( if a meshio or PyVista mesh was given. See: https://github.com/nschloe/meshio - See: https://github.com/pyvista/pyvista - See: :any:`Field.__call__` + See: https://github.com/GeoStat-Framework/ogs5py + + See: https://github.com/pyvista/pyvista """ has_pyvista = False has_ogs5py = False @@ -174,6 +175,8 @@ def mesh_call( out = f_cls.unstructured(pos=pnts, **kwargs) # Deal with the output fields = [out] if isinstance(out, np.ndarray) else out + if f_cls.value_type == "vector": + fields = [f.T for f in fields] for f_name, field in zip(_names(name, len(fields)), fields): mesh[f_name] = field # convert ogs5py mesh @@ -202,6 +205,8 @@ def mesh_call( pnts = pnts.T[select] out = f_cls.unstructured(pos=pnts, **kwargs) fields = [out] if isinstance(out, np.ndarray) else out + if f_cls.value_type == "vector": + fields = [f.T for f in fields] f_lists = [] for field in fields: f_list = [] @@ -213,6 +218,8 @@ def mesh_call( else: out = f_cls.unstructured(pos=mesh.points.T[select], **kwargs) fields = [out] if isinstance(out, np.ndarray) else out + if f_cls.value_type == "vector": + fields = [f.T for f in fields] for f_name, field in zip(_names(name, len(fields)), fields): mesh.point_data[f_name] = field else: From 619f24270a9633c85910de7cd8dde77b79e40ead Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Fri, 12 Feb 2021 23:18:45 +0100 Subject: [PATCH 73/97] Examples: update 3D vector example to use mesh method --- examples/04_vector_field/01_3d_vector_field.py | 17 +++-------------- 1 file changed, 3 insertions(+), 14 deletions(-) diff --git a/examples/04_vector_field/01_3d_vector_field.py b/examples/04_vector_field/01_3d_vector_field.py index 89b22b11b..3089ba05e 100755 --- a/examples/04_vector_field/01_3d_vector_field.py +++ b/examples/04_vector_field/01_3d_vector_field.py @@ -7,7 +7,6 @@ externally defined and it will be generated by PyVista. """ # sphinx_gallery_thumbnail_path = 'pics/GS_3d_vector_field.png' -import numpy as np import gstools as gs import pyvista as pv @@ -16,25 +15,15 @@ ############################################################################### # create a uniform grid with PyVista -nx, ny, nz = 40, 30, 10 dim, spacing, origin = (40, 30, 10), (1, 1, 1), (-10, 0, 0) mesh = pv.UniformGrid(dim, spacing, origin) -x = mesh.points[:, 0] -y = mesh.points[:, 1] -z = mesh.points[:, 2] ############################################################################### # create an incompressible random 3d velocity field on the given mesh +# with added mean velocity in x-direction model = gs.Gaussian(dim=3, var=3, len_scale=1.5) -srf = gs.SRF(model, generator='VectorField', seed=198412031) -srf((x, y, z), mesh_type='unstructured') - -# add a mean velocity in x-direction -srf.field[0, :] += 0.5 - -############################################################################### -# add the velocity field to the mesh object -mesh["Velocity"] = srf.field.T +srf = gs.SRF(model, mean=(0.5, 0, 0), generator="VectorField", seed=198412031) +srf.mesh(mesh, points="points", name="Velocity") ############################################################################### # Now, we can do the plotting From 2ba9e432010b02697def9251c5ce7b331d8d1e4f Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sat, 13 Feb 2021 10:10:46 +0100 Subject: [PATCH 74/97] Field: fix repr of vector mean --- gstools/field/tools.py | 3 +++ gstools/tools/misc.py | 2 +- 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/gstools/field/tools.py b/gstools/field/tools.py index 074a44b2e..596f58948 100644 --- a/gstools/field/tools.py +++ b/gstools/field/tools.py @@ -15,6 +15,7 @@ from gstools.normalizer import Normalizer from gstools.tools.export import to_vtk, vtk_export +from gstools.tools.misc import list_format __all__ = ["fmt_mean_norm_trend", "to_vtk_helper", "mesh_call"] @@ -25,6 +26,8 @@ def _fmt_func_val(f_cls, func_val): return str(None) if callable(func_val): return "" # or format(func_val.__name__) + if np.size(func_val) > 1: + return list_format(func_val, prec=f_cls.model._prec) return "{0:.{p}}".format(float(func_val), p=f_cls.model._prec) diff --git a/gstools/tools/misc.py b/gstools/tools/misc.py index 6edacd688..f8a8d5123 100755 --- a/gstools/tools/misc.py +++ b/gstools/tools/misc.py @@ -20,7 +20,7 @@ def list_format(lst, prec): """Format a list of floats.""" return "[{}]".format( - ", ".join("{x:.{p}}".format(x=x, p=prec) for x in lst) + ", ".join("{x:.{p}}".format(x=float(x), p=prec) for x in lst) ) From 36ac370deb12fee05dbb35bbca9992e4dc417e19 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 22 Feb 2021 13:31:50 +0100 Subject: [PATCH 75/97] CondSRF: simplify field saving; remove unneeded generator value type check --- gstools/field/cond_srf.py | 11 +++-------- 1 file changed, 3 insertions(+), 8 deletions(-) diff --git a/gstools/field/cond_srf.py b/gstools/field/cond_srf.py index 33b5b73a6..34243b904 100644 --- a/gstools/field/cond_srf.py +++ b/gstools/field/cond_srf.py @@ -94,12 +94,9 @@ def __call__(self, pos, seed=np.nan, mesh_type="unstructured", **kwargs): ) field, krige_var = self.krige(pos, **kwargs) var_scale, nugget = self.get_scaling(krige_var, shape) - field = self.krige.post_field( - field=field + var_scale * self.raw_field + nugget, save=False - ) - self.krige.post_field(self.krige.field) - # processing is done in the kriging class - return self.post_field(field, process=False) + # need to use a copy to not alter "field" by reference + self.krige.post_field(self.krige.field.copy()) + return self.post_field(field + var_scale * self.raw_field + nugget) def get_scaling(self, krige_var, shape): """ @@ -146,8 +143,6 @@ def set_generator(self, generator, **generator_kwargs): self.value_type = self.generator.value_type else: raise ValueError("gstools.CondSRF: Unknown generator " + generator) - if self.value_type != "scalar": - raise ValueError("CondSRF: only scalar field value type allowed.") @property def krige(self): From d8059ff2662d7d22ce0fe8ad90df73357e58fcc3 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 22 Feb 2021 13:33:10 +0100 Subject: [PATCH 76/97] Field: skip some string formatting routines in coverage --- gstools/field/tools.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/gstools/field/tools.py b/gstools/field/tools.py index 596f58948..7f8757b92 100644 --- a/gstools/field/tools.py +++ b/gstools/field/tools.py @@ -21,7 +21,7 @@ __all__ = ["fmt_mean_norm_trend", "to_vtk_helper", "mesh_call"] -def _fmt_func_val(f_cls, func_val): +def _fmt_func_val(f_cls, func_val): # pragma: no cover if func_val is None: return str(None) if callable(func_val): @@ -31,12 +31,12 @@ def _fmt_func_val(f_cls, func_val): return "{0:.{p}}".format(float(func_val), p=f_cls.model._prec) -def _fmt_normalizer(f_cls): +def _fmt_normalizer(f_cls): # pragma: no cover norm = f_cls.normalizer return str(None) if norm.__class__ is Normalizer else norm.name -def fmt_mean_norm_trend(f_cls): +def fmt_mean_norm_trend(f_cls): # pragma: no cover """Format string repr. for mean, normalizer and trend of a field.""" args = [ "mean=" + _fmt_func_val(f_cls, f_cls.mean), From 79f7de48b749ca51d4148a97f4c93c700a5a465f Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 22 Feb 2021 13:33:52 +0100 Subject: [PATCH 77/97] Tools: skip corner case for mean from vector in coverage --- gstools/tools/misc.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gstools/tools/misc.py b/gstools/tools/misc.py index f8a8d5123..58006f347 100755 --- a/gstools/tools/misc.py +++ b/gstools/tools/misc.py @@ -84,7 +84,7 @@ def eval_func( def _func_from_single_val(value, dim=None, value_type="scalar"): # care about broadcasting vector values for each dim v_d = dim if value_type == "vector" else 1 # value dim - if v_d is None: + if v_d is None: # pragma: no cover raise ValueError("_func_from_single_val: dim needed for vector value.") value = np.array(value, dtype=np.double).ravel()[:v_d] # fill up vector valued output to dimension with last value From cf01e0db19a787f25e35b2f3d3ba5a79146ccf3d Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 22 Feb 2021 13:34:44 +0100 Subject: [PATCH 78/97] VectorField: add test for vector mean --- tests/test_incomprrandmeth.py | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) diff --git a/tests/test_incomprrandmeth.py b/tests/test_incomprrandmeth.py index 8d1a67bc6..6760fd72e 100644 --- a/tests/test_incomprrandmeth.py +++ b/tests/test_incomprrandmeth.py @@ -6,13 +6,13 @@ import copy import unittest import numpy as np -from gstools import Gaussian +import gstools as gs from gstools.field.generator import IncomprRandMeth class TestIncomprRandMeth(unittest.TestCase): def setUp(self): - self.cov_model_2d = Gaussian(dim=2, var=1.5, len_scale=2.5) + self.cov_model_2d = gs.Gaussian(dim=2, var=1.5, len_scale=2.5) self.cov_model_3d = copy.deepcopy(self.cov_model_2d) self.cov_model_3d.dim = 3 self.seed = 19031977 @@ -43,9 +43,20 @@ def test_unstruct_3d(self): self.assertAlmostEqual(modes[1, 0], -0.28049855754819514) def test_assertions(self): - cov_model_1d = Gaussian(dim=1, var=1.5, len_scale=2.5) + cov_model_1d = gs.Gaussian(dim=1, var=1.5, len_scale=2.5) self.assertRaises(ValueError, IncomprRandMeth, cov_model_1d) + def test_vector_mean(self): + srf = gs.SRF( + self.cov_model_2d, + mean=(0.5, 0), + generator="VectorField", + seed=198412031, + ) + srf.structured((self.x_grid, self.y_grid)) + self.assertAlmostEqual(np.mean(srf.field[0]), 1.3025621393180298) + self.assertAlmostEqual(np.mean(srf.field[1]), -0.04729596839446052) + if __name__ == "__main__": unittest.main() From 12bfaa0788c43a0f8d053b81ea51b159a6e349f7 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 22 Feb 2021 15:16:57 +0100 Subject: [PATCH 79/97] CovModel: remove redundant __ne__ method --- gstools/covmodel/base.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/gstools/covmodel/base.py b/gstools/covmodel/base.py index 563be3036..ec34f9997 100644 --- a/gstools/covmodel/base.py +++ b/gstools/covmodel/base.py @@ -1131,10 +1131,6 @@ def __eq__(self, other): return False return compare(self, other) - def __ne__(self, other): - """Compare CovModels.""" - return not self.__eq__(other) - def __repr__(self): """Return String representation.""" return model_repr(self) From 05eba49036baf994ebb13d3f3847572ee70e86ad Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 22 Feb 2021 15:17:52 +0100 Subject: [PATCH 80/97] Normalizer: add __eq__ magic method --- gstools/normalizer/base.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/gstools/normalizer/base.py b/gstools/normalizer/base.py index d4ce96613..debe4ae23 100644 --- a/gstools/normalizer/base.py +++ b/gstools/normalizer/base.py @@ -232,6 +232,17 @@ def _neg_kllf(par, dat): setattr(self, name, val) return {name: getattr(self, name) for name in all_names} + def __eq__(self, other): + """Compare Normalizers.""" + # check for correct base class + if type(self) is not type(other): + return False + # if base class is same, this is save + for val in self.default_parameter: + if not np.isclose(getattr(self, val), getattr(other, val)): + return False + return True + @property def name(self): """:class:`str`: The name of the normalizer class.""" From 801b6c9ffb789cbbbcbb88f26288daca84eaf620 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 22 Feb 2021 15:18:16 +0100 Subject: [PATCH 81/97] CondSRF: minor fix --- gstools/field/cond_srf.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/gstools/field/cond_srf.py b/gstools/field/cond_srf.py index 34243b904..0514c23e1 100644 --- a/gstools/field/cond_srf.py +++ b/gstools/field/cond_srf.py @@ -187,8 +187,8 @@ def trend(self): return self.krige.trend @trend.setter - def trend(self, tren): - self.krige.trend = tren + def trend(self, trend): + self.krige.trend = trend @property def value_type(self): From db02681a3a2cf3d595a5753d992655b09aab2126 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 22 Feb 2021 15:18:55 +0100 Subject: [PATCH 82/97] Tests: more CondSRF testing --- tests/test_condition.py | 52 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 52 insertions(+) diff --git a/tests/test_condition.py b/tests/test_condition.py index 32d4b079a..830681d48 100644 --- a/tests/test_condition.py +++ b/tests/test_condition.py @@ -97,6 +97,58 @@ def test_ordinary(self): self.assertAlmostEqual(val, field_1[i], places=2) self.assertAlmostEqual(val, field_2[dim * (i,)], places=2) + def test_raise_error(self): + self.assertRaises(ValueError, gs.CondSRF, gs.Gaussian()) + krige = gs.krige.Ordinary(gs.Stable(), self.cond_pos, self.cond_val) + self.assertRaises(ValueError, gs.CondSRF, krige, generator="unknown") + + def test_nugget(self): + model = gs.Gaussian( + nugget=0.01, + var=0.5, + len_scale=2, + anis=[0.1, 1], + angles=[0.5, 0, 0], + ) + krige = gs.krige.Ordinary( + model, self.cond_pos, self.cond_val, exact=True + ) + crf = gs.CondSRF(krige, seed=19970221) + field_1 = crf.unstructured(self.pos) + field_2 = crf.structured(self.pos) + for i, val in enumerate(self.cond_val): + self.assertAlmostEqual(val, field_1[i], places=2) + self.assertAlmostEqual(val, field_2[3 * (i,)], places=2) + + def test_setter(self): + krige1 = gs.krige.Krige(gs.Exponential(), self.cond_pos, self.cond_val) + krige2 = gs.krige.Krige( + gs.Gaussian(var=2), + self.cond_pos, + self.cond_val, + mean=-1, + trend=-2, + normalizer=gs.normalizer.YeoJohnson(), + ) + crf1 = gs.CondSRF(krige1) + crf2 = gs.CondSRF(krige2, seed=19970221) + # update settings + crf1.model = gs.Gaussian(var=2) + crf1.mean = -1 + crf1.trend = -2 + crf1.normalizer = gs.normalizer.YeoJohnson() + # check if setting went right + self.assertTrue(crf1.model == crf2.model) + self.assertTrue(crf1.normalizer == crf2.normalizer) + self.assertAlmostEqual(crf1.mean, crf2.mean) + self.assertAlmostEqual(crf1.trend, crf2.trend) + # reset kriging + crf1.krige.set_condition() + # compare fields + field1 = crf1(self.pos, seed=19970221) + field2 = crf2(self.pos) + self.assertTrue(np.all(np.isclose(field1, field2))) + if __name__ == "__main__": unittest.main() From 7f4af3cea143ee62d1c743400b88d239b260c64e Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 22 Feb 2021 15:19:19 +0100 Subject: [PATCH 83/97] Tests: test comparing normalizers --- tests/test_normalize.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/tests/test_normalize.py b/tests/test_normalize.py index a722f64be..49ce8070c 100644 --- a/tests/test_normalize.py +++ b/tests/test_normalize.py @@ -158,6 +158,17 @@ def test_parameterless(self): ) ) + def test_compare(self): + norm1 = gs.normalizer.BoxCox() + norm2 = gs.normalizer.BoxCox(lmbda=0.5) + norm3 = gs.normalizer.YeoJohnson() + norm4 = "this is not a normalizer" + # check campare + self.assertTrue(norm1 == norm1) + self.assertTrue(norm1 != norm2) + self.assertTrue(norm1 != norm3) + self.assertTrue(norm1 != norm4) + if __name__ == "__main__": unittest.main() From 79687ffa4616da967cc91f4f9d9fc57771484085 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 22 Feb 2021 19:53:28 +0100 Subject: [PATCH 84/97] Tests: check autofitting normalizer --- tests/test_normalize.py | 46 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 46 insertions(+) diff --git a/tests/test_normalize.py b/tests/test_normalize.py index 49ce8070c..413bce4b8 100644 --- a/tests/test_normalize.py +++ b/tests/test_normalize.py @@ -169,6 +169,52 @@ def test_compare(self): self.assertTrue(norm1 != norm3) self.assertTrue(norm1 != norm4) + def test_check(self): + self.assertRaises(ValueError, gs.field.Field, gs.Cubic(), normalizer=5) + + def test_auto_fit(self): + x = y = range(60) + pos = gs.tools.geometric.gen_mesh([x, y]) + model = gs.Gaussian(dim=2, var=1, len_scale=10) + srf = gs.SRF( + model, seed=20170519, normalizer=gs.normalizer.LogNormal() + ) + srf(pos) + ids = np.arange(srf.field.size) + samples = np.random.RandomState(20210201).choice( + ids, size=60, replace=False + ) + # sample conditioning points from generated field + cond_pos = pos[:, samples] + cond_val = srf.field[samples] + krige = gs.krige.Ordinary( + model=gs.Stable(dim=2), + cond_pos=cond_pos, + cond_val=cond_val, + normalizer=gs.normalizer.BoxCox(), + fit_normalizer=True, + fit_variogram=True, + ) + # test fitting during kriging + self.assertTrue(np.abs(krige.normalizer.lmbda - 0.0) < 1e-1) + self.assertAlmostEqual(krige.model.len_scale, 10.267877910391935) + self.assertAlmostEqual( + krige.model.sill, krige.normalizer.normalize(cond_val).var() + ) + # test fitting during vario estimate + emp_vario = gs.vario_estimate( + cond_pos, + cond_val, + normalizer=gs.normalizer.BoxCox, + fit_normalizer=True, + ) + model = gs.Stable(dim=2) + model.fit_variogram(*emp_vario) + self.assertAlmostEqual(model.var, 0.6426670183) + self.assertAlmostEqual(model.len_scale, 9.635193952) + self.assertAlmostEqual(model.nugget, 0.001617908408) + self.assertAlmostEqual(model.alpha, 2.0) + if __name__ == "__main__": unittest.main() From 63df1c03a2adc1db261e5814f622c51a3c5008a3 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 22 Feb 2021 19:54:29 +0100 Subject: [PATCH 85/97] Tests: check auto-init normalizers in krige --- tests/test_condition.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tests/test_condition.py b/tests/test_condition.py index 830681d48..587469d13 100644 --- a/tests/test_condition.py +++ b/tests/test_condition.py @@ -136,7 +136,8 @@ def test_setter(self): crf1.model = gs.Gaussian(var=2) crf1.mean = -1 crf1.trend = -2 - crf1.normalizer = gs.normalizer.YeoJohnson() + # also checking correctly setting uninitialized normalizer + crf1.normalizer = gs.normalizer.YeoJohnson # check if setting went right self.assertTrue(crf1.model == crf2.model) self.assertTrue(crf1.normalizer == crf2.normalizer) From 3a3fd47e4756f59dd654d7b5a2586e57a84d40fc Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 22 Feb 2021 20:00:22 +0100 Subject: [PATCH 86/97] Tests: reduce places in comparison (normalizer autofit) --- tests/test_normalize.py | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/tests/test_normalize.py b/tests/test_normalize.py index 413bce4b8..e4d1051a4 100644 --- a/tests/test_normalize.py +++ b/tests/test_normalize.py @@ -197,9 +197,11 @@ def test_auto_fit(self): ) # test fitting during kriging self.assertTrue(np.abs(krige.normalizer.lmbda - 0.0) < 1e-1) - self.assertAlmostEqual(krige.model.len_scale, 10.267877910391935) + self.assertAlmostEqual(krige.model.len_scale, 10.267877, places=4) self.assertAlmostEqual( - krige.model.sill, krige.normalizer.normalize(cond_val).var() + krige.model.sill, + krige.normalizer.normalize(cond_val).var(), + places=4, ) # test fitting during vario estimate emp_vario = gs.vario_estimate( @@ -210,10 +212,10 @@ def test_auto_fit(self): ) model = gs.Stable(dim=2) model.fit_variogram(*emp_vario) - self.assertAlmostEqual(model.var, 0.6426670183) - self.assertAlmostEqual(model.len_scale, 9.635193952) - self.assertAlmostEqual(model.nugget, 0.001617908408) - self.assertAlmostEqual(model.alpha, 2.0) + self.assertAlmostEqual(model.var, 0.6426670183, places=4) + self.assertAlmostEqual(model.len_scale, 9.635193952, places=4) + self.assertAlmostEqual(model.nugget, 0.001617908408, places=4) + self.assertAlmostEqual(model.alpha, 2.0, places=4) if __name__ == "__main__": From 3af20e0f7d14071bea503fa8a236e3705021584a Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 22 Feb 2021 20:44:46 +0100 Subject: [PATCH 87/97] Krige: better setting of drift_functions --- gstools/krige/base.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/gstools/krige/base.py b/gstools/krige/base.py index e60196fdb..e924a260d 100755 --- a/gstools/krige/base.py +++ b/gstools/krige/base.py @@ -148,9 +148,8 @@ def __init__( self._krige_pos = None self._cond_trend = None self._cond_ext_drift = np.array([]) - self._drift_functions = [] - if drift_functions is not None: - self.set_drift_functions(drift_functions) + self._drift_functions = None + self.set_drift_functions(drift_functions) self.set_condition( cond_pos, cond_val, From 7dce30574d4e3a604f15dcf35cf1f9790c4307f1 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 22 Feb 2021 20:46:14 +0100 Subject: [PATCH 88/97] Tests: test Error raise for wrong latlon input in vario estimate --- tests/test_variogram_unstructured.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/tests/test_variogram_unstructured.py b/tests/test_variogram_unstructured.py index bf0264548..6ca3ca7c1 100644 --- a/tests/test_variogram_unstructured.py +++ b/tests/test_variogram_unstructured.py @@ -404,6 +404,12 @@ def test_standard_bins(self): # no pos given self.assertRaises(ValueError, gs.standard_bins) + def test_raise(self): + # 1d field given for latlon estimation -> needs 2d + self.assertRaises( + ValueError, gs.vario_estimate, [[1, 2]], [1, 2], latlon=True + ) + if __name__ == "__main__": unittest.main() From b1be3322601e68fe6c02231b9160438c4ce5eeed Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 22 Feb 2021 20:46:48 +0100 Subject: [PATCH 89/97] Tests: krige: raise check; better check kriging the mean --- tests/test_krige.py | 64 ++++++++++++++++++++++++++++++++++++--------- 1 file changed, 51 insertions(+), 13 deletions(-) diff --git a/tests/test_krige.py b/tests/test_krige.py index 6d4d935d5..af7984b41 100644 --- a/tests/test_krige.py +++ b/tests/test_krige.py @@ -5,16 +5,20 @@ import numpy as np import unittest -from gstools import Gaussian, Exponential, Spherical, krige, SRF +import gstools as gs def trend(*xyz): return xyz[0] +def mean_func(*xyz): + return 2 * xyz[0] + + class TestKrige(unittest.TestCase): def setUp(self): - self.cov_models = [Gaussian, Exponential, Spherical] + self.cov_models = [gs.Gaussian, gs.Exponential, gs.Spherical] self.dims = range(1, 4) self.data = np.array( [ @@ -58,7 +62,7 @@ def test_simple(self): anis=[0.9, 0.8], angles=[2, 1, 0.5], ) - simple = krige.Simple( + simple = gs.krige.Simple( model, self.cond_pos[:dim], self.cond_val, self.mean ) field_1, __ = simple.unstructured(self.grids[dim - 1]) @@ -83,7 +87,7 @@ def test_ordinary(self): anis=[0.9, 0.8], angles=[2, 1, 0.5], ) - ordinary = krige.Ordinary( + ordinary = gs.krige.Ordinary( model, self.cond_pos[:dim], self.cond_val, @@ -112,7 +116,7 @@ def test_universal(self): anis=[0.9, 0.8], angles=[2, 1, 0.5], ) - universal = krige.Universal( + universal = gs.krige.Universal( model, self.cond_pos[:dim], self.cond_val, drift ) field_1, __ = universal.unstructured(self.grids[dim - 1]) @@ -137,7 +141,7 @@ def test_detrended(self): anis=[0.5, 0.2], angles=[0.4, 0.2, 0.1], ) - detrended = krige.Detrended( + detrended = gs.krige.Detrended( model, self.cond_pos[:dim], self.cond_val, trend ) field_1, __ = detrended.unstructured(self.grids[dim - 1]) @@ -157,14 +161,14 @@ def test_extdrift(self): cond_drift = [] for i, grid in enumerate(self.grids): dim = i + 1 - model = Exponential( + model = gs.Exponential( dim=dim, var=2, len_scale=10, anis=[0.9, 0.8], angles=[2, 1, 0.5], ) - srf = SRF(model) + srf = gs.SRF(model) field = srf(grid) ext_drift.append(field) field = field.reshape(self.grid_shape[:dim]) @@ -179,7 +183,7 @@ def test_extdrift(self): anis=[0.5, 0.2], angles=[0.4, 0.2, 0.1], ) - extdrift = krige.ExtDrift( + extdrift = gs.krige.ExtDrift( model, self.cond_pos[:dim], self.cond_val, @@ -202,7 +206,6 @@ def test_extdrift(self): ) def test_pseudo(self): - for Model in self.cov_models: for dim in self.dims: model = Model( @@ -213,7 +216,7 @@ def test_pseudo(self): angles=[0.4, 0.2, 0.1], ) for meth in self.p_meth: - krig = krige.Krige( + krig = gs.krige.Krige( model, self.p_data[:dim], self.p_vals, unbiased=False ) field, __ = krig([0, 0, 0][:dim]) @@ -224,7 +227,6 @@ def test_pseudo(self): ) def test_error(self): - for Model in self.cov_models: for dim in self.dims: model = Model( @@ -235,7 +237,7 @@ def test_error(self): anis=[0.9, 0.8], angles=[2, 1, 0.5], ) - ordinary = krige.Ordinary( + ordinary = gs.krige.Ordinary( model, self.cond_pos[:dim], self.cond_val, @@ -248,6 +250,42 @@ def test_error(self): self.assertAlmostEqual(err[1], model.nugget, places=2) self.assertAlmostEqual(err[4], model.nugget, places=2) + def test_raise(self): + # no cond_pos/cond_val given + self.assertRaises(ValueError, gs.krige.Krige, gs.Stable(), None, None) + + def test_krige_mean(self): + # check for constant mean (simple kriging) + krige = gs.krige.Simple(gs.Gaussian(), self.cond_pos, self.cond_val) + mean_f = krige.structured(self.pos, only_mean=True) + self.assertTrue(np.all(np.isclose(mean_f, 0))) + krige = gs.krige.Simple( + gs.Gaussian(), + self.cond_pos, + self.cond_val, + mean=mean_func, + normalizer=gs.normalizer.YeoJohnson, + trend=trend, + ) + # check applying mean, norm, trend + mean_f1 = krige.structured(self.pos, only_mean=True) + mean_f2 = gs.normalizer.tools.apply_mean_norm_trend( + self.pos, + np.zeros(tuple(map(len, self.pos))), + mean=mean_func, + normalizer=gs.normalizer.YeoJohnson, + trend=trend, + mesh_type="structured", + ) + self.assertTrue(np.all(np.isclose(mean_f1, mean_f2))) + krige = gs.krige.Simple(gs.Gaussian(), self.cond_pos, self.cond_val) + mean_f = krige.structured(self.pos, only_mean=True) + self.assertTrue(np.all(np.isclose(mean_f, 0))) + # check for constant mean (ordinary kriging) + krige = gs.krige.Ordinary(gs.Gaussian(), self.cond_pos, self.cond_val) + mean_f = krige.structured(self.pos, only_mean=True) + self.assertTrue(np.all(np.isclose(mean_f, krige.get_mean()))) + if __name__ == "__main__": unittest.main() From c6d711fa682229aed639cdd8cd4498371995848b Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 22 Feb 2021 20:57:54 +0100 Subject: [PATCH 90/97] Krige: remove ext-drift deleter (confusing) --- gstools/krige/base.py | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/gstools/krige/base.py b/gstools/krige/base.py index e924a260d..1124c7c1e 100755 --- a/gstools/krige/base.py +++ b/gstools/krige/base.py @@ -471,12 +471,6 @@ def set_condition( assuming the sill to be the data variance and with the standard bins provided by the :any:`standard_bins` routine. Default: False - - Notes - ----- - In order to remove an existing external drift use: - - >>> del Krige.cond_ext_drift """ # only use existing external drift, if no new positions are given ext_drift = ( @@ -618,10 +612,6 @@ def cond_ext_drift(self): """:class:`numpy.ndarray`: The ext. drift at the conditions.""" return self._cond_ext_drift - @cond_ext_drift.deleter - def cond_ext_drift(self): - self._cond_ext_drift = np.array([]) - @property def cond_mean(self): """:class:`numpy.ndarray`: Trend at the conditions.""" From 5d16fd5d7a2af69bb7b67a9e5efd45174a5640ff Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 22 Feb 2021 21:14:11 +0100 Subject: [PATCH 91/97] Tests: minor improvements --- gstools/tools/misc.py | 2 +- tests/test_covmodel.py | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/gstools/tools/misc.py b/gstools/tools/misc.py index 58006f347..72818f8c7 100755 --- a/gstools/tools/misc.py +++ b/gstools/tools/misc.py @@ -17,7 +17,7 @@ __all__ = ["list_format", "eval_func"] -def list_format(lst, prec): +def list_format(lst, prec): # pragma: no cover """Format a list of floats.""" return "[{}]".format( ", ".join("{x:.{p}}".format(x=float(x), p=prec) for x in lst) diff --git a/tests/test_covmodel.py b/tests/test_covmodel.py index 727ea64ae..15123af20 100644 --- a/tests/test_covmodel.py +++ b/tests/test_covmodel.py @@ -304,6 +304,7 @@ def test_covmodel_class(self): # checking some properties model_par = Stable() + self.assertFalse(model_par.do_rotation) self.assertEqual(len(model_par.arg), len(model_par.arg_list)) self.assertEqual(len(model_par.iso_arg), len(model_par.iso_arg_list)) self.assertEqual(len(model_par.arg), len(model_par.iso_arg) + 2) @@ -344,6 +345,7 @@ def test_covmodel_class(self): with self.assertWarns(AttributeWarning): Gau_fix(dim=1) self.assertRaises(ValueError, Gaussian, dim=0) + self.assertRaises(ValueError, Gau_fix, latlon=True) # check inputs self.assertRaises(ValueError, model_std.percentile_scale, per=-1.0) self.assertRaises(ValueError, Gaussian, anis=-1.0) From a470040b23e361251b360bb7a1f89fc0f5cefbe8 Mon Sep 17 00:00:00 2001 From: LSchueler Date: Wed, 10 Mar 2021 21:34:20 +0100 Subject: [PATCH 92/97] Fix typos in docstrings and comments --- examples/05_kriging/README.rst | 4 ++-- examples/10_normalizer/01_auto_fit.py | 8 ++++---- gstools/krige/base.py | 2 +- gstools/normalizer/base.py | 2 +- gstools/tools/misc.py | 2 +- gstools/variogram/variogram.py | 4 ++-- 6 files changed, 11 insertions(+), 11 deletions(-) diff --git a/examples/05_kriging/README.rst b/examples/05_kriging/README.rst index bfbc7e461..ef92e425f 100644 --- a/examples/05_kriging/README.rst +++ b/examples/05_kriging/README.rst @@ -34,14 +34,14 @@ the features you want: the given functions, the external drift has to given for each point form an "external" source. This results in :any:`ExtDrift` kriging. * `trend`, `mean`, `normalizer`: These are used to pre- and post-process data. - If you already have fitted a trend model, that is provided as a callable function, + If you already have fitted a trend model that is provided as a callable function, you can give it to the kriging routine. Normalizer are power-transformations to gain normality. `mean` behaves similar to `trend` but is applied at another position: 1. conditioning data is de-trended (substracting trend) 2. detrended conditioning data is then normalized (in order to follow a normal distribution) - 3. normalized conditioning data is set to zero mean (subsracting mean) + 3. normalized conditioning data is set to zero mean (subtracting mean) Cosequently, when there is no normalizer given, trend and mean are the same thing and only one should be used. diff --git a/examples/10_normalizer/01_auto_fit.py b/examples/10_normalizer/01_auto_fit.py index 78138565d..e95534056 100644 --- a/examples/10_normalizer/01_auto_fit.py +++ b/examples/10_normalizer/01_auto_fit.py @@ -3,7 +3,7 @@ ----------------- In order to demonstrate how to automatically fit normalizer and variograms, -we generate synthetic lognormal data, that should be interpolated with +we generate synthetic log-normal data, that should be interpolated with ordinary kriging. Normalizers are fitted by minimizing the likelihood function and variograms @@ -14,7 +14,7 @@ Artificial data ^^^^^^^^^^^^^^^ -Here we generate log-normal data following a gaussian covariance model. +Here we generate log-normal data following a Gaussian covariance model. We will generate the "original" field on a 60x60 mesh, from which we will take samples in order to pretend a situation of data-scarcity. """ @@ -44,7 +44,7 @@ # Fitting and Interpolation # ^^^^^^^^^^^^^^^^^^^^^^^^^ # -# Now we want to interpolate the "meassured" samples +# Now we want to interpolate the "measured" samples # and we want to normalize the given data with the BoxCox transformation. # # Here we set up the kriging routine and use a :any:`Stable` model, that should @@ -76,7 +76,7 @@ # As we see, it went quite well. Variance is a bit underestimated, but # length scale and nugget are good. The shape parameter of the stable model # is correctly estimated to be close to `2`, -# so we result in a gaussian like model. +# so we result in a Gaussian like model. # # The BoxCox parameter `lmbda` was estimated to be almost 0, which means, # the log-normal distribution was correctly fitted. diff --git a/gstools/krige/base.py b/gstools/krige/base.py index 1124c7c1e..56a1bd493 100755 --- a/gstools/krige/base.py +++ b/gstools/krige/base.py @@ -221,7 +221,7 @@ def __call__( chunk_size = pnt_cnt if chunk_size is None else int(chunk_size) chunk_no = int(np.ceil(pnt_cnt / chunk_size)) ext_drift = self._pre_ext_drift(pnt_cnt, ext_drift) - # iterate of chunks + # iterate chunks for i in range(chunk_no): # get chunk slice for actual chunk chunk_slice = ( diff --git a/gstools/normalizer/base.py b/gstools/normalizer/base.py index debe4ae23..8832fc595 100644 --- a/gstools/normalizer/base.py +++ b/gstools/normalizer/base.py @@ -237,7 +237,7 @@ def __eq__(self, other): # check for correct base class if type(self) is not type(other): return False - # if base class is same, this is save + # if base class is same, this is safe for val in self.default_parameter: if not np.isclose(getattr(self, val), getattr(other, val)): return False diff --git a/gstools/tools/misc.py b/gstools/tools/misc.py index 72818f8c7..aca1d276c 100755 --- a/gstools/tools/misc.py +++ b/gstools/tools/misc.py @@ -93,7 +93,7 @@ def _func_from_single_val(value, dim=None, value_type="scalar"): ) def _f(*pos): - # zip uses shortes len of iterables given (correct for scalar value) + # zip uses shortest len of iterables given (correct for scalar value) return np.concatenate( [ np.full_like(p, val, dtype=np.double) diff --git a/gstools/variogram/variogram.py b/gstools/variogram/variogram.py index a06b6be5b..417d590a5 100644 --- a/gstools/variogram/variogram.py +++ b/gstools/variogram/variogram.py @@ -133,7 +133,7 @@ def vario_estimate( same points, with the same statistical properties. bin_edges : :class:`numpy.ndarray`, optional the bins on which the variogram will be calculated. - If None are given, standard bins provided by the :any:`standard_bins` + If :any:`None` are given, standard bins provided by the :any:`standard_bins` routine will be used. Default: :any:`None` sampling_size : :class:`int` or :any:`None`, optional for large input data, this method can take a long @@ -205,7 +205,7 @@ def vario_estimate( Default: False mean : :class:`float`, optional mean value used to shift normalized input data. - Could also be a callable. The default is None. + Can also be a callable. The default is None. normalizer : :any:`None` or :any:`Normalizer`, optional Normalizer to be applied to the input data to gain normality. The default is None. From 9a875fc8276eae68b285a30f7a937edbf948b33c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Sat, 13 Mar 2021 15:48:41 +0100 Subject: [PATCH 93/97] Normalizer: add _dx as class attribute --- gstools/normalizer/base.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/gstools/normalizer/base.py b/gstools/normalizer/base.py index 8832fc595..40e2006f2 100644 --- a/gstools/normalizer/base.py +++ b/gstools/normalizer/base.py @@ -34,6 +34,7 @@ class Normalizer: """:class:`tuple`: Valid range for input data.""" denormalize_range = (-np.inf, np.inf) """:class:`tuple`: Valid range for output/normal data.""" + _dx = 1e-6 # dx for numerical derivative def __init__(self, data=None, **parameter): # only use parameter, that have a provided default value @@ -54,7 +55,7 @@ def _normalize(self, data): return data def _derivative(self, data): - return spm.derivative(self._normalize, data, dx=1e-6) + return spm.derivative(self._normalize, data, dx=self._dx) def _loglikelihood(self, data): add = -0.5 * np.size(data) * (np.log(2 * np.pi) + 1) From 8ac2194a73105a7dbe6e75b83991f56fa9bc9f90 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Sat, 13 Mar 2021 15:54:24 +0100 Subject: [PATCH 94/97] Field: rename mesh_call to generate_on_mesh; minor doc fixes --- gstools/field/base.py | 9 +++++++-- gstools/field/tools.py | 4 +++- 2 files changed, 10 insertions(+), 3 deletions(-) diff --git a/gstools/field/base.py b/gstools/field/base.py index 12a6e5915..30224a794 100755 --- a/gstools/field/base.py +++ b/gstools/field/base.py @@ -14,7 +14,11 @@ import numpy as np from gstools.covmodel.base import CovModel from gstools.tools.geometric import format_struct_pos_dim, gen_mesh -from gstools.field.tools import mesh_call, to_vtk_helper, fmt_mean_norm_trend +from gstools.field.tools import ( + generate_on_mesh, + to_vtk_helper, + fmt_mean_norm_trend, +) from gstools.normalizer.tools import apply_mean_norm_trend, _check_normalizer __all__ = ["Field"] @@ -137,7 +141,7 @@ def mesh( See: :any:`Field.__call__` """ - return mesh_call(self, mesh, points, direction, name, **kwargs) + return generate_on_mesh(self, mesh, points, direction, name, **kwargs) def pre_pos(self, pos, mesh_type="unstructured"): """ @@ -369,6 +373,7 @@ def name(self): return self.__class__.__name__ def _fmt_mean_norm_trend(self): + # fmt_mean_norm_trend for all child classes return fmt_mean_norm_trend(self) def __repr__(self): diff --git a/gstools/field/tools.py b/gstools/field/tools.py index 7f8757b92..0e2f5424e 100644 --- a/gstools/field/tools.py +++ b/gstools/field/tools.py @@ -102,13 +102,15 @@ def to_vtk_helper( ) -def mesh_call( +def generate_on_mesh( f_cls, mesh, points="centroids", direction="all", name="field", **kwargs ): """Generate a field on a given meshio, ogs5py or pyvista mesh. Parameters ---------- + f_cls : :any:`Field` + The field class in use. mesh : meshio.Mesh or ogs5py.MSH or PyVista mesh The given meshio, ogs5py, or PyVista mesh points : :class:`str`, optional From 92600ab2d63904ece1c5aebc69443ffc91faf909 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Sat, 13 Mar 2021 17:03:32 +0100 Subject: [PATCH 95/97] BoxCox: (bugfix) use correct data range for lmbda<0 --- gstools/normalizer/methods.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gstools/normalizer/methods.py b/gstools/normalizer/methods.py index 1a08bbf7c..7ff19d03c 100644 --- a/gstools/normalizer/methods.py +++ b/gstools/normalizer/methods.py @@ -83,7 +83,7 @@ def denormalize_range(self): if np.isclose(self.lmbda, 0): return (-np.inf, np.inf) if self.lmbda < 0: - return (-np.inf, np.divide(1, self.lmbda)) + return (-np.inf, -np.divide(1, self.lmbda)) return (-np.divide(1, self.lmbda), np.inf) def _denormalize(self, data): From bf374718d528aaabdf4a4423a38f0e2a7c537582 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Sat, 13 Mar 2021 17:08:05 +0100 Subject: [PATCH 96/97] BoxCoxShift: (bugfix) use correct data range for lmbda<0 --- gstools/normalizer/methods.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gstools/normalizer/methods.py b/gstools/normalizer/methods.py index 7ff19d03c..544f7c71a 100644 --- a/gstools/normalizer/methods.py +++ b/gstools/normalizer/methods.py @@ -157,7 +157,7 @@ def denormalize_range(self): if np.isclose(self.lmbda, 0): return (-np.inf, np.inf) if self.lmbda < 0: - return (-np.inf, np.divide(1, self.lmbda)) + return (-np.inf, -np.divide(1, self.lmbda)) return (-np.divide(1, self.lmbda), np.inf) def _denormalize(self, data): From 4b12aed20136241a4ea43c8b977769d273b02cee Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Sat, 13 Mar 2021 17:08:44 +0100 Subject: [PATCH 97/97] Examples: change field size in Auto-fit example to prevent confusion with sample-size --- examples/10_normalizer/01_auto_fit.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/10_normalizer/01_auto_fit.py b/examples/10_normalizer/01_auto_fit.py index e95534056..7654f5a26 100644 --- a/examples/10_normalizer/01_auto_fit.py +++ b/examples/10_normalizer/01_auto_fit.py @@ -22,8 +22,8 @@ import gstools as gs import matplotlib.pyplot as plt -# structured field with a size of 60x60 -x = y = range(60) +# structured field with edge length of 50 +x = y = range(51) pos = gs.tools.geometric.gen_mesh([x, y]) model = gs.Gaussian(dim=2, var=1, len_scale=10) srf = gs.SRF(model, seed=20170519, normalizer=gs.normalizer.LogNormal()) @@ -94,9 +94,9 @@ # of the measurement samples and the field is reconstructed quite accurately. fig, ax = plt.subplots(1, 3, figsize=[8, 3]) -ax[0].imshow(srf.field.reshape(60, 60).T, origin="lower") +ax[0].imshow(srf.field.reshape(len(x), len(y)).T, origin="lower") ax[1].scatter(*cond_pos, c=cond_val) -ax[2].imshow(krige.field.reshape(60, 60).T, origin="lower") +ax[2].imshow(krige.field.reshape(len(x), len(y)).T, origin="lower") # titles ax[0].set_title("original field") ax[1].set_title("sampled field")