Reference and Disclaimer: For this exercise, I have used explanation and major chunk of codes from the book Spatial Ecology and Conservation Modeling by Robert Fletcher & Marie-Josée Fortin.

Load the necessary packages

rm(list = ls())

## First specify the packages of interest
packages = c("rgdal","spatstat","pgirmess","ncf","spdep", "geoR","gstat","raster", "waveslim","fields","vegan","reshape2","adespatial","ggplot2")

## Now load or install&load all
package.check <- lapply(
  packages,
  FUN = function(x) {
    if (!require(x, character.only = TRUE)) {
      install.packages(x, dependencies = TRUE)
      library(x, character.only = TRUE)
    }
  }
)

Load the raster data and the point data to be analysed

lc_wt41 <- raster("./SpaStatAssign/ccapsubset_wetland_wi41.tif")
point.data = readOGR("./SpaStatAssign/SpaStatData.gdb","WhiteIbis41pointsAugust2021")

Checking the projection for the landcover data: wetland

proj4string(lc_wt41)
## [1] "+proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs"

Checking the projection for the point locations of White Ibis-14

proj4string(point.data)
## [1] "+proj=longlat +datum=WGS84 +no_defs"

Since both are in different CRS we have to transform projection for point data from the raster crs.

ibis.crs <- CRS("+proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs") #defines new crs
point.data <- spTransform(point.data, ibis.crs)

On checking the file we see that the last two observations are erroneous and that there are many fields that contains no information, hence, sub-setting our dataframe.

ibis41 <- data.frame(point.data)
ibis41 <- ibis41[1:128,c(3,9:11)]

Vizualizing the raster and the ibis point location

plot(lc_wt41)
points(ibis41[ibis41$SymbolID==0, c("coords.x1","coords.x2")], col="black")

We now use extract() from raster package to extract covariate values from layers corresponding to the locations from the survey data of IBIS and merge them with the point data

coords <- cbind(ibis41$coords.x1, ibis41$coords.x2)
colnames(coords) <- c("x", "y")
distmat <- as.matrix(dist(coords))
#Maximum distance to consider in correlogram (~1/2 to 2/3 total dist)
maxdist <- 2/3*max(distmat)
land.cov <- extract(x=lc_wt41, y=coords)
#combine
ibis41_df <- cbind(ibis41,land.cov)

Creating distance-band weights

knn1 <- knearneigh(coords)
k1 <- knn2nb(knn1)
critical.threshold <- max(unlist(nbdists(k1,coords)))
#Computing distance-band weights
nb.dist.band <- dnearneigh(coords, 0, critical.threshold)
# #Weight summary
# summary(nb.dist.band)
#Connectivity histogram
dist.band.card <- card(nb.dist.band)
# dist.band.card

ggplot() +
  geom_histogram(aes(x=dist.band.card),bins = 15) +
  xlab("Number of Neighbors")

#Connectivity graph
plot(nb.dist.band, coords, lwd=.2, col="blue", cex = .5)

K nearest neighbour analysis

k6 <- knn2nb(knearneigh(coords, k = 6))
k6.card <- card(k6)
ggplot() +
  geom_histogram(aes(x=k6.card), binwidth = .01) +
  xlab("Number of Neighbors")

plot(k6, coords, lwd=.2, col="blue", cex = .5)

Thiessen polygon

theissen <- dismo::voronoi(coords)
theissen <- terra::voronoi(coords)

#plot
spplot(theissen, "id")

#or plot
plot(theissen)
points(coords, cex=0.5,col="red")

We use correlog() from the pgirmess package specifying to use method = “Moran”

#Correlogram with pgirmess
correlog.pgirmess <- pgirmess::correlog(coords, ibis41_df$land.cov, method="Moran",
                                        nbclass=14, alternative = "two.sided")

#Moran and P values for each distance class
round(correlog.pgirmess,2)
## Moran I statistic 
##       dist.class  coef p.value    n
##  [1,]    4308.58  0.80    0.00 6014
##  [2,]   12925.75 -0.07    0.49  634
##  [3,]   21542.91  0.47    0.00  608
##  [4,]   30160.07 -0.06    0.08   74
##  [5,]   38777.24 -0.25    0.00  814
##  [6,]   47394.40 -0.46    0.00  334
##  [7,]   56011.56 -0.57    0.06   30
##  [8,]   64628.73 -0.15    0.62   62
##  [9,]   73245.89 -0.38    0.00  348
## [10,]   81863.05  0.64    0.00  904
## [11,]   90480.22  0.16    0.00 1000
## [12,]   99097.38 -0.13    0.04  296
## [13,]  107714.55 -0.45    0.00 2938
## [14,]  116331.71 -0.87    0.00 2200

Correlog() creates a matrix that contains information like distance class, the Moran coefficient for that distance, the p-value and the sample size used for that distance. We now plot the correlogram

plot(correlog.pgirmess[,1], correlog.pgirmess[,2],
     xlab="Distance (m)", ylab="Moran's I", main="Correlogram with pgirmess", pch=21,col="red", bg="lightgreen")
abline(h=0, col="blue")

The plot suggests that there exist positive spatial dependence till 20.5 km followed by negative spatial dependence and the cycle goes on. Something is not right here.

Interpret spatial dependence with Moran’s I using ncf package

plot(correlog.ncf)
abline(h=0)

WOW!!! Things are bad… very bad

#Spline correlogram with 95% pointwise bootstrap confidence intervals
spline.corr <- spline.correlog(x = ibis41_df$coords.x1, y = ibis41_df$coords.x2, z = ibis41_df$land.cov,
                               xmax = maxdist, resamp=99, type="boot")
## 10  of  99 
20  of  99 
30  of  99 
40  of  99 
50  of  99 
60  of  99 
70  of  99 
80  of  99 
90  of  99

Interpret spatial dependence using spline.correlog()

plot(spline.corr, main="Spline correlogram")

Can’t be any worst

May be the error is in choosing the right method or the covariate…

Thanks!!