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.

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.

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?

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

@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

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

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

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

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

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!