Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
def setUp(self):
self.test_dir = tempfile.mkdtemp()
# structured field with a size 100x100x100 and a grid-size of 1x1x1
x = y = z = range(50)
model = Gaussian(dim=3, var=0.6, len_scale=20)
self.srf_structured = SRF(model)
self.srf_structured((x, y, z), mesh_type="structured")
# unstrucutred field
seed = MasterRNG(19970221)
rng = np.random.RandomState(seed())
x = rng.randint(0, 100, size=1000)
y = rng.randint(0, 100, size=1000)
model = Exponential(
dim=2, var=1, len_scale=[12.0, 3.0], angles=np.pi / 8.0
)
self.srf_unstructured = SRF(model, seed=20170519)
self.srf_unstructured([x, y])
def test_sample_dist(self):
model = Gaussian(dim=1, var=3.5, len_scale=8.0)
pdf, cdf, ppf = model.dist_func
rad = self.rng.sample_dist(
size=self.few_modes, pdf=pdf, cdf=cdf, ppf=ppf, a=0
)
self.assertEqual(rad.shape[0], self.few_modes)
model = Gaussian(dim=2, var=3.5, len_scale=8.0)
pdf, cdf, ppf = model.dist_func
rad = self.rng.sample_dist(
size=self.few_modes, pdf=pdf, cdf=cdf, ppf=ppf, a=0
)
self.assertEqual(rad.shape[0], self.few_modes)
model = Gaussian(dim=3, var=3.5, len_scale=8.0)
pdf, cdf, ppf = model.dist_func
rad = self.rng.sample_dist(
"""
Zinn & Harvey transformation
----------------------------
Here, we transform a field with the so called "Zinn & Harvey" transformation presented in
`Zinn & Harvey (2003) `__.
With this transformation, one could overcome the restriction that in ordinary
Gaussian random fields the mean values are the ones being the most connected.
"""
import gstools as gs
# structured field with a size of 100x100 and a grid-size of 1x1
x = y = range(100)
model = gs.Gaussian(dim=2, var=1, len_scale=10)
srf = gs.SRF(model, seed=20170519)
srf.structured([x, y])
gs.transform.zinnharvey(srf, conn="high")
srf.plot()
"""
Generating a Random Vector Field
--------------------------------
As a first example we are going to generate a vector field with a Gaussian
covariance model on a structured grid:
"""
import numpy as np
import gstools as gs
# the grid
x = np.arange(100)
y = np.arange(100)
# a smooth Gaussian covariance model
model = gs.Gaussian(dim=2, var=1, len_scale=10)
srf = gs.SRF(model, generator="VectorField", seed=19841203)
srf((x, y), mesh_type="structured")
srf.plot()
###############################################################################
# Let us have a look at the influence of the covariance model. Choosing the
# exponential model and keeping all other parameters the same
# a rougher exponential covariance model
model2 = gs.Exponential(dim=2, var=1, len_scale=10)
srf.model = model2
srf((x, y), mesh_type="structured", seed=19841203)
srf.plot()
.. math:: S(k) = \left(\frac{1}{2\pi}\right)^n \cdot
\frac{(2\pi)^{n/2}}{(bk)^{n/2-1}}
\int_0^\infty r^{n/2-1} C(r) J_{n/2-1}(bkr) r dr
Where :math:`k=\left\Vert\mathbf{k}\right\Vert`.
Depending on the spectrum, the spectral-density is defined by:
.. math:: \tilde{S}(k) = \frac{S(k)}{\sigma^2}
You can access these methods by:
"""
import gstools as gs
model = gs.Gaussian(dim=3, var=2.0, len_scale=10)
ax = model.plot("spectrum")
model.plot("spectral_density", ax=ax)
"""
Creating an Ensemble of Fields
------------------------------
Creating an ensemble of random fields would also be
a great idea. Let's reuse most of the previous code.
"""
import numpy as np
import matplotlib.pyplot as pt
import gstools as gs
x = y = np.arange(100)
model = gs.Gaussian(dim=2, var=1, len_scale=10)
srf = gs.SRF(model)
###############################################################################
# This time, we did not provide a seed to :any:`SRF`, as the seeds will used
# during the actual computation of the fields. We will create four ensemble
# members, for better visualisation and save them in a list and in a first
# step, we will be using the loop counter as the seeds.
ens_no = 4
field = []
for i in range(ens_no):
field.append(srf.structured([x, y], seed=i))
###############################################################################
# Now let's have a look at the results:
Compare Kriging
---------------
"""
import numpy as np
from gstools import Gaussian, krige
import matplotlib.pyplot as plt
# 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)
###############################################################################
# spatial random field class
model = Gaussian(dim=1, var=0.5, len_scale=2)
###############################################################################
kr1 = krige.Simple(model=model, mean=1, cond_pos=cond_pos, cond_val=cond_val)
kr2 = krige.Ordinary(model=model, cond_pos=cond_pos, cond_val=cond_val)
kr1(gridx)
kr2(gridx)
###############################################################################
plt.plot(gridx, kr1.field, label="simple kriged field")
plt.plot(gridx, kr2.field, label="ordinary kriged field")
plt.scatter(cond_pos, cond_val, color="k", zorder=10, label="Conditions")
plt.legend()
plt.show()
In the following we will show the influence of the nugget and
measurement errors.
"""
import numpy as np
import gstools as gs
# condtions
cond_pos = [0.3, 1.1, 1.9, 3.3, 4.7]
cond_val = [0.47, 0.74, 0.56, 1.47, 1.74]
cond_err = [0.01, 0.0, 0.1, 0.05, 0]
# resulting grid
gridx = np.linspace(0.0, 15.0, 151)
# spatial random field class
model = gs.Gaussian(dim=1, var=0.9, len_scale=1, nugget=0.1)
###############################################################################
# Here we will use Simple kriging (`unbiased=False`) to interpolate the given
# conditions.
krig = gs.krige.Krige(
model=model,
cond_pos=cond_pos,
cond_val=cond_val,
mean=1,
unbiased=False,
exact=False,
cond_err=cond_err,
)
krig(gridx)
"""
Anisotropy and Rotation
=======================
The internally used (semi-) variogram
represents the isotropic case for the model.
Nevertheless, you can provide anisotropy ratios by:
"""
import gstools as gs
model = gs.Gaussian(dim=3, var=2.0, len_scale=10, anis=0.5)
print(model.anis)
print(model.len_scale_vec)
###############################################################################
# As you can see, we defined just one anisotropy-ratio
# and the second transversal direction was filled up with ``1.``.
# You can get the length-scales in each direction by
# the attribute :any:`CovModel.len_scale_vec`. For full control you can set
# a list of anistropy ratios: ``anis=[0.5, 0.4]``.
#
# Alternatively you can provide a list of length-scales:
model = gs.Gaussian(dim=3, var=2.0, len_scale=[10, 5, 4])
model.plot("vario_spatial")
print("Anisotropy representations:")
"""
binary fields
-------------
Here we transform a field to a binary field with only two values.
The dividing value is the mean by default and the upper and lower values
are derived to preserve the variance.
"""
import gstools as gs
# structured field with a size of 100x100 and a grid-size of 1x1
x = y = range(100)
model = gs.Gaussian(dim=2, var=1, len_scale=10)
srf = gs.SRF(model, seed=20170519)
srf.structured([x, y])
gs.transform.binary(srf)
srf.plot()