Is it possible to fit a variogram on 2D image and generate an equivalent 3D field using gstools?

[rkenko]

Hi,

I am new in geostats and i wonder if it is possible to fit a variogram on a 2D field then generate a 3D equivalent field.

This is the code i am using:

#importing the image from which we seek the variogram 
image       = plt.imread("field.png")[:,:,0]
# spreading values of imported image between ]-inf;+inf[
image       = np.tan(((10/np.pi)*(image))+(np.pi/2)) 


# estimating empirical variogram
gamma = gs.vario_estimate_structured(image)
bin_center = np.arange(len(gamma))


# fiiting the empirical variogram at small range, 128 px
bins        = np.arange(128)
fit_model   = gs.Gaussian(dim=2)
fit_model.fit_variogram(bin_center[:128], gamma[:128], nugget=False)
# printing all the parameters of the covariance model fitted
print(fit_model)


# generating 3D field from fitted variogram with size 256
x  = y  = z = range(256)
model       = gs.Gaussian(dim=3, var=fit_model.var, len_scale=fit_model.len_scale)
# generating the spatial random field
srf         = gs.SRF(model)
srf((x, y, z), mesh_type = 'structured')
srf         = srf.field


# passing the field through atan in order to spread distribution between [0;1]
srf         = (np.pi/10)*(np.arctan(srf)-(np.pi/2))+1
print(np.amin(srf))
print(np.amax(srf))
plt.imsave("x_middle_slice.png",srf[:,:,128],cmap="gray",vmin=0,vmax=1)
# saving to a VTK file for visualization
srf         = pv.wrap(srf)
srf.save('gs.vtk')

My question is: am i doing right ? the purpose of the tanh is to ensure values stays in the range [0,1] since they represents pixels values in a picture. without the tanh, the field isn't realistic in terms of colors.

I observed the slices of the generated 3D volume dont have the same variogram than the 2D picture of the begining. I expected them to be the same... Is there a way to correct my code please ?

[MuellerSeb]

Hi there,

This looks right to me. I would check for normality of the data after the tan transformation since we need a normal distribution to apply this workflow.

You could adopt two things:

1. I would use vario_estimate and not vario_estimate_structured to get an omnidirectional estimation for the empirical variogram:

   x = y = np.arange(image.shape[0]), np.arange(image.shape[1])
   bin_center, gamma = gs.vario_estimate((x,y), image, mesh_type="structured")

2. You can use the 3D model directly to fit the omnidirectional variogram

   model = gs.Gaussian(dim=3)
   model.fit_variogram(bin_center, gamma, nugget=False)

The rest looks good to me 😉

Cheers,
Sebastian

[MuellerSeb]

For the technical question: the mesh_type is crucial:

bin_center, gamma = gs.vario_estimate((x,y), image, mesh_type="structured")

1. theoretical question:
   Just look at the histogram. It should just look rather normal (not skewed or something). It doesn't need to be standardized.

2. theoretical question: Have a look at the JBessel covariance model. That is a holemodel and should be what you are looking for: https://geostat-framework.readthedocs.io/projects/gstools/en/v1.4.0/generated/gstools.covmodel.JBessel.html

Cheers, Sebastian

[rkenko]

Hi @MuellerSeb, thank you for yout time.

I have one last remark. The code is not working using:

x = y = image.shape[0], image.shape[1]
bin_center, gamma = gs.vario_estimate((x,y), image, mesh_type="structured")

Here is the error:

Traceback (most recent call last):
  File "/home/generate_grf.py", line 96, in <module>
    bin_center, gamma = gs.vario_estimate((x, y), image, mesh_type="structured")
  File "/home/anaconda3/lib/python3.9/site-packages/gstools/variogram/variogram.py", line 270, in vario_estimate
    pos, __, dim = format_struct_pos_shape(
  File "/home/anaconda3/lib/python3.9/site-packages/gstools/tools/geometric.py", line 513, in format_struct_pos_shape
    raise ValueError("Formatting: position tuple doesn't match dimension.")
ValueError: Formatting: position tuple doesn't match dimension.

The code works using this line. But doing this i wonder if am only giving the values of the diagonal pixels.

x , y = np.arange(image.shape[0]), np.arange(image.shape[1])
bin_center, gamma = gs.vario_estimate((x, y), image,bins)

I want to do something like this, but it don't work...

nb_points = 10000
maxy, maxx  = image.shape
x, y = np.random.randint(0, maxx, nb_points), np.random.randint(0, maxy, nb_points)
bin_center, gamma = gs.vario_estimate((x,y), image, mesh_type="structured")

[MuellerSeb]

Sorry, there was a mistake on my side. You correctly put np.arange around the image shapes. But you still need to put mesh_type="structured" in the call signature. You don't need bins (but you can give bin_edges of course), since they are calculated automatically.
If you want to have a sampled variogram (if I interpret your np.random.randint(0, maxx, nb_points) correctly), you can do that directly in the estimation routine. So the call could look like this:

x , y = np.arange(image.shape[0]), np.arange(image.shape[1])
bin_center, gamma = gs.vario_estimate((x, y), image, mesh_type="structured", sampling_size=10000)

[rkenko]

yes, i copy-pasted the wrong code, i used mesh_type="structured" for this test. Thanks a lot

[MuellerSeb]

One side-note: take a look at the image, since I think the image is not given in xy coordinates but in matrix style ij. So it could be, you need to transpose the image after loading it

image = image.T