Introduction

This is an illustration of generation of sampling sites using the sp package. Much of the text and R code has been taken from California Soil Resource Lab at UC Davis.

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.

Load objects of interest

Setup the environment, define working directory, load a shapefile of polygons representing objects in the study area, and visualize them.

if(!require(sp)){
    install.packages("sp")
    library(sp)
}
## Loading required package: sp
if(!require(rgdal)){
    install.packages("rgdal")
    library(rgdal)
}
## Loading required package: rgdal
## rgdal: version: 1.1-10, (SVN revision 622)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 1.11.4, released 2016/01/25
##  Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/gdal
##  Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
##  Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/proj
##  Linking to sp version: 1.2-3
# set working directory:
setwd("/Users/laura/documents/unal/ame")

# get data:
# if you need the shapefile used here, please go to https://goo.gl/xSKLj2
# 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('datos/oregon/orotl.shp', layer='orotl')
## OGR data source with driver: ESRI Shapefile 
## Source: "datos/oregon/orotl.shp", layer: "orotl"
## with 36 features
## It has 1 fields
plot(x, main="Objects of interest")

Random sampling schema

Let’s try a random sampling scheme using the spsample() function. Note that spsample will not sample each polygon, rather it works on the union of polys. This means that a random sampling scheme may lead to unsampled polygons.

# our sample size will have 100 sites
# try it:
plot(x) ; points(spsample(x, n=100, type='random'), col='red', pch=3, cex=0.5)

Systematic sampling schema

Now, let’s try a systematic sampling scheme. Take note of results in your notebook.

# our sample size will have 100 sites
# try it:
plot(x) ; points(spsample(x, n=100, type='regular'), col='red', pch=3, cex=0.5)

Stratified sampling schema

What about a stratified sampling scheme? Let’s see:

# our sample size will have 100 sites
# try it:
plot(x) ; points(spsample(x, n=100, type='stratified'), col='red', pch=3, cex=0.5)

Another sampling technique

In this approach, we need to iterate through all polygons in our original dataset, generating approximately 10 sample points within each polygon. Note that we use sapply() to iterate through the list of polygons, and do.call(‘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=10, type='hexagonal', offset=c(0,0)))

# what the hell is s?
typeof(s)
## [1] "list"
# a list of what?
class(s[[1]])
## [1] "SpatialPoints"
## attr(,"package")
## [1] "sp"
# I see. Let's summarize such a list:
summary(s)
##       Length Class         Mode
##  [1,]  9     SpatialPoints S4  
##  [2,]  8     SpatialPoints S4  
##  [3,] 10     SpatialPoints S4  
##  [4,]  9     SpatialPoints S4  
##  [5,] 10     SpatialPoints S4  
##  [6,] 10     SpatialPoints S4  
##  [7,]  9     SpatialPoints S4  
##  [8,] 12     SpatialPoints S4  
##  [9,]  9     SpatialPoints S4  
## [10,]  9     SpatialPoints S4  
## [11,] 12     SpatialPoints S4  
## [12,] 10     SpatialPoints S4  
## [13,] 11     SpatialPoints S4  
## [14,] 11     SpatialPoints S4  
## [15,] 11     SpatialPoints S4  
## [16,]  8     SpatialPoints S4  
## [17,] 10     SpatialPoints S4  
## [18,] 11     SpatialPoints S4  
## [19,]  9     SpatialPoints S4  
## [20,] 10     SpatialPoints S4  
## [21,]  9     SpatialPoints S4  
## [22,]  9     SpatialPoints S4  
## [23,]  6     SpatialPoints S4  
## [24,] 11     SpatialPoints S4  
## [25,]  9     SpatialPoints S4  
## [26,]  9     SpatialPoints S4  
## [27,] 10     SpatialPoints S4  
## [28,]  9     SpatialPoints S4  
## [29,] 11     SpatialPoints S4  
## [30,] 10     SpatialPoints S4  
## [31,] 10     SpatialPoints S4  
## [32,] 12     SpatialPoints S4  
## [33,]  9     SpatialPoints S4  
## [34,] 10     SpatialPoints S4  
## [35,] 10     SpatialPoints S4  
## [36,] 10     SpatialPoints S4
# we now have a list of SpatialPoints objects, one entry per polygon of original data
plot(x, main="Sampling points whithin a given polygon") ; points(s[[35]], col='red', pch=3, cex=.5)

# stack into a single SpatialPoints object
s.merged <- do.call('rbind', s)

As 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
s.final <- SpatialPointsDataFrame(s.merged, data=data.frame(poly_id=pt_id))
 
# check:
plot(x, main="Sampling points associated to each polygon") ; points(s.final, col=s.final$poly_id, pch=3, cex=0.5)

The last task is to 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
s.final@proj4string <- x@proj4string
 
# write out to new file
writeOGR(s.final, dsn='datos/oregon/', layer='sample_pts', driver='ESRI Shapefile')

Do not forget to try this code by yourself. You might want to use your own data.

REFERENCES

Pebesma, E.J., R.S. Bivand, 2005. Classes and methods for spatial data in R. R News 5 (2), http://cran.r-project.org/doc/Rnews/.

Roger S. Bivand, Edzer Pebesma, Virgilio Gomez-Rubio, 2013. Applied spatial data analysis with R, Second edition. Springer, NY. http://www.asdar-book.org/