Strange simulations generation

[Martin20494]

Hi team,

I am using gstools to generate simulations of random points based on variogram using Gaussian method and unconditional simulation. Below is what I did:

# Packages
import pandas as pd
import gstools as gs
from gstools.random import MasterRNG

# Read file. The file has only 3 columns - z, x, y
river_points = pd.read_csv(<path to the file>, index_col=0)

# Estimate variogram
bin_center, gamma = gs.vario_estimate((river_points.x, river_points.y), river_points.z)

# Fit variogram
var = gs.Gaussian(dim=2)

# Calculate variogram
var.fit_variogram(bin_center, gamma, nugget=False) 

This is how the original/emperical variogram looks like

import matplotlib.pyplot as plt

# Plot original variogram
ax = var.plot(x_max=max(bin_center))
ax.scatter(bin_center, gamma)

Then I generate unconditional simulations

# Create model for simulations
srf = gs.SRF(var)
srf.set_pos([river_points.x, river_points.y], 'unstructured')
ens_no = 100

# Generate simulations
seed = MasterRNG(1)

for i in range(ens_no):
     srf(seed=seed(), store=f'simulated_rivers_{i}')

And this is how I check simulated variograms

# Get lists of bin_center and gamma
bin_center_list = []
gamma_list = []

# Calculate variogram
for i in range(ens_no):
     bin_each, gamma_each = gs.vario_estimate((river_points.x, river_points.y), srf['simulated_rivers_{i}'])
     bin_center_list.append(bin_each)
     gamma_list.append(gamma_each)

# Model variogram
model_list = []

# Fit variogram
for i in range(ens_no):
    fit = gs.Gaussian(dim=2)
    fit.fit_variogram(bin_center_list[i], gamma_list[i])
    model_list.append(fit)

This is how simulated variograms look like. I thought they would be slightly different, but in fact they are extremely different (in terms of range, sill, nugget, etc.). Did I do anything wrong or is there any better way to do this? Please help me with this.

The data I used is private but I can modify a bit and send you through email. Or I think if you use other data, you can also reproduce the issue.

fig, ax = plt.subplots(figsize=(5, 5))

for i in range(ens_no):
    model_list[i].plot(ax=ax)
    
var.plot(ax=ax, linestyle='--')
ax.get_legend().remove()

[MuellerSeb]

Hi there! That looks indeed strange. But it is not surprising, that at least the sill changes, since the mean sill of your ensemble looks like it is close to the given sill of around 3.

One option could be to force the moments of the generated fields after the srf(seed=seed(), store=f'simulated_rivers_{i}') call:

srf.transform("force_moments", field=f'simulated_rivers_{i}')

Then you could also set nugget=False in the fitting as well.

How many points does your data set have?

[Martin20494]

Hi Mueller,

Thank you so much for your reply.

I have tried srf.transform("force_moments", field=f'simulated_rivers_{i}') and I actually set the nugget=False as my previous code (posted above). However, the results were slightly changed and there were also some 'vertical-line' variograms (as below). My data have 668 points.

Do you know what happened and have any suggestion on what I should do?