"""
A First and Simple Example
--------------------------

As a first example, we will create a two dimensional plurigaussian field
(PGS). Thus, we need two spatial random fields(SRF) and on top of that, we
need a field describing the categorical data and its spatial relation.
We will start off by creating the two SRFs with a Gaussian variogram, which
makes the fields nice and smooth. But before that, we will import all
necessary libraries and define a few variables, like the number of grid
cells in each dimension.
"""

import matplotlib.pyplot as plt
import numpy as np

import gstools as gs

dim = 2
# no. of cells in both dimensions
N = [180, 140]

x = np.arange(N[0])
y = np.arange(N[1])

###############################################################################
# In this first example we will use the same geostatistical parameters for
# both fields for simplicity. Thus, we can use the same SRF instance for the
# two fields.

model = gs.Gaussian(dim=dim, var=1, len_scale=10)
srf = gs.SRF(model)
field1 = srf.structured([x, y], seed=20170519)
field2 = srf.structured([x, y], seed=19970221)

###############################################################################
# Now, we will create the lithotypes field describing the categorical data. For
# now, we will only have two categories and we will address them by the
# integers 0 and 1. We start off by creating a matrix of 0s from which we will
# select a rectangle and fill that with 1s. This field does not have to match
# the shape of the SRFs.

centroid = [200, 160]

# size of the rectangle
rect = [40, 32]

lithotypes = np.zeros(centroid)
lithotypes[
    centroid[0] // 2 - rect[0] // 2 : centroid[0] // 2 + rect[0] // 2,
    centroid[1] // 2 - rect[1] // 2 : centroid[1] // 2 + rect[1] // 2,
] = 1

###############################################################################
# With the two SRFs and the L-field ready, we can create our first PGS. First, we
# set up an instance of the PGS class and then we are ready to calculate the
# field by calling the instance and handing over the L-field.

pgs = gs.PGS(dim, [field1, field2])
P = pgs(lithotypes)

###############################################################################
# Finally, we can plot the PGS, but we will also show the L-field and the two
# original Gaussian fields.

fig, axs = plt.subplots(2, 2)

axs[0, 0].imshow(field1, cmap="copper", origin="lower")
axs[0, 1].imshow(field2, cmap="copper", origin="lower")
axs[1, 0].imshow(lithotypes, cmap="copper", origin="lower")
axs[1, 1].imshow(P, cmap="copper", origin="lower")

# For more information on Plurigaussian fields and how they naturally extend
# truncated Gaussian fields, we can recommend the book
# [Plurigaussian Simulations in Geosciences](https://doi.org/10.1007/978-3-642-19607-2)