Is it possible with the current implementation to conduct Kriging for unstructured grid or a set of (x, y) points?

[flydream0428]

In the example/tutorial notebook, Kriging is always done for a structured regular-spaced grid. I am wondering if it is possible to only conduct kriging for a list of (x, y) points? Thanks!

[PauloCarvalhoRJ]

Hello,

Depending on your application, you can model the variogram and perform kriging in depositional regular space (UVW) instead of Cartesian space (XYZ). The figure below illustrates an irregular grid representing some folded geology (A) and the corresponding regular grid in UVW space (B):

It's easy to follow but there's a caveat: you have to unfold the primary data (e.g. well or borehole data), as illustrated in the figure below:

So, it would be necessary for GeostatsPy to include an unfolding routine as described here: https://github.com/PauloCarvalhoRJ/gammaray/blob/master/docs/GammaRayManual.docx (Section 15.3) and implemented in C++ here: https://github.com/PauloCarvalhoRJ/gammaray/blob/master/domain/geogrid.cpp (look for function GeoGrid::unfold()). Of course, there are other unfolding algorithms out there. It's likely there is one already implemented in Python.

regards,

Paulo

[PauloCarvalhoRJ]

The figure below illustrates how to compute 2D stratigraphic coordinates of a sample point P with respect to an irregular cell of a 2D unstructured surface:

source: https://math.stackexchange.com/questions/13404/mapping-irregular-quadrilateral-to-a-rectangle

This computation is used to compute texture coordinates in computer graphics. So, I extended it to 3D to compute stratigraphic coordinates in irregular meshes (reference to the code posted above).

[flydream0428]

@PauloCarvalhoRJ Thanks a lot for your reply, Paulo. I may be wrong here, but isn't easier to just add a function to allow Kriging on a list of points based on their locations rather than a structured grid? When we supply a structured grid to GSLIB, my understanding is it will calculate the coordinates internally anyway. The function can return the results as a list rather than a matrix.

[PauloCarvalhoRJ]

One of the reasons behind using a regular grid is efficiency. But if you wish to krige over a point set, you can study the kb2d.py's kb2d() function code to create a new function intended for point sets. Here's the loop over the grid blocks:

GeostatsPy/working/kb2d.py

Lines 109 to 239 in a256e1b

	# MAIN LOOP OVER ALL THE BLOCKS IN THE GRID:
	nk = 0
	ak = 0.0
	vk = 0.0
	for iy in range(0,ny):
	yloc = ymn + (iy-0)*ysiz
	for ix in range(0,nx):
	xloc = xmn + (ix-0)*xsiz
	current_node = (yloc,xloc)
	# Find the nearest samples within each octant: First initialize
	# the counter arrays:
	na = -1 # accounting for 0 as first index
	dist.fill(1.0e+20)
	nums.fill(-1)
	dist, nums = tree.query(current_node,ndmax) # use kd tree for fast nearest data search
	na = len(dist) - 1
	# Is there enough samples?
	if na + 1 < ndmin: # accounting for min index of 0
	est = UNEST
	estv = UNEST
	print('UNEST at ' + str(ix) + ',' + str(iy))
	else:
	# Put coordinates and values of neighborhood samples into xa,ya,vra:
	for ia in range(0,na+1):
	jj = int(nums[ia])
	xa[ia] = x[jj]
	ya[ia] = y[jj]
	vra[ia] = vr[jj]
	# Handle the situation of only one sample:
	if na == 0: # accounting for min index of 0 - one sample case na = 0
	cb1 = cova2(xa[0],ya[0],xa[0],ya[0],nst,c0,PMX,cc,aa,it,ang,anis,rotmat,maxcov)
	xx = xa[0] - xloc
	yy = ya[0] - yloc
	# Establish Right Hand Side Covariance:
	if ndb <= 1:
	cb = cova2(xx,yy,xdb[0],ydb[0],nst,c0,PMX,cc,aa,it,ang,anis,rotmat,maxcov)
	else:
	cb = 0.0
	for i in range(0,ndb):
	cb = cb + cova2(xx,yy,xdb[i],ydb[i],nst,c0,PMX,cc,aa,it,ang,anis,rotmat,maxcov)
	dx = xx - xdb(i)
	dy = yy - ydb(i)
	if (dx*dx+dy*dy) < EPSLON:
	cb = cb - c0
	cb = cb / real(ndb)
	if ktype == 0:
	s[0] = cb/cbb
	est = s[0]*vra[0] + (1.0-s[0])*skmean
	estv = cbb - s[0] * cb
	else:
	est = vra[0]
	estv = cbb - 2.0*cb + cb1
	else:
	# Solve the Kriging System with more than one sample:
	neq = na + 1 + ktype # accounting for first index of 0
	nn = (neq + 1)*neq/2
	# Set up kriging matrices:
	iin=-1 # accounting for first index of 0
	for j in range(0,na+1):
	# Establish Left Hand Side Covariance Matrix:
	for i in range(0,na+1): # was j - want full matrix
	iin = iin + 1
	a[iin] = cova2(xa[i],ya[i],xa[j],ya[j],nst,c0,PMX,cc,aa,it,ang,anis,rotmat,maxcov)
	xx = xa[j] - xloc
	yy = ya[j] - yloc
	# Establish Right Hand Side Covariance:
	if ndb <= 1:
	cb = cova2(xx,yy,xdb[0],ydb[0],nst,c0,PMX,cc,aa,it,ang,anis,rotmat,maxcov)
	else:
	cb = 0.0
	for j1 in range(0,ndb):
	cb = cb + cova2(xx,yy,xdb[j1],ydb[j1],nst,c0,PMX,cc,aa,it,ang,anis,rotmat,maxcov)
	dx = xx - xdb[j1]
	dy = yy - ydb[j1]
	if (dx*dx+dy*dy) < EPSLON:
	cb = cb - c0
	cb = cb / real(ndb)
	r[j] = cb
	rr[j] = r[j]
	# Set the unbiasedness constraint:
	if ktype == 1:
	for i in range(0,na+1):
	iin = iin + 1
	a[iin] = unbias
	iin = iin + 1
	a[iin] = 0.0
	r[neq] = unbias
	rr[neq] = r[neq]
	# Solve the Kriging System:
	s = ksol_numpy(neq,a,r)
	ising = 0 # need to figure this out
	# Write a warning if the matrix is singular:
	if ising != 0:
	print('WARNING KB2D: singular matrix')
	print(' for block' + str(ix) + ',' + str(iy)+ ' ')
	est = UNEST
	estv = UNEST
	else:
	# Compute the estimate and the kriging variance:
	est = 0.0
	estv = cbb
	sumw = 0.0
	if ktype == 1:
	estv = estv - real(s[na+1])*unbias
	for i in range(0,na+1):
	sumw = sumw + s[i]
	est = est + s[i]*vra[i]
	estv = estv - s[i]*rr[i]
	if ktype == 0:
	est = est + (1.0-sumw)*skmean
	kmap[ny-iy-1,ix] = est
	vmap[ny-iy-1,ix] = estv
	if est > UNEST:
	nk = nk + 1
	ak = ak + est
	vk = vk + est*est
	# END OF MAIN LOOP OVER ALL THE BLOCKS:

. You'd need to modify that loop to read the explicit XYZ coordinates coming from the point set instead of computing them from the grid specs. You also have to modify the part dealing with the block variance, afterall you'd be kriging on points.