GeoStat-Framework / GSTools

GSTools - A geostatistical toolbox: random fields, variogram estimation, covariance models, kriging and much more
https://geostat-framework.org
GNU Lesser General Public License v3.0
544 stars 71 forks source link

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

Closed rkenko closed 1 year ago

rkenko commented 1 year ago

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 ?