"""
Creating PGS with periodic boundaries
-------------------------------------

Plurigaussian fields with periodic boundaries (P-PGS) are used in various
applications, including the simulation of interactions between the landsurface
and the atmosphere, as well as the application of homogenisation theory to
porous media, e.g. [Ricketts 2024](https://doi.org/10.1007/s11242-024-02074-z).

In this example we will use GSTools's Fourier generator to create periodic
random fields, which can in turn be used to generate P-PGS.
"""

import matplotlib.pyplot as plt
import numpy as np

import gstools as gs

dim = 2
# define the spatial grid, see the periodic random field [examples](https://geostat-framework.readthedocs.io/projects/gstools/en/latest/examples/01_random_field/08_fourier.html)
# for details.

# domain size and periodicity
lithotypes = 200
# no. of cells in both dimensions
N = [170, 153]

x = np.linspace(0, lithotypes, N[0], endpoint=False)
y = np.linspace(0, lithotypes, N[1], endpoint=False)

###############################################################################
# The parameters of the covariance model are very similar to previous examples.
# The interesting part is the SRF class. We set the `generator` to `"Fourier"`,
# which inherently generates periodic SRFs. The Fourier generator needs an
# extra parameter `period` which defines the periodicity.

model = gs.Gaussian(dim=dim, var=0.8, len_scale=40)
srf = gs.SRF(model, generator="Fourier", period=lithotypes)
field1 = srf.structured([x, y], seed=19770319)
field2 = srf.structured([x, y], seed=19860912)

###############################################################################
# Very similar to previous examples, we create a simple lithotypes field.

M = [200, 160]

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

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

###############################################################################
# With the two SRFs and the lithotypes ready, we can create our first P-PGS.

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

###############################################################################
# Finally, we can plot the PGS, but we will also show the lithotypes and the
# two original periodic Gaussian fields. Especially with `field1` you can
# nicely see the periodic structures in the black structure in the upper right
# corner. This transfers to the P-PGS, where you can see that the structures
# seemlessly match the opposite boundaries.

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")

plt.show()