library(spatstat)
library(sf)

Read R code for reading Urkiola shape file and converting to ppp object

# Tree locations
urkiola.sf <- st_read("datafiles/urkiola/urkiola.shp")
## Reading layer `urkiola' from data source 
##   `/Users/jessieeastburn/Documents/Fall 2021/GEOG 6000/datafiles/urkiola/urkiola.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1245 features and 5 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 2.1 ymin: 0.1 xmax: 219.1 ymax: 149.9
## CRS:           NA
# Park boundary
urkiola.win <- st_read("datafiles/urkiola/urkiolaWindow.shp")
## Reading layer `urkiolaWindow' from data source 
##   `/Users/jessieeastburn/Documents/Fall 2021/GEOG 6000/datafiles/urkiola/urkiolaWindow.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 3 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 0.05 ymin: 0.05 xmax: 219.95 ymax: 149.95
## CRS:           NA
urkiola.win <- as.owin(urkiola.win)
# First extract coordinates
urkiola.crds <- st_coordinates(urkiola.sf)
## Convert to ppp
urkiola.ppp <- ppp(x = urkiola.crds[, 1], y = urkiola.crds[,2], 
                   marks = as.factor(urkiola.sf$tree), window = urkiola.win)
plot(urkiola.ppp)

Give the intensity (spatial density) of the combined set of two tree species by using the summary() function

summary(urkiola.ppp)
## Marked planar point pattern:  1245 points
## Average intensity 0.06564029 points per square unit
## 
## Coordinates are given to 1 decimal place
## i.e. rounded to the nearest multiple of 0.1 units
## 
## Multitype:
##       frequency proportion  intensity
## birch       886  0.7116466 0.04671269
## oak         359  0.2883534 0.01892760
## 
## Window: polygonal boundary
## single connected closed polygon with 44 vertices
## enclosing rectangle: [0.05, 219.95] x [0.05, 149.95] units
##                      (219.9 x 149.9 units)
## Window area = 18967 square units
## Fraction of frame area: 0.575

The average intensity 0.06564029 points per square unit

Plot the distribution of the two tree species as separate figures.

oak = plot(split(urkiola.ppp)$oak, main = "Oak Tree Distribution")

birch = plot(split(urkiola.ppp)$birch, main = "Birch Tree Distribution")

Make plots of the kernel density of each species, and compare these plots.

oak.k = split(urkiola.ppp)$oak
birch.k = split(urkiola.ppp)$birch
plot(density(oak.k, sigma = 45))

plot(density(birch.k, sigma = 45))

oak.k.bw = bw.diggle(oak.k)
plot(density(oak.k, sigma = oak.k.bw))

birch.k.bw = bw.diggle(birch.k)
plot(density(birch.k, sigma = birch.k.bw))

This does not suggest that the two species co-occur within the park.

Using Ripley’s K function, examine the individual distributions of the two species for non-random spatial organization.

oak.kest = Kest(oak.k)
plot(oak.kest)

birch.kest = Kest(birch.k)
plot(birch.kest)

Use the envelope() function to test the obtained Ripley’s K against a hypothesis of Complete Spatial Randomness.

oak.kest.mc <- envelope(oak.k, fun = 'Kest', nsim = 99, verbose = FALSE, global = TRUE)
plot(oak.kest.mc, shade = c("hi", "lo"))

summary(oak.kest.mc)
## Simultaneous critical envelopes for K(r)
## and observed value for 'oak.k'
## Obtained from 99 simulations of CSR
## Alternative: two.sided
## Envelopes computed as theoretical curve plus/minus maximum simulated value of 
## maximum absolute deviation
## Significance level of Monte Carlo test: 1/100 = 0.01
## Data: oak.k
birch.kest.mc <- envelope(birch.k, fun = 'Kest', nsim = 99, verbose = FALSE, golbal = TRUE)
plot(birch.kest.mc, shade = c("hi", "lo"))

summary(birch.kest.mc)
## Pointwise critical envelopes for K(r)
## and observed value for 'birch.k'
## Obtained from 99 simulations of CSR
## Alternative: two.sided
## Upper envelope: pointwise maximum of simulated curves
## Lower envelope: pointwise minimum of simulated curves
## Significance level of Monte Carlo test: 2/100 = 0.02
## Data: birch.k

Since the observed line falls above the theoretical line, both the birch and oak species are clustered.

Test the spatial dependence between the two distributions to test if the two species truly co-occur or not.

urkiola.kc <- envelope(urkiola.ppp, 
                       Kcross, 
                       i = "oak", 
                       j = "birch",
                       nsim = 99, 
                       verbose = FALSE)
plot(urkiola.kc)

The combined distribution is random because the observed curve is within the envelope.