library(spatstat)
###The spatstat package offers a lot of super handy, super simple functions to quickly and effectively analyze spatial data
###This code is designed to help users become acquainted with basic functions for analyzing planar point patterns
###Although these plots may look not very polished, they are intended to offer initial analysis of spatial data
######How to make a spatial point pattern -- ppp()
#####Showing this merely as a demonstration so that others can bring in their own ppp
#####A ppp is nothing more than a georeferenced event
####There is no additional information other than an event occurred here
x <- rnorm(100, 0.5, 0.25)
y <- rnorm(100, 0.5, 0.25)
ppp <- ppp(x, y, c(0,1), c(0,1))
plot.ppp(ppp) ####Notice how points falling outside the window of view are marked differently than those within
####Let's work with some built in data....
####This is already a ppp, detailing where swedish pine trees are found
attach(swedishpines)
swedes <- swedishpines ##Copy the data
str(swedes) ###Pay attention that the class of the data is "ppp" so we're good to go
plot.ppp(swedes)
summary(swedes)
#####Find the distance between nearest neighbors
neighbors <- nndist(swedes)
summary(neighbors)
boxplot(neighbors)
#####Let's say we wanna make a plot with arrows pointing to nearest neighbors
near <- nnwhich(swedes)
ref <- swedes[near]
plot.ppp(swedes)
arrows(swedes$x, swedes$y, ref$x, ref$y, angle = 15, length = 0.15, code = 2, col = "red")
###Let's say we want to make a down and dirty density plot (a DDDP of sorts) of our spatial data
plot(density(swedes, 1)) ###the function takes the formula density(ppp, bw)
plot(density(swedes, 5)) ###bw refers to the smoothing bandwith
plot(density(swedes, 10)) ###think of it as a smoothing factor of sorts
plot(density(swedes, 25)) ###see how the density plots compare with different smoothing bandwiths
####Or perhaps a contour plot
contour(density(swedes), axes = FALSE)
contour(density(swedes), axes = TRUE)
contour(density(swedes, 5), axes = FALSE) ###You can also add a smoothing bandwith for higher resolution
###Let's say you would like to split the Swedish pines up into quadrats
quadrats <- quadratcount(swedes, nx = 4, ny = 4) ###This splits our ppp into a 4x4 matrix -- 16 quadrats
quadrats
table(quadrats) ###Table enumerates the number of quadrats with a given number of pines
###The top row is the number of swedish pines in a quadrat.
###The bottom row is the number of quadrats that number of swedish pines
###Let's say you would like to put those enumerations on the plot.ppp image
plot(swedes)
plot(quadrats, add = TRUE, cex = 2)
###Let's say you'd like to test for spatial clustering or regular spacing
###Note that there are several parameters that I did not include
###These perform a quadrat test using a Chi-Squared distribution for a 2x2 matrix (4 quadrats)
###Note that the null hypothesis is that the data are completely spatially random (CSR)
###We can easily set the alternative hypothesis as two.sided, regular, or clustered
quadrat.test(swedes, nx = 2, ny = 2, alternative = "two.sided")
quadrat.test(swedes, nx = 2, ny = 2, alternative = "regular")
quadrat.test(swedes, nx = 2, ny = 2, alternative = "clustered")
###Nearest neighbor functions
###The math behind these functions can be quite complicated.
###I offer a brief explanation, but I recommend any number of supplemental info for deeper understanding
###The G function compares radial distances to the nearest neighbor for each event
###Where G < Gest we should see clustering
###Where G > Gest we should see regular spacing
G <- Gest(swedes)
plot(G)
###The F function compares radial distances to the nearest neighbor from a point in a regularly spaced grid o
###Where F < Fest we should see regular spacing
###Where F > Fest we should see clustering
F <- Fest(swedes)
plot(F)
###The J function compares the F and G functions by looking at (1-G)/(1-F)
###For J = 1 --> CSR
###For J < 1 --> Clustering
###For J > 1 --> Regular spacing
###CAREFUL!! Tend to become unstable for larger radii
J <- Jest(swedes)
plot(J)
###Super robust test
###Performs a pairwise comparison for the radii of all points from each other
###Since we have 71 points, it will make 70 comparison 71 times.
###Needless to say, this take A LOT of computing power. BUT it is incredibly robust because of the pairwise comparisons
L <- Lest(swedes)
plot(L)
###Let's say we want to make some confidence intervals, using 499 simulations of each function
envelope.G <- envelope(swedes, fun = Gest, nsim = 499)
plot(envelope.G)
envelope.F <- envelope(swedes, fun = Fest, nsim = 499)
plot(envelope.F)
envelope.J <- envelope(swedes, fun = Jest, nsim = 499)
plot(envelope.J)
envelope.L <- envelope(swedes, fun = Lest, nsim = 499)
plot(envelope.L)
###Pretty convenient, no?
par(mfrow=c(2,2))
plot(envelope.G)
plot(envelope.F)
plot(envelope.J)
plot(envelope.L)
####How to make a down and dirty Voronoi plot
plot(dirichlet(swedes))
plot(swedes, add = TRUE)
####For the brave of heart....3D plots made quickly
threeDPpp <- pp3(runif(50), runif(50), runif(50), box3(c(0,1)))
plot(threeDPpp)