Overlay Operations, Mixing spatstat and sp Operations

In spatial statistics, we often handle data measured on different spatial units e.g., point-process data (where the underlying stochastic process governs where the events happen), geostatistical data (where the spatial location of data are pre-determined) and aggregated data (where data pertain to areal regions in space). Overlay operations are one way of moving information from one spatial unit to another.

Suppose we have access to a dataset consisting of points in an obervation window (either point-process or geostatistical data):

library(spatstat)
library(spatstat.utils)
library(sp)
win <- owin()
patt <- rpoispp(lambda = 500, win = win)

Suppose further that we also have some environmental covariate data:

xseq <- seq(0.05, 0.95, length.out = 10)
yseq <- xseq
spix <- SpatialPixels(SpatialPoints(expand.grid(xseq,
    yseq)))
spoly <- as(spix, "SpatialPolygons")
df <- data.frame(deprivation = runif(100, 0, 100),
    elevation = runif(100, 100, 300))
spdf <- SpatialPolygonsDataFrame(spoly, data = df,
    match.ID = FALSE)

spplot(spdf, "deprivation")
spplot(spdf, "elevation")

We often need to find out the values of the covariates at each of the point locations. This can be done using an overlay operation.

We first convert the ppp object into a SpatialPoints object:

pts <- as(patt, "SpatialPoints")

then perform the overlay

olay <- over(pts, geometry(spdf))
olay

so the relevant covariate data for each point is given by,

cvdata <- spdf@data[olay, ]
head(cvdata)

we can do a visual check to make sure that this is correct as follows:

spplot(spdf, "deprivation")
pts <- SpatialPointsDataFrame(pts, data = cvdata)
spplot(pts, "deprivation")

The patterns look similar.

If we wanted to count the number of events in each polygon in spdf, we could do the following:

spdf$counts <- sapply(1:length(spdf), function(i) {
    sum(olay == i)
})
spplot(spdf, "counts")
plot(pts, pch = 1)
plot(spdf, add = TRUE)