3D non-stationary lognormal distribution

[YuanZhongYan]

Hello,

This is a wonderful Python module; thank you so much to the developer.

I am new to geostatistics and want to generate a 3D non-stationary lognormal distribution field. However, I have encountered a problem applying the trend function (e.g., a linear function). I am unsure whether and how to transform the trend function in the lognormal field into the form of a normal distribution field. Below is my code:

import numpy as np
import gstools as gs

# Mean and standard deviation values in lognormal distribution
mu_log = 1*10**(-5)
sigma_log = 5*10**(-6)
cov = sigma_log/mu_log

# Mean and standard deviation values in normal distribution
mu_nor = np.log(mu_log)-0.5*sigma_log**2
sigma_nor = np.log(1+cov**2)
dk = 0.5

def trend(x, y, z):
    return (1 + dk)*mu_log +2*dk*mu_log*(z-30)/30 # This is a linear trend function in the lognormal field.

x = np.arange(-7.8, 7.8, 0.3)
y = np.arange(-7.8, 7.8, 0.3)
z = np.arange(0, 30.3, 0.3)


pos = gs.generate_grid([x, y, z])

for i in range(5):
    x_coords = pos[0, :]
    y_coords = pos[1, :]
    z_coords = pos[2, :]
    model = gs.Exponential(dim=3, var=sigma_nor, len_scale=[4,4,2])
    srf = gs.SRF(model, mean=mu_nor, normalizer=gs.normalizer.LogNormal)
    field = srf.structured([x,y,z])
    field_log_norm = srf.field

[LSchueler]

Hi YuanZhongYan,

thank you for the kind words.

I cannot really follow your example. For example, what is it with that loop near the bottom of your skript?
But in case I understood your question correctly, you could try something like this

x_u, y_u, z_u = np.meshgrid(x, y, z, indexing='ij')
tr = trend(x_u, y_u, z_u)

And after generating the field do

srf.field += tr

[YuanZhongYan]

The loop generates 5 lognormal distribution fields. Sorry for the ambiguity.

Anyway, thank you very much! The problem is solved with your reminder.

[MuellerSeb]

Hey there,

GSTools provides trends as a built-in feature. You can pass your function to the SRF class and it will be applied inplace.
Also, to generate ensembles of fields, you can reuse the same srf with a new seed, and store the fields with a specific name to access them afterwads. I modified your example to demonstrate this. (Note, that I corrected an error in your calculations for the mean and variance. sigma_nor as you calculated it, is actually the variance and should be used to calculate, mu_nor. I renamed both to be clear):

import numpy as np
import gstools as gs

# Mean and standard deviation values in lognormal distribution
mu_log = 1 * 10 ** (-5)
sigma_log = 5 * 10 ** (-6)
cov = sigma_log / mu_log

# Mean and variance values in normal distribution
var_nor = np.log(1 + cov**2)
mean_nor = np.log(mu_log) - 0.5 * var_nor
dk = 0.5


def trend(x, y, z):
    """This is a linear trend function in the lognormal field."""
    return (1 + dk) * mu_log + 2 * dk * mu_log * (z - 30) / 30


x = np.arange(-7.8, 7.8, 0.3)
y = np.arange(-7.8, 7.8, 0.3)
z = np.arange(0, 30.3, 0.3)

model = gs.Exponential(dim=3, var=var_nor, len_scale=[4, 4, 2])
srf = gs.SRF(model, mean=mean_nor, normalizer=gs.normalizer.LogNormal, trend=trend)
srf.set_pos([x, y, z], "structured")

# generate ensemble
for i in range(5):
    srf(seed=i, store=f"field_{i}")

# access fields by name
field = srf["field_0"]

# plot fields by name
srf.plot(field="field_0")

• See this example for details about ensemble generation.
• See this section for details about mean, normalization and trend.

Hope this helps!
Sebastian

[YuanZhongYan]

Thank you for your kind reply. I still have a question:

The trend function (1 + dk) * mu_log + 2 * dk * mu_log * (z - 30) / 30 is a trend function in the lognormal field. Does it need to be transformed to the form in the normal field?

[MuellerSeb]

Hey there,

we could have been a bit more detailed in the documentation about mean, normalization and trend to clarify questions like this.

Here a short explanation how these three things interact:

Field generation in GSTools starts with a raw field that is generated as a normal distributed field with zero mean. Now, you can apply a mean, a normalizer and/or a trend when specified.

The mean is added to the raw field and can be non-stationary, which means, it can be a callable function depending on the spatial location.

After that the field is denormalized to generate the target distribution if a normalizer is given (in your case the log-normal normalizer, where the denormalizer function is the classic exponential function).

If a trend is given (either a fix value or again a callable function), this trend is then added (!) to the denormalized field:

field = trend + denormalizer(raw_field + mean)
# or with position dependent trend and mean
field = trend(pos) + denormalizer(raw_field + mean(pos))

I don't know exactly what you want to achieve with your trend function, but I hope you now understand that this trend function is added to the log-normal field. Since you use mu_log within the function, I am not sure if the function really does what you expect, since, as I just wrote, the log-normal field already applied the given mean.

If you want a spatial varying geometric mean, you could try using a callable mean function, that is added to the raw normal field.
My gut-feeling tells me, this could be your aim, but you would need to replace mu_log with mean_nor.

Hope this helps.
Sebastian

[YuanZhongYan]

Thank you very much for your explicit explanation of my questions.
Thanks again!