Example command files

Following are the example command files, as distributed with gstat. The data (zinc concentrations of the top soil) are collected in a flood plain of the river Maas, not far from where the Maas entered the Netherlands (Borgharen, Itteren, about 3 to 5 km North of Maastricht). All coordinates are in metres, using the standard coordinates of Dutch topographical maps. Moving from the river, zinc concentrations tend to decrease.

For the universal kriging examples, let the function $D(x)$ be the function that is for every location $x$ the (normalised) square root distance to the river. This function is physically stored for the observation locations in column 4 of zinc_map.eas (as obtained in example 9), and for the prediction locations in the map sqrtdist.map . The prediction area is split in two separate sub-regions, $A$ and $B$ . Let the function $I_A (x)$ be 1 if $x \in A$ or else 0; and let the function $I_B (x)$ be 1 if $x \in B$ or else 0.

The functions $I_A (x)$ and $I_B (x)$ are physically stored for the observation locations in columns 5 and 6 of zinc_map.eas, and for the prediction locations in the maps part_a.map and part_b.map. (Note that the partitioning in $A$ and $B$ is arbitrary, it serves only illustrational purposes.)

The following command files, data and maps are distributed with gstat. These command files are purely for illustrational purposes and are by no means intended as some form of a `correct' analysis. Remind that everything from a # to the end of the line is comment.

<

Example 1. Variogram modelling

#
# One variable definition:
# to start the variogram modelling user interface.
#
data(zinc): 'zinc.eas', x=1, y=2, v=3;

Example 2. Variogram modelling of two variables

#
# Two variables with (initial estimates of) variograms,
# start the variogram modelling user interface
#
data(zinc): 'zinc.eas', x=1, y=2, v=3;
data(ln_zinc): 'zinc.eas', x=1, y=2, v=3, log;

variogram(zinc): 10000 Nug() + 140000 Sph(800);
variogram(ln_zinc): 1 Nug() + 1 Sph(800);

Example 3. Inverse distance interpolation

#
# Inverse distance interpolation on a mask map
#
data(zinc): 'zinc.eas', x=1, y=2, v=3;
mask: 'mask_map'; # the prediction locations
predictions(zinc): 'id_pr'; # result map

Example 4. Ordinary block kriging

#
# Local ordinary block kriging at non-gridded locations
#
data(zinc): 'zinc.eas', x=1, y=2, v=3,
min=20, max=40, radius=1000; # local neighbourhood
variogram(zinc): 2.42e+04 Nug(0) + 1.34e+05 Sph(800);
data(): 'locs.eas', x=1, y=2; # prediction locations
blocksize: dx=40, dy=40; # 40 $\times$ 40 block averages
set output = 'zinc_ok.out'; # ascii output file

Example 5. Simple kriging on a mask map

#
# Local simple point kriging on a mask map
#
data(ln_zinc): 'zinc.eas', x=1, y=2, v=3, log,
min=20, max=40, radius=1000, sk_mean=5.9;
variogram(ln_zinc): 0.0554 Nug(0) + 0.581 Sph(900);
mask: 'mask_map';
predictions(ln_zinc): 'lzn_skpr';
variances(ln_zinc): 'lzn_skvr';

Example 6. Unconditional simulation

#
# Unconditional Gaussian simulation on a mask
# (local neigbourhoods, simple kriging)
#
# defines empty variable:
data(ln_zn_dummy): dummy, sk_mean=5.9, max=20;
variogram(ln_zn_dummy): 0.0554 Nug(0) + 0.581 Sph(900);
mask: 'mask_map';
method: gs; # Gaussian simulation instead of kriging
predictions(ln_zn_dummy): 'lzn_uspr';

Example 7. Conditional simulation

#
# Gaussian simulation, conditional upon data
# (local neighbourhoods, simple and ordinary kriging)
#
data(ln_zinc): 'zinc.eas', x=1, y=2, v=3, log,
sk_mean=5.9, max=20;
variogram(ln_zinc): 0.0554 Nug(0) + 0.581 Sph(900);
mask: 'mask_map';
method: gs;
predictions(ln_zinc): 'lzn_cspr';
set n_uk = 20;
# use ordinary kriging when $\ge$ 20 data in neighbouhood
# set nsim=10;

Example 8. Ordinary block kriging on a mask map

#
# Change of support: local ordinary block kriging on a mask
#
data(ln_zinc): 'zinc.eas', x=1, y=2, v=3, log,
min=20, max=40, radius=1000;
variogram(ln_zinc): 0.0554 Nug(0) + 0.581 Sph(900);
mask: 'mask_map';
predictions(ln_zinc): 'lzn_okbp';
variances(ln_zinc): 'lzn_okbv';
blocksize: dx=40, dy=40; # define block dimensions

Example 9. Map values at point locations

#
# Obtain map values at data() locations
# (Point-map overlay)
#
data(a): dummy; # define n dummy data variable (n=n_masks/2)
data(b): dummy;
data(): 'zinc.eas', x=1, y=2, v=3; # prediction locations
method: map; # mapvalues as `predictions'
masks: 'sqrtdist', 'part_a', 'part_b'; # the maps
set output = 'zinc_map.eas'; # ascii output file.

