Merge pull request #124 from GeoStat-Framework/normalizer

Mean, Normalizer, Trend and Fitting

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
 
 
@@ -83,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.
@@ -219,15 +223,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()
@@ -363,6 +367,7 @@ You can contact us via <[email protected]>.
 [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
+[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/

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",

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

docs/source/index.rst

@@ -6,9 +6,18 @@ 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,10 +108,9 @@ showing the most important use cases of GSTools, which are
 - `Field transformations <examples/07_transformations/index.html>`__
 - `Geographic Coordinates <examples/08_geo_coordinates/index.html>`__
 - `Spatio-Temporal Modelling <examples/09_spatio_temporal/index.html>`__
+- `Normalizing Data <examples/10_normalizer/index.html>`__
 - `Miscellaneous examples <examples/00_misc/index.html>`__
 
-Some more examples are provided in the examples folder.
-
 
 Spatial Random Field Generation
 ===============================
@@ -253,15 +261,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()

docs/source/normalizer.rst

@@ -0,0 +1,8 @@
+gstools.normalizer
+==================
+
+.. automodule:: gstools.normalizer
+
+.. raw:: latex
+
+    \clearpage

docs/source/package.rst

@@ -18,3 +18,4 @@ GSTools API
    random.rst
    tools.rst
    transform.rst
+   normalizer.rst

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

docs/source/variogram.rst

@@ -2,8 +2,6 @@ gstools.variogram
 =================
 
 .. automodule:: gstools.variogram
-   :members:
-   :undoc-members:
 
 .. raw:: latex
 

examples/03_variogram/06_auto_bin_latlon.py

@@ -72,5 +72,17 @@
 ###############################################################################
 # 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!

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

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")

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 (subtracting 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.

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
+)

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()