Geog 380: Point Pattern Analysis

Today we’ll start looking at our first topic in spatial statistics: point pattern analysis! Often, we might have a dataset of point locations. These can take many forms: you might want point data of business locations, sightings of endangered species, coordinates of fires or other disasters, or even the exact locations of spectators in a sports arena. For all of these examples, it would be helpful to know some information about how the points cluster in space

library(sf)

Linking to GEOS 3.9.1, GDAL 3.4.3, PROJ 7.2.1; sf_use_s2() is TRUE

library(maptools)

Loading required package: sp
Checking rgeos availability: FALSE
Please note that 'maptools' will be retired by the end of 2023,
plan transition at your earliest convenience;
some functionality will be moved to 'sp'.
    Note: when rgeos is not available, polygon geometry     computations in maptools depend on gpclib,
    which has a restricted licence. It is disabled by default;
    to enable gpclib, type gpclibPermit()

library(spatstat)

Loading required package: spatstat.data
Loading required package: spatstat.geom
spatstat.geom 3.0-3
Loading required package: spatstat.random
spatstat.random 2.2-0
Loading required package: spatstat.core
Loading required package: nlme
Loading required package: rpart
spatstat.core 2.4-4
Loading required package: spatstat.linnet
spatstat.linnet 2.3-2

spatstat 2.3-4       (nickname: ‘Watch this space’) 
For an introduction to spatstat, type ‘beginner’ 

library(raster)   

Attaching package: ‘raster’

The following object is masked from ‘package:nlme’:

    getData

library(dplyr)   

Warning: package ‘dplyr’ was built under R version 4.2.2
Attaching package: ‘dplyr’

The following objects are masked from ‘package:raster’:

    intersect, select, union

The following object is masked from ‘package:nlme’:

    collapse

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

library(spData) 

To access larger datasets in this package, install the
spDataLarge package with: `install.packages('spDataLarge',
repos='https://nowosad.github.io/drat/', type='source')`

library(urbnmapr)
library(ggplot2)

Warning: package ‘ggplot2’ was built under R version 4.2.3

library(tidyr)

Attaching package: ‘tidyr’

The following object is masked from ‘package:raster’:

    extract

load(url("https://github.com/mgimond/Spatial/raw/main/Data/ppa.RData"))


#In this tutorial, we need our data to be in spatstat format
#the spatstat package requires a ppp type object, not an sf object
#don't worry too much about this now, but take a look at the following printout to verify that it's correct
starbucks

Planar point pattern: 171 points
window: rectangle = [648032.3, 917740.7] x [4609785, 4748107] units

We’ll get back to the analysis in a second, but if my object isn’t already a ppp object? We’ll have to convert it, like this:

*Note that “marks” are attributes associated with a ppp object. We won’t need any attributes in this example, so I’ve just deleted them to avoid confusion.

#Note that we will only use this data to illustrate the ppp conversion - we can delete afterwards

setwd("~/Binghamton/geog380/NYS_Place_Points_SHP")

Warning: The working directory was changed to C:/Users/mhaller/Documents/Binghamton/geog380/NYS_Place_Points_SHP inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

points <- st_read("NYS_Place_Points.shp")

Reading layer `NYS_Place_Points' from data source 
  `C:\Users\mhaller\Documents\Binghamton\geog380\NYS_Place_Points_SHP\NYS_Place_Points.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 4571 features and 8 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 108752.5 ymin: 4484510 xmax: 756007.4 ymax: 4984062
Projected CRS: NAD83 / UTM zone 18N

points_ppp <- as.ppp(points)

Warning: only first attribute column is used for marks

#marks are attributes - we need any for this example, so set to null
marks(points_ppp) <- NULL

plot(points_ppp, main=NULL, cols=rgb(0,0,0,.2), pch=20)

Returning to the starbucks data:

#We just want the points, not the attributes attached to them
#this code removes the attributes (variables)

library(spatstat)
marks(starbucks)  <- NULL

#This line allows us to bind the points to the MA boundary file
Window(starbucks) <- ma

#For simplicity, we will just use base R plotting for now
plot(starbucks, main=NULL, cols=rgb(0,0,0,.2), pch=20)

First, we’ll calculate some basic descriptives statistics using our point data. Let’s start with the mean center of the points:

xy <- coords(starbucks)

#This allows us to calculate the mean center
mc <- apply(xy, 2, mean)
#select just the mean center
mean <- cbind(mc[1], mc[2])
#plot it
plot(starbucks, main=NULL, cols=rgb(0,0,0,.2), pch=20)
points(mean, col='red', pch = "*", cex = 3)

The red * is the central most point on the plot, based on the distribution of points. It is calculated based on the average of the x and y values.

Next we’ll calculate the standard distance! This is a measure of the variance of the distance from each point to the mean center - think of it like a standard deviation for point data.

#don't worry too much about the complicated code here
sd <- sqrt(sum((xy[,1] - mc[1])^2 + (xy[,2] - mc[2])^2) / nrow(xy))

#now we'll plot the standard distance as a circle
#don't worry about this code, I won't ask you to reproduce it
bearing <- 1:360 * pi/180
cx <- mc[1] + sd * cos(bearing)
cy <- mc[2] + sd * sin(bearing)
circle <- cbind(cx, cy)

#now we can plot it
plot(starbucks, main=NULL, cols=rgb(0,0,0,.2), pch=20)
points(mean, col='red', pch = "*", cex = 3)

lines(circle, col='red', lwd=2)

Those are two very basic descriptive statistics for point data. Next, we’ll move on to measures of point density.

Because we’ll use the population data in a bit, we’ll first start by transforming the data.

#Let's take a look at the population density raster, pop
hist(pop, main=NULL, las=1)

This data looks very skewed! Data that is approximately normally distributed would work better here - let’s transform our variable

pop.lg <- log(pop)
hist(pop.lg, main=NULL, las=1)

#One way to calculate the intensity of the point process is with the summary() function
summary(starbucks)

Planar point pattern:  171 points
Average intensity 8.268627e-09 points per square unit

Coordinates are given to 1 decimal place
i.e. rounded to the nearest multiple of 0.1 units

Window: polygonal boundary
single connected closed polygon with 128 vertices
enclosing rectangle: [623157.2, 921923.8] x [4602709, 4756659] units
                     (298800 x 154000 units)
Window area = 20680600000 square units
Fraction of frame area: 0.45

Remember, point pattern intensity is a global or first-order characteristic

Next, we will calculate a measure of density using a quadrat grid. In this way, we can look at the number of points in each quadrat.

Q <- quadratcount(starbucks, nx= 6, ny=3)

plot(starbucks, pch=20, cols="grey70", main=NULL)  # Plot points
plot(Q, add=TRUE)  # Add quadrat grid

Here, we have quadrat point counts. Let’s calculate the density of points in each quadrat!

# Compute the density for each quadrat
Q.d <- intensity(Q)

# Plot the density
plot(intensity(Q, image=TRUE), main=NULL, las=1)  # Plot density raster
plot(starbucks, pch=20, cex=0.6, col=rgb(0,0,0,.5), add=TRUE)  # Add points

In this plot, blue areas are low density, pink/purple are intermediate, and yellow/orange are high density regions. The units here are points per square meter. Those units aren’t very nice, are they? We can rescale them to kilometers like this:

#rescale all of our data:

starbucks.km <- rescale(starbucks, 1000, "km")
ma.km <- rescale(ma, 1000, "km")
pop.km    <- rescale(pop, 1000, "km")
pop.lg.km <- rescale(pop.lg, 1000, "km")

Now we can plot it again!

# Compute the density for each quadrat (in counts per km2)
Q   <- quadratcount(starbucks.km, nx= 6, ny=3)
Q.d <- intensity(Q)

# Plot the density
plot(intensity(Q, image=TRUE), main=NULL, las=1)  # Plot density raster
plot(starbucks.km, pch=20, cex=0.6, col=rgb(0,0,0,.5), add=TRUE)  # Add points

Now our scale looks much nicer. Another thing to note is that we do not have to use uniformly square quadrats. Here, I’ll show an example of using the population density object to calculate new quadrats. The regions defined by population density are called “tessellated surfaces”.

brk  <- c( -Inf, 4, 6, 8 , Inf)  # Define the breaks
Zcut <- cut(pop.lg.km, breaks=brk, labels=1:4)  # Classify the raster
E    <- tess(image=Zcut)  # Create a tesselated surface

#Now we'll plot it
#note that las=1 tells R to show the legend labels upright (as opposed to flipped 90 degrees to the left - that is the default)
plot(E, main="", las=1)

#Repeat the same density analysis as before, but with the new quadrats
#we just need to tell R which tessellated object to use for quadrats
Q   <- quadratcount(starbucks.km, tess = E)  # Tally counts
Q.d <- intensity(Q)  # Compute density

# Plot the density
plot(intensity(Q, image=TRUE), main=NULL, las=1)  # Plot density raster
plot(starbucks.km, pch=20, cex=0.6, col=rgb(0,0,0,.5), add=TRUE)  # Add points

Advanced: Kernel Density

Another measure of density is called Kernel Density. Rather than use quadrats, this form of density is calculated using a moving sub-region window, called a Kernel window. This is quite similar to the focal operations that we discussed in our raster unit. This measure of density is nice because it gives us a more detailed picture of point density across the study unit. We can also convert the output of the kernel density function into a raster object, which would allow us to use raster methods on it.

library(raster)
#This is really easy - calculate kernel density
dens <- density(starbucks.km)

#convert into a raster
rast <- raster(dens)

#use the native raster plotting method to plot it
plot(rast)

#Another cool feature of kernel density: we can plot it using contour lines

plot(dens, main=NULL, las=1)
contour(dens, add=TRUE)

We can also change the bandwidth of the kernel density estimation. By default, R chooses the bandwidth for us, based on the size of the plotting window. We can re-specify the bandwidth using the sigma = argument. Here, I show an example where I change the bandwidth to 50 km.

dens2 <- density(starbucks.km, sigma=50) # Using a 50km bandwidth
plot(dens2, main=NULL, las=1)
contour(dens2, add=TRUE)

Let’s return to complete spatial randomness for a moment. If we want to test whether our point data was generated by complete spatial randomness, there is a very simple function that will allow us to test this. Remember that, whenever we set up statistical tests, we always want to set up a null and alternative hypothesis. Ours might look like this:

Null: The point data was generated by a random spatial process Alternative: The point data was not generated by a random spatial process

Now we can perform our test. Since this is not a statistics class, I won’t go into the technical details here, but it should be easy enough to interpret the results of the test. If you want to learn more about what is going on under the hood of this test, you can always look it up later.

quadrat.test(starbucks.km)

Warning: Some expected counts are small; chi^2 approximation may be inaccurate

    Chi-squared test of CSR using quadrat counts

data:  starbucks.km
X2 = 695.5, df = 19, p-value < 2.2e-16
alternative hypothesis: two.sided

Quadrats: 20 tiles (irregular windows)

Our test statistic is a Chi-squared value - we know that we have sufficient evidence to reject the null hypothesis if our test statistic is large enough. With a P-value well under 0.05, we can safely reject the null hypothesis and conclude that our point pattern was not generated by a random spatial process.

Distance-Based Methods

Let’s dig deeper into the pattern of clustering. The first method we’ll look at is called the Average Nearest Neighbor method (or ANN). ANN is essentially the average distance between all points within a dataset and their individual nearest point (known as Nearest-Neighbor Distance (NND)). We can compare the actual ANN pattern in our data with the expected ANN pattern under CSR as a global indicator of the overall clustering in our data.

This is similar to quadrat analysis, but it uses distance instead.

r #Let's calculate the average nearest neighbor distance! mean(nndist(starbucks.km, k = 1))

[1] 3.275492

This tells us that each point is, on average, 3.28 km away from its nearest neighbor

```r #Next we’ll plot the ANN for different numbers of neighbors (K)

ANN <- apply(nndist(starbucks.km, k=1:100),2,FUN=mean) plot(ANN ~ eval(1:100), type=“b”, main=NULL, las=1) ```

This shows us that the ANN increases as the number of neighbors increases, but it doesn’t give us too much info. ANN is particularly helpful if we want to compare clustering in one point dataset to another.

Our next statistic will give us a lot more details about the pattern of clustering in our data.

Essentially, The black line (K-pois) represents the theoretical
K-function under the null hypothesis that the points are completely randomly distributed (CSR/IRP). Where K falls under the theoretical K-pois line, the points are deemed more dispersed than expected at distance r. Where K falls above the theoretical k-pois line, the points are deemed more clustered than expected at distance r.

Calculation: Ripley’s K function essentially summarises the distance between points for all distances using radial distance bands. The calculation is relatively straightforward (Wilkin, 2015):

-For point event A, count the number of points inside a buffer (radius) of a certain size. Then count the number of points inside a slightly larger buffer (radius). -Repeat this for every point event in the dataset. -Compute the average number of points in each buffer (radius) and divide this by the overall point density. -Repeat this using points drawn from a Poisson random model for the same set of buffers (don’t worry too much about this model - you just need to be able to interpret the plot). -Compare the observed distribution with the distribution with the Poisson distribution.

Note that this function will be very slow for large datasets.

plot(Kest(starbucks.km, correction = "border"))

At a radius of about 0-30 km, the points are much more clustered than we would expect, given a random spatial process. Greater than 30km, the pattern is a bit less clear, and may even become more dispersed than expected. Overall, this gives us a very clear picture that our data is unlikely to have been generated randomly.

One caveat with this analysis: it assumes that our point process is stationary (i.e. that it doesn’t change across space). This may pose some validity concerns for our K-function (the same is also true for ANN analysis).

Another issue is that it can also be challenging to see small differences between the actual and predicted values. The L-function is a corrected version of the K-function that transforms the K function into a straight line - this should make the difference easier to visualize.

plot(Lest(starbucks.km, correction = "border"))

Resources

Dorman, M. (2022). Chapter 13 Point pattern analysis (in preparation). Retrieved from: https://geobgu.xyz/r/point-pattern-analysis.html.

Gimond, M. (2022). Chapter 11: Point Pattern Analysis. Retrieved from: https://mgimond.github.io/Spatial/chp11_0.html.

Peeples, M. (N. D.). Point Pattern Analysis. Retrieved from: https://www.mattpeeples.net/modules/PointPattern.html

R Spatial, (2021). Point Pattern Analysis. Retrieved from: https://rspatial.org/terra/analysis/7-pointpat.html.

Wilkin, J. (2015). 8 Analysing Spatial Patterns III: Point Pattern Analysis. Retrieved from: https://jo-wilkin.github.io/GEOG0030/coursebook/analysing-spatial-patterns-iii-point-pattern-analysis.html.