Skip to content

Commit

Permalink
Merge 84bc379 into 21c97fc
Browse files Browse the repository at this point in the history
  • Loading branch information
MuellerSeb committed Aug 2, 2021
2 parents 21c97fc + 84bc379 commit f71821b
Show file tree
Hide file tree
Showing 21 changed files with 1,143 additions and 314 deletions.
2 changes: 0 additions & 2 deletions docs/source/transform.rst
Original file line number Diff line number Diff line change
Expand Up @@ -2,8 +2,6 @@ gstools.transform
=================

.. automodule:: gstools.transform
:members:
:undoc-members:

.. raw:: latex

Expand Down
17 changes: 9 additions & 8 deletions examples/01_random_field/01_srf_ensemble.py
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,8 @@
Creating an ensemble of random fields would also be
a great idea. Let's reuse most of the previous code.
We will set the position tuple `pos` before generation to reuse it afterwards.
"""

import numpy as np
Expand All @@ -14,26 +16,26 @@

model = gs.Gaussian(dim=2, var=1, len_scale=10)
srf = gs.SRF(model)
srf.set_pos([x, y], "structured")

###############################################################################
# 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
# members, for better visualisation, save them in to srf class 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))
srf(seed=i, store=f"field{i}")

###############################################################################
# Now let's have a look at the results:
# Now let's have a look at the results. We can access the fields by name or
# index:

fig, ax = pt.subplots(2, 2, sharex=True, sharey=True)
ax = ax.flatten()
for i in range(ens_no):
ax[i].imshow(field[i].T, origin="lower")
ax[i].imshow(srf[i].T, origin="lower")
pt.show()

###############################################################################
Expand All @@ -44,9 +46,8 @@
# provides a seed generator :any:`MasterRNG`. The loop, in which the fields are
# generated would then look like


from gstools.random import MasterRNG

seed = MasterRNG(20170519)
for i in range(ens_no):
field.append(srf.structured([x, y], seed=seed()))
srf(seed=seed(), store=f"better_field{i}")
11 changes: 7 additions & 4 deletions examples/06_conditioned_fields/00_condition_ensemble.py
Original file line number Diff line number Diff line change
Expand Up @@ -26,15 +26,18 @@
model = gs.Gaussian(dim=1, var=0.5, len_scale=1.5)
krige = gs.krige.Ordinary(model, cond_pos, cond_val)
cond_srf = gs.CondSRF(krige)
cond_srf.set_pos(gridx)

###############################################################################
# We can specify individual names for each field by the keyword `store`:

fields = []
for i in range(100):
fields.append(cond_srf(gridx, seed=i))
cond_srf(seed=i, store=f"f{i}")
label = "Conditioned ensemble" if i == 0 else None
plt.plot(gridx, fields[i], color="k", alpha=0.1, label=label)
plt.plot(gridx, cond_srf.krige(gridx, only_mean=True), label="estimated mean")
plt.plot(gridx, cond_srf[f"f{i}"], color="k", alpha=0.1, label=label)

fields = [cond_srf[f"f{i}"] for i in range(100)]
plt.plot(gridx, cond_srf.krige(only_mean=True), label="estimated mean")
plt.plot(gridx, np.mean(fields, axis=0), linestyle=":", label="Ensemble mean")
plt.plot(gridx, cond_srf.krige.field, linestyle="dashed", label="kriged field")
plt.scatter(cond_pos, cond_val, color="k", zorder=10, label="Conditions")
Expand Down
33 changes: 22 additions & 11 deletions examples/06_conditioned_fields/01_2D_condition_ensemble.py
Original file line number Diff line number Diff line change
Expand Up @@ -20,34 +20,39 @@
model = gs.Gaussian(dim=2, var=0.5, len_scale=5, anis=0.5, angles=-0.5)
krige = gs.Krige(model, cond_pos=cond_pos, cond_val=cond_val)
cond_srf = gs.CondSRF(krige)
cond_srf.set_pos([x, y], "structured")

###############################################################################
# We create a list containing the generated conditioned fields.
# By specifying ``store=[f"fld{i}", False, False]``, only the conditioned field
# is stored with the specified name. The raw random field and the raw kriging
# field is not stored. This way, we can access each conditioned field by index
# ``cond_srf[i]``:

ens_no = 4
field = []
for i in range(ens_no):
field.append(cond_srf.structured([x, y], seed=i))
cond_srf(seed=i, store=[f"fld{i}", False, False])

###############################################################################
# Now let's have a look at the pairwise differences between the generated
# fields. We will see, that they coincide at the given conditions.

fig, ax = plt.subplots(ens_no + 1, ens_no + 1, figsize=(8, 8))
# plotting kwargs for scatter and image
sc_kwargs = dict(c=cond_val, edgecolors="k", vmin=0, vmax=np.max(field))
im_kwargs = dict(extent=2 * [0, 5], origin="lower", vmin=0, vmax=np.max(field))
vmax = np.max(cond_srf.all_fields)
sc_kw = dict(c=cond_val, edgecolors="k", vmin=0, vmax=vmax)
im_kw = dict(extent=2 * [0, 5], origin="lower", vmin=0, vmax=vmax)
for i in range(ens_no):
# conditioned fields and conditions
ax[i + 1, 0].imshow(field[i].T, **im_kwargs)
ax[i + 1, 0].scatter(*cond_pos, **sc_kwargs)
ax[i + 1, 0].set_ylabel(f"Field {i+1}", fontsize=10)
ax[0, i + 1].imshow(field[i].T, **im_kwargs)
ax[0, i + 1].scatter(*cond_pos, **sc_kwargs)
ax[0, i + 1].set_title(f"Field {i+1}", fontsize=10)
ax[i + 1, 0].imshow(cond_srf[i].T, **im_kw)
ax[i + 1, 0].scatter(*cond_pos, **sc_kw)
ax[i + 1, 0].set_ylabel(f"Field {i}", fontsize=10)
ax[0, i + 1].imshow(cond_srf[i].T, **im_kw)
ax[0, i + 1].scatter(*cond_pos, **sc_kw)
ax[0, i + 1].set_title(f"Field {i}", fontsize=10)
# absolute differences
for j in range(ens_no):
ax[i + 1, j + 1].imshow(np.abs(field[i] - field[j]).T, **im_kwargs)
ax[i + 1, j + 1].imshow(np.abs(cond_srf[i] - cond_srf[j]).T, **im_kw)

# beautify plots
ax[0, 0].axis("off")
Expand All @@ -56,3 +61,9 @@
a.set_xticks([]), a.set_yticks([])
fig.subplots_adjust(wspace=0, hspace=0)
fig.show()

###############################################################################
# To check if the generated fields are correct, we can have a look at their
# names:

print(cond_srf.field_names)
4 changes: 3 additions & 1 deletion examples/07_transformations/00_log_normal.py
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,8 @@
-----------------
Here we transform a field to a log-normal distribution:
See :any:`transform.normal_to_lognormal`
"""
import gstools as gs

