Kriging within an irregularly shaped orebody

Any question regarding the Interpolation method using Kriging

Kriging within an irregularly shaped orebody

Postby Gabriel_W-H » Sat Aug 25, 2018 8:03 pm


As a relatively inexperienced user of R and RGeostats* I have been reluctant to post this problem which I assume quite simple and more down to ignorance.

I am trying to krige [2D] within an orebody, the extent of which is delineated by the drill holes plotted as points on my grid in the below kriging interpolation.


Code: Select all
neigh=neigh.create(ndim = 2, type = 0),nodes=c(200,150))


# Display Kriging Estimates #
quilt.plot(OK_Zn[,2], OK_Zn[,3], OK_Zn[,4],nx =200, ny=150,col=tim.colors(256),xlab = "Easting(m)", ylab = "Northing(m)",main="Ordinary Kriging Estimates(Isotropic Model) [Zn %]",asp=1)

# Display Kriging Standard Deviations#
quilt.plot(OK_Zn[,2], OK_Zn[,3], OK_Zn[,5],nx =200, ny=150,col=tim.colors(256),xlab = "Easting(m)", ylab = "Northing(m)",main="Ordinary Kriging Standard Deviations (Isotropic Model) [Zn]",asp=1)

As you can see I have been using the function db.grid.init to build a grid from my drill holes on which to krige. However I cannot force the function, using selections or any other modification to only krige within the dimensions of the drill holes, rather I krige to the extent of my grid which is a rectangle.

I learned RGeoStats with the Meuse dataset, for which we relied on the meuse_grid.csv file to krige within the desired area by setting up a target db with the code:

To krige in this manner I first have to create a grid csv from the datapoints. Unfortunately I have failed to do this using other R functions so I am back with RGeostats wondering if there is a simple solution to this problem using existing functionality to which I am ignorant.

Thanks for taking the time to read this, and for the other posts on these forums that have helped me troubleshoot plenty of other problems on the way :)
Last edited by Gabriel_W-H on Tue Aug 28, 2018 3:54 pm, edited 3 times in total.
Posts: 3
Joined: Wed Aug 15, 2018 1:13 pm

Re: Kriging within an irregularly shaped orebody

Postby Gabriel_W-H » Sat Aug 25, 2018 8:09 pm

Sorry I'm not sure I used the Img function correctly, instead I've attached it.
Kriging estimate.png
Kriging estimate.png (227.49 KiB) Viewed 715 times
Posts: 3
Joined: Wed Aug 15, 2018 1:13 pm

Re: Kriging within an irregularly shaped orebody

Postby Didier Renard » Tue Aug 28, 2018 1:35 pm


We just received your post and did not have time to prepare the answer earlier...
Your problem has a (if not many) simple solution(s).
However, before starting, I want to make sure that when you mention RGeos (twice in your post), you mean "RGeostats". As a matter of fact, another package exists, called RGeoS, for which I won't be able to answer.

Within RGeostats, you want to perform an interpolation / mapping (by kriging) at the nodes of a regular grid, in an area around your boreholes.

First you read your boreholes correctly and put them in a RGeostats data base (let us call it db.well). Then you defined a grid by using db.grid.init(): this function is meant to let the system create a grid automatically overlaying the db.well (provided as argument). There is another way to create it (db.create) where you clearly define the origin, the mesh and the number of meshes along each space direction. A third solution is obviously (as you mentioned) to import it from an already existing CSV file (but this CSV file is usually the result of a dump from RGeostats as its format is quite particular).

This grid stands as a data base (or the Db class again) which will serve in collecting the results of the interpolation procedure. The grids are always considered as rectangular (or square). Then, if you are only interested in a limited portion of it, you have two possibilities:
1) create a selection on the grid. This selection is created once for all, stored as a standard variable in the grid db (although with a specific locator named "sel" which tells the system about its status), and remains active until canceled or deleted. To create this secltion, you have several possibilities again:
1.a) Graphically: by digitizing a polygon on the figure showing your boreholes for example (polygon.digit) and applying it to select the set of active grid nodes (db.polygon)
1.b) Calculating the distance of each grid node to all data and applying a threshold
1.c) to calculate the convex hull (or its dilated version) around the data (db.selhull)
1.d) and probably many others

2) to choose the nodes that will be kriged on the fly. This can be done by designing a (moving) neighborhood. The neighborhood is centered it on each target grid node in turn. You can dimension its maximum radius so that only the boreholes located within this circle will be used. Furthermore, you can tell the system that if there are not at least XXX samples in this circle, the kriging must not take place. This will naturally create a selection of the only estimated grid nodes located close to the borehole locations.

Hope this will help.
PS: The way you posted your question was absolutely correct. Unfortunately, it may take time before getting a response.
Didier Renard
Posts: 294
Joined: Thu Sep 20, 2012 4:22 pm

Re: Kriging within an irregularly shaped orebody

Postby Gabriel_W-H » Thu Aug 30, 2018 10:12 am


Thank you for the detailed response.
I was indeed referring to RGeostats, sorry for any confusion and I've edited my original post.

I couldn't get the selection method to work, but option #2 kriging with a moving neighbourhood was relatively successful.
My follow on question is more resource evaluation applicable: how to extract the average grade over an area of the kriged domain.

E.g. I could use the polygon.digit function to manually draw the polygon, or the threshold function (although I couldn't code the syntax properly). How would one apply the selection to the kriged surface and which function is use to extract summary information from these datapoints?

Thanks again,
Posts: 3
Joined: Wed Aug 15, 2018 1:13 pm

Re: Kriging within an irregularly shaped orebody

Postby Didier Renard » Sun Apr 07, 2019 7:26 pm

I answer to your question by constructing a small test. The data set is constituted of 200 samples generated randomly in a square area (with edge equal to 1 unit) as presented in the next graphic:

Code: Select all
nech = 200
data = db.create(x1=runif(nech),x2=runif(nech),z1=rnorm(nech))

data.png (14.3 KiB) Viewed 437 times

The next step consists in digitizing a polygon (from the previous display) and in saving the resulting polygon as a R object (called poly):

Code: Select all
poly = polygon.digit()

The next operation consists in using this polygon to select the samples located within the polygon (called active). The next figure displays the active samples together with the polygon:

Code: Select all
data = db.polygon(data,poly)
plot(data,pch=19,title="Data with Polygon")

Data with Polygon
data-poly.png (11.53 KiB) Viewed 437 times

We create a fine grid covering the area of interest, with a 100 x 100 meshes (sqaure with a side equal to 0.01). We use the polygon to restrain the grid to the only active cells located within the polygon. The restricted grid is represented together with the data and the polygon.

Code: Select all
grid = db.create(nx=c(100,100),dx=c(0.01,0.01))
grid = db.polygon(grid,poly)

Grid and Data and Polygon
grid-data-poly.png (40.81 KiB) Viewed 437 times

In order to continue, we must introduce:
- a spatial model: it is created as a cubic structure with a range of 0.5 and a sill of 3
- a neighborhood: it is created as a Unique neighborhood (all samples are involved when kriging each target node)

Code: Select all
neigh = neigh.create(ndim=2,type=0)
model = model.create(vartype="Cubic",range=0.5,sill=3)

Finally, we can perform a local estimation using Kriging, considering the only data located within the polygon (the polygon is used as a selection for the input data base) and performing the estimation at the only grid nodes located within the polygon (the polygon is used as a selection for the output data base). The result is presented next:

Code: Select all
grid = kriging(data,grid,model,neigh)

Local estimation
local-estimate.png (54.44 KiB) Viewed 437 times

Another concern may be to estimate the average value of the information contained in the polygon area. This can be obtained by using the kriging procedure through the global function:

Code: Select all
result = global(data,model,polygon=poly,dbout=grid)

The following result is produced:

Code: Select all
Global estimation kriging
Total number of data                    = 200
Number of active data                  = 60
Number of variables                     = 1
Cvv                                                  = 0.773599
Estimation by kriging                    = 0.248588
Lagrange Parameter #1               = -0.003314
Estimation St. Dev. of the mean  = 0.027993
CVgeo                                             = 0.112608
Surface                                           = 0.322828
Q (Estimation * Surface)              = 0.080251

Hoping this will help.
Didier Renard
Posts: 294
Joined: Thu Sep 20, 2012 4:22 pm

Return to Kriging

Who is online

Users browsing this forum: No registered users and 0 guests