Interpolate Gaussian random field to be used as initial condition

Hi everyone,
I’m currently working on a diffusion-reaction problem where I would like to use a Gaussian random field (GRF) as the initial condition. I’ve been using the GSTools package to generate GRFs quite easily in Python:
(GSTools Quickstart — GSTools 1.1.1 documentation).
I can generate a 2D grid that visualizes a GRF, and I’d like to interpolate this data so that the GRF can be used as an initial condition for my PDE system. But I’m not sure how to do this, as I don’t have access to a formula to do something like

u = fem.Function(V)
u.interpolate(lambda x: GRF(x)) //interpolate GRF here

Does anyone know how I can accomplish this? I tried to do some digging online and was little surprised at not being able to find examples of others using GRFs as initial conditions. Any help is greatly appreciated, thank you!

The question should probably be asked to the GSTools team, rather than the FEniCS one: “given a np.array x which contains the list of coordinates, how can I evaluate the GRF at those points?”

If I browse to the source code
https://geostat-framework.readthedocs.io/projects/gstools/en/v1.1.1/_modules/gstools/field/srf.html#SRF
I see at some point I line

        self.raw_field = self.generator.__call__(x, y, z, mesh_type)

which may be an indication, but it’s hard to give guidance on an external library I have never seen before.