Generation of Sample Site Locations [sp package for R]

Posted: Friday, August 15th, 2008

Setting up sampling designs is a non-trivial aspect to any field experiment that includes a spatial component. The sp package for R provides a simple framework for generating point sampling schemes based on region-defining features (lines or polygons) and a sampling type (regular spacing, non-aligned, random, random-stratified, hexagonal grid, etc.). The rgdal package provides a set of functions for importing/exporting common vector data formats. This example demonstrates simple import/export, iterating over sp objects, and reconstituting new objects from lists of objects. A more complex sampling scheme is demonstrated here.

1. Setup the environment, load some sample polygons, and tryout the spsample() function.

# load required packages
# read data:
# note the funky syntax
# note that this should have a REAL projection defined
# an incorrect definition may result in an error from readOGR
x <- readOGR('polys/polys.shp', layer='polys')
# spsample will not sample each polygon, rather it works on the union of polys
# try it:
plot(x) ; points(spsample(x, n=100, type='random'), col='red', pch=3, cex=0.5)

Sampling with sample example 1

2. Iterate through all polygons in our original dataset, generating approximately 100 sample points within each polygon. Note that we use sapply() it iterate through the list of polygons, and'rbind', ...) to 'stack' the list elements back into a single SpatialPoints object.

# hexagonal grid from lower-left corner
s <- sapply(slot(x, 'polygons'), function(i) spsample(i, n=100, type='hexagonal', offset=c(0,0)))
# we now have a list of SpatialPoints objects, one entry per polygon of original data
plot(x) ; points(s[[4]], col='red', pch=3, cex=.5)
# stack into a single SpatialPoints object
s.merged <-'rbind', s)

Sampling with sample example 2

3. Now that the sample points for each polygon have been merged into a single SpatialPoints object, we need to attach a dataframe with the ID associating each sample point with its parent polygon. Attaching this data will "promote" the SpatialPoints object to a SpatialPointsDataFrame object.

# add an id, based on source polygon:
# extract the original IDs
ids <- sapply(slot(x, 'polygons'), function(i) slot(i, 'ID'))
# determine the number of ACTUAL sample points generated for each poly
npts <- sapply(s, function(i) nrow(i@coords))
# generate a reconstituted vector of point IDs
pt_id <- rep(ids, npts)
# promote to SpatialPointsDataFrame <- SpatialPointsDataFrame(s.merged, data=data.frame(poly_id=pt_id))
# check:
plot(x) ; points(,$poly_id, pch=3, cex=0.5)

Sampling with spsample example 3

4. Copy over the spatial reference system data from the polygons object, and save sample points to a new shapefile. Note that you can only write to a shapefile if the object in question is a SpatialPointsDataFrame object.

# copy source data spatial reference system to new object <- x@proj4string
# write out to new file
writeOGR(, dsn='polys/', layer='rnd_pts', driver='ESRI Shapefile')


Customizing Maps in R: spplot() and latticeExtra functions

Working with Spatial Data

Ordinary Kriging Example: GRASS-R Bindings