Submitted by dylan on Thu, 2006-10-26 04:38.

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.