Spatial Coverage Sampling with spcosa

Spatial coverage sampling without prior points

In the first example, spatial coverage sampling will be applied to map clay content in a study area by means of ordinary point kriging.

First, a vector or raster representation of the study area has to be loaded. As an example, the ASCII-grid demoGrid.asc residing in the maps directory of the spcosa-package will be read by means of the readGDAL-function of the rgdal-package (Bivand et al., 2015)

library(rgdal)
filename <- system.file("maps", "demoGrid.asc", package = "spcosa")
ascii_grid <- readGDAL(fname = filename, silent = TRUE)

To obtain a uniform distribution of sampling points over the study area, the sampling points will be selected at the centroids of compact spatial strata. Compact strata can be constructed by invoking the stratify method:

stratification <- stratify(ascii_grid, nStrata = 75, nTry = 10)

In this example, the study area has been partitioned into 75 compact strata. The resulting stratification can be plotted by means of:

plot(stratification)

Each plot can be modified by adding ggplot2-functions to the plot-method:

plot(stratification) +
    scale_x_continuous(name = "Easting (km)") +
    scale_y_continuous(name = "Northing (km)")

Scale for 'x' is already present. Adding another scale for 'x', which will
replace the existing scale.

Scale for 'y' is already present. Adding another scale for 'y', which will
replace the existing scale.

The spsample-method, an overloaded method from the sp-package, can be used to select the centroid of each stratum:

sampling_pattern <- spsample(stratification)

The plot-method can be used to visualise the resulting sampling pattern:

plot(sampling_pattern)

The sampling pattern can also be plotted on top of the stratification:

plot(stratification, sampling_pattern)

Sampling locations can be extracted by means of a simple type cast to class data.frame:

sampling_points <- as(sampling_pattern, "data.frame")
head(sampling_points)

          x         y
1 79.971831  5.253521
2 13.703704 37.802469
3 11.338710 10.435484
4 26.808824 26.323529
5  4.090909 24.038961
6 25.298246 17.298246

Next, field work needs to be done to acquire data at these locations.

Loading required package: gstat

Registered S3 method overwritten by 'xts':
  method     from
  as.zoo.xts zoo 

These data will be stored in a data.frame:

head(my_data)

  clay SOM
1 11.2 5.8
2 10.6 4.8
3 10.3 5.3
4 12.8 6.2
5 10.9 6.0
6  9.1 4.9

Model-based inference will be applied to predict clay contents. Model-based inference is not implemented in the spcosa-package, because other add-on packages already provide excellent support for it. See for instance, gstat (Pebesma, 2004), spatial (Venables & Ripley, 2002), and geoR (Diggle & Ribeiro, 2007).

As an example, we will use the gstat package for predicting clay contents by means of ordinary point kriging. First, the observations and sampling pattern will be combined in a single object of class SpatialPointsDataFrame:

#sp_data <- data.frame(sampling_points, my_data)
#coordinates(sp_data) <- ~ x * y
sp_data <- SpatialPointsDataFrame(coords = sampling_points, data = my_data)

Next, the sample variogram for clay has to be estimated, and needs to be fitted by a permissable model:

sample_variogram <- variogram(clay ~ 1, sp_data)
variogram_model <- fit.variogram(sample_variogram, 
    model = vgm(psill = 1, model = "Exp", range = 5, nugget = 0))
plot(sample_variogram, model = variogram_model)

Finally, ordinary point kriging predictions can be obtained by means of:

g <- gstat(id = "clay", formula = clay ~ 1, data = sp_data, model = variogram_model)
y_hat <- predict(g, newdata = ascii_grid)

[using ordinary kriging]

The resulting predictions and variances of the prediction error are given in the figures below:

Spatial coverage sampling with prior points

Sometimes, samples from previous sampling campaigns are available. In these situations, spatial infill sampling may be performed. This type of spatial coverage sampling aims to distribute new sampling points evenly over the study area, while taking the locations of existing sampling points into account. Suppose a data.frame is available containing the coordinates of 50 existing sampling points:

head(prior_points)

     x    y
1 98.2 38.1
2 47.3 38.2
3 23.5 40.0
4 19.9 25.4
5  3.2 32.4
6 28.4  1.3

Twenty-five new points can be assigned to sparsely sampled regions by means of:

coordinates(prior_points) <- ~ x * y
stratification <- stratify(ascii_grid, priorPoints = prior_points, nStrata = 75, nTry = 100)
sampling_pattern <- spsample(stratification)

Note that the total number of strata, and therefore the total number of points, equals 50+25=75. In addition, also note that the nTry argument has been set to 100. The algorithm will now use 100 random starting configurations and keeps the best solution to reduce the risk of ending up in an unfavourable local optimum.

plot(stratification, sampling_pattern)

Note that prior points and new points are represented by different symbols.

A map of clay contents can be obtained by means of ordinary point kriging. Suppose again, that data have been collected in the field and stored in data.frame sp_data:

Again, we will use gstat for model-based inference:

sample_variogram <- variogram(clay~1, sp_data)
variogram_model <- fit.variogram(sample_variogram,
    model = vgm(psill = 1, model = "Exp", range = 1, nugget = 0))
plot(sample_variogram, model = variogram_model)

g <- gstat(id = "clay", formula = clay ~ 1, data = sp_data, model = variogram_model)
y_hat <- predict(g, newdata = ascii_grid)

[using ordinary kriging]

The resulting predictions and variances of the prediction error are given in the figures below:

Random sampling from compact geographical strata

In this section, the global mean clay and organic matter contents of the study area will be estimated by means of stratified simple random sampling. Again, we will use ascii_grid as a representation of the study area. The study area will be partitioned into 25 compact strata:

stratification <- stratify(ascii_grid, nStrata = 25, nTry = 10)

plot(stratification)

The spsample-method can be used to randomly sample two sampling units per stratum:

sampling_pattern <- spsample(stratification, n = 2)

plot(stratification, sampling_pattern)

Sampling points can be extracted by means of a type cast to class data.frame:

sampling_points <- as(sampling_pattern, "data.frame")
head(sampling_points)

         x        y
1 33.83659 12.15027
2 28.81836 18.51965
3 84.69064 36.65839
4 79.47751 50.19131
5 27.47371 18.48114
6 25.62352 17.01825

Next, some field work has to be done. Again, the results will be stored in data.frame my_data:

The spatial mean clay and soil organic matter contents can be estimated by (de Gruijter et al., 2006):

estimate("spatial mean", stratification, sampling_pattern, my_data)

   clay     SOM 
9.55672 4.81517 

In estimating the spatial mean, differences in surface area of the strata are taken into account. Note, that the spatial mean is estimated for all columns in my_data. The standard error can be estimated in a similar way:

estimate("standard error", stratification, sampling_pattern, my_data)

     clay       SOM 
0.2509970 0.1152259 

The spatial cumulative distribution function (SCDF) (see deGruijter et al., 2006) can be estimated by means of

scdf <- estimate("scdf", stratification, sampling_pattern, my_data)

The SCDFs are returned as a list of matrices, i.e. one matrix for each property:

lapply(X = scdf, FUN = head, n = 4)

$clay
     value cumFreq
[1,]   5.0  0.0000
[2,]   5.2  0.0205
[3,]   6.2  0.0412
[4,]   6.6  0.0612

$SOM
     value cumFreq
[1,]   2.3  0.0000
[2,]   2.9  0.0205
[3,]   3.3  0.0612
[4,]   3.5  0.1025

lapply(X = scdf, FUN = tail, n = 4)

$clay
      value cumFreq
[33,]  13.9  0.9204
[34,]  14.7  0.9415
[35,]  14.9  0.9590
[36,]  15.0  0.9796

$SOM
      value cumFreq
[26,]   7.1  0.9204
[27,]   7.2  0.9379
[28,]   7.4  0.9590
[29,]   7.6  0.9794

The SCDFs for clay and SOM are visualised below.

Stratified simple random sampling for composites

In this example, the aim is to estimate the global mean clay and organic matter contents of a field. To reduce laboratory costs, the soil aliquots collected at the sampling locations will be bulked into composite samples. The study area is a field near the village of Farmsum, in the North-East of the Netherlands. An ESRI shape file of this field is available in the maps directory of the spcosa-package. It can be loaded by means of readOGR, a function in the rgdal package (Bivand et al., 2015):

directory <- system.file("maps", package = "spcosa")
shp_farmsum <- readOGR(dsn = directory, layer = "farmsum", verbose = FALSE)

First, the field will be stratified into, say, 20 compact strata of equal size. Strata of equal size are desirable to simplify fieldwork, i.e. equal supports of soil can be collected at the sampling locations (Brus et al., 1999, de Gruijter et al., 2006).

stratification <- stratify(shp_farmsum, nStrata = 20, equalArea = TRUE, nTry = 10)
plot(stratification)

Next, two sampling units will be selected at random in each stratum. At least two sampling units per stratum are required to estimate the sampling variance of the estimated mean.

sampling_pattern <- spsample(stratification, n = 2, type = "composite")
plot(stratification, sampling_pattern)

Sampling points can be extracted means of:

sampling_points <- as(sampling_pattern, "data.frame")
head(sampling_points)

  composite       x1       x2
1         1 259338.4 587256.1
2         2 259317.3 587245.8
3         1 259353.0 587297.2
4         2 259334.0 587289.3
5         1 259340.4 587376.9
6         2 259340.6 587362.0

Note that an extra column has been added specifying the sampling units to be bulked into each composite. A composite sample is formed by bulking one aliquot (sampling unit) per stratum. Field work now results in a composite sample of size two. These data will be stored in a data.frame called my_data:

The spatial mean and its standard error can be estimated by means of:

estimate("spatial mean", stratification, sampling_pattern, my_data)

 clay   SOM 
10.70  5.35 

estimate("standard error", stratification, sampling_pattern, my_data)

  clay    SOM 
0.0100 0.0025 

If we do not want to bulk soil aliquots, the same stratification can be used to select a sample of 20 \(\times\) 2 sampling locations:

sampling_pattern <- spsample(stratification, n = 2)
sampling_points <- as(sampling_pattern, "data.frame")
head(sampling_points)

        x1       x2
1 259333.2 587266.0
2 259341.3 587258.7
3 259347.4 587303.3
4 259360.9 587285.2
5 259330.6 587387.2
6 259341.5 587359.6

plot(stratification, sampling_pattern)

Map projections

In the examples above, the maps didn’t have projection attributes. At field scale, projection attributes are not really necessary. However, at continental and global scales, for example, projection attributes can’t be ignored. The spcosa-package is capable of handling projection attributes of class CRS. More information on projections is available at the `PROJ.4’ project website and in the rgdal-package documentation (Bivand et al., 2015).

To illustrate the effect of stratification on smaller spatial scales, consider two (relatively coarse) grids covering the surface of the earth:

grd <- expand.grid(
    longitude = seq(from = -176, to = 180, by = 4),
    latitude  = seq(from =  -86, to =  86, by = 4)
)
gridded(grd) <-  ~ longitude * latitude

grd_crs <- grd
proj4string(grd_crs) <- CRS("+proj=longlat +ellps=WGS84")

Note that grd is identical to grd_crs, except that grd_crs has projection attributes. If no projection attributes are available, the algorithms in the spcosa-package use squared Euclidean distances in the k-means algorithms. However, if projection attributes have been set, the coordinates will be transformed to lat/long format (latitude/longitude) and squared great circle distances will be used instead.

Both grids will be partitioned into 50 compact geographical strata:

strata     <- stratify(grd,     nStrata = 50)
strata_crs <- stratify(grd_crs, nStrata = 50)

plot(strata)

plot(strata_crs)

Note that grd seems to have more compact strata near the geographic poles than grd_crs. However, the contrary is true. This becomes evident when both stratifications are projected on a sphere (use your mouse to interact with these spheres):

You must enable Javascript to view this page properly.

The left figure is based on grd, and the right figure on grd_crs. The strata of grd_crs are clearly more compact than those of grd. In addition, grd suffers from pronounced edge effects near the poles and at 180 degrees longitude. The strata are discontinuous at this meridian, i.e., two points on opposite sides of the meridian are treated as very distant when squared Euclidean distances are used. The strata of grd_crs, on the other hand, have been optimised on a sphere by using squared great circle distances, and don’t suffer from edge effects, i.e., the great circle distance between two nearby points on opposite sides of this meridian is very small.

Simple random sampling

Although this package is about sampling from compact strata, it can also be used for simple random sampling, by setting nStrata = 1:

shp_mijdrecht <- readOGR(
    dsn = system.file("maps", package = "spcosa"), 
    layer = "mijdrecht", verbose = FALSE)
stratification <- stratify(shp_mijdrecht, nStrata = 1, nGridCells = 5000)
sampling_pattern <- spsample(stratification, n = 30)

plot(stratification, sampling_pattern)

Sampling of non-convex areas

In case of spatial coverage sampling, sampling the centroid of each cluster may become problematic in case of non-convex areas. A centroid may be situated well outside the area of interest. If this happens, the sampling point will be relocated to the nearest grid cell that is part of the target universe. This pragmatic solution usually gives reasonable results. However, in some extreme situations the solution may be less desirable. As an example, consider the ‘doughnut’-shaped field below.

doughnut <- expand.grid(s1 = -25:25, s2 = -25:25)
d <- with(doughnut, sqrt(s1^2 + s2^2))
doughnut <- doughnut[(d < 25) & (d > 15), ]
coordinates(doughnut) <- ~ s1 * s2
gridded(doughnut) <- TRUE

stratification <- stratify(doughnut, nStrata = 2, nTry = 100)
sampling_pattern <- spsample(stratification)
plot(stratification, sampling_pattern)

Note that this problem does not arise in random sampling from compact geographical strata.