13.1. GEOSTATISTICS WITH GRASS AND GSTAT

The gstat4 package is Free Software for geostatistical modeling, prediction and simulation in one, two or three dimensions (Pebesma and Wesseling, 1998 and Pebesma, 2001). It requires the gnuplot5 graphical plotting software for the display of empirical variograms and variogram models.

With gstat you can perform geostatistical modeling in terms of generating empirical (sample) variograms and cross variograms (or covariograms). The software calculates sample (co-)variograms from ordinary, weighted or generalized least squares residuals. Models can be fitted to these variograms to predict data distributions. Using weighted least squares, nested models are fitted to sample (co-)variograms. Restricted maximum likelihood estimation of partial sills is also implemented. Variograms are plotted using the plotting program gnuplot, when working in interactive variogram modeling user interface.

The gstat software provides prediction and estimation using a model that is the sum of a trend modeled as a linear function of polynomials of the coordinates or of user-defined base functions, and an independent or dependent, geostatistically modeled residual. This allows for simple, ordinary and univer-

north-east direction. Burrough and McDonnell, 1998, use a subset of the same data set in their book (reduced area). The data set provided with the R/GRASS interface was re-projected from the Dutch standard coordinate system (TDN coordinates in stereographic projection) to UTM coordinate system zone 32, on WGS84 ellipsoid. A GRASS LOCATION was defined with following parameters: projection UTM, ellipsoid WGS84, zone 32, north 5652930, south 5650610, west 269870, east 272460, nsres 10, ewres 10, rows 232, cols 259. A pre-defined LOCATION including the data stored column-wise in sites lists can be downloaded from the GRASS Web site.3 The data sets are stored as sites lists plus one raster map; they can be used to experiment with interpolation or other methods.

This data set contains the following columns (topsoil data were collected as bulk samples during fieldwork in 1990 within a radius of 5 m according to Burrough and McDonnell, 1998:102, 309):

East, north (UTM zone 32 coordinates in meters); x, y (local coordinates in meters); elev (elevation above local reference level in meters); d.river (distance from main river Maas channel in meters); Cd (cadmium in ppm); Cu (copper in ppm); Pb (lead in ppm); Zn (zinc in ppm); LOI (percentage organic matter loss on ignition); flfd (flood frequency class, 1: annual, 2: 2-5 years, 3: every 5 years); soil (3 unnamed soil types).

sal kriging, simple, ordinary and universal cokriging, standardized cokriging, kriging with external drift, block kriging and kriging the trend , as well as uncorrelated, ordinary or weighted least squares regression prediction. Simulation in gstat comprises uni- or multivariable conditional or unconditional multi-Gaussian sequential simulation of point values or block averages, or (multi-) indicator sequential simulation (features cited after Pebesma, 2001).

The gstat/GRASS interface allows the user to read point data from site lists and raster maps. This requires to have the GRASS support compiled into gstat. You can check your version with flag -v:

gstat -V

The line with libraries must list grass besides other supported formats (e.g. grass gdal netcdf ).

Output of gstat (prediction or simulation results) is written to raster maps and also to site lists. You need to run gstat from inside GRASS as the program requires the GRASS environment to internally set up the LOCATION definitions. When a subregion is set in GRASS, gstat will only interpolate or simulate the raster cells according to the current region. The variables of interest need to be floating point numbers (DOUBLE) in sites list or stored in a raster map. The instructions for gstat are stored in an ASCII file. When using GRASS sites lists as input maps, following column order conventions have to be followed:

Easting I Northing I#site no %FP data [%FP data]

These are the same conventions as for standard GRASS sites lists. The program gstat reads GRASS site data from the current MAPSET with the data() function. Variable positions are defined as:

x=l: coordinate column 1 contains the x-coordinate y=2: coordinate column 2 contains the y-coordinate z=3: coordinate column 3 contains the z-coordinate (optional) v=l: data column 1 contains the first data (measurement) variable, when 0, a grid map is read

To illustrate how it works, we run a sample session based on the Maas river bank data set. First start GRASS with the Maas UTM LOCATION, then copy the Zn (zinc) concentrations sites map to the current MAPSET:

grass53 /usr/local/share/grassdata/maas/userl/ g.copy sites=Zn,zinc s.info zinc

The following example is based on the manual of gstat6. Store the following commands to the file gstat.maas1.zn in your home-directory:

# Two variables with (initial estimates of) variograms,

# start the variogram modeling user interface data (zinc): Zn, x=l, у=2, v=l;

data(ln zinc): Zn, x=l, y=2, v=l, log; variogram(zinc): 10000 NugO + 140000 Sph(800); variogram(ln zinc): 1 NugO + 1 Sph(800);

As the zinc concentrations are stored as first DOUBLE attribute in the sites list (in ppm, reported by s.info zinc), we select this data column through v=1. Run the analysis by:

gstat gstat.maasl.zn

The program starts to analyze the data and subsequently displays univariate statistics:

gstat: Linux version 2.4.3 (04 January 2004) Copyright (C) 1992, 2004 Edzer J. Pebesma using Marsaglias random number generator data(zinc): gisrc: [/home/neteler/.grassrcS] GRASS site list zinc: 0 cat, 2 dim, 0 str, 1 dbl. gstat/grass: 98 sites read successfully.

(GRASS site list)

[x:] x l : [ 269800, 272500]

[y:] y 2 : [5.6506e+06, 5.653e+06]

sample std.: 398.808

0 cat, 2 dim, 0 str, 1 dbl. data(ln zinc): gstat/grass: 98 sites read successfully, zinc (GRASS site list)

zinc

attribute: col[1] n: 98

sample mean: 481.031 GRASS site list zinc

attribute: log (col[1]) [x:] x l n: 98 [y:] y 2

sample mean: 5.87065 sample std. [starting interactive mode] press return to continue...

[ 269800, 272500] [5.6506e+06, 5.653e+06] 0.778309

After pressing <ENTER> we reach the main menu, which allows us to interactively analyze the loaded data set:

gstat 2.4.3 (04 January 2004), gstat.maas1.zn

enter/modify data choose variable calculate what cutoff, width direction variogram model fit method >show plot <Tab>

zxnc

semivariogram 1204.15, 80.2765

total

10000 Nug(O) + 140000 Sph(800) no fit

Command: