Skip to content

Commit

Permalink
add multiple examples
Browse files Browse the repository at this point in the history
  • Loading branch information
MuellerSeb committed Aug 28, 2019
1 parent d147af5 commit 38e48cc
Show file tree
Hide file tree
Showing 10 changed files with 182 additions and 0 deletions.
14 changes: 14 additions & 0 deletions examples/10_simple_kriging.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,14 @@
import numpy as np
from gstools import Gaussian, krige
# 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)
krig = krige.Simple(model, mean=1, cond_pos=[cond_pos], cond_val=cond_val)
krig([gridx])
ax = krig.plot()
ax.scatter(cond_pos, cond_val, color="k", zorder=10, label="Conditions")
ax.legend()
14 changes: 14 additions & 0 deletions examples/11_ordinary_kriging.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,14 @@
import numpy as np
from gstools import Gaussian, krige
# 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)
krig = krige.Ordinary(model, cond_pos=[cond_pos], cond_val=cond_val)
krig([gridx])
ax = krig.plot()
ax.scatter(cond_pos, cond_val, color="k", zorder=10, label="Conditions")
ax.legend()
19 changes: 19 additions & 0 deletions examples/12_compare_kriging.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,19 @@
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()
23 changes: 23 additions & 0 deletions examples/13_condition_ensemble.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,23 @@
import numpy as np
from gstools import Gaussian, SRF
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]
gridx = np.linspace(0.0, 15.0, 151)
# spatial random field class
model = Gaussian(dim=1, var=0.5, len_scale=2)
srf = SRF(model)
srf.set_condition([cond_pos], cond_val, "ordinary")
fields = []
for i in range(100):
if i % 10 == 0: print(i)
fields.append(srf([gridx], seed=i))
label = "Conditioned ensemble" if i == 0 else None
plt.plot(gridx, fields[i], color="k", alpha=0.1, label=label)
plt.plot(gridx, np.full_like(gridx, srf.mean), label="estimated mean")
plt.plot(gridx, np.mean(fields, axis=0), linestyle=':', label="Ensemble mean")
plt.plot(gridx, srf.krige_field, linestyle='dashed', label="kriged field")
plt.scatter(cond_pos, cond_val, color="k", zorder=10, label="Conditions")
plt.legend()
plt.show()
9 changes: 9 additions & 0 deletions examples/14_transform_01.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,9 @@
from gstools import SRF, Gaussian
from gstools import transform as tf
# structured field with a size of 100x100 and a grid-size of 1x1
x = y = range(100)
model = Gaussian(dim=2, var=1, len_scale=10)
srf = SRF(model, seed=20170519)
srf.structured([x, y])
tf.normal_to_lognormal(srf)
srf.plot()
9 changes: 9 additions & 0 deletions examples/15_transform_02.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,9 @@
from gstools import SRF, Gaussian
from gstools import transform as tf
# structured field with a size of 100x100 and a grid-size of 1x1
x = y = range(100)
model = Gaussian(dim=2, var=1, len_scale=10)
srf = SRF(model, seed=20170519)
srf.structured([x, y])
tf.binary(srf)
srf.plot()
9 changes: 9 additions & 0 deletions examples/16_transform_03.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,9 @@
from gstools import SRF, Gaussian
from gstools import transform as tf
# structured field with a size of 100x100 and a grid-size of 1x1
x = y = range(100)
model = Gaussian(dim=2, var=1, len_scale=10)
srf = SRF(model, seed=20170519)
srf.structured([x, y])
tf.zinnharvey(srf, conn="high")
srf.plot()
9 changes: 9 additions & 0 deletions examples/17_transform_04.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,9 @@
from gstools import SRF, Gaussian
from gstools import transform as tf
# structured field with a size of 100x100 and a grid-size of 1x1
x = y = range(100)
model = Gaussian(dim=2, var=1, len_scale=10)
srf = SRF(model, seed=20170519)
field = srf.structured([x, y])
tf.normal_to_arcsin(srf)
srf.plot()
12 changes: 12 additions & 0 deletions examples/18_transform_05.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,12 @@
from gstools import SRF, Gaussian
from gstools import transform as tf
# structured field with a size of 100x100 and a grid-size of 1x1
x = y = range(100)
model = Gaussian(dim=2, var=1, len_scale=10)
srf = SRF(model, mean=-9, seed=20170519)
srf.structured([x, y])
tf.normal_force_moments(srf)
tf.zinnharvey(srf, conn="low")
tf.binary(srf)
tf.normal_to_lognormal(srf)
srf.plot()
64 changes: 64 additions & 0 deletions examples/19_check_rand_meth_sampling.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,64 @@
import numpy as np
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

from gstools import SRF, Stable


def norm_rad(vec):
"""Direction on the unit sphere."""
vec = np.array(vec, ndmin=2)
norm = np.zeros(vec.shape[1])
for i in range(vec.shape[0]):
norm += vec[i]**2
norm = np.sqrt(norm)
return np.einsum("j,ij->ij", 1 / norm, vec), norm


def plot_rand_meth_samples(generator):
"""Plot the samples of the rand meth class."""
norm, rad = norm_rad(generator._cov_sample)

fig = plt.figure(figsize=(10, 4))

if generator.model.dim == 3:
ax = fig.add_subplot(121, projection=Axes3D.name)
u = np.linspace(0, 2 * np.pi, 100)
v = np.linspace(0, np.pi, 100)
x = np.outer(np.cos(u), np.sin(v))
y = np.outer(np.sin(u), np.sin(v))
z = np.outer(np.ones(np.size(u)), np.cos(v))
ax.plot_surface(x, y, z, rstride=4, cstride=4, color='b', alpha=0.1)
ax.scatter(norm[0], norm[1], norm[2])
ax.set_aspect('equal')
elif generator.model.dim == 2:
ax = fig.add_subplot(121)
u = np.linspace(0, 2 * np.pi, 100)
x = np.cos(u)
y = np.sin(u)
ax.plot(x, y, color='b', alpha=0.1)
ax.scatter(norm[0], norm[1])
ax.set_aspect('equal')
else:
ax = fig.add_subplot(121)
ax.bar(-1, np.sum(np.isclose(norm, -1)), color="C0")
ax.bar(1, np.sum(np.isclose(norm, 1)), color="C0")
ax.set_xticks([-1, 1])
ax.set_xticklabels(('-1', '1'))
ax.set_title("Direction sampling")

ax = fig.add_subplot(122)
# x = np.linspace(0, np.max(rad))
x = np.linspace(0, 10 / generator.model.integral_scale)
y = generator.model.spectral_rad_pdf(x)
ax.plot(x, y)
sample_in = np.sum(rad <= np.max(x))
ax.hist(rad[rad <= np.max(x)], bins=sample_in // 50, density=True)
ax.set_xlim([0, np.max(x)])
ax.set_title("Radius samples shown {}/{}".format(sample_in, len(rad)))
fig.show()


model = Stable(dim=3, alpha=1)
srf = SRF(model)
plot_rand_meth_samples(srf.generator)

0 comments on commit 38e48cc

Please sign in to comment.