Geosci. Model Dev.

, 15, 3161–3182, 2022

https://doi.org/10.5194/gmd-15-3161-2022
© Author(s) 2022. This work is distributed under
the Creative Commons Attribution 4.0 License.

GSTools v1.3: a toolbox for geostatistical modelling in Python

Sebastian Müller1,2 , Lennart Schüler2,1,3 , Alraune Zech4,1 , and Falk Heße2,1
1 Department of Computational Hydrosystems, UFZ – Helmholtz Centre for Environmental Research, Leipzig, Germany
2 Instituteof Earth and Environmental Sciences, University of Potsdam, Potsdam, Germany
3 Center for Advanced Systems Understanding (CASUS), Görlitz, Germany
4 Department of Earth Sciences, Utrecht University, Utrecht, the Netherlands

Correspondence: Sebastian Müller (sebastian.mueller@ufz.de)

Received: 1 September 2021 – Discussion started: 19 October 2021

Revised: 27 January 2022 – Accepted: 3 February 2022 – Published: 12 April 2022

Abstract. Geostatistics as a subfield of statistics accounts for available to practitioners (Pyrcz and Deutsch, 2014; Rubin,
the spatial correlations encountered in many applications of, 2003; Diggle and Ribeiro, 2007; Kitanidis, 2008; Banerjee
for example, earth sciences. Valuable information can be ex- et al., 2014).
tracted from these correlations, also helping to address the Yet, the rate of adoption of geostatistics by practitioners
often encountered burden of data scarcity. Despite the value has been slow and uneven (Zhang and Zhang, 2004; Ra-
of additional data, the use of geostatistics still falls short of its jaram, 2016). One reason is the perceived lack of ready-made
potential. This problem is often connected to the lack of user- geostatistical software (Zhang and Zhang, 2004; Neuman,
friendly software hampering the use and application of geo- 2004; Winter, 2004; Rajaram, 2016; Cirpka and Valocchi,
statistics. We therefore present GSTools, a Python-based 2016; Fiori et al., 2016). Although a decent number of geo-
software suite for solving a wide range of geostatistical prob- statistical software solutions are available (Bellin and Rubin,
lems. We chose Python due to its unique balance between 1996; Deutsch and Journel, 1997; Brouste et al., 2008; Rubin
usability, flexibility, and efficiency and due to its adoption et al., 2010; Pebesma, 2004; Savoy et al., 2017; Heße et al.,
in the scientific community. GSTools provides methods for 2014; Vrugt, 2016), user-friendliness and licensing can ham-
generating random fields; it can perform kriging, variogram per their adoption as pointed out by Rubin et al. (2018).
estimation and much more. We demonstrate its abilities by Addressing these challenges, we present GSTools – an
virtue of a series of example applications detailing their use. extensive Python package for geostatistical analysis (Müller
and Schüler, 2021). To the best of our knowledge, no open-
source Python package that provides such a comprehensive
collection of random field generation, forward modelling,
1 Introduction kriging and data analysis currently exists.
We believe that the choice of Python has the potential to
Geostatistics emerged as a distinct branch of statistics in the address several of the challenges for geostatistical applica-
early 1950s through the pioneering work by Krige (1951). tions. First, a script language like Python allows for striking
Krige’s goal of estimating the abundance of mineral re- a balance between ease of use (as provided by graphical user
sources led him to develop some of the first methods, but it interfaces) and flexibility (as provided by command-line-
was the French mathematician Georges Matheron who devel- based tools). The presence of a graphical user interface (GUI)
oped the mathematical foundations (Matheron, 1962). Today, is sometimes seen as an indication of usability (Remy, 2005;
geostatistics is applied in fields like geology (Hohn, 1999), Rubin et al., 2018). However, a GUI does not necessarily
hydrogeology (Kitanidis, 2008), hydrology or soil sciences make software more user-friendly, almost always limits flex-
(Goovaerts, 1999), meteorology (Cecinati et al., 2017), ecol- ibility, increases the programming effort, and makes repro-
ogy (Rossi et al., 1992; Sales et al., 2007), oceanography ducibility of results and workflows harder or even unfeasible
(Monestiez et al., 2004), and epidemiology (Schüler et al., (Queiroz et al., 2017). Also note that in recent years pow-
2021), and a large number of textbooks make the theory

Published by Copernicus Publications on behalf of the European Geosciences Union.

3162 S. Müller et al.: GeoStatTools

erful tools emerged in data science to represent data along is the nugget or sub-scale variance, and cor(h) is the model-
with code snippets, documentation and results in interactive defining, normalized correlation function depending on the
computational notebooks provided by Jupyter (Perkel, 2018; non-dimensional distance h = s · r` .
Beg et al., 2021). Second, Python is known as a glue lan- The associated covariance and correlation functions are
guage, being able to combine independent software solutions given by the following:
to achieve complex workflows. This is particularly impor-  r
tant since geostatistics often relies on ready-made solvers C(r) = σ 2 · cor s · , (2)
for data generation or partial differential equation (PDE)-  r `
based model solvers. Third, Python is a simple yet powerful ρ(r) = cor s · . (3)
`
language with an increasing user base and community sup-
port for scientific computing and data analysis. It thus has Note that covariance and correlation are neglecting the
a wide appeal and excellent prospects for the foreseeable fu- nugget effect at the origin. Thus, the variance is interpreted as
ture. This guarantees that engineers and scientists with only a the variation above the nugget, which is sometimes referred
moderate background in computer science are able to apply to as the partial sill of the semi-variogram or the correlated
the toolbox and to make the necessary application-specific variability (Rubin, 2003). Consequently, the sill or limit of
adjustments. Finally, the licensing should be as permissive as the semi-variogram is calculated as the sum of variance and
possible to guarantee adoption and even further development nugget.
by interested users. The (semi-)variogram, covariance and correla-
We introduce GSTools and present its main features with tion functions of a model are accessible through
a general overview of its functionality and abilities in Sect. 2. model.variogram, model.covariance and
We focus on the covariance model, field generation kriging model.correlation, respectively. Every covari-
and variogram estimation. In Sect. 3, we discuss the wider ance model is defined by at least six parameters: dimension
context of GSTools, namely a number of Python packages dim, variance var, main length scale len_scale, rescale
connected with GSTools which can be used to seamlessly factor rescale, anisotropy ratios anis and rotation angles
model geostatistical workflows. Section 4 presents a num- angles, with the last two being dimension dependent. Fig-
ber of workflows to showcase the abilities of GSTools and ure 1 shows an example code for instantiating an exponential
demonstrate its usage. We close out with a short summary of model and the resulting model functions exemplifying the
the main advantages of GSTools and concluding remarks. parameters. Table 1 provides an overview of the predefined
models in GSTools.
In addition to the predefined covariance models, users can
2 GSTools features specify their own model functions by providing a normalized
correlation function. Figure 2 shows a reimplementation of
2.1 Covariance models and variography
the exponential model in only three lines of code.
The powerful CovModel class represents covariance and The dimension-dependent spectrum of an isotropic covari-
variogram models. Methods provided by this class are the ance model can be called with model.spectrum. It is di-
basis for most of the functionality of GSTools, such as var- rectly calculated from the covariance function by the follow-
iography, spatial random field generation and kriging. ing:
 d Z
2.1.1 Covariance models 1
S (k) = · C(|r|) · ei·k·r dr
2π
GSTools implements a CovModel class to define covari- Rd
ance models of weakly stationary (spatial) processes. Here, |k| n d o
= · H d −1 r 2 −1 C (·) (|k|). (4)
weak stationarity means that the associated semi-variogram (2π |k|) 2
d 2
is bounded, since we assume a constant mean and a finite
variance. To approximate unbounded variograms such as Here, r = |r| and k = |k| are the norms of the correspond-
the power-law model (Webster and Oliver, 2007), we pro- ing vectors, and H is the Hankel transform, which provides
vide a set of truncated power-law models following Di Fed- a mathematically self-contained and numerically robust for-
erico and Neuman (1997). The internal representation of a mulation of the radially symmetric Fourier transformation.
(semi-)variogram γ is given by GSTools makes use of an implementation of H provided
  r  by the Python package hankel (Murray and Poulin, 2019;
γ (r) = σ 2 · 1 − cor s · +n , (1) Ogata, 2005). For models with a known analytical solution,
`
GSTools uses them for improved computations.
where r is the (isotropic) lag distance, ` is the (main) cor- A prerequisite for kriging or random field genera-
relation length, s is a rescaling factor to adjust model repre- tion is that the applied covariance function is positive
sentation (default is 1), σ 2 is the variance or partial sill, n (semi-)definite. That can be checked through the spectral

Geosci. Model Dev., 15, 3161–3182, 2022 https://doi.org/10.5194/gmd-15-3161-2022

S. Müller et al.: GeoStatTools 3163

Figure 1. Initialization of an exponential covariance model given by cor(h) = exp(−h) (Rubin, 2003). Note that the rescaling factor is 1
by default. The right panel shows the plot of the variogram, covariance and correlation functions of the model, which can be created with
convenience plotting methods.

Table 1. Predefined covariance models in GSTools.

Model cor(h) Source

 
Gaussian exp −h2 Webster and Oliver (2007)

Exponential exp (−h) Webster and Oliver (2007)

exp −hα


Stable Wackernagel (2003)

Matern 21−ν · √ν · h
ν √
· Kν ν · h


Rasmussen and Williams (2005)
0(ν)
2 −α
 
Rational 1 + hα Rasmussen and Williams (2005)

Cubic (1 − 7h2 + 35 3 7 5 3 7
4 h − 2h + 4h ) (h<1) Chilès and Delfiner (2012)
Linear (1 − h) (h<1) Webster and Oliver (2007)
 p 
2 −1 2
Circular π · cos (h) − h · 1 − h (h<1) Webster and Oliver (2007)

Spherical (1 − 32 · h + 12 · h3 ) (h<1) Webster and Oliver (2007)

 !
1 d−1 3 2
2 F1 2 ,− 2 , 2 ,h
HyperSpherical 1−h· 
1 d−1 3
 (h<1) Matérn (1960)
2 F1 2 ,− 2 , 2 ,1
 !
1 3 2
2 F1 2 ,−ν, 2 ,h
SuperSpherical 1−h· 
1 3
 (h<1) Matérn (1960)
2 F1 2 ,−ν, 2 ,1
 −ν
JBessel 0(ν + 1) · h2 · Jν (h) Chilès and Delfiner (2012)

TPLSimple (1 − h)ν (h<1) Wendland (1995)

 
TPLGaussian H · E1+H h2 Di Federico and Neuman (1997)

TPLExponential 2H · E1+2H (h) Di Federico and Neuman (1997)

TPLStable 2H · E hα


Müller et al. (2021a)
α 1+ 2H
α

Formulas including the subscript (h < 1) are piecewise-defined functions being constantly zero for h ≥ 1. Here, h is the
non-dimensional distance, d is the dimension, 0(x) is the gamma function, Kν (x) is the modified Bessel function of the second kind,
Jν (x) is the Bessel function of the first kind, 2 F1 (a, b, c, x) is the ordinary hyper-geometric function and Eν (x) is the exponential
integral function (Abramowitz et al., 1972). All other variables are shape parameters of the respective models.

https://doi.org/10.5194/gmd-15-3161-2022 Geosci. Model Dev., 15, 3161–3182, 2022

3164 S. Müller et al.: GeoStatTools

Figure 2. Initialization of a user-defined exponential covariance model. The only thing that needs to be defined is the normalized correlation
function cor.

Figure 3. Spatial covariance structure of an anisotropic exponential model in 3D plotted with the built-in interactive routines of GSTools.
The example shows an eighth turn on the xy plane with anisotropy factors (1/2, 1/4). Rotation angles are given in radians.

density which is derived by the following (note that S only 2003):

depends on the norm of k): v v
u d  2 u d  2
uX r i su X ri r̃
S(k) n d o h=t = t =s· , (6)
− d2 2 −1 ρ (·) (k) . `i ` i=1 ei `
E(k) = = k(2π k) · H d r (5) i=1
σ2 2 −1

where ` = s · `1 is the main length scale incorporating the

From Bochner’s theorem (Rudin, 1990), it follows that rescale factor s, ei = ``1i represents the anisotropy ratios and
the spectral density is a probability density function if r = (r1 , r2 , . . .) represents the distances along the main axis
and only if the underlying covariance function is positive of correlation resulting in the isotropic distance r̃. Con-
(semi-)definite, which all predefined models in GSTools sequently, GSTools uses a main length scale, a set of
satisfy. As a consequence, the error variance during kriging anisotropy ratios and a set of rotation angles to fully describe
is always non-negative. an anisotropic model.
In practice, the main directions of correlation do not neces-
2.1.2 Anisotropy and rotation sarily follow the principal axis. The CovModel accounts for
rotation through rotation angles, where their number m de-
Variograms are typically defined based on the lag distance r, pends on the dimension d: m = d·(d−1) 2 . In two dimensions,
resulting in an isotropic model. However, many natural pro- rotation is fully described by a single angle for rotation in the
cesses involve anisotropy with varying correlation ranges in xy plane, and in three dimensions three angles are applied to
different (orthogonal) directions. An example is hydraulic the xy plane, xz plane and yz plane, respectively. In three
conductivity, where anisotropy typically arises from the ge- dimensions these are often referred to as Tait–Bryan angles
ologic stratification. The implementation of anisotropy in yaw, pitch and roll (Goldstein, 1980); see Fig. 3 for an exam-
GSTools is based on the non-dimensional distance (Rubin, ple.

Geosci. Model Dev., 15, 3161–3182, 2022 https://doi.org/10.5194/gmd-15-3161-2022

S. Müller et al.: GeoStatTools 3165

Figure 4. Initialization of a Yadrenko covariance model. We use the Earth’s radius as the rescaling factor to have a meaningful length scale.
The routine vario_yadrenko still depends on the central angle given in radians.

Figure 5. Estimating an empirical variogram of synthetic unstructured data and fitting an exponential model. The number of bins was
calculated to be 21 with a maximum bin distance of ca. 45.

One unique feature of GSTools is the support of arbi- I3 = [2, 3] (yz plane), etc. These values define a rotation ma-
trary dimensionality in all provided routines. For rotation trix Rot to transform principal axes to the directions of cor-
in higher dimensions, we apply the following scheme: the relation and the back-rotation matrix DeRot = Rot−1 for the
first angles coincide with those of the next lower dimension, inverse:
and the added d − 1 angles describe rotations in the planes x
of the added dimension (in 3D: xz and yz). Thus, there are m
Y
6 rotation angles in 4D, 10 in 5D, etc. Rotation in higher Rot = G((−1)i−1 αi , Ii )
dimensions is only relevant for spatio-temporal modelling i=1
with three spatial dimensions and application to other fields = G((−1)m−1 αm , Im ) · . . . · G(α1 , I1 ). (8)
of research with high-dimension data. The scheme was cho-
sen for metric spatio-temporal models to account for spatial The alternating signs of the rotation angles (−1)i−1 αi were
anisotropy in a similar way as a simple spatial model. chosen to match Tait–Bryan angles in 3D.
Rotation in the xi xj plane is described by the matrix For applying or removing anisotropy, we define the
G(α, [i, j ]) ∈ Rd×d . isotropify matrix Iso = diag(e1−1 , e2−1 , . . .) and anisotropify
 matrix AnIso = Iso−1 . Combining these two types of matri-

 cos θ k = l = i, j ces allows us to isometrize (i.e. isotropify and derotate) and


sin θ anisometrize (i.e. rotate and anisotropify) spatial points via
 k = i, l = j
 the following:
G(α, [i, j ])kl = − sin θ k = j, l = i (7)

1 k = l 6 = i, j Isom = Iso · DeRot, (9)





0 else AnIsom = Rot · AnIso. (10)


The order of rotating planes is determined by the described GSTools provides the routine CovModel.isometrize
scheme, i.e. I1 = [1, 2] (xy plane), I2 = [1, 3], (xz plane) to convert spatial positions to their derotated and isotropic

https://doi.org/10.5194/gmd-15-3161-2022 Geosci. Model Dev., 15, 3161–3182, 2022

3166 S. Müller et al.: GeoStatTools

Figure 6. Estimation of directional variograms for given main axes: the code snippet shows the setup for estimating and fitting the variogram
to an anisotropic field. The graphs show the main axes of the rotated model and the fitting results. Plotting commands have been omitted.

Figure 7. Estimating an empirical variogram (bottom left) of synthetic unstructured data (top left) after Box–Cox normalization of skewed
input values. Panels on the right show the histogram of the data values before (top) and after (bottom) the normalization. For demonstration
purposes, a Matern model was fitted to the empirical variogram. Plotting commands have been omitted.

counterparts as required by Eq. (6) and the routine p1 = (φ1 , λ1 ) and p2 = (φ2 , λ2 ) is given by their latitude (φ)
CovModel.anisometrize to invert this: and longitude (λ) and can be described by a central angle
calculated from the great circle distance:
x isom = Isom · x, (11)
x anisom = AnIsom · x. (12) ζ (p1 , p2 ) = arccos (sin φ1 sin φ2 + cos φ1 cos φ2 cos(1λ)) (13)
A huge family of valid models on the sphere can be derived
from 3D models by inserting the chordal distance which re-
2.1.3 Geographical coordinates sults in the associated Yadrenko covariance model CY (Lan-
tuéjoul et al., 2019):
Earth’s surface is a non-Euclidean manifold, and all large-   
ζ
scale, geographically referenced data will necessarily reflect CY (ζ ) = C3D 2 · sin . (14)
that. We deal with the non-Euclidean nature of this kind 2
of data by assuming the Earth to be a perfect sphere and The underlying manifold introduces new restrictions for
then using the fact that the distance between two points covariance models to be positive definite. The manifold

Geosci. Model Dev., 15, 3161–3182, 2022 https://doi.org/10.5194/gmd-15-3161-2022

S. Müller et al.: GeoStatTools 3167

Both estimators require predefined bins M(r) to group the

pairwise point distances of the given field. GSTools pro-
vides a standard binning routine, where the maximal bin
width is set to one-third of the diameter of the containing
box of the field, the number of bins is determined by Sturges’
rule (Sturges, 1926) and all bins have equal width. Figure 5
provides an example of the variogram estimation of an un-
structured spatial random field with automatic binning.
GSTools accounts for anisotropy by providing estimat-
ing routines for directional variograms along a specified di-
rection with a certain angle tolerance and bandwidth. When
providing orthogonal axes, it is possible to fit a theoretical
model and its anisotropy ratios as shown in Fig. 6. Determin-
ing the main rotation axes from given data, however, is up
to the user and beyond the scope of the presented GSTools
version.
Field data often do not follow a normal distribution,
which is a crucial assumption for variogram estimation.
For example, transmissivity is usually assumed to be log-
normally distributed (Dagan, 1989), while rainfall data are
normalized by applying the Box–Cox transformation (Ce-
cinati et al., 2017). GSTools provides a set of normal-
izers based on power transforms, which can be fitted to a
given data set using a maximum likelihood approach (Elia-
Figure 8. Comparison of parametric normalizers in GSTools.
son, 1993): LogNormal, BoxCox (Box and Cox, 1964),
YeoJohnson (Yeo and Johnson, 2000), Modulus (John
structure of the sphere only allows for isotropic models. For and Draper, 1980) and Manly (Manly, 1976). An example
small-scale applications, it is valid to assume anisotropy. An application is shown in Fig. 7, and a comparison of all pro-
appropriate adaption is the use of a 2D projection like Gauss– vided normalizers can be seen in Fig. 8.
Krüger coordinates. We provide Yadrenko models as a uni- GSTools also provides routines to detrend data. For ex-
fied representation for non-Euclidian coordinates since they ample, temperature could decrease with elevation or con-
facilitate all presented models to be used with geographical ductivity could decrease with depth. Another application is
coordinates as demonstrated in Fig. 4. analysing spatial correlation of residuals after application of
a regression model to the data. All routines dealing with data
2.1.4 Empirical variogram, data preparation and have the keywords trend, normalizer and mean, where
model fitting the last keyword describes the mean of the normalized data.

The empirical variogram is an important tool for analysing 2.2 Kriging, random fields and conditioned random
spatially correlated data. It is estimated with the subpackage fields
gstools.variogram, which provides two estimators for
the empirical variogram: Matheron and Cressie (Webster and 2.2.1 Kriging
Oliver, 2007). The default Matheron’s estimator for a vari-
ogram γ of a spatial random field U is given by The subpackage gstools.krige provides routines for
Gaussian process regression, also known as kriging, which
1 X
2
γ (r) = · |M(r)|−1 U (x i ) − U (x j ) , (15) is a method of data interpolation based on predefined covari-
2 M(r) ance models (Wackernagel, 2003). Kriging aims to derive the
value of a field z at some point x 0 , when there are fixed
where M(r) is the set of all pairwise spatial random field conditioning values z(x 1 ). . .z(x n ) at given points x 1 . . .x n .
points, separated by the distance r and a certain tolerance The resulting value z0 at x 0 is calculated as a weighted mean
ε > 0. z0 = ni=1 wi · zi , where the weights, w = (w1 , . . ., wn ), are
P
Cressie’s estimator, which is more robust to outliers, is determined by the specific kriging routine.
given by We provide multiple kriging routines derived from the
 4 Krige class. (i) Simple: the data are interpolated with a
1 −1
P p
2 · |M(r)| M(r) |U (x i ) − U (x j )| given mean value for the kriging field. (ii) Ordinary: the
γ (r) = . (16) mean of the resulting field is unknown and estimated along-
0.457 + 0.494/|M(r)| + 0.045/|M(r)|2

https://doi.org/10.5194/gmd-15-3161-2022 Geosci. Model Dev., 15, 3161–3182, 2022

3168 S. Müller et al.: GeoStatTools

Figure 9. Comparison of simple, ordinary and universal kriging. All three routines have a similar setup, where simple kriging needs an
estimated mean and where universal kriging needs additional drift functions. Plotting commands have been omitted.

Figure 10. A simple setup for linear regression kriging. Although the interpolation coincides with a piecewise linear function, we gain
information about the error variances between the conditioning points as shown in the right plot.

side the interpolation (unbiasedness). (iii) Universal: in ciated covariance function:

addition to ordinary kriging, one can provide drift functions
f1 , . . ., fk . (iv) ExtDrift: like universal kriging but the
drift is provided by an external source.
The advantage of using the general Krige class is the
combination of all described features, such as for instance
using universal kriging with a functional drift together with (17)
additional external drifts. A typical scenario is a temperature
interpolation with an assumed north–south drift (functional with C = (C(x i , x j ))ij =1...n being the covariance matrix, de-
drift) and a linear correlation to altitude (external drift). pending on the conditioning points and the given model.
Since all variogram models in GSTools assume weak C0 = (C(x i , x 0 ))Ti=1...n is the covariance vector for the target
stationarity, the kriging system is always built on the asso- point x 0 . F = (fj (x i ))i=1...n,j =1...k is a submatrix contain-
ing the functional drift values at the conditioning points and
F0 = (fi (x 0 ))Ti=1...k at the target point, where k is the num-
ber of functional drifts. E = (eij )i=1...n,j =1...l is a submatrix
containing the external drift values at the conditioning points

Geosci. Model Dev., 15, 3161–3182, 2022 https://doi.org/10.5194/gmd-15-3161-2022

S. Müller et al.: GeoStatTools 3169

Figure 11. Generation of a structured random field following a Gaussian variogram.

Figure 12. An example for an ensemble of 1D random fields conditioned to five measurements (dots). Plotting commands have been omitted.

and E0 = (ei0 )Ti=1...l at the target point, where l is the num- surements at the same location are averaged out in the result-
ber of external drifts. The parameters µ, φ = (φ1 , . . ., φk )T ing field (Mohammadi et al., 2017).
and ψ = (ψ1 , . . ., ψl )T are Lagrange multipliers for the un- One last feature is the capability of kriging the mean
biased condition, the functional drifts and the external drifts, (Wackernagel, 2003), which allows for deriving the mean
respectively. The vector 1 and its Lagrange multiplier µ are value estimated during ordinary kriging or estimating the
given in brackets since their appearance depends on whether mean drift determined from given functional and/or exter-
the system should be unbiased or not (ordinary vs. simple nal drift terms as shown in Fig. 9. A minimal example for
kriging). Note that the number of functional drifts k and ex- regression kriging is shown in Fig. 10.
ternal drifts l can be zero, depending on whether they are
given or not. 2.2.2 Random fields
GSTools also provides the possibility to incorporate
measurement errors variances σi2 for each conditioning point A core element of GSTools is the spatial random field gen-
by adjusting the covariance matrix (Wackernagel, 2003): erator class SRF. A covariance model (Sect. 2.1) is needed to
C̃ = C + diag(σ12 , . . ., σn2 ) instantiate a spatial random field. We provide two ways for
field generation: structured or unstructured. In both cases, the
C(x 1 , x 1 ) + σ12 . . .
 
C(x 1 , x n ) positions at which the field will be evaluated are given by a
.. .. .. (18)
= . pos argument. In the structured case, pos contains one tu-
 
. . .
C(x n , x 1 ) ... C(x n , x n ) + σn2 ple per dimension, each defining the subdivision of the cor-
responding axis resulting in a rectilinear grid. For unstruc-
By default, the measurement error variances, σi2 , are set to tured grids, the pos tuple contains the x, y and z coordi-
the model nugget. In order to get numerically stable results, nates of every evaluation point. SRF allows for controlling
we solve the kriging system with the pseudo-inverse matrix, the underlying pseudo-random number generation by a seed
which has the advantage that redundant data or multiple mea- to reproduce field generation. A code example is given in

https://doi.org/10.5194/gmd-15-3161-2022 Geosci. Model Dev., 15, 3161–3182, 2022

3170 S. Müller et al.: GeoStatTools

Figure 13. Generation of a structured incompressible random vector field with exponential variogram.

Figure 14. Example of a log-transformed binary field with the low values being connected by applying the zinnharvey, binary and
lognormal transformations successively.

Fig. 11. Field generation is performed through the random- flow allows users to generate a random field only from a
ization method (Kraichnan, 1970; Heße et al., 2014), which given correlation, covariance or variogram function.
utilizes the spectral density (Eq. 5) of the variogram model A key advantage of the randomization method implemen-
to approximate a Wiener process in Fourier space by tation is the possibility to extend a generated SRF seamlessly,
which not only preserves its statistical properties but also
s
N preserves the actual realization of the SRF. Potential appli-
σ2 X

cations are the following. (i) Particle simulations, where ran-
U (x) = · Z1,i · cos (k i · x) + Z2,i · sin (k i · x) , (19)
N i=1 dom incompressible velocity fields can be generated exactly
at the location of the individual particles (see workflow in
where N is the number of Fourier modes of the approxima- Sect. 2.3.1). It avoids interpolation errors, arising from grid-
tion. The random variables Z1,i , Z2,i ∼ N (0, 1) are mutually based velocity fields. (ii) If concentration plumes are sim-
independent and are drawn from a standard normal distribu- ulated on a large domain, the SRF can be calculated on de-
tion. The k i values are mutually independent random sam- mand only for the time dependent spatial extent of the plume.
ples, which are drawn from the spectral density with the aid And (iii) for high-performance computing applications, the
of emcee, a Python package for Markov chain Monte Carlo field generation can be directly coupled to the domain de-
sampling (Foreman-Mackey et al., 2013). composition, and each task only generates the SRF for its
The randomization method is implemented in the part of the domain. There are two main classes of alternative
RandMeth class and used by default. The RandMeth rou- methods to the randomization method (Heße et al., 2014).
tines create isotropic random fields. Thus, the corresponding By decomposing the covariance function, small spatial ran-
covariance is radially symmetric, and the spectral density can dom fields can be computed very fast. But the computational
be calculated by the Hankel transformation. Anisotropy is re- cost quickly becomes unfeasible as the field grows in size.
alized by rescaling and rotating the input points. The work- A second and quite popular class is the sequential Gaussian

Geosci. Model Dev., 15, 3161–3182, 2022 https://doi.org/10.5194/gmd-15-3161-2022

S. Müller et al.: GeoStatTools 3171

Figure 15. A workflow to generate a spatio-temporal random field with one spatial dimension.

method, which can also create conditioned spatial random method to generate a single random field, and (ii) we only
fields. However, numerical problems can arise if the underly- need to solve the kriging problem once and not sequentially.
ing correlation function is smooth at the origin, and the com- Figure 12 shows an example of an ensemble of condi-
putational costs also become unfeasible for highly resolved tioned random fields in one dimension. Where measurements
random fields (Emery, 2004). of the target variables are available, all realizations satisfy
Just like the kriging routines, the spatial random field gen- them. However, random fields behave as unconditional fields
erator allows for incorporating predefined trend, normalizer (i.e. of an ensemble with identical parameters, like mean,
and mean for a greater variety of distributions. A special SRF variance and correlation length) where no point measure-
class feature is the capability to perform variance upscaling ments are available (x > 6). Characteristics, such as the en-
to respect generation of random fields on mesh cells with a semble variance, significantly change given the distribution
certain volume. We hereby use the upscaling method coarse of measurements and conditioning. The ensemble mean and
graining (Attinger, 2003) to rescale the variance in Eq. (19) the kriging field coincide, proving that the kriging field is the
at each target point based on a given filter volume size λ: best linear unbiased predictor for the given data.
!d/2
`2 2.3 Additional features
2 2
σ (λ) = σ · , (20)
λ 2


`2 + 2 2.3.1 Incompressible random vector field generation
√
where ` is the correlation length, and λ = d V is the filter Kraichnan (1970) was the first to suggest a randomization
size derived from the cell volume V depending on the field method for studying the diffusion of single particles in a ran-
dimension, assuming the cell element is a hypercube. This dom incompressible velocity field. He came up with a ran-
approach was derived from the groundwater flow equation, domization method which includes a projector, ensuring the
assuming a Gaussian covariance model and should therefore incompressibility of the vector field.
be used with caution in differing scenarios. An example is When U is the mean velocity (oriented towards the first
provided in the workflow in Sect. 4.3. basis vector e1 ), we generate random fluctuations with a
given covariance model around U . And at the same time,
2.2.3 Conditioned random fields we ensure that the velocity field remains incompressible; that
is, ∇ ·U = 0 by using the randomization method (Eq. 19) and
When point measurements of a target variable are available, adding a projector p(k i ) to every mode being summed:
they need to be considered when generating random fields.
GSTools provides a class CondSRF combining kriging and s
N
random field generation, where we first derive the kriged σ2 X
field and the error variance and then generate a random field U(x) = U e1 − p(k i ) · Z1,i cos (k i · x)
N i=1
with zero mean where the variance in Eq. (19) is replaced

with the estimated error variance. This procedure is advan- +Z2,i sin (k i · x) , (21)
tageous to classical sequential Gaussian simulation (Webster ki1
and Oliver, 2007), as (i) we make use of the randomization p(k i ) = e1 − · ki . (22)
|k i |2

https://doi.org/10.5194/gmd-15-3161-2022 Geosci. Model Dev., 15, 3161–3182, 2022

3172 S. Müller et al.: GeoStatTools

Calculating ∇ ·U = 0 verifies that the resulting field is indeed GSTools provides an interface for a number of mesh stan-
incompressible. An example is shown in Fig. 13. Things like dards, such as meshio (Schlömer et al., 2021), PyVista
boundary conditions cannot be modelled with this method, (Sullivan and Kaszynski, 2019) and ogs5py (Müller et al.,
but it can be used, for example, in transport simulations of 2020). When using meshio or PyVista, the generated
the saturated subsurface (Schüler et al., 2016) or for studying fields will be stored immediately in the mesh container. There
turbulent open water (Kraichnan, 1970). are two options to generate a field on a given mesh: either on
the points (points="points") or on the cell centroids
2.3.2 Field transformations (points="centroids"), which is important depending
on the specification of the variable in the numerical scheme.
GSTools generates Gaussian random fields while real Figure 16 provides an example.
data often do not follow a Gaussian distribution. This is
typically addressed through data transformation. GSTools
provides a number of appropriate transformations beyond 3 GSTools within the ecosystem of the GeoStat
power transformations provided by the normalizer sub- Framework
module (Sect. 2.1.4): (i) binary, (ii) discrete, (iii)
boxcox (Box and Cox, 1964), (iv) zinnharvey (Zinn GSTools is part of a larger suite of Python packages,
and Harvey, 2003), (v) normal_force_moments, collectively hosted on GitHub under https://github.com/
(vi) normal_to_lognormal, (vii) GeoStat-Framework (last access: 31 March 2022). The other
normal_to_uniform, (viii) normal_to_arcsin and packages in the GeoStat Framework complement some of the
(ix) normal_to_uquad. abilities of GSTools and form a comprehensive framework
Transformations can be combined sequentially to create for geostatistical applications. We introduce some packages
more complex scenarios as in Fig. 14. Note that, in contrast and demonstrate how they interact with, enhance and lever-
to normalizers, transformations cannot be fitted to given data, age the abilities of GSTools.
which leaves the choice of the best transformation to the user.
3.1 ogs5py
2.3.3 Spatio-temporal modelling
ogs5py (Müller et al., 2020) provides a Python API for the
Spatio-temporal modelling provides insights into time- FEM-based OpenGeoSys 5 (Kolditz et al., 2012) scientific
dependent stochastic processes like rainfall, air temperature software suite for hydrogeological processes like ground-
or crop yield, which is of high practical relevance. GSTools water flow and transport modelling where data scarcity is
provides the metric spatio-temporal model (Cressie and a typical shortcoming. Examples are point measurements
Wikle, 2011) for all covariance models by enhancing the spa- of hydraulic head from observation wells and breakthrough
tial with a time dimension, resulting in the spatio-temporal curves from tracer experiments. Inferring hydraulic conduc-
dimension dst : tivity from these data requires a modelling framework with
integrated stochastic data generation. The combination of
vu d  2  2 
uX r i GSTools with ogs5py enables a user to integrate the geo-
1t 
Cm (r, 1t) = C t + statistical modelling of an aquifer with hydrogeological sim-
i=1
ei κ ulations. Such an example for pumping test simulations is
p  provided in Sect. 4.3.
=C r̃ 2 + 1t˜2 , (23)
3.2 welltestpy and AnaFlow
where r̃ is the isotropified spatial distance as defined in
Eq. (6), 1t is the temporal distance and 1t˜ the isotropified A pumping test is a cost-effective subsurface observation
temporal distance. The parameter κ accounts for a spatio- method typically used by hydrogeologists for aquifer charac-
temporal anisotropy ratio and is the last entry of the anis terization. The package welltestpy (Müller et al., 2021b)
array appended to the spatial anisotropy ratios. The imple- provides tools to handle, process, plot and analyse data from
mentation in GSTools enables the direct incorporation of pumping test campaigns. It assists practitioners in identifying
spatial anisotropy and rotation in a spatio-temporal model. hydrogeological parameters by fitting measured drawdowns
It further supports the use of arbitrary spatial dimensions in to some conceptual flow model. The package contains a num-
spatio-temporal models. Figure 15 shows the generation of a ber of examples that illustrate these abilities.
spatio-temporal random field with one spatial dimension. AnaFlow (Müller et al., 2021a) provides a wide range
of analytical expressions for pumping tests under various
2.3.4 Working on meshes conditions. Classical examples are Thiem’s and Theis’ so-
lution assuming homogeneous aquifer properties. In addi-
For improved handling of spatial random fields as input tion, AnaFlow provides extended versions of both solutions,
for PDE solvers like the finite element method (FEM), which account for aquifer heterogeneity and allow for esti-

Geosci. Model Dev., 15, 3161–3182, 2022 https://doi.org/10.5194/gmd-15-3161-2022

S. Müller et al.: GeoStatTools 3173

Figure 16. Generating spatial random fields on finite element method (FEM) meshes: either on cell centroids (middle) or mesh points (right).
Plotting commands have been omitted.

Figure 17. Estimating the temperature variogram with geographic coordinates using the spherical Yadrenko model. Estimated length scale
is ca. 0.9 (radians) and sill is around 13.

Figure 18. Setting up universal kriging with a drift function. Figure 19. Workflow for regression kriging with a linear regression
model.

mating higher-order geostatistical parameters like variance become the kriging toolbox for the GeoStat Framework, pro-
and correlation length (Zech et al., 2012, 2016). viding advanced methods.

3.3 PyKrige 3.4 Development, documentation and installation

GSTools provides an interface to the stand-alone package GSTools is compatible with Python versions ≥ 3.6
PyKrige (Murphy et al., 2021) for more specialized krig- , although previous releases support older versions of
ing applications. After 10 years of independent development, Python. Performance-critical parts, like variogram estima-
PyKrige has recently been migrated to the GeoStat Frame- tion (Sect. 2.1.4), kriging summation (Sect. 2.2.1) and the
work, and its functionality is currently integrated with the summation of the randomization method (Sect. 2.2.2), are
other packages. So far, the covariance model can be ex- implemented in Cython (Behnel et al., 2011). GSTools
changed between the packages. In the future, PyKrige will mainly depends on the SciPy ecosystem with its mandatory

https://doi.org/10.5194/gmd-15-3161-2022 Geosci. Model Dev., 15, 3161–3182, 2022

3174 S. Müller et al.: GeoStatTools

Figure 20. Plot of temperature measurements (a), universal kriging interpolation (b) and regression kriging results (c).

by Sphinx. Continuous integration is established through

GitHub actions where Python wheels are pre-built for the
most common operating systems (Windows, Linux, MacOS)
and Python versions to enable simple installation. Each re-
lease on GitHub is directly deployed to the Python package
index (https://www.pypi.org; last access: 31 March 2022) as
well as conda-forge (conda-forge community, 2015). An ex-
tensive set of unit tests are performed automatically and con-
tinuously through GitHub actions.

3.5 Interoperability

To integrate GSTools into the Scientific Python Stack,

we provide a set of interfaces to other packages.
Figure 21. Scatterplot of latitude–temperature values (grey dots), These include the above-mentioned packages ogs5py,
the linear regression result (dashed orange line), universal kriging
meshio, PyVista and pyevtk (https://github.com/
mean drift (dashed blue line) and the cross-sections of the respective
pyscience-projects/pyevtk, last access: 31 March 2022
kriging interpolation (solid lines).
) for mesh operations. Other packages for geostatistics
are also supported, such as PyKrige (Sect. 3.3) and
scikit-gstat (Mälicke, 2022), the latter having a fo-
dependencies numpy (Harris et al., 2020) and scipy (Vir- cus on variography and can be used for more detailed var-
tanen et al., 2020). The source code is maintained under a iogram estimation. For both packages, interfaces are pro-
GitHub organization for optimizing team efforts. Users have vided to convert covariance models of GSTools to or from
the opportunity to communicate with developers by asking their counterparts in the respective package. Another pack-
questions in a discussions forum, raising issues or improv- age worth mentioning is verde (Uieda, 2018), a Python li-
ing code by making pull requests. All packages come with brary for processing and gridding spatial data. Some of the
detailed documentation via https://readthedocs.org (last ac- features provided there can be easily combined with capa-
cess: 31 March 2022), which contains a range of tutorials ex- bilities of GSTools such as detrending data to preprocess
plaining the features and a full API documentation created inputs.

Geosci. Model Dev., 15, 3161–3182, 2022 https://doi.org/10.5194/gmd-15-3161-2022

S. Müller et al.: GeoStatTools 3175

Figure 22. Initialization of an OGS model with mesh generation.

Figure 23. Generating correlated log-normal SRFs adapted to the

mesh settings of the numerical model for three connectivity struc-
tures following the Zinn and Harvey (2003) transformation.

Figure 24. Tracer transport simulation results: spatially distributed

4 Workflows concentration plumes after 15 d with transmissivity distributions in
the background.
Having explained the core features of GSTools, we now
provide a couple of example applications covering the topic
of kriging, variogram estimation, random field generation
and coupling with other tools to achieve more elaborate trend model is used as the internal drift in the kriging sys-
workflows. The examples illustrate the abilities of GSTools tem. The methods differ in use of the covariance model. The
and serve as a starting point for a user’s project development. linear RK does not incorporate spatial correlation informa-
All shown code snippets are taken from the actual workflow tion, while UK does through the drift function for calculat-
scripts and are not self-contained. ing the kriging weights. Both methods are often considered to
provide mathematically equal results, but we show that there
4.1 Regression kriging vs. universal kriging: finding a are sensitive differences. The resources for this workflow are
north–south temperature trend provided in Müller (2021).
As a data basis, we use measured temperature of the
Kriging is a well-established interpolation method applied German weather service retrieved with the Python package
in many fields of natural science. We compare two options wetterdienst (Gutzmann et al., 2021), which we exam-
of incorporating auxiliary variables to calculate the kriging ine for a linear north–south trend. We use the established
weights: (i) regression kriging (RK), where the trend of input spherical covariance model in its Yadrenko variant suitable
data is estimated by regression and simple kriging is applied for geographical coordinates. Variogram estimation and fit-
to the residuals, and (ii) universal kriging (UK), where the ting results are shown in Fig. 17.

https://doi.org/10.5194/gmd-15-3161-2022 Geosci. Model Dev., 15, 3161–3182, 2022

3176 S. Müller et al.: GeoStatTools

Figure 25. Workflow to generate an ensemble of transmissivity fields on a given mesh (a). A single realization in shown in the right plot (b).

Figure 26. Comparison the ensemble mean drawdown hh(r, t)i (a, d) with the effective head solution hCG (r, t) (c, f) for two parameter sets.
The vanishing absolute difference between both (b, e) shows that they perfectly agree.

Figure 18 shows how to set up the UK estimator, including Figure 21 shows the estimated mean trends for both UK
the drift function, and Fig. 19 shows the setup of the RK es- and RK, revealing completely contrary results. The RK re-
timator. RK requires the preceding step of fitting the regres- sult indicates an increase in mean temperature with increas-
sion model for the trend of the Detrended kriging routine. ing latitude, which seems reasonable given a raising terrain
The interpolation results are shown in Fig. 20, indicating that elevation from the Baltic Sea in the north towards the Alps
both methods provide equally good results. in the south. The estimated mean of UK shows the opposite

Geosci. Model Dev., 15, 3161–3182, 2022 https://doi.org/10.5194/gmd-15-3161-2022

S. Müller et al.: GeoStatTools 3177

Figure 27. Transmissivity of the Herten aquifer analogue and locations of virtual observations marked as black dots.

Figure 28. Variogram estimation and resulting experimental (dots) and fitted variogram γ (line) of the Herten aquifer analogue.

with temperature decreasing with latitude. A potential expla- model.msh, mesh details are transferred for generating dis-
nation here is the general temperature increase towards the tributed values at particular mesh locations, even for unstruc-
Equator. While the UK mean fits better with the cross-section tured grids. The subroutine transform.zinnharvey al-
at 10◦ longitude (Fig. 21), the RK mean fits the scatter dia- lows for generating Gaussian structures where the mean val-
gram better, as expected. ues of the field are not connected but the low- or high-
conductivity areas are connected, using the transformation by
4.2 Heterogeneous transport simulation: the impact of Zinn and Harvey (2003). Note that the correlation lengths un-
connectivity dergo rescaling (Gong et al., 2013). The concept of connec-
tivity follows the paper by Zinn and Harvey (2003), where
The combination of ogs5py and GSTools makes it possi- connectedness refers to connected paths of extreme or spe-
ble to quickly set up and run subsurface flow and transport cial values in the conductivity field.
simulations in a heterogeneous aquifer setting. The critical Simulated tracer plumes in Fig. 24 show the particular ef-
step is the generation of a spatially distributed hydraulic con- fects of connectivity: the plume remains relatively compact
ductivity distribution, adapted to the numerical simulation for classical Gaussian fields, where mean values are con-
grid. We further demonstrate GSTools’ ability to generate nected. Transformed fields lead to more disrupted, dynamic
different connectivity structures, and we discuss their impact plumes, which is mostly caused by trapping in areas of con-
on transport results. The resources for this workflow are pro- nected low conductivity and preferential flow in connected
vided in Müller and Zech (2021a). high-conductivity areas.
A flow and transport model is initialized through an in-
stance of the OGS class from ogs5py, with simple mesh
generation (Fig. 22) and specification of model parame-
ters and boundary conditions. Random fields are initialized
through the SRF class (Fig. 23). By passing the subclass

https://doi.org/10.5194/gmd-15-3161-2022 Geosci. Model Dev., 15, 3161–3182, 2022

3178 S. Müller et al.: GeoStatTools

Figure 29. One realization of the conditioned SRF (a) and absolute difference (b) between the “true” (Fig. 27) and conditioned transmissivity,
showing increasing differences with distance from the conditioning area (rectangle).

Figure 30. Generation of an ensemble of 20 conditional realizations and transmissivity transects T (x) at y = 4 m. The thick blue line is the
true transmissivity (Fig. 27), and the shaded area shows the range of 1 standard deviation calculated from 20 realizations of conditioned
fields. Black points indicate the observations.

4.3 Characterizing mean drawdowns of a pumping test the coarse graining procedure for Gaussian variograms
ensemble according to Eq. (20).
Calculated ensemble means can be compared to analyt-
Combining flow simulations in ogs5py with random fields ical solutions (Fig. 26), such as Theis’ solution for ho-
of GSTools allows for performing Monte Carlo studies to mogeneous media or the effective drawdown solution by
identify ensemble mean behaviour. Zech et al. (2016) made Zech et al. (2016), making use of their implementations in
use of this workflow to prove the applicability of an effec- welltestpy and AnaFlow.
tive drawdown solution for pumping tests in random conduc-
tivity. We present a short form of their workflow, which is 4.4 Geostatistical exercises with the Herten aquifer
accessible in Müller and Zech (2021b).
The flow model is initialized through the OGS class We demonstrate how to estimate variograms and how to con-
with model parameters and boundary conditions creat- dition spatial random fields on observations using data from
ing the convergent flow setting of a pumping test. The the Herten aquifer analogue (Bayer et al., 2011). The aquifer
mesh generation and time stepping can be specifically analogue was created from surveying multiple outcrop faces
adapted to the non-uniform flow conditions. Ensembles at a gravel pit, situated in the Rhine valley in southern Ger-
of heterogeneous transmissivity fields are generated with many. The 2D information was interpolated to a 3D data set,
the SRF class where reproducibility is controlled by the including hydraulic, thermal and chemical information. The
seed and where normal fields are converted in place with workflow files are provided in Schüler and Müller (2021).
normalizer.LogNormal as shown in Fig. 25. We determine spatial correlations through variogram es-
The implementation of the randomization method timation using gstools.variogram. First, we identify
(Sect. 2.2.2) allows for the adaption of random fields to the the full transmissivity structure. The aquifer analogue data
non-uniform grid. The associated variance upscaling follows are given in a facies structure with one hydraulic conduc-

Geosci. Model Dev., 15, 3161–3182, 2022 https://doi.org/10.5194/gmd-15-3161-2022

S. Müller et al.: GeoStatTools 3179

tivity value K per facies. We calculate transmissivity by in- access: 31 March 2022). The workflows can be found in
tegrating the separate repositories https://github.com/GeoStat-Examples/
R hydraulic conductivity over the vertical axis
T (x, y) = K(x, y, z) dz. The structured transmissivity is gstools-temperature-trend, (Müller, 2021), https://github.
shown in Fig. 27, which we consider as true distribution for com/GeoStat-Examples/gstools-connectivity-and-transport,
the following exercises. (Müller and Zech, 2021a), https://github.com/GeoStat-Examples/
gstools-pumping-test-ensemble, (Müller and Zech, 2021b),
We select 13×13 virtual observations on a rectangular grid
https://github.com/GeoStat-Examples/gstools-herten-example, and
(Fig. 27), covering a subarea of about 42 m2 . These observa- (Schüler and Müller, 2021).
tions are used to determine the empirical variogram shown
in Fig. 28. We fit an exponential covariance model to the
data, which suits well with a coefficient of determination of Data availability. All data can be accessed by the given
R 2 = 0.913. DOIs (https://doi.org/10.5281/zenodo.5159728, Müller, 2021,
We use the fitted exponential variogram model and ordi- https://doi.org/10.5281/zenodo.5159578, Müller and Zech,
nary kriging to create conditioned spatial random fields with 2021a, https://doi.org/10.5281/zenodo.4891875, Müller and
CondSRF. Figure 29 shows one realization and the absolute Zech, 2021b, https://doi.org/10.5281/zenodo.5159658, Schüler
difference to the true transmissivity (Fig. 27). Differences and Müller, 2021) and the related repositories are hosted under
grow with increasing distance from observations. This trend https://github.com/GeoStat-Examples (last access: 31 March
can be even better seen in a transmissivity transect shown 2022).
in Fig. 30. The standard deviation calculated from 20 real-
izations of conditioned SRFs shows that deviations from the
reference field are significantly lower close to observation Author contributions. SM and LS are the main authors of the
GSTools package, with contributions by Falk Heße to an older
points.
version of the package. SM, LS and AZ contributed to the imple-
mentation of the workflows. AZ acted as supervisor for SM with
respect to the scientific applications of GSTools (workflows). The
5 Conclusions
article was collectively written by SM, FH, AZ and LS with major
contributions by SM.
The GSTools package provides a Python-based platform
for geostatistical applications. It is similar to software
packages like gstat for R or stand-alone packages like
Competing interests. The contact author has declared that neither
TPROGS (Carle, 1999), GSLIB (Deutsch and Journel, 1997) they nor their co-authors have any competing interests.
and S-GeMS (Remy, 2005). However, we believe that a com-
prehensive and ready-made geostatistical software package
for Python has advantages, simply through the choice of the Disclaimer. Publisher’s note: Copernicus Publications remains
programming language, as it has a gentle learning curve, is neutral with regard to jurisdictional claims in published maps and
often used as a glue language and is widely adopted by the institutional affiliations.
scientific community. Salient features of GSTools are its
random field generation and its versatile covariance model.
It is furthermore integrated with other Python packages, like Acknowledgements. We would like to thank everyone who con-
PyKrige (Murphy et al., 2021), ogs5py (Müller et al., tributed to GSTools and all the people using it and asking questions
2020) or scikit-gstat (Mälicke, 2022), and provides in- about it. We want to especially thank Mirko Mälicke for the good
terfaces to meshio (Schlömer et al., 2021) and PyVista cooperation to make SciKit-GStat and GSTools work together har-
(Sullivan and Kaszynski, 2019). GeoStat Examples (https: moniously and complement each other.
//github.com/GeoStat-Examples, last access: 31 March 2022
) provides a number of applications, including the four pre-
Financial support. Sebastian Müller and Falk Heße were finan-
sented workflows. They showcase the abilities of GSTools
cially supported by the Deutsche Forschungsgemeinschaft via grant
and can serve as a starting point for practitioners to develop
number HE-7028-1/2. Sebastian Müller was also funded by the
their own solutions for the geostatistical problems they face. German Federal Environmental Foundation (grant no. 20016/432).
This work was partially funded by the Center of Advanced Systems
Understanding (CASUS), which is financed by Germany’s Federal
Code availability. As part of the GeoStat Framework, the Ministry of Education and Research (BMBF) and by the Saxon
code of GSTools is developed at https://github.com/ Ministry for Science, Culture and Tourism (SMWK) with tax funds
GeoStat-Framework/GSTools and available via Zenodo at on the basis of the budget approved by the Saxon State Parliament.
https://doi.org/10.5281/zenodo.5883346 (Müller and Schüler,
2021). It is distributed under the GNU LGPL v3.0 licence. The article processing charges for this open-access
The documentation, which includes a quick-start guide, some publication were covered by the Helmholtz Centre for
more in-depth tutorials and a complete overview over the Environmental Research – UFZ.
API, can be accessed via https://gstools.readthedocs.io/ (last

https://doi.org/10.5194/gmd-15-3161-2022 Geosci. Model Dev., 15, 3161–3182, 2022

3180 S. Müller et al.: GeoStatTools

Review statement. This paper was edited by Fabien Maussion and conda-forge community: The conda-forge Project: Community-
reviewed by two anonymous referees. based Software Distribution Built on the conda
Package Format and Ecosystem, Zenodo [software],
https://doi.org/10.5281/zenodo.4774217, 2015.
Cressie, N. and Wikle, C. K.: Statistics for spatio-temporal data,
Wiley Series in Probability and Statistics, 1st edn., John Wiley &
References Sons, Hoboken, New Jersey, ISBN 978-0-471-69274-4, 2011.
Dagan, G.: Flow and Transport in Porous Formations, 1st edn.,
Abramowitz, M., and Stegun, I. A.: Handbook of mathematical Springer, Berlin, Heidelberg, https://doi.org/10.1007/978-3-642-
functions, 10th edn., Dover Publications, New York, ISBN 978- 75015-1, 1989.
0-486-61272-0, 1972. Deutsch, C. V. and Journel, A. G.: GSLIB: Geostatistical software
Attinger, S.: Generalized coarse graining procedures for library and user’s guide, Applied geostatistics series, 2. edn., Ox-
flow in porous media, Computat. Geosci., 7, 253–273, ford University Press, ISBN 9780195100150, 1997.
https://doi.org/10.1023/B:COMG.0000005243.73381.e3, 2003. Di Federico, V. and Neuman, S. P.: Scaling of random
Banerjee, S., Carlin, B. P., and Gelfand, A. E.: Hierarchical Mod- fields by means of truncated power variograms and as-
eling and Analysis for Spatial Data, 2 edn., Chapman and Hal- sociated spectra, Water Resour. Res., 33, 1075–1085,
l/CRC, Boca Raton, https://doi.org/10.1201/b17115, 2014. https://doi.org/10.1029/97WR00299, 1997.
Bayer, P., Huggenberger, P., Renard, P., and Comunian, Diggle, P. and Ribeiro, P. J.: Model-based Geostatistics, Springer
A.: Three-dimensional high resolution fluvio-glacial Series in Statistics, 1st edn., Springer-Verlag, New York,
aquifer analog: Part 1: Field study, J. Hydrol., 405, 1–9, https://doi.org/10.1007/978-0-387-48536-2, 2007.
https://doi.org/10.1016/j.jhydrol.2011.03.038, 2011. Eliason, S. R.: Maximum likelihood estimation: Logic and prac-
Beg, M., Taka, J., Kluyver, T., Konovalov, A., Ragan-Kelley, tice, Sage Publications, 1st edn., Thousand Oaks, CA, US, ISBN
M., Thiéry, N. M., and Fangohr, H.: Using Jupyter for Re- 9781506315904, 1993.
producible Scientific Workflows, in: Computing in Science Emery, X.: Testing the correctness of the sequential algorithm for
Engineering, Computing in Science Engineering, 23, 36–46, simulating Gaussian random fields, Stoch. Env. Res. Risk A., 18,
https://doi.org/10.1109/MCSE.2021.3052101, 2021. 401–413, https://doi.org/10.1007/s00477-004-0211-7, 2004.
Behnel, S., Bradshaw, R., Citro, C., Dalcin, L., Seljebotn, D. S., Fiori, A., Cvetkovic, V., Dagan, G., Attinger, S., Bellin, A., Diet-
and Smith, K.: Cython: The Best of Both Worlds, in: Computing rich, P., Zech, A., and Teutsch, G.: Debates – Stochastic sub-
in Science Engineering, Computing in Science Engineering, 13, surface hydrology from theory to practice: The relevance of
31–39, https://doi.org/10.1109/MCSE.2010.118, 2011. stochastic subsurface hydrology to practical problems of con-
Bellin, A. and Rubin, Y.: HYDRO_GEN: A spatially distributed taminant transport and remediation. What is characterization and
random field generator for correlated properties, Stoch. Hydrol. stochastic theory good for?, Water Resour. Res., 52, 9228–9234,
Hydraul., 10, 253–278, https://doi.org/10.1007/BF01581869, https://doi.org/10.1002/2015WR017525, 2016.
1996. Foreman-Mackey, D., Hogg, D. W., Lang, D., and Goodman, J.:
Box, G. E. P. and Cox, D. R.: An Analysis of Transformations, J. emcee: The MCMC Hammer, Publ. Astron. Soc. Pac., 125, 306–
Roy. Stat. Soc. B, 26, 211–243, https://doi.org/10.1111/j.2517- 312, https://doi.org/10.1086/670067, 2013.
6161.1964.tb00553.x, 1964. Goldstein, H.: Classical mechanics, 2nd edn., Addison-Wesley,
Brouste, A., Istas, J., and Lambert-Lacroix, S.: On Fractional ISBN 9780201029185, 1980.
Gaussian Random Fields Simulations, J. Stat. Softw., 23, 1–23 Gong, R., Haslauer, C. P., Chen, Y., and Luo, J.: Analytical re-
, https://doi.org/10.18637/jss.v023.i01, 2008. lationship between Gaussian and transformed-Gaussian spa-
Carle, S. F.: T-PROGS: Transition probability geostatistical soft- tially distributed fields, Water Resour. Res., 49, 1735–1740,
ware, version 2.1, Tech. Rep., University of California, Davis, https://doi.org/10.1002/wrcr.20143, 2013.
http://gmsdocs.aquaveo.com/t-progs.pdf (last access: 31 March Goovaerts, P.: Geostatistics in soil science: state-of-the-art and per-
2022), 1999. spectives, Geoderma, 89, 1–45, https://doi.org/10.1016/S0016-
Cecinati, F., Wani, O., and Rico-Ramirez, M. A.: Comparing Ap- 7061(98)00078-0, 1999.
proaches to Deal With Non-Gaussianity of Rainfall Data in Gutzmann, B., Motl, A., Lassahn, D., Kamenshchikov,
Kriging-Based Radar-Gauge Rainfall Merging, Water Resour. I., Bachmann, M., and Schrammel, M.: earthob-
Res., 53, 8999–9018, https://doi.org/10.1002/2016WR020330, servations/wetterdienst: v0.18.0, Zenodo [software],
2017. https://doi.org/10.5281/zenodo.4737739, 2021.
Chilès, J.-P. and Delfiner, P.: Geostatistics: Modeling Spatial Un- Harris, C. R., Millman, K. J., van der Walt, S. J., Gommers,
certainty, second edn., Wiley Series in Probability and Statistics, R., Virtanen, P., Cournapeau, D., Wieser, E., Taylor, J., Berg,
edited by: Balding, D. J., Cressie, N. A. C., Fitzmaurice, G. M., S., Smith, N. J., Kern, R., Picus, M., Hoyer, S., van Kerk-
Goldstein, H., Johnstone, I. M., Molenberghs, G., Scott, D. W., wijk, M. H., Brett, M., Haldane, A., Fernández del Río, J.,
Smith, A. F. M., Tsay, R. S., and Weisberg, S., John Wiley & Wiebe, M., Peterson, P., Gérard-Marchant, P., Sheppard, K.,
Sons, https://doi.org/10.1002/9781118136188, 2012. Reddy, T., Weckesser, W., Abbasi, H., Gohlke, C., and Oliphant,
Cirpka, O. A. and Valocchi, A. J.: Debates – Stochastic sub- T. E.: Array programming with NumPy, Nature, 585, 357–362,
surface hydrology from theory to practice: Does stochastic https://doi.org/10.1038/s41586-020-2649-2, 2020.
subsurface hydrology help solving practical problems of con- Heße, F., Prykhodko, V., Schlüter, S., and Attinger, S.: Generating
taminant hydrogeology?, Water Resour. Res., 52, 9218–9227, random fields with a truncated power-law variogram: A compar-
https://doi.org/10.1002/2016WR019087, 2016.

Geosci. Model Dev., 15, 3161–3182, 2022 https://doi.org/10.5194/gmd-15-3161-2022

S. Müller et al.: GeoStatTools 3181

ison of several numerical methods, Environ. Modell. Softw., 55, Müller, S. and Zech, A.: GeoStat – Examples/gstools-
32–48, https://doi.org/10.1016/j.envsoft.2014.01.013, 2014. connectivity-and-transport: v1.0, Zenodo [data set],
Hohn, M.: Geostatistics and Petroleum Geology, Computer https://doi.org/10.5281/zenodo.5159578, 2021a.
Methods in the Geosciences, 2 edn., Springer Netherlands, Müller, S. and Zech, A.: GeoStat – Examples/gstools-
https://doi.org/10.1007/978-94-011-4425-4, 1999. pumping-test-ensemble: v1.0, Zenodo [data set],
John, J. A. and Draper, N. R.: An Alternative Family of https://doi.org/10.5281/zenodo.4891875, 2021b.
Transformations, J. Roy. Stat. Soc. C-App., 29, 190–197, Müller, S., Zech, A., and Heße, F.: ogs5py: A Python – API for the
https://doi.org/10.2307/2986305, 1980. OpenGeoSys 5 Scientific Modeling Package, Groundwater, 59,
Kitanidis, P.: Introduction to Geostatistics: Applications in Hydro- 117–122, https://doi.org/10.1111/gwat.13017, 2020.
geology, 1st edn., Cambridge University Press, Cambridge, New Müller, S., Heße, F., Attinger, S., and Zech, A.: The extended gen-
York, ISBN 9780521587471, 2008. eralized radial flow model and effective conductivity for trun-
Kolditz, O., Bauer, S., Bilke, L., Böttcher, N., Delfs, J.-O., Fis- cated power law variograms, Adv. Water Resour., 156, 104027,
cher, T., Görke, U. J., Kalbacher, T., Kosakowski, G., McDer- https://doi.org/10.1016/j.advwatres.2021.104027, 2021a.
mott, C. I., Park, C. H., Radu, F., Rink, K., Shao, H., Shao, Müller, S., Leven, C., Dietrich, P., Attinger, S., and Zech, A.:
H. B., Sun, F., Sun, Y. Y., Singh, A. K., Taron, J., Walther, How to Find Aquifer Statistics Utilizing Pumping Tests? Two
M., Wang, W., Watanabe, N., Wu, Y., Xie, M., Xu, W., and Field Studies Using welltestpy, Groundwater, 60, 137–144,
Zehner, B.: OpenGeoSys: An open-source initiative for numer- https://doi.org/10.1111/gwat.13121, 2021b.
ical simulation of thermo-hydro-mechanical/chemical (THM/C) Murphy, B., Müller, S., and Yurchak, R.: GeoStat-
processes in porous media, Environ. Earth Sci., 67, 589–599, Framework/PyKrige: v1.6.0, Zenodo [code],
https://doi.org/10.1007/s12665-012-1546-x, 2012. https://doi.org/10.5281/zenodo.4661732, 2021.
Kraichnan, R.: Diffusion by a Random Velocity Field, Phys. Fluids, Murray, S. G. and Poulin, F. J.: hankel: A Python li-
13, 22–31, https://doi.org/10.1063/1.1692799, 1970. brary for performing simple and accurate Hankel transfor-
Krige, D. G.: A statistical approach to some basic mine valuation mations, The Journal of Open Source Software, 4, 1397,
problems on the Witwatersrand, J. S. Afr. I. Min. Metall., 52, https://doi.org/10.21105/joss.01397, 2019.
119–139, 1951. Neuman, S. P.: Stochastic groundwater models in practice, Stoch.
Lantuéjoul, C., Freulon, X., and Renard, D.: Spectral Simu- Env. Res. Risk A., 18, 268–270, https://doi.org/10.1007/s00477-
lation of Isotropic Gaussian Random Fields on a Sphere, 004-0192-6, 2004.
Math. Geosci., 51, 999–1020, https://doi.org/10.1007/s11004- Ogata, H.: A Numerical Integration Formula Based on the
019-09799-4, 2019. Bessel Functions, Publ. Res. I. Math. Sci., 41, 949–970,
Mälicke, M.: SciKit-GStat 1.0: a SciPy-flavored geostatisti- https://doi.org/10.2977/prims/1145474602, 2005.
cal variogram estimation toolbox written in Python, Geosci. Pebesma, E. J.: Multivariable geostatistics in S: the
Model Dev., 15, 2505–2532, https://doi.org/10.5194/gmd-15- gstat package, Comput. Geosci., 30, 683–691,
2505-2022, 2022. https://doi.org/10.1016/j.cageo.2004.03.012, 2004.
Manly, B. F. J.: Exponential Data Transformations, J. Roy. Stat. Soc. Perkel, J. M.: Why Jupyter is data scientists’ compu-
D-Sta., 25, 37–42, https://doi.org/10.2307/2988129, 1976. tational notebook of choice, Nature, 563, 145–146,
Matheron, G.: Traité de géostatistique appliquée, no. 14 in Mé- https://doi.org/10.1038/d41586-018-07196-1, 2018.
moires du BRGM, Editions Technip, Tome II: le krigeage, no Pyrcz, M. J. and Deutsch, C. V.: Geostatistical Reservoir Modeling,
24, Editions BRGM, Paris, 1962. 2 edn., Oxford University Press, Oxford, ISBN 978-0199731442,
Matern, B.: Spatial variation – stochastic models and their applica- 2014.
tions to some problems in forest survey sampling investigations, Queiroz, F., Silva, R., Miller, J., Brockhauser, S., and
Report of the Forest Research Institute of Sweden 49, 1–144, Fangohr, H.: Track 1 Paper: Good Usability Prac-
1960 (in English, Swedish summary). tices in Scientific Software Development, Figshare,
Mohammadi, H., Riche, R. L., Durrande, N., Touboul, E., and Bay, https://doi.org/10.6084/m9.figshare.5331814.v3, 2017.
X.: An analytic comparison of regularization methods for Gaus- Rajaram, H.: Debates – Stochastic subsurface hydrology from the-
sian Processes, arXiv [preprint], arXiv:1602.00853, 5 May 2017. ory to practice: Introduction, Water Resour. Res., 52, 9215–9217,
Monestiez, P., Petrenko, A., Leredde, Y., and Ongari, B.: Geosta- https://doi.org/10.1002/2016WR020066, 2016.
tistical analysis of three dimensional current patterns in coastal Rasmussen, C. E. and Williams, C. K. I.: Gaus-
oceanography: Application to the gulf of lions (NW mediter- sian Processes for Machine Learning, 1st
ranean sea), in: geoENV IV – Geostatistics for environmental ap- edn., The MIT Press, ISBN 9780262256834,
plications, edited by: Sanchez-Vila, X., Carrera, J., and Gómez- https://doi.org/10.7551/mitpress/3206.001.0001, 2005.
Hernández, J. J., pp. 367–378, Springer Netherlands, Dordrecht, Remy, N.: S-GeMS: The Stanford Geostatistical Modeling Soft-
https://doi.org/10.1007/1-4020-2115-1_31, 2004. ware: A Tool for New Algorithms Development, in: Quantitative
Müller, S.: GeoStat – Examples/gstools-temperature-trend: v1.0, Geology and Geostatistics, Geostatistics Banff 2004, 14, 865–
Zenodo [data set], https://doi.org/10.5281/zenodo.5159728, 871, https://doi.org/10.1007/978-1-4020-3610-1_89, 2005.
2021. Rossi, R. E., Mulla, D. J., Journel, A. G., and Franz,
Müller, S. and Schüler, L.: GeoStat – Frame- E. H.: Geostatistical tools for modeling and interpreting
work/GSTools: v1.3.5 “Pure Pink”, Zenodo [code], ecological spatial dependence, Ecol. Monogr., 62, 277–314,
https://doi.org/10.5281/zenodo.5883346, 2021. https://doi.org/10.2307/2937096, 1992.

https://doi.org/10.5194/gmd-15-3161-2022 Geosci. Model Dev., 15, 3161–3182, 2022

3182 S. Müller et al.: GeoStatTools

Rubin, Y.: Applied Stochastic Hydrogeology, 1st edn., Ox- Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy,
ford University Press, New York, ISBN 9780195138047, T., Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W.,
https://doi.org/10.1093/oso/9780195138047.001.0001, 2003. Bright, J., van der Walt, S. J., Brett, M., Wilson, J., Millman,
Rubin, Y., Chen, X., Murakami, H., and Hahn, M.: A Bayesian K. J., Mayorov, N., Nelson, A. R. J., Jones, E., Kern, R., Larson,
approach for inverse modeling, data assimilation, and condi- E., Carey, C. J., Polat, i., Feng, Y., Moore, E. W., VanderPlas,
tional simulation of spatial random fields, Water Resour. Res., J., Laxalde, D., Perktold, J., Cimrman, R., Henriksen, I., Quin-
46, W10523, https://doi.org/10.1029/2009WR008799, 2010. tero, E. A., Harris, C. R., Archibald, A. M., Ribeiro, A. H., Pe-
Rubin, Y., Chang, C.-F., Chen, J., Cucchi, K., Harken, B., Heße, F., dregosa, F., van Mulbregt, P., and SciPy 1.0 Contributors: SciPy
and Savoy, H.: Stochastic hydrogeology’s biggest hurdles ana- 1.0: Fundamental algorithms for scientific computing in python,
lyzed and its big blind spot, Hydrol. Earth Syst. Sci., 22, 5675– Nat. Methods, 17, 261–272, https://doi.org/10.1038/s41592-019-
5695, https://doi.org/10.5194/hess-22-5675-2018, 2018. 0686-2, 2020.
Rudin, W.: Fourier Analysis on Groups, 1st edn., Wiley- Vrugt, J. A.: Markov chain Monte Carlo simulation using
Interscience, John Wiley & Sons, ISBN 9780470744819, the DREAM software package: Theory, concepts, and MAT-
https://doi.org/10.1002/9781118165621, 1990. LAB implementation, Environ. Modell. Softw., 75, 273–316,
Sales, M. H., Souza, C. M., Kyriakidis, P. C., Roberts, https://doi.org/10.1016/j.envsoft.2015.08.013, 2016.
D. A., and Vidal, E.: Improving spatial distribution es- Wackernagel, H.: Multivariate Geostatistics: An Introduction
timation of forest biomass with geostatistics: A case with Applications, 3 edn., Springer-Verlag, Berlin Heidel-
study for Rondônia, Brazil, Ecol. Model., 205, 221–230, berg, ISBN 978-3-540-44142-7, https://doi.org/10.1007/978-3-
https://doi.org/10.1016/j.ecolmodel.2007.02.033, 2007. 662-05294-5, 2003.
Savoy, H., Heße, F., and Rubin, Y.: anchoredDistr: a Pack- Webster, R. and Oliver, M. A.: Geostatistics for Environmental Sci-
age for the Bayesian Inversion of Geostatistical Parameters entists, 2 edn., John Wiley & Sons, ISBN 978-0-470-02858-2,
with Multi-type and Multi-scale Data, R Journal, 9, 6–17, 2007.
https://doi.org/10.32614/RJ-2017-034, 2017. Wendland, H.: Piecewise polynomial, positive definite and com-
Schlömer, N., McBain, G. D., Luu, K., christos, Li, T., Hochste- pactly supported radial functions of minimal degree, Adv. Com-
ger, M., Keilegavlen, E., Ferrándiz, V. M., Barnes, C., put. Math., 4, 389–396, https://doi.org/10.1007/BF02123482,
Lukeš, V., Dalcin, L., Jansen, M., Wagner, N., Gupta, A., 1995.
Müller, S., Woodsend, B., Andersen, K., Schwarz, L., Blechta, Winter, C. L.: Stochastic hydrology: practical alterna-
J., Christovasilis, I. P., Coutinho, C., Beurle, D., ffilotto, tives exist, Stoch. Env. Res. Risk A., 18, 271–273,
Dokken, J. S., blacheref, so1291, Cervone, A., Shrimali, B., https://doi.org/10.1007/s00477-004-0198-0, 2004.
Bill, and Jones, D.: nschloe/meshio: None, Zenodo [code], Yeo, I. and Johnson, R. A.: A new family of power transforma-
https://doi.org/10.5281/zenodo.4900671, 2021. tions to improve normality or symmetry, Biometrika, 87, 954–
Schüler, L. and Müller, S.: GeoStat – Examples/gstools- 959, https://doi.org/10.1093/biomet/87.4.954, 2000.
herten-example: v1.0, Zenodo [data set], Zech, A., Schneider, C. L., and Attinger, S.: The extended Thiem’s
https://doi.org/10.5281/zenodo.5159658, 2021. solution – Including the impact of heterogeneity, Water Resour.
Schüler, L., Suciu, N., Knabner, P., and Attinger, S.: A time Res., 48, W10535, https://doi.org/10.1029/2012WR011852,
dependent mixing model to close PDF equations for trans- 2012.
port in heterogeneous aquifers, Adv. Water Resour., 96, 55–67, Zech, A., Müller, S., Mai, J., Heße, F., and Attinger, S.: Extending
https://doi.org/10.1016/j.advwatres.2016.06.012, 2016. Theis’ solution: Using transient pumping tests to estimate pa-
Schüler, L., Calabrese, J. M., and Attinger, S.: Data driven rameters of aquifer heterogeneity, Water Resour. Res., 52, 6156–
high resolution modeling and spatial analyses of the 6170, https://doi.org/10.1002/2015WR018509, 2016.
COVID-19 pandemic in Germany, PLOS ONE, 16, 1–14, Zhang, Y.-K. and Zhang, D.: Forum: The state of stochas-
https://doi.org/10.1371/journal.pone.0254660, 2021. tic hydrology, Stoch. Env. Res. Risk A., 18, 265,
Sturges, H. A.: The Choice of a Class Interval, J. Am. Stat. Assoc., https://doi.org/10.1007/s00477-004-0190-8, 2004.
21, 65–66, https://doi.org/10.1080/01621459.1926.10502161, Zinn, B. and Harvey, C. F.: When good statistical models of
1926. aquifer heterogeneity go bad: A comparison of flow, disper-
Sullivan, C. B. and Kaszynski, A. A.: PyVista: 3D plotting and sion, and mass transfer in connected and multivariate Gaus-
mesh analysis through a streamlined interface for the Visualiza- sian hydraulic conductivity fields, Water Resour. Res., 39, 1051,
tion Toolkit (VTK), Journal of Open Source Software, 4, 1450, https://doi.org/10.1029/2001WR001146, 2003.
https://doi.org/10.21105/joss.01450, 2019.
Uieda, L.: Verde: Processing and gridding spatial data using
Green’s functions, Journal of Open Source Software, 3, 957,
https://doi.org/10.21105/joss.00957, 2018.

Geosci. Model Dev., 15, 3161–3182, 2022 https://doi.org/10.5194/gmd-15-3161-2022