This function completes what I have started with the functions variogram and variogramfit. It is not intended to be a highly optimized code for performing kriging but may have educational value. Note that for using kriging, you should download the latest version of variogramfit.
kriging uses ordinary kriging to interpolate a variable z measured at locations with the coordinates x and y at unsampled locations xi, yi. The function requires the variable vstruct that contains all necessary information on the variogram. vstruct is the forth output argument of the function variogramfit.
This is a rudimentary, but easy to use function to perform an ordinary kriging interpolation. I call it rudimentary since it always includes ALL observations to estimate values at unsampled locations. This may not be necessary when sample locations are not within the autocorrelation range but would require something like a k nearest neighbor search algorithm or something similar. Thus, the algorithms works best for relatively small numbers of observations (100-500). For larger numbers of observations I recommend the use of GSTAT.
Inspired by: Experimental (Semi-) Variogram, variogramfit
Christian Brembilla (view profile)
Dear Wolfgang,
I was trying the example in your code. I just cut and past in the command windod but I get an error. Matlab does not recognize the variable variogram. Do I have to built up another function variogram? How can I get of rid of this problem?
Ravindar Reddy (view profile)
how to run the code??2016b
Eric Roots (view profile)
Great code, works well for the small(ish) number of data points I'm working with.
One question, is there a way to make sure that the interpolated values stay within a certain range? Specifically, can I ensure the interpolated points are all positive? For now I've just been taking the absolute values of the interpolated image, but that sometimes leads to odd results.
Ricardo Arévalo (view profile)
I really like the code. I have tried Ordinary Kriging in gstat and geoR. My grid has 8,392 rows and 12,345 columns with a cellsize of 250 m. It take really a long time. I believe your code is faster even and the result is quite good. Will be good to add the error image?
Dorothy (view profile)
yuanyuan ding (view profile)
I really appreciate your files,but I hope you can adopt your file into 3D data ,it is too difficult to adopt it for me.
Ahmed Zidan (view profile)
Hello,
I have a data in a regular grid (dx=dy) and x and y have different length.
I am facing a problem to create a covariance matrix, because of the different length and regular grid.
I tried to meshgrid, but the matrix will be very large that matlab can't accept it.
I will be grateful for any suggesting.
Thank you all.
Jonathan Fan (view profile)
@Wolfgrang: I really appreciate your Kriging.m codes.
Below are the inputs and outputs for the "Kriging.m".
function [zi,s2zi] = kriging(vstruct,x,y,z,xi,yi,chunksize)
I was wondering can I ask you a quick question?
My question is:
If we have more than two observations at one location or very close to each other with very different value, it would result in 'singularities' in the 2D/3D surface.
I was wondering how would you solve this issue, and incorporate it into the "Kriging.m"function.
meo (view profile)
Willy (view profile)
@Wolfgang: I really appreciate your interactivity and your pertinence in all your comments. I just coded the Deutsch method to solve the problem of the negative weights. Thanks to your doi link. Thus, the gaussian variogram is now well working :)
Here are the lines code to insert within the "for" loop, between the paragraphs "solve system" and "estimate zi".
However, I'm not an excellent developper. I'll let you find some better/faster way to code this.
% solve system
lambda = A*b;
% correct negative weights (Deutsch 1996)
if sum(lambda(:)<0) % if negative weigths are found
% compute the "average absolute magnitude of the negative weights"
lambdaLim=lambda(1:end-1,:);
PosInd=lambdaLim>=0;
lambdaLim(PosInd)=NaN;
lambdaLim=nanmean(abs(lambdaLim),1);
% compute the mean covariance between "negative" locations and the IX ones
bNoOnes=b(1:end-1,:);
bLim=bNoOnes;
bLim(PosInd)=NaN;
bLim=nanmean(bLim,1);
% correct the weights
lambdaC=lambda(1:end-1,:); % step 1 of Deutsch 1996
lambdaC(lambdaC<0)=0; % step 2 of Deutsch 1996
lambdaC(PosInd & bNoOnes-repmat(bLim,numobs,1)<0 & lambda(1:end-1,:)-repmat(lambdaLim,numobs,1)<0)=0; % step 3 of Deutsch 1996
% restandardization of the weights
lambda(1:end-1,:)=lambdaC./repmat(sum(lambdaC,1),numobs,1);
end
% estimate zi
zi(IX) = lambda'*z;
Wolfgang Schwanghart (view profile)
@Sébastien: Thanks for pointing at the error when using a nugget variogram.
I am currently abroad, so I probably won't find much time to update the function. So, in the meantime please adopt the changes proposed by Sébastien.
The problem with the gaussian variogram is most likely due to negative weights (see http://dx.doi.org/10.1016/0098-3004(96)00005-2 ). I will try to implement the procedure suggested by Deutsch (1996) in a future update.
Sébastien (view profile)
Very nice submission. I was amazed by its speed.
However, two things to say:
1) there is a mistake in the code when expanding b with ones. It should be the following:
% expand b with ones
b = vstruct.func([vstruct.range vstruct.sill],b);
if ~isempty(vstruct.nugget)
b = b+vstruct.nugget;
end
b=[b;ones(1,chunksize)];
Indeed, the ones should be added AFTER adding the nugget. Otherwise, the nugget is also applied to the ones: an offset is then unfortunately observed on the results...
2) I have the same problem than Kirsten: it's working well with spherical variograms. However, when I apply a gaussian one, huge values are found. I'm actually investigating the problem. If I find something, I'll let you know. But I would appreciate if you find the solution first ;)
Wolfgang Schwanghart (view profile)
@Imran Zafar: Fig. 4 shows the kriging variance which quantifies the uncertainty of the kriging estimates. You will note, that the kriging variance increases with increasing distance from observations.
Imran Zafar (view profile)
@Wolfgang Schwanghart: can you please illustrate the figure showing 4 plots? No description given. Please get back
Kirsten (view profile)
I sent this dataset to Wolfgang and Yesid, but in case anyone else wants to take a look here's my issue... I've used this dataset with an exponential and spherical model and get good results, but when I use the gaussian model it goes crazy at the domain edges. It will even predict min and max values smaller/larger than those in the original dataset, which should not be possible (and doesn't happen with the other variogram models). The problem seems to be worst at the edge where x=23.45 (I'm using leave-one-out validation and getting strange answers). Thanks for your help!
Kirsten
Here's my data:
x,y,z=[0.000000 0.000000 0.143917
5.860000 0.000000 0.159720
11.730000 0.000000 0.147167
17.590000 0.000000 0.181833
23.450000 0.000000 0.094583
0.000000 7.400000 0.170333
5.860000 7.400000 0.161250
11.730000 7.400000 0.160000
17.590000 7.400000 0.130750
23.450000 7.400000 0.119333
5.860000 8.750000 0.132833
11.730000 8.750000 0.109583
17.590000 8.750000 0.113000
23.450000 8.750000 0.090167
0.000000 14.800000 0.131231
5.860000 14.800000 0.137417
11.730000 14.800000 0.175333
17.590000 14.800000 0.191167
23.450000 14.800000 0.125333
5.860000 17.490000 0.037083
11.730000 17.490000 0.138750
17.590000 17.490000 0.093500
23.450000 17.490000 0.076917
0.000000 22.200000 0.159923
5.860000 22.200000 0.112000
11.730000 22.200000 0.397667
17.590000 22.200000 0.189750
5.860000 26.240000 0.146250
11.730000 26.240000 0.113083
17.590000 26.240000 0.116083
23.450000 26.240000 0.108750
0.000000 29.610000 0.126833
5.860000 29.610000 0.122667
11.730000 29.610000 0.200167
17.590000 29.610000 0.128333
23.450000 29.610000 0.118500
5.860000 34.990000 0.077750
11.730000 34.990000 0.084167
17.590000 34.990000 0.240167
23.450000 34.990000 0.100167
0.000000 37.010000 0.122583
5.860000 37.010000 0.139417
11.730000 37.010000 0.145083
17.590000 37.010000 0.142250
23.450000 37.010000 0.085833
5.860000 43.740000 0.046538
11.730000 43.740000 0.120583
17.590000 43.740000 0.105833
23.450000 43.740000 0.079500
0.000000 44.410000 0.133500
5.860000 44.410000 0.148250
11.730000 44.410000 0.136182
17.590000 44.410000 0.151083
23.450000 44.410000 0.116917
0.000000 51.810000 0.080750
5.860000 51.810000 0.127500
11.730000 51.810000 0.124000
17.590000 51.810000 0.108917
23.450000 51.810000 0.144250
0.000000 52.480000 0.079583
5.860000 52.480000 0.038750
11.730000 52.480000 0.142545
17.590000 52.480000 0.081917
23.450000 52.480000 0.066417
0.000000 59.210000 0.120917
5.860000 59.210000 0.106333
11.730000 59.210000 0.118750
17.590000 59.210000 0.108500
23.450000 59.210000 0.091667
0.000000 61.230000 0.191750
5.860000 61.230000 0.086667
11.730000 61.230000 0.102417
17.590000 61.230000 0.083455
23.450000 61.230000 0.058250
0.000000 66.610000 0.080083
5.860000 66.610000 0.161750
11.730000 66.610000 0.092833
17.590000 66.610000 0.047500
23.450000 66.610000 0.138083
0.000000 69.980000 0.104667
5.860000 69.980000 0.059000
11.730000 69.980000 0.048833
17.590000 69.980000 0.038167
23.450000 69.980000 0.052917
0.000000 74.020000 0.056500
5.860000 74.020000 0.104833
11.730000 74.020000 0.084250
17.590000 74.020000 0.080083
23.450000 74.020000 0.065583
0.000000 78.730000 0.051333
5.860000 78.730000 0.076333
11.730000 78.730000 0.069083
17.590000 78.730000 0.067250
23.450000 78.730000 0.047167
0.000000 81.420000 0.099417
5.860000 81.420000 0.082000
11.730000 81.420000 0.064917
17.590000 81.420000 0.054000
23.450000 81.420000 0.026333
0.000000 87.470000 0.036333
5.860000 87.470000 0.032083
11.730000 87.470000 0.040250
17.590000 87.470000 0.025000
23.450000 87.470000 0.014000
0.000000 88.820000 0.070750
5.860000 88.820000 0.054000
11.730000 88.820000 0.058083
17.590000 88.820000 0.049917
0.000000 96.220000 0.035667
5.860000 96.220000 0.028750
11.730000 96.220000 0.031125
17.590000 96.220000 0.030217
23.450000 96.220000 0.005333];
Yesid Goyes (view profile)
Send me , i need to run kriging with example. (g.es@hotmail.es)
Wolfgang Schwanghart (view profile)
@Kirsten: I have just ran kriging with the example I provided in the function help and modified it to take a gaussian variogram. So far, I cannot reconstruct your experience. Can you please provide sample data?
Kirsten (view profile)
Has anyone used this function with the gaussian variogram model? I get good results using a spherical or exponential mode, but when I use the gaussian model it is giving me crazy values (larger or smaller than any of the observations) at the edges. I don't think it should be able to give any values out of the range of the observations, so I believe something isn't be calculated properly. I'd appreciate any suggestions!
banana li (view profile)
very good! Thanks!
Roya O (view profile)
Yadong (view profile)
Mohammad Javad Tourian (view profile)