Example 10. Multiple kriging

#
# Multiple kriging: prediction on more than one variable
# (ordinary kriging of two variables)
# (note that zinc_map.eas wass obtained through ex09.gst)
#
data(ln_zinc): 'zinc_map.eas', x=1, y=2, v=3, log,
min=20, max=40, radius=1000;
data(sq_dist): 'zinc_map.eas', x=1, y=2, v=4,
min=20, max=40, radius=1000;
variogram(ln_zinc): 0.0554 Nug(0) + 0.581 Sph(900);
variogram(sq_dist): 0.0631 Sph(900);
mask: 'mask_map';
predictions(ln_zinc): 'lzn_okpr';
variances(ln_zinc): 'lzn_okvr';
predictions(sq_dist): 'sqd_okpr';
variances(sq_dist): 'sqd_okvr';

Example 11. Multivariable kriging (cokriging)

#
# Multivariable kriging: ordinary local cokriging of two variables
#
data(ln_zinc): 'zinc_map.eas', x=1, y=2, v=3, log,
min=20, max=40, radius=1000;
data(sq_dist): 'zinc_map.eas', x=1, y=2, v=4,
min=20, max=40, radius=1000;
variogram(ln_zinc): 0.0554 Nug(0) + 0.581 Sph(900);
variogram(sq_dist): 0.0001 Nug(0) + 0.0631 Sph(900);
variogram(ln_zinc, sq_dist): 0 Nug(0)-0.156 Sph(900);
# NOTE: the 0 Nug(0)'s are added to make gstat recognize
# the Linear Model of Coregionalization
mask: 'mask_map';
predictions(ln_zinc): 'lzn_ckpr';
variances(ln_zinc): 'lzn_ckvr';
predictions(sq_dist): 'sqd_ckpr';
variances(sq_dist): 'sqd_ckvr';
# the next map holds the prediction error covariances:
covariances(sq_dist,ln_zinc): 'znsqdcov';

Example 12. Stratified ordinary kriging

#
# Stratified ordinary kriging (within-categorie ordinary kriging)
#
data(zinc.at.0): 'zinc_at0.eas', x=1, y=2, v=3, log,
min=20, max=40, radius=1000; # where $\htmladdnormallink{part\_a}{gif/part\_a.gif} = 0$
data(zinc.at.1): 'zinc_at1.eas', x=1, y=2, v=3, log,
min=20, max=40, radius=1000; # where $\htmladdnormallink{part\_a}{gif/part\_a.gif} = 1$
variogram(zinc.at.0): 0.0654 Nug(0) + 0.548 Sph(900);
variogram(zinc.at.1): 0.716 Sph(900);
# the mask map is 0 for zinc.at.0 locations, 1 for zinc.at.1
mask: 'part_a';
# stratified mode: one map holds predictions for all vars:
predictions: 'lzn_stpr';
# another the prediction variances for all vars:
variances: 'lzn_stvr';

Example 13. Universal kriging (a)

using the model $Z(x)=\beta_1 + D(x) \beta_2 +e(x)$:
#
# Local universal kriging, using one continuous variable
#
data(ln_zinc): 'zinc_map.eas', x=1, y=2, v=3, log,
X=4, # sqrtdist values at data locations
min=20, max=40, radius=1000; # apply model locally
# the variogram should be that of the residual:
variogram(ln_zinc): 0.0674 Nug(0) + 0.149 Sph(700);
mask: 'sqrtdist'; # sqrtdist values at prediction locations
predictions(ln_zinc): 'lzn_ukpr';
variances(ln_zinc): 'lzn_ukvr';

Example 14. Universal kriging (b)

