3.875

3.9 | 8 ratings Rate this file 71 Downloads (last 30 days) File Size: 3.4 KB File ID: #29025
image thumbnail

Ordinary Kriging

by

 

2D-interpolation using geostatistics

| Watch this File

File Information
Description

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.

Acknowledgements

Experimental (Semi ) Variogram and Variogramfit inspired this file.

MATLAB release MATLAB 7.10 (R2010a)
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (14)
24 Apr 2013 meo  
04 Apr 2013 Willy

@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;

03 Apr 2013 Wolfgang Schwanghart

@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.

03 Apr 2013 Sébastien

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 ;)

15 Feb 2013 Wolfgang Schwanghart

@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.

11 Feb 2013 Imran Zafar

@Wolfgang Schwanghart: can you please illustrate the figure showing 4 plots? No description given. Please get back

16 Oct 2012 Kirsten

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];

16 Oct 2012 Yesid Goyes

Send me , i need to run kriging with example. (g.es@hotmail.es)

16 Oct 2012 Wolfgang Schwanghart

@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?

15 Oct 2012 Kirsten

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!

12 Oct 2012 banana li

very good! Thanks!

08 Sep 2012 Roya O  
11 Sep 2011 Yadong  
21 Apr 2011 Mohammad Javad Tourian  

Contact us