Code covered by the BSD License  

Highlights from
Kriging and Inverse Distance Interpolation using GSTAT

2.5

2.5 | 2 ratings Rate this file 36 Downloads (last 30 days) File Size: 5.7 KB File ID: #31055
image thumbnail

Kriging and Inverse Distance Interpolation using GSTAT

by

 

15 Apr 2011 (Updated )

InterpolationGstat performs inverse distance or kriging interpolation using gstat.

| Watch this File

File Information
Description

Kriging and inverse distance are popular interpolation methods, especially in earth sciences. There are some routines already available on matlab but are severely limited by matlabs memory constraints. By using gstat to handle interpolation and variogram modelling, very large datasets are easily handled.

USAGE: Requires gstat.exe (freeware) to run (www.gstat.org).

% INVERSE DISTANCE:
% INPUTS:
% x :: sample locations in x
% y :: sample locations in y
% z :: sample data
% res :: resolution of interpolation grid
% radius :: maximum search radius for interpolation
% method = 'invdist'
% uselog :: 1 if z values should be log transformed, 0 if not
% name :: name of output file without extension. e.g. 'predictions'.
% (file is given extension '.grd')
% useloc (OPTIONAL) :: if 1, will use file 'locations.dat' if it exist for
% interpolation locations. If 0, will prompt user to choose whether to use
% locations.dat (if it exists) or rewrite. Useful when calling the
% interpolation routine in a loop.
% S is not required for inverse distance
%
% KRIGING:
% INPUTS:
% x,y,z,res,radius :: as for inverse distance
% method = 'kriging'
% uselog :: as for inverse distance
% name :: as for inverse distance
% S :: Structure containing variogram information as output
% by the variogramfit function.
% useloc (OPTIONAL):: as for inverse distance
%
% OUTPUTS:
% xi :: the prediction locations in x
% yi :: the prediction locations in y
% zi :: the predicted result
% pcolor can be used to visualise. e.g. pcolor(xi,yi,zi); shading interp;

For variogram modelling, use SampleVarioGstat.m to create the experimental variogram. The excellent variogramfit by Wolfgram Schwanghart should be used to fit the experimental variogram. File ID: #25948.

Look at the readme for tips on kriging and inverse distance interpolation, and help InterpolationGstat and help SampleVarioGstat for correct usage.

NOTE:
This routine is intended to make it easy to perform kriging or inverse distance interpolation. It is not intended as a full interface to all the capabilities of gstat. As such only 2D ordinary kriging, inverse distance and experimental variogram modelling are supported, with only key parameters being offered as options.
 

Acknowledgements

Variogramfit inspired this file.

MATLAB release MATLAB 7.9 (R2009b)
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (10)
22 Mar 2014 Shruti

Sir, I founf the functions very helpful. But when I'm running the function InterpolationGstat.m I'm getting an error gstat.exe has stopped working whereas I could able to run the function sampleVarioGstat.m . Can you please let me know what todo. Its really urgent.

15 Jun 2012 Cynthia Gerlein

James,
I would like to get the semi-variogram of a very large dataset (230 000 datapoints). Gstat seems to be the perfect solution but i get this error:

Error using fgets
Invalid file identifier. Use fopen to generate
a valid file identifier.

Error in fgetl (line 34)
[tline,lt] = fgets(fid);

Error in SampleVarioGstat (line 68)
fgetl(fid);

Error in AutocorrGPR (line 361)
[av_dist,gamma,h,np]=SampleVarioGstat(x,y,z)

I used Wolfgang's example for x, y and z : x = rand(1000,1)*4-2; y = rand(1000,1)*4-2; z = 3*sin(x*15)+ randn(size(x)); and I am using Matlab 2012a on MacOS 10.7.4.

Any suggestion on what the problem might be ?

15 Jun 2012 Cynthia Gerlein

I just noticed that gstat doesn't seem to be working and I'm guessing this might be the source of the error above. On top of the error message that I posted previously, I also get:

Trying to run GSTAT on gstat.cmd
GSTAT FAILED ............

I tried to download the version for Mac (http://www.gstat.org/download.html) and to put in my file, but it's still not working...

15 Jun 2012 Cynthia Gerlein

PS: I am using Matlab 2012a on MacOS 10.7.4.

27 Apr 2011 Wolfgang Schwanghart

Hi James,
thanks for the update. There is still a minor error, that could be easily avoided. If you call InterpolationGstat for the first time and set useloc to true, than an error message tells me that File does not exist. Calling it with useloc = false is fine then. Calling it with three output arguments and useloc = true again, returns an error again.

Error in ==> InterpolationGstat at 61
if nargin<8

??? Output argument "xi" (and maybe others) not assigned during call to
"U:\WDmatlab\toolboxes\InterpGstat\InterpolationGstat.m>InterpolationGstat".

Yet, otherwise, the function works well now. The only thing I would like to see in the future is that InterpolationGstat returns the kriging variance as well. I have included it in my kriging function on the FEX because I think it is one of the most important features of kriging.

Best regards, Wolfgang

21 Apr 2011 James Ramm

Hi Wolfgang,
Thanks for letting me know the errors.
I just uploaded a corrected version (should appear tomorrow, april 22nd)
Tested with 'peaks' data sets for both inverse distance and kriging.
Please let me know about any other errors.
Gstat still needs to be in your working directory. I will make this more flexible in the near future.

19 Apr 2011 Wolfgang Schwanghart

Hi James, thanks for submitting this function and for the credits! However, I think the functions still need some polishing, since they throw errors.

E.g.
****************
[xi,yi,zi] = InterpolationGstat(x,y,z,0.01,0.5,'kriging',0,'predictions',S,0);
??? Undefined function or variable "use".

Error in ==> InterpolationGstat at 96
if strcmp(use,'n')

*******************
[xi,yi,zi] = InterpolationGstat(x,y,z,0.01,0.5,'kriging',0,'predictions',S,1);
Writing Data File
Invoking Gstat...
gstat: Win32/MinGW version 2.4.1 (12 March 2003)
Copyright (C) 1992, 2003 Edzer J. Pebesma
using Marsaglia's random number generator
data(Krig): SampleData.dat (Ascii table file)
attribute: col[3] [x:] x_1 : [ -1.998, 1.998]
n: 1000 [y:] y_2 : [ -1.999, 1.998]
sample mean: 0.073122 sample std.: 2.40105
[using ordinary kriging]
gstat: read failed on file `locations'
ans =

7

??? Error using ==> fgetl at 44
Invalid file identifier. Use fopen to generate a valid file identifier.

Error in ==> ImportSurferGrd at 15
code=fgetl(grdfile); % Reads surfer code 'DSAA'

Error in ==> InterpolationGstat at 178
zi =ImportSurferGrd([name,'.grd']);

****************

In addition, the interpolationgstat seems to work only, if your current folder is the folder where gstat.exe is located. Can this restriction be avoided?

I think, it is a great idea to call gstat from within matlab. At their current state, however, your functions still need some editing to work. Looking forward to an update, Wolfgang

18 Apr 2011 James Ramm

Gstat is freely available from gstat.org. Its a small executable - just save it to the same directory as InterpolationGstat and you are set.

18 Apr 2011 Sean de

I don't have gstat; by looking at the code this appears to be a file worthy of much more than 1 star.

15 Apr 2011 Ian Howat  
Updates
21 Apr 2011

Fixed bugs in 'InterpolationGstat.m'
- errors involving use of 'uselog' and pre-existing locations files resolved
- if gstat fails, it no longer attempts to open and read the resulting predictions file (as it doesnt exist)

Contact us