using the model $Z(x)= D(x) \beta_1 + I_A (x) \beta_2 + I_B (x) \beta_3 +e(x)$:
#
# Universal kriging, using one continuous and
# two binary variables.
#
data(ln_zinc): 'zinc_map.eas', x=1, y=2, v=3, log,
X=-1&4&5&6;
# -1: no default intercept (col. 5 and 6 form an intercept)
# use global kriging: local kriging would lead to a singularity
# the variogram of e is:
variogram(ln_zinc): 0.0698 Nug(0) + 0.147 Sph(709);
# mask maps holding the independent variable values
# at prediction locations, their order corresponding
# that of the X-columns:
mask: 'sqrtdist', 'part_a', 'part_b';
predictions(ln_zinc): 'lzn_vkpr';
variances(ln_zinc): 'lzn_vkvr';

Example 14a. Stratified universal kriging

using the model $Z_i (x)=\beta_{i,1} + D(x) \beta_{i,2} +e_i (x)$, $i$ being the category number:
#
# Stratified universal kriging (within-categorie universal kriging)
#
data(zinc.at.0): 'zinc_at0.eas', x=1, y=2, v=3, X=4, log,
min=20, max=40, radius=1000; # where $\htmladdnormallink{part\_a}{gif/part\_a.gif} = 0$
data(zinc.at.1): 'zinc_at1.eas', x=1, y=2, v=3, X=4, log,
min=20, max=40, radius=1000; # where $\htmladdnormallink{part\_a}{gif/part\_a.gif} = 1$

# residual variograms:
variogram(zinc.at.0): 0.096572 Nug(0) + 0.226367 Sph(1069.33);
variogram(zinc.at.1): 0.115766 Sph(237.257);

mask: 'part_a', # 0 for zinc.at.0 locations, 1 for zinc.at.1 locs.
'sqrtdist', # predictor values corresp. to col. 4 for zinc.at.0
'sqrtdist'; # predictor values corresp. to col. 4 for zinc.at.1
predictions: 'lzn_stup';
variances: 'lzn_stuv';

Example 15. Linear model prediction (OLS)

using the model $Z(x)=\beta_1 + D(x) \beta_2 +e(x)$:
#
# Local linear model, using one continuous variable
#
data(ln_zinc): 'zinc_map.eas', x=1, y=2, v=3, X=4, log,
min=20, max=40, radius=1000; # apply linear model locally
# no variogram definition: assume residual to be IID.
mask: 'sqrtdist';
predictions(ln_zinc): 'lzn_trpr';
# prediction variance for point locations:
variances(ln_zinc): 'lzn_trvr';

Example 16. Multivariable (cumulative) indicator simulation

#
# Multivariable indicator cosimulation
#
data(i200): 'zinc.eas', x=1, y=2, v=3, max=20, radius=1000,
I=200, sk_mean = 0.28;
data(i400): 'zinc.eas', x=1, y=2, v=3, max=20, radius=1000,
I=400, sk_mean = 0.56;
data(i800): 'zinc.eas', x=1, y=2, v=3, max=20, radius=1000,
I=800, sk_mean = 0.85;

# define an LMC:
variogram(i200): 0.0490637 Nug(0) + 0.182814 Exp(300);
variogram(i400): 0.0608225 Nug(0) + 0.21216 Exp(300);
variogram(i200, i400): 0 Nug() + 0.14806 Exp(300);
variogram(i800): 0.0550284 Nug(0) + 0.0842966 Exp(300);
variogram(i200, i800): 0 Nug() + 0.0525584 Exp(300);
variogram(i400, i800): 0 Nug() + 0.102852 Exp(300);

method: is;
mask: 'mask_map';
# apply order corrections for cumulative indicators:
set order = 4;

predictions(i200): 'i200pr';
predictions(i400): 'i400pr';
predictions(i800): 'i800pr';
# uncomment next line to get 5 simulations:
# set nsim = 5;

Example 17. Trend surface prediction

#
# global coordinate polynomial trend surfaces
# trend orders 0-3.
#
data(zinc.0): 'zinc.eas', x=1, y=2, v=3, d=0, log;
data(zinc.1): 'zinc.eas', x=1, y=2, v=3, d=1, log;
data(zinc.2): 'zinc.eas', x=1, y=2, v=3, d=2, log;
data(zinc.3): 'zinc.eas', x=1, y=2, v=3, d=3, log;
mask: 'mask_map';
# predict block averages for very small blocks:
blocksize: dx=1, dy=1;
# variances apply to mean values,
# not for single observations
predictions(zinc.0): 'lzn_tr0';
variances (zinc.0): 'lzn_vr0';
predictions(zinc.1): 'lzn_tr1';
variances (zinc.1): 'lzn_vr1';
predictions(zinc.2): 'lzn_tr2';
variances (zinc.2): 'lzn_vr2';
predictions(zinc.3): 'lzn_tr3';
variances (zinc.3): 'lzn_vr3';