Point-process modelling with the sp and spatstat packages
Some simple examples of importing spatial data from text files, converting between R datatype, creation of a point process model and evaluating the model. Input data sources are: soil pit locations with mollic and argillic diagnostic horizons (mollic-pedons.txt and argillic-pedons.txt), and a simplified outline of Pinnacles National Monument (pinn.txt). The outline polygon is used to define a window in which all operations are conducted.
The 'sp' package for R contains the function spsample(), can be used to create a sampling plan for a given region of interest: i.e. the creation of n points within that region based on several algorithms. This example illustrates the creation of 50 sampling points within Pinnacles, according to the following criteria: regular (points are located on a regular grid), nonaligned (points are located on a non-aligned grid-like pattern), random (points are located at random), stratified (collectively exhaustive, see details here).
The 'spatstat' package for R contains several functions for creating point-process models: models describing the distribution of point events: i.e. the distribution of tree species within a forest. If covariate data is present (i.e. gridded precipitation, soil moisture, aspect, etc.) these covariates can be incorporated into the point-process model. Without covariate data, the model is based on an spatial distribution estimator function. Note that the development of such models is complicated by factors such as edge-effects, degree of stochasticity, spatial connectivity, and stationarity. These are rather contrived examples, so please remember to read up on any functions you plan to use for your own research. An excellent tutorial on Analyzing spatial point patterns in R was recently published.
Helpful links
Spatstat Quick Reference
Print Version with Links
![]() Four sampling designs |
![]() Spatial density of each pedon type |
![]() Spatial density of the four sampling designs |
![]() Example point-process model of mollic soils |
![]() Diagnostics of a simple point-process model |
Note: This code should be updated to use the slot() function instead of the '@' syntax for accessing slots!
Load required packages and input data (see attached files at bottom of this page)
# load required packages library(sp) library(spatstat) # read in pinnacles boundary polygon: x,y coordinates # use the GRASS vector, as it should be topologically correct # v.out.ascii format=standard in=pinn_bnd > pinn.txt ... edit out extra crud pinn <- read.table('pinn.txt') # read in mollic and argillic pedons # see ogrinfo hack m <- read.table('mollic-pedons.txt', col.names=c('x','y')) a <- read.table('argillic-pedons.txt', col.names=c('x','y')) # add a flag for the type of pedon m$type <- factor('M') a$type <- factor('A') #combine into a single dataframe pedons <- rbind.data.frame(a,m)
Using the functions from the 'sp' package create a polygon object from the pinn.txt coordinates
# create a spatial polygon class object pinn.poly <- SpatialPolygons(list(Polygons(list(Polygon( pinn )), "x"))) # inspect this new object with str() # total area of all polygons pinn.poly@polygons[[1]]@area # coordinates of first polygon: this is rough syntax! pinn.poly@polygons[[1]]@Polygons[[1]]@coords
Generate a sampling plan for 50 sites using regular grid, non-aligned grid, random, and random stratified approaches
# generate random points within the pinnacled boundary p.regular <- spsample(pinn.poly, n = 50, "regular") p.nalign <- spsample(pinn.poly, n = 50, "nonaligned") p.random <- spsample(pinn.poly, n = 50, "random") p.stratified <- spsample(pinn.poly, n = 50, "stratified") # setup plot environment par(mfrow=c(2,2)) # each of the sampling designs plot(pinn.poly) title("Regular") points(p.regular, pch=16, col='red', cex=0.3) plot(pinn.poly) title("Nonaligned") points(p.nalign, pch=16, col='red', cex=0.3) plot(pinn.poly) title("Random") points(p.random, pch=16, col='red', cex=0.3) plot(pinn.poly) title("Stratified") points(p.stratified, pch=16, col='red', cex=0.3)
Convert 'sp' class objects to 'spatstat' analogues note the use of 'slots'
# pinn boundary: # extract coordinates: and get a length - 1 value p.temp <- pinn.poly@polygons[[1]]@Polygons[[1]]@coords n <- length(p.temp[,1]) - 1 # create two vectors: x and y # these will contain the reversed vertices, minus the last point # in order to adhere to the spatstat specs: no repeating points, in counter-clockwise order x <- rev(p.temp[,1][1:n]) y <- rev(p.temp[,2][1:n]) # make a list of coordinates: note that we are removing the last vertex p.list <- list(x=x,y=y) # make a spatstat window object from the polygon vertices W <- owin(poly=p.list) # pedons with their 'marks' i.e. pedon type, and the pinn boundary as the 'window' pedons.ppp <- ppp(pedons$x, pedons$y, xrange=c(min(pedons$x), max(pedons$x)), yrange=c(min(pedons$y), max(pedons$y)) , window=W, marks=pedons$type)
Plot and summarize the new combined set of pedon data
# plot and summarize the pedons data: # note the method used to subset the two 'marks' par(mfrow=c(1,2)) plot(density.ppp(pedons.ppp[pedons.ppp$marks == 'M']), main="Mollic Point Density") points(pedons.ppp[pedons.ppp$marks == 'M'], cex=0.2, pch=16) plot(density.ppp(pedons.ppp[pedons.ppp$marks == 'A']), main="Argillic Point Density") points(pedons.ppp[pedons.ppp$marks == 'A'], cex=0.2, pch=16) summary(pedons.ppp) Marked planar point pattern: 151 points Average intensity 1.38e-06 points per unit area Marks: frequency proportion intensity A 62 0.411 5.67e-07 M 89 0.589 8.14e-07 Window: polygonal boundary single connected closed polygon with 309 vertices enclosing rectangle: [ 657228.3125 , 670093.8125 ] x [ 4030772.75 , 4047986.25 ] Window area = 109337135.585938
Convert the sampling design points (from above) to 'spatstat' objects and plot their density
# convert the random datasets: using the same window: ppp.regular <- ppp(p.regular@coords[,1], p.regular@coords[,2], window=W) ppp.nalign <- ppp(p.nalign@coords[,1], p.nalign@coords[,2], window=W) ppp.random <- ppp(p.random@coords[,1], p.random@coords[,2], window=W) ppp.stratified <- ppp(p.stratified@coords[,1], p.stratified@coords[,2], window=W) # visually check density of random points: par(mfrow=c(2,2)) plot(density.ppp(ppp.regular), main="Regular Sampling") points(ppp.regular, pch=16, cex=0.2) plot(density.ppp(ppp.nalign), main="Non-Aligned Sampling") points(ppp.nalign, pch=16, cex=0.2) plot(density.ppp(ppp.random), main="Random Sampling") points(ppp.random, pch=16, cex=0.2) plot(density.ppp(ppp.stratified), main="Stratified Sampling") points(ppp.stratified, pch=16, cex=0.2)
Simple, and probably flawed attempt to use a point-process model for the pedon data Third order polynomial model for the distribution of pedons with a mollic epipedon. See manula page for ppm() for detailed examples.
# model the spatial occurance of Mollic epipedons with a 3rd-order polynomial, using the Poisson Process Theory fit <- ppm( unmark(pedons.ppp[pedons.ppp$marks == 'M']), ~polynom(x,y,3), Poisson()) # view the fitted model par(mfcol=c(2,2)) plot(unmark(pedons.ppp[pedons.ppp$marks == 'M']), main="Mollic Pedons") plot(fit) # plot some diagnostics on the fitted model: Pearson residuals (see references) diagnose.ppm(fit, type="pearson") #another example using a buil-in dataset: the Lansing Forest # fit nonstationary marked Poisson process # with different log-cubic trend for each species data(lansing) fit <- ppm(lansing, ~ marks * polynom(x,y,3), Poisson()) plot(fit)
Point-process model diagnostic references from the spatstat manual
Baddeley, A., Turner, R., Moller, J. and Hazelton, M. (2005) Residual analysis for spatial point processes. Journal of the Royal Statistical Society, Series B 67, 617–666.
Stoyan, D. and Grabarnik, P. (1991) Second-order characteristics for stochastic structures connected with Gibbs point processes. Mathematische Nachrichten, 151:95–100.
Attachments:
Software
- General Purpose Programming with Scripting Languages
- LaTeX Tips and Tricks
- PostGIS: Spatially enabled Relational Database Sytem
- PROJ: forward and reverse geographic projections
- GDAL and OGR: geodata conversion and re-projection tools
- R: advanced statistical package
- Access Data Stored in a Postgresql Database
- Additive Time Series Decomposition in R: Soil Moisture and Temperature Data
- Aggregating SSURGO Data in R
- Cluster Analysis 1: finding groups in a randomly generated 2-dimensional dataset
- Color Functions
- Comparison of Slope and Intercept Terms for Multi-Level Model
- Comparison of Slope and Intercept Terms for Multi-Level Model II: Using Contrasts
- Creating a Custom Panel Function (R - Lattice Graphics)
- Customized Scatterplot Ideas
- Estimating Missing Data with aregImpute() {R}
- Exploration of Multivariate Data
- Interactive 3D plots with the rgl package
- Making Soil Property vs. Depth Plots
- Numerical Integration/Differentiation in R: FTIR Spectra
- Plotting XRD (X-Ray Diffraction) Data
- Using lm() and predict() to apply a standard curve to Analytical Data
- Working with Spatial Data
- Customizing Maps in R: spplot() and latticeExtra functions
- Converting Alpha-Shapes into SP Objects
- Some Ideas on Interpolation of Categorical Data
- Visual Interpretation of Principal Coordinates (of) Neighbor Matrices (PCNM)
- Visualizing Random Fields and Select Components of Spatial Autocorrelation
- Generation of Sample Site Locations [sp package for R]
- Target Practice and Spatial Point Process Models
- Ordinary Kriging Example: GRASS-R Bindings
- Point-process modelling with the sp and spatstat packages
- Simple Map Creation
- Comparison of PSA Results: Pipette vs. Laser Granulometer
- GRASS GIS: raster, vector, and imagery analysis
- Generic Mapping Tools: high quality map production