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)