How to implement kriging in Matlab?

I need to perform kriging in Matlab based on three borehole samples.

Numbers for soil layers indicate depth from ground zero, boreholes - position on model.

The outcome should look something like this (already done that in Zsoil):

Any ideas how to perform this?

[Image Analyst]

Are you using commas for decimal points instead of periods? Like (x,y) = 7,5 -- okay that means commas are commas. But SoilLevel1 = (3, 27) - what the heck does that mean? Is the value 3.27? Or do you have two values: 3 and 37? And if two values, what does each value represent?

[KarolN]

No, these are commas all right. 3,27 is one value. It means soil layer 1 ENDS at depth of 3,27 meters, as measured from the surface level.

Similarly soil layer 3 ends at 8,95 meters and so on.

Answer 1

[the cyclist]

Googling "MATLAB" and "kriging" turned up this submission from the File Exchange.

[KarolN]

Looks OK but I wonder how to write data from my table into these inputs:

"vstruct structure array with variogram information as returnedvariogramfit (forth output argument)

x,y coordinates of observations

z values of observations

xi,yi coordinates of locations for predictions

chunksize nr of elements in zi that are processed at one time.

The default is 100, but this depends largely on your

available main memory and numel(x)."

Do you have any hint how to start with that? Thanks

[the cyclist]

@KarolN wrote the following as an "answer", but I am moving it here, because it was actually a comment on my answer:

I'd like to point out that https://www.mathworks.com/matlabcentral/fileexchange/29025-ordinary-kriging is out of date, as currently neither nargin or nargout aren't supported in the script. As for my original question, I hoped for some hints how to implement my samples, I know there are many kriging tools but I simply don't have time to rewrite their code or figure how it works.

[the cyclist]

Sorry, but I am completely ignorant of this method, so I can't really give you a better "hint" than providing you with a pointer to code where someone else solved the problem. I would have hoped that this would give you a solid headstart on writing your own.

I hadn't noticed how old the submission I mentioned was. There are other kriging examples in the File Exchange, that look more up-to-date. Maybe they will be more helpful.

[KarolN]

Thanks anyway for your effort. Maybe the thread will attract someone who had similar exercise.

Answer 2

[Image Analyst]

I believe I saw somewhere in the documentation that scatteredInterpolant did kriging. Attached is a demo.

[KarolN]

I studied the code but it's for image analysis, while I have numerical samples to implement. Anyway, thanks for effort.

[Image Analyst]

It should work. As well as griddedInterpolant(). Not sure why you couldn't adapt it to your data. Can you attach your data so we can read it in (I don't want to type in stuff from your image.)

Also, I don't understand your desired output image. What is displayed on the vertical axis and horizontal axis and how does that relate to the boreholes and soil layers?

[KarolN]

The problem is that scatteredinterpolant requries a png image, I do not see any options for inputting numerical data, but here you go my data from other kriging tool:

sample.name = 'Borehole nr 1';
sample.location = 7.5;
sample.depth = [3.27 0 8.95 12.2]

sample(2).name = 'Borehole nr 2';
sample(2).location = 17.2;
sample(2).depth = [2.52 0 8.68 12.2]

sample(3).name = 'Borehole nr 3';
sample(3).location = 30.2;
sample(3).depth = [1.53 4.58 8.47 12.2]

vstruct = sample

[KarolN]

as for my output image it is from Zsoil software and what you have there is a grid of finite elements and borehole samples and then the program interpolates data from samples onto the grid.

[Image Analyst]

It doesn't. Read the documentation - nothing in there about requiring an image or a PNG image. It just deals with 2-D data -- ANY kind of 2-D data. But I don't know what your two dimensions in your image represent because the axes are not labeled and you didn't answer my question. Again:

What is displayed on the vertical axis and horizontal axis and how does that relate to the boreholes and soil layers?

[KarolN]

Sorry about that. On the y axis you have depth of each soil layer and on the x axis you have horizontal position of them.

Ex. Layer 1 has y: 3,27 2,52 and 1,53 meters counting from ground level. X is: 7,5 17,2 and 30,2. So x,y in first sample, 1st soil layer is: (7,5 ; 3,27). Layer 2 has 0 (means is absent here) 0 and 4,58 and same x coordinates. And so on.

The meaning of the exercise is an interpolation of data, so you place our 3 boreholes in a fragment of the ground. Ground model is 31x13 meters. I am supposed to interpolate data from samples, so I get an approximate picture of how soil layers may be placed in our ground model.

Basically we are speaking of an 31x13 image, which is visible for us only in 3 thin strips - and we have to reconstruct the whole of it based on those 3 strips.

[Image Analyst]

So is this what you have:

sample.name = 'Borehole nr 1';
sample.location = 7.5;
sample.depth = [3.27 0 8.95 12.2]
sample(2).name = 'Borehole nr 2';
sample(2).location = 17.2;
sample(2).depth = [2.52 0 8.68 12.2]
sample(3).name = 'Borehole nr 3';
sample(3).location = 30.2;
sample(3).depth = [1.53 4.58 8.47 12.2]
x = [sample.location]
soilLevels = vertcat(sample.depth)'
for k = 1 : size(soilLevels, 1)
    thisLevel = soilLevels(k, :);
    plot(x, thisLevel, '.-', 'LineWidth', 2, 'MarkerSize', 20)
    hold on
end
grid on
xlabel('x distance', 'fontSize', fontSize)
ylabel('Depth of top of layer', 'fontSize', fontSize)
title('Soil depths as function of location', 'fontSize', fontSize)
legend('Soil 1', 'Soil 2', 'Soil 3', 'Soil 4', 'location', 'west')

I just don't see how it makes sense to interpolate. Like what would be in the graph for x=17.2 between 2.52 and 8.68? Everything between the blue line and the yellow line will be soil 1 (the blue) right? Because from the top of soil 1 to the top of the next soil is all soil 1, right?

[KarolN]

Yes, between 2.52 and 8.68 there will be layer 1. Well from samples you can approximate the soil composition. I don't see the sense in interpolating it from just 3 samples, for area that size you would need 20-30 to get it faithfully, but it is a student exercise so I guess they tried to make it easy for us.

Ideally, kriging in matlab should get us result like that:

But I as a beginner I am not versed enough to code that properly.

If I got you interested in kriging models, you can download ooDace Toolbox and try it for yourself. My skills are too low for it.

[Wolfgang Schwanghart]

Actually, I don't understand why you need kriging here. The figure created by Zsoil shows a linear interpolation of the individual depths above which soils are classified into the respective soil layers. Kriging won't help here, inparticular since you would need to have a variogram.

[KarolN]

I do not deny what you say is right. But student exercises do not always make sense. I am supposed to use kriging in matlab, not only zsoil...