Expand All @@ -11,5 +13,5 @@
model = gs.Gaussian(dim=2, var=1, len_scale=10)
srf = gs.SRF(model, seed=20170519)
srf.structured([x, y])
gs.transform.normal_to_lognormal(srf)
srf.transform("normal_to_lognormal") # also "lognormal" works
srf.plot()
4 changes: 3 additions & 1 deletion examples/07_transformations/01_binary.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,8 @@
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.
See :any:`transform.binary`
"""
import gstools as gs

Expand All @@ -13,5 +15,5 @@
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.transform("binary")
srf.plot()
40 changes: 23 additions & 17 deletions examples/07_transformations/02_discrete.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,32 +6,38 @@
If we do not give thresholds, the pairwise means of the given
values are taken as thresholds.
If thresholds are given, arbitrary values can be applied to the field.
See :any:`transform.discrete`
"""
import numpy as np
import gstools as gs

# structured field with a size of 100x100 and a grid-size of 0.5x0.5
# Structured field with a size of 100x100 and a grid-size of 0.5x0.5
x = y = np.arange(200) * 0.5
model = gs.Gaussian(dim=2, var=1, len_scale=5)
srf = gs.SRF(model, seed=20170519)

# create 5 equidistanly spaced values, thresholds are the arithmetic means
srf.structured([x, y])
discrete_values = np.linspace(np.min(srf.field), np.max(srf.field), 5)
gs.transform.discrete(srf, discrete_values)
srf.plot()

# calculate thresholds for equal shares
###############################################################################
# Create 5 equidistanly spaced values, thresholds are the arithmetic means

values1 = np.linspace(np.min(srf.field), np.max(srf.field), 5)
srf.transform("discrete", store="f1", values=values1)
srf.plot("f1")

###############################################################################
# Calculate thresholds for equal shares
# but apply different values to the separated classes
discrete_values2 = [0, -1, 2, -3, 4]
srf.structured([x, y])
gs.transform.discrete(srf, discrete_values2, thresholds="equal")
srf.plot()

# user defined thresholds
values2 = [0, -1, 2, -3, 4]
srf.transform("discrete", store="f2", values=values2, thresholds="equal")
srf.plot("f2")

###############################################################################
# Create user defined thresholds
# and apply different values to the separated classes

values3 = [0, 1, 10]
thresholds = [-1, 1]
# apply different values to the separated classes
discrete_values3 = [0, 1, 10]
srf.structured([x, y])
gs.transform.discrete(srf, discrete_values3, thresholds=thresholds)
srf.plot()
srf.transform("discrete", store="f3", values=values3, thresholds=thresholds)
srf.plot("f3")
4 changes: 3 additions & 1 deletion examples/07_transformations/03_zinn_harvey.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,8 @@
`Zinn & Harvey (2003) <https://www.researchgate.net/publication/282442995_zinnharvey2003>`__.
With this transformation, one could overcome the restriction that in ordinary
Gaussian random fields the mean values are the ones being the most connected.
See :any:`transform.zinnharvey`
"""
import gstools as gs

Expand All @@ -14,5 +16,5 @@
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.transform("zinnharvey", conn="high")
srf.plot()
8 changes: 5 additions & 3 deletions examples/07_transformations/04_bimodal.py
Original file line number Diff line number Diff line change
@@ -1,13 +1,15 @@
"""
bimodal fields
---------------
Bimodal fields
--------------
We provide two transformations to obtain bimodal distributions:
* `arcsin <https://en.wikipedia.org/wiki/Arcsine_distribution>`__.
* `uquad <https://en.wikipedia.org/wiki/U-quadratic_distribution>`__.
Both transformations will preserve the mean and variance of the given field by default.
See: :any:`transform.normal_to_arcsin` and :any:`transform.normal_to_uquad`
"""
import gstools as gs

Expand All @@ -16,5 +18,5 @@
model = gs.Gaussian(dim=2, var=1, len_scale=10)
srf = gs.SRF(model, seed=20170519)
field = srf.structured([x, y])
gs.transform.normal_to_arcsin(srf)
srf.transform("normal_to_arcsin") # also "arcsin" works
srf.plot()
21 changes: 16 additions & 5 deletions examples/07_transformations/05_combinations.py
Original file line number Diff line number Diff line change
Expand Up @@ -9,21 +9,32 @@
Then we apply the Zinn & Harvey transformation to connect the low values.
Afterwards the field is transformed to a binary field and last but not least,
we transform it to log-values.
We can select the desired field by its name and we can define an output name
to store the field.
If you don't specify `field` and `store` everything happens inplace.
"""
# sphinx_gallery_thumbnail_number = 1
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, mean=-9, seed=20170519)
srf.structured([x, y])
gs.transform.normal_force_moments(srf)
gs.transform.zinnharvey(srf, conn="low")
gs.transform.binary(srf)
gs.transform.normal_to_lognormal(srf)
srf.plot()
srf.transform("force_moments", field="field", store="f_forced")
srf.transform("zinnharvey", field="f_forced", store="f_zinnharvey", conn="low")
srf.transform("binary", field="f_zinnharvey", store="f_binary")
srf.transform("lognormal", field="f_binary", store="f_result")
srf.plot(field="f_result")

###############################################################################
# The resulting field could be interpreted as a transmissivity field, where
# the values of low permeability are the ones being the most connected
# and only two kinds of soil exist.
#
# All stored fields can be accessed and plotted by name:

print("Max binary value:", srf.f_binary.max())
srf.plot(field="f_zinnharvey")
21 changes: 17 additions & 4 deletions examples/07_transformations/README.rst
Original file line number Diff line number Diff line change
Expand Up @@ -20,18 +20,31 @@ common transformations:
normal_to_uniform
normal_to_arcsin
normal_to_uquad
apply_function


All the transformations take a field class, that holds a generated field,
as input and will manipulate this field inplace.
as input and will manipulate this field inplace or store it with a given name.

Simply import the transform submodule and apply a transformation to the srf class:
Simply apply a transformation to a field class:

.. code-block:: python
from gstools import transform as tf
import gstools as gs
...
tf.normal_to_lognormal(srf)
srf = gs.SRF(model)
srf(...)
gs.transform.normal_to_lognormal(srf)
Or use the provided wrapper:

.. code-block:: python
import gstools as gs
...
srf = gs.SRF(model)
srf(...)
srf.transform("lognormal")
Examples
--------
13 changes: 7 additions & 6 deletions examples/08_geo_coordinates/01_dwd_krige.py
Original file line number Diff line number Diff line change
Expand Up @@ -130,8 +130,9 @@ def north_south_drift(lat, lon):
g_lat = np.arange(47, 56.1, 0.1)
g_lon = np.arange(5, 16.1, 0.1)

field, k_var = uk((g_lat, g_lon), mesh_type="structured")
mean = uk((g_lat, g_lon), mesh_type="structured", only_mean=True)
uk.set_pos((g_lat, g_lon), mesh_type="structured")
uk(return_var=False, store="temp_field")
uk(only_mean=True, store="mean_field")

###############################################################################
# And that's it. Now let's have a look at the generated field and the input
Expand All @@ -140,8 +141,8 @@ def north_south_drift(lat, lon):
levels = np.linspace(5, 23, 64)
fig, ax = plt.subplots(1, 3, figsize=[10, 5], sharey=True)
sca = ax[0].scatter(lon, lat, c=temp, vmin=5, vmax=23, cmap="coolwarm")
co1 = ax[1].contourf(g_lon, g_lat, field, levels, cmap="coolwarm")
co2 = ax[2].contourf(g_lon, g_lat, mean, levels, cmap="coolwarm")
co1 = ax[1].contourf(g_lon, g_lat, uk["temp_field"], levels, cmap="coolwarm")
co2 = ax[2].contourf(g_lon, g_lat, uk["mean_field"], levels, cmap="coolwarm")

[ax[i].plot(border[:, 0], border[:, 1], color="k") for i in range(3)]
[ax[i].set_xlim([5, 16]) for i in range(3)]
Expand All @@ -160,8 +161,8 @@ def north_south_drift(lat, lon):
# a look at a cross-section at a longitude of 10 degree:

fig, ax = plt.subplots()
ax.plot(g_lat, field[:, 50], label="Interpolated temperature")
ax.plot(g_lat, mean[:, 50], label="North-South mean drift")
ax.plot(g_lat, uk["temp_field"][:, 50], label="Interpolated temperature")
ax.plot(g_lat, uk["mean_field"][:, 50], label="North-South mean drift")
ax.set_xlabel("Lat in deg")
ax.set_ylabel("T in [°C]")
ax.set_title("North-South cross-section at 10°")
Expand Down
Loading

0 comments on commit f71821b

Please sign in to comment.