Introduction to Spatial Autocorrelation

Waldo Tobler is known for many things, one of them being Tobler’s First Law of Geography. According to Tobler, “…everything is related to everything else, but near things are more related than distant things” (Tobler, 1970). By assuming that this law is correct, spatial analysis concepts like that of spatial autocorrelation can be explored with greater certainty. Spatial autocorrelation is a term used to describe how a variable is correlated with itself through space. Positive spatial autocorrelation occurs when objects with similar values are closer together (i.e. clustered) while negative spatial autocorrelation occurs when objects with dissimilar values are closer together (i.e. dispersed).

Source: https://bit.ly/topic7notes

The Moran’s I statistic will be used to quantify the magnitude of spatial autocorrelation experienced. By using R to perform spatial autocorrelation analysis, an array of statistical, cartographic, and plotting packages will be readily available to the user. Moreover, the user may contribute towards the development of R as it is an open-source programming language.

Using Spatial Autocorrelation for Vegetation Resources Inventory (VRI)

The VRI is a provincial forest inventory program designed to help secure economic prosperity and environmental sustainability. (Province of British Columbia, n.d.). Spatially, the VRI dataset is split into multiple polygons representing homogenous stands within the ownership of the Province of British Columbia. Stand attributes such as stand height and stand age are estimated via photo interpretation and ground sampling.

In order to determine whether spatial patterns exist within BC’s managed forests, techniques to determine spatial autocorrelation must be employed. For example, if there were recent catastrophic events (such as mountain pine beetle) that occurred in the forest, it would be beneficial to study which areas were most affected why. Spatial autocorrelation can help answer part of questions like these; therefore, helping the Province of British Columbia make durable decisions on the land that support a healthy and competitive forest sector.

How to format and map VRI data

This tutorial will be focusing on the Greater Victoria Water Supply Area (GVWSA). The spatial autocorrelation will be calculated for two variables: estimated stand height and estimated stand age (projected to 2020 from estimated establishment date). It’s important to study stand height and stand age because they are both valuable metrics for estimating stand health and value (Hassegawa et al., 2015). By using stand height and stand age as variables, the fundamentals of how spatial autocorrelation can be used to study VRI data will be more effectively communicated. Using more complex variables will require the need for increased interpretation of the results. If future studies in the GVWSA include either of these variables, it would be highly valuable to know if spatially autocorrelation is present. If it is, then it may be an indication that there is something of interest in the distribution of the variable being tested (Haining, 2001). Moreover, further investigation would be required to understand the reasons behind the observed spatial variation.  

Install and load the proper libraries.

#Install packages
#install.packages("spdep")
#install.packages("raster")
#install.packages("rgdal")
#install.packages("tmap")
#install.packages("shinyjs")

#Load libraries
library(rgdal)
library(tmap)
library(spdep)
library(raster)
library(shinyjs)

The next step is to set the working directory and read in the desired data.

#Set working directory
dir <- "C:/Users/Lucas Aubert/OneDrive/University/Year 4/Fall 2020/GEOG 418/Labs/Lab_3"
setwd(dir)

#Read in shapefile
VRI <- readOGR(dsn = "./WatershedVRI/WatershedVRI.shp", layer = "WatershedVRI")
#Project data to NAD83 UTM Zone 10N
VRI <- spTransform(VRI, CRS("+init=epsg:26910"))

Change the field names to increase readability.

#View field names
head(VRI@data)

#Select desired fields
vriCleanCols <- c("FID_VEG_CO", "POLYGON_ID", "PROJ_AGE_1",
                  "SITE_INDEX", "SPECIES__4", "SPECIES__5",
                  "PROJ_HEI_1", "SPECIES_PC", "SPECIES__6",
                  "VRI_LIVE_S", "BASAL_AREA", "WHOLE_STEM",
                  "CROWN_CL_1")

#Create subset of data
vriClean <- VRI[,vriCleanCols]

#Rename columns to make them more readable
newNames <- c("FID", "PolyID", "Stand_Age", "Site_Index",
              "CoDom_Sp", "Dom_Sp", "Stand_HT", "DomSP_Perc", 
              "CDomSP_Perc", "Stand_Dens", "Stand_BA", "Stand_StemBio", "Stand_CrownCl")

colnames(vriClean@data) <- newNames
head(vriClean@data)

#Meta Data (https://www2.gov.bc.ca/assets/gov/farming-natural-resources-and-industry/forestry/stewardship/forest-analysis-inventory/data-management/standards/vegcomp_poly_rank1_data_dictionaryv5_2019.pdf)
# FID = Field ID
# PolyID = VRI Polygon ID
# Stand_Age = Estimated stand age projected to 2020 from estimated establishment date
# Site_Index = A value to estimate site quality. This describes the height that the stand could grow to by age 50 in meters.
# CoDom_Sp = The species code for the co-dominant tree species. Full list of codes: https://www.for.gov.bc.ca/hfp/publications/00026/fs708-14-appendix_d.htm
# Dom_Sp = The species code for the dominant tree species. Full list of codes: https://www.for.gov.bc.ca/hfp/publications/00026/fs708-14-appendix_d.htm
# Stand_HT = The estimated height for the stand
# DomSP_Perc = The estimated percentage of the dominent species
# CDomSP_Perc = The estimated percentage of the co-dominent species
# Stand_Dens = Estimated density of stand (Stems per hectare)
# Stand_BA = Estimated Basal area of the stand (square meters)
# Stand_StemBio = Estimated stand level stem biomass (tonnes per hectare)
# Stand_CrownCl = The percentage of ground area covered by tree crowns

Use this portion of code to remove any NA values from the stand height and stand age variables.

vriClean <- vriClean[!is.na(vriClean@data$Stand_HT) & !is.na(vriClean@data$Stand_Age), ]

Visualize the stand height and dominant species variables to gain a preliminary understanding of the data.

#Changes the map view to interactive
tmap_mode("plot")

#Tool for selecting colour palettes
tmaptools::palette_explorer() 
#Create a map of the study area using the stand height attribute
#n = number of desired classes
#type ?tm_polygons to see all the possible arguments for customizing the polygons
map_StdHT <- tm_shape(vriClean) + 
  tm_polygons(col = "Stand_HT", 
              title = "Stand Height (meters)", 
              style = "jenks", 
              palette = "viridis", n = 6)
#View map
map_StdHT

There appears to be relatively taller stands in the northeast and southeast regions of the GVWSA. Note: The large area in the east-central portion of the GVWSA containing a ‘short’ stand is the Sooke Lake Reservoir.

#Create a map of the study area using the stand age attribute
#Create stand age map
map_StdAGE <- tm_shape(vriClean) + 
  tm_polygons(col = "Stand_Age", 
              title = "Stand Age (years)", 
              style = "jenks", 
              palette = "viridis", n = 6)
#View map
map_StdAGE

Older stands appear to be clustered in the north-central region of the GVWSA.

How to Quantify Nearby Stands?

In order to conclude whether the chosen stand attributes are spatially autocorrelated or not, the neighbours of each stand need to be defined. Each neighbour is then given a weight that distinguishes it from the rest of the stands in the study area. By weighting the neighbours, the degree to which a given stand attribute is spatially correlated with itself across the study area can then be quantified (i.e. global spatial autocorrelation). In addition, the degree to which the neighbours’ attributes are related to a given stand or not can be computed (i.e. local spatial autocorrelation).

Methods for Defining Neighbours

Two methods for defining neighbourhoods will be discussed: queen weights and rook weights. Both methods have their strengths and weaknesses. By implementing both methods in this analysis, the qualities and drawbacks of these methods will be clear.

Queen Weights

This method defines a given area’s neighbours (nearest areas) by distributing a weight of one to the areas that are in the path of which a queen might take on a chess board. Unlike the image below, the method used in this report will actually set each one of the neighbours to weigh 0.25 since there are four neighbours.

Source: https://bit.ly/topic8notes 
#Defining neighbourhood using queen's case:
#This function builds a neighbours list based on regions with contiguous boundaries (i.e. sharing one or more boundary points)
vri.nb <- poly2nb(vriClean)
#Define a network grid
vri.net <- nb2lines(vri.nb, coords=coordinates(vriClean))
#Apply a coordinate system
crs(vri.net) <- crs(vriClean)

In order to get a visual representation of the neighbours that each stand in the GVWSA is connected to, a map can be created.

#This map will show the neighbours that each polygon is connected to
#It utilizes the centroid of the polygons as the nodes.
#The lines show the connections to other stands.
tm_shape(vriClean) + tm_borders(col='lightgrey') + 
  tm_shape(vri.net) + tm_lines(col='red')

Rook Weights

This method defines a given area’s neighbours (nearest areas) by distributing a weight of one to the areas that are in the path of which a Rook might take on a chess board. Unlike the image below, the method used in this report will actually set each one of the neighbours to weigh 0.125 since there are eight neighbours.

Source: https://bit.ly/topic8notes 
#Defining neighbourhood using rook's case
vri.nb2 <- poly2nb(vriClean, queen = FALSE)
#Define a network grid
vri.net2 <- nb2lines(vri.nb2, coords=coordinates(vriClean))
#Apply a coordinate system
crs(vri.net2) <- crs(vriClean)

Create a map to see the neighbourhoods using rook’s case. It can be assumed that where there is a yellow line, there is a blue line below it. The blue lines are the connections that were included in the queen’s case and not the rook’s case.

tm_shape(vriClean) + tm_borders(col='lightgrey') + 
  tm_shape(vri.net) + tm_lines(col='blue', lwd = 2) +
  tm_shape(vri.net2) + tm_lines(col='yellow', lwd = 2)

Creating a Weights Matrix

The weights for the queen and rook matrices are created using the formulas below (Smith, n.d.). Regarding queen’s case, polygons (j) that share a boundary (bnd) with a given area (i) are assigned a weight of 1. A stronger condition is applied to rook’s case to remove the neighbours that only share a single boundary point (e.g. a shared corner point). In doing so, the rook’s weight matrix applies a condition defining a neighbour as a polygon that shares a positive portion of their boundary with the given area. In other words, neighbours must have a shared boundary of length greater than zero.

Queen weights matrix formula. Source: https://bit.ly/weightsUPenn
Rook weights matrix formula. Source: https://bit.ly/weightsUPenn 
#Create a matrix using queen’s case
#If there are two neighbours, the weight of 1 will be distributed over those two - resulting in a weight of 0.5 for each neighbour
#Non-neighbours will receive a weight of zero
#By setting zero.policy = TRUE, the number zero will be assigned to the lagged value of zones without #neighbours
#Setting zero.policy = FALSE will result in an error for any empty neighbour sets
#The style describes how the weights list with values will be added
#W is row standardized (sums over all links to n)
vri.lw1 <- nb2listw(vri.nb, zero.policy = TRUE, style = "W")
#The following prints a summary of the weights and links used for the weights matrix. 
print.listw(vri.lw1, zero.policy = TRUE)
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 6162 
## Number of nonzero links: 37350 
## Percentage nonzero weights: 0.0983665 
## Average number of links: 6.061344 
## 29 regions with no links:
## 5958 5962 5966 5975 5985 5986 5991 6006 6013 6020 6036 6052 6072 6075 6084 6090 6091 6093 6102 6106 6113 6131 6134 6135 6136 6139 6141 6142 6156
## 
## Weights style: W 
## Weights constants summary:
##      n       nn   S0       S1       S2
## W 6133 37613689 6133 2471.906 28801.38

Repeat for rook’s case.

#Create a matrix using rook’s case
vri.lw2 <- nb2listw(vri.nb2, zero.policy = TRUE, style = "W")
print.listw(vri.lw2, zero.policy = TRUE)
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 6162 
## Number of nonzero links: 32576 
## Percentage nonzero weights: 0.0857935 
## Average number of links: 5.286595 
## 30 regions with no links:
## 5958 5962 5966 5975 5985 5986 5991 6006 6013 6020 6036 6052 6072 6075 6084 6090 6091 6092 6093 6102 6106 6113 6131 6134 6135 6136 6139 6141 6142 6156
## 
## Weights style: W 
## Weights constants summary:
##      n       nn   S0       S1       S2
## W 6132 37601424 6132 2798.685 29199.94

As can be seen by the matrix and links summaries, queen weighting resulted in 37350 nonzero links (i.e. neighbours) while rook weighting resulted in 32576 nonzero links. This is to be expected because the weighting method for queen’s case allows for more flexibility when defining neighbours. The queen weights method yielded approximately 6 average links per polygon while the rook weights method yielded approximately 5 links per polygon. This is unexpected since queen’s weighting is theoretically supposed to yield twice the amount of links (8) as rook’s weighting (4). Furthermore, both weighting methods resulted in similar amounts of regions with no neighbours (29: queen, 30: rook). For consistency, queen weights will be used for the remainder of the project.

Using Global Moran’s I Statistic

The most common method for testing spatial autocorrelation is by using the Global Moran’s I test. This method compares how similar each polygon is with its neighbours and averages all these comparisons. The result is one value that will yield the overall impression of how a given variable is correlated with itself across the study area.

How Does Global Moran’s I Test Work?

In order to quantify the similarity between a given polygon (i) and its neighbours (j), the variation from the mean is calculated for both i and j. If the values of both i and j are either greater or less than the mean, positive spatial autocorrelation exists. When the values of i are above the mean and the values of j are below the mean (or visa versa), negative spatial autocorrelation exists.

Given the formula for the Moran’s I statistic (see below), the weights (wi,j) created in the last step can be applied to distinguish the neighbours (xj) of a given polygon (xi) from the rest of the study area. Polygons given a weight of zero will effectively be removed from the calculation, resulting in only the polygons required for Moran’s I to be calculated.

Global Moran’s I formula. Source: https://bit.ly/topic8notes

Calculating Global Moran’s I

In order to test whether the results produced with the Global Moran’s I statistic are statistically significant, the following needs to be calculated: Global Moran’s I statistic, expected Moran’s I for a random distribution, and the variance of values in the dataset.

#Compute Moran's I Test for stand height
#Uses the data from weights matrix, polygon values, and zero.policy settings
#Remember, setting zero.policy = True will bypass the calculation of the polygons with no neighbours
#Setting zero.policy = False will result in a crashed program if there are any polygons that have no neighbours in the data (which there are in this case)
mi_HT <- moran.test(vriClean$Stand_HT, vri.lw1, zero.policy = TRUE)
mi_HT
## 
##  Moran I test under randomisation
## 
## data:  vriClean$Stand_HT  
## weights: vri.lw1  n reduced by no-neighbour observations
##   
## 
## Moran I statistic standard deviate = 34.953, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      2.830355e-01     -1.630789e-04      6.564808e-05
#Repeat for stand age
mi_AGE <- moran.test(vriClean$Stand_Age, vri.lw2, zero.policy = TRUE)
mi_AGE
## 
##  Moran I test under randomisation
## 
## data:  vriClean$Stand_Age  
## weights: vri.lw2  n reduced by no-neighbour observations
##   
## 
## Moran I statistic standard deviate = 44.569, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      3.840214e-01     -1.631055e-04      7.430493e-05

The following R code will calculate the range of Moran’s I values that represents a random distribution. In other words, a Moran’s I value in between this range suggests that the variable is experiencing random spatial autocorrelation. Moreover, values above the maximum threshold experience positive spatial autocorrelation and values below the minimum threshold experience negative spatial autocorrelation.

If this function were to be called, the ranges for negative (dispersed distribution) and positive (clustered distribution) spatial autocorrelation would be -2.022889 and 2.238584, respectively.

#WARNING: This function will take approximately 20 minutes to run (or longer depending on your
#machine)
#This calculates the expected range of the spatial autocorrelation given connectivity of the polygons #(derived from the weights matrix)
moran.range <- function(lw) {
  wmat <- listw2mat(lw)
  return(range(eigen((wmat + t(wmat))/2)$values))
}
moran.range(vri.lw1)

Upon completion, a Z-test can be computed to determine whether the Global Moran’s I value is significantly different than the expected range of values for a random distribution.

#This extracts the values needed to calculate the z-score
#Global Moran's I
mI_HT <- mi_HT$estimate[[1]]
#Expected Moran's I for a random distribution
eI_HT <- mi_HT$estimate[[2]]
#Variance of values in the dataset
var_HT <- mi_HT$estimate[[3]]
#Calculate Z-score
z_HT <- (mI_HT-eI_HT)/var_HT**0.5

#Repeat for stand age  
mI_AGE <- mi_AGE$estimate[[1]]
eI_AGE <- mi_AGE$estimate[[2]]
var_AGE <- mi_AGE$estimate[[3]]
#Calculate Z-score
z_AGE <- (mI_AGE-eI_AGE)/var_AGE**0.5

With the null hypothesis stating that the observed Global Moran’s I statistic is no different than the expected random statistic and the alternate hypothesis stating the converse, a confidence statement can be made. For this analysis a significance level of 0.05 will be used; therefore, if the observed Moran’s I values are above +1.96 or below -1.96, the null hypothesis can be rejected. Given the z-scores of stand height (35.22) and stand age (44.57), we may say that we are 95% confident that both variables experience spatial autocorrelation. In particular, since the observed values are above the expected range for a random Moran’s I (2.238584), it can be said that significant positive spatial autocorrelation exists in both stand height and stand age.

Using Local Moran’s I Statistic

Like Global Moran’s I, the Local Moran’s I statistic can be used to quantify spatial patterns occurring across a given area. However, instead of producing one value that encompasses the entire study area, the Local Moran’s I method produces a statistic for each polygon in the study area. The advantage of this localized method is that it shows exactly which polygons are similar or different to the objects in their neighbourhood. In other words, it identifies where spatial autocorrelation occurs as well as the level of significance that is experienced. Using Local Moran’s I Statistic Like Global Moran’s I, the Local Moran’s I statistic can be used to quantify spatial patterns occurring across a given area. However, instead of producing one value that encompasses the entire study area, the Local Moran’s I method produces a statistic for each polygon in the study area. The advantage of this localized method is that it shows exactly which polygons are similar or different to the objects in their neighbourhood. In other words, it identifies where spatial autocorrelation occurs as well as the level of significance that is experienced.

How Does Local Moran’s I Test Work?

The Local Moran’s I statistic is a measure of the difference between a given observation (i) and the mean multiplied by the sum of differences of its neighbours (j) and the mean (see below). The difference between i and the mean is also divided by the standard deviation of i (i.e. a measure of how much the values at neighbourhood locations vary from the mean). When the attributes of i and j are either above (i.e. hot spot) or below (i.e. cold spot) the mean, it means positive spatial autocorrelation is present at location i (i.e. spatial clustering). Spatial outliers occur where the attribute at i is above than the mean and the attributes of j are higher than the mean, visa versa.

Local Moran’s I formula. Source: https://bit.ly/topic8notes

Calculating Local Moran’s I

#This will calculate the observed Local Moran's I, expected Local Moran’s I, variance, z-score, and p
#value for each polygon
#Stand Height
lisa.test_HT <- localmoran(vriClean$Stand_HT, vri.lw1)
#Type ‘lisa.test_HT’ in the console to see the output
#Stand Age
lisa.test_AGE <- localmoran(vriClean$Stand_Age, vri.lw1)
#Type ‘lisa.test_Age’ in the console to see the output

#Extracting the results from lisa.test (Local Moran’s I test) and appending it to the spatial dataset
#Stand Height
vriClean$Ii_HT <- lisa.test_HT[,1]
vriClean$E.Ii_HT <- lisa.test_HT[,2]
vriClean$Var.Ii_HT <- lisa.test_HT[,3]
vriClean$Z.Ii_HT <- lisa.test_HT[,4]
vriClean$P_HT <- lisa.test_HT[,5]

#Repeat for Stand Age
vriClean$Ii_AGE <- lisa.test_AGE[,1]
vriClean$E.Ii_AGE <- lisa.test_AGE[,2]
vriClean$Var.Ii_AGE <- lisa.test_AGE[,3]
vriClean$Z.Ii_AGE <- lisa.test_AGE[,4]
vriClean$P_AGE <- lisa.test_AGE[,5]

Once the calculations for Local Moran’s I are complete, a map of the z-scores from each polygon can be plotted to display the polygons that are significantly different than their neighbours. Given a 95% confidence level, polygons with z-scores above 1.96 experience significant positive spatial autocorrelation while polygons with z-scores below -1.96 experience significant negative spatial autocorrelation.

#Map individual results
#Polygons with z-scores > 1.96 experience significant positive spatial autocorrelation (SAC)
#Polygons with z-scores < -1.96 experience significant negative spatial autocorrelation (SAC)
#Polygons with z-scores in between -1.96 and 1.96 experience a SAC no different than a polygon exhibiting random SAC
#Use > tmaptools::palette_explorer() to explore colour palettes
#Map Stand Height
map_LISA_HT <- tm_shape(vriClean) + 
  tm_polygons(col = "Z.Ii_HT", 
              title = "Local Moran's I - Stand Height", 
              style = "fixed",
              breaks = c(-Inf, -1.96, 1.96, Inf),
              labels = c("Negative SAC", "Random SAC", "Positive SAC"),
              palette = "RdBu", n = 3,
              midpoint = NA,
              border.alpha = 0.3) 
map_LISA_HT

As can be seen, clusters of stands in the southeastern and central parts of the GVWSA experienced positive spatial autocorrelation. In other words, these stands, along with their neighbours, were either significantly taller or shorter than the mean stand height of the study area. There were few stands experiencing negative spatial autocorrelation; those that did are primarily located in the southeast portion of the GVWSA.

#Map Stand Age
map_LISA_AGE <- tm_shape(vriClean) + 
  tm_polygons(col = "Z.Ii_AGE", 
              title = "Local Moran's I - Stand Age", 
              style = "fixed",
              breaks = c(-Inf, -1.96, 1.96, Inf),
              labels = c("Negative SAC", "Random SAC", "Positive SAC"),
              palette = "RdBu", n = 3,
              midpoint = NA,
              border.alpha = 0.3) 
map_LISA_AGE

As can be seen, clusters of stands in the central part of the GVWSA experienced positive spatial autocorrelation. In other words, these stands, along with their neighbours, were either significantly older or younger than the mean stand age of the study area. There were few stands experiencing negative spatial autocorrelation; those that did are primarily located in the southeast portion of the GVWSA.

Plot Local and Global Spatial Autocorrelation

To summarize both the local and global spatial autocorrelation observed for both stand height and stand age in the GVWSA, two scatter plots can be made. Each point in the plot represents a polygon’s relationship with its neighbours. The points in quadrants 1 and 3 represent the stands that experience positive spatial autocorrelation. This is because both the values of i and j (think back to the Local Moran’s I formula) are either above or below the mean. Points in quadrants 2 and 4 represents the stands that experience negative spatial autocorrelation. These are the stands in which their values are above the mean while their neighbours are below the mean, visa versa. Points represented with a diamond shape are values that are deemed to be significant at a specific level of confidence (set by the p-values calculated using the lisa.test); the numbers next to these points represent their ID number. The trend lines in each plot indicates that positive spatial autocorrelation is present throughout the entire dataset. In other words, the trend line is representing the global spatial autocorrelation. The grey circles are simply the polygons with no defined neighbours.

#Plot stand height
local_moran_plot_HT <- moran.plot(vriClean$Stand_HT, vri.lw1, zero.policy=TRUE, spChk=NULL, labels=NULL, xlab="Stand Height", 
                               ylab="Nieghbour's Height", quiet=NULL)

#Plot stand age
local_moran_plot_AGE <- moran.plot(vriClean$Stand_Age, vri.lw1, zero.policy=TRUE, spChk=NULL, labels=NULL, xlab="Stand Age",
                                   ylab="Nieghbour's Age", quiet=NULL)

Summary

In this tutorial, you have learned a few techniques for analyzing spatial patterns in a VRI dataset. One of them is Global Moran’s I, where I can be used to quantify whether spatial autocorrelation exists in the dataset or not. Another is Local Moran’s I, where z-scores can be used to visualize where the clusters of spatially autocorrelated stands are. Along with mapping and plotting the data, the Global and Local Moran’s I statistics provide valuable information regarding the spatial autocorrelation of values in your dataset.

Regarding the results produced from this project, several takeaways can be made regarding the presence of spatial autocorrelation in stand height and stand age within the GVWSA. The Province of BC may be interested in logging the areas where stands are old and tall to maximize profits. The local spatial autocorrelation maps provide useful information for supporting the decisions relating to this. However, these maps do not show whether the stands experiencing spatial autocorrelation are below or above the mean.

To build on the methods discussed in this project, one may choose to perform a hot spot analysis and outlier analysis to determine the directionality of the stands experiencing positive spatial autocorrelation. In other words, are i and j experiencing values below or above the mean? Moreover, the areas experiencing negative spatial autocorrelation can be further inspected to see which of the stands represent significant outliers in the dataset. By implementing these analyses, the Province of BC will be able to gain an in-depth analysis to the attributes within the VRI dataset they are interested in. Furthermore, further analysis of other variables may point to reasons for why the observed spatial patterns exits in stand height and stand age in the GVWSA.

References

Haining, R. P. (2001). Spatial Autocorrelation. In International Encyclopedia of the Social & Behavioral Sciences (pp. 14763–14768). https://doi.org/10.1016/b0-08-043076-7/02511-0 Hassegawa, M., Havreljuk, F., Ouimet, R., Auty, D., Pothier, D., & Achim, A. (2015). Large-scale variations in lumber value recovery of yellow birch and sugar maple in Quebec, Canada. PLoS ONE, 10(8), e0136674. https://doi.org/10.1371/journal.pone.0136674 Province of British Columbia. (n.d.). Vegetation Resources Inventory. Retrieved October 17, 2020, from https://www2.gov.bc.ca/gov/content/industry/forestry/managing-our-forest-resources/forest-inventory Smith, T. E. (n.d.). Spatial Weight Matrices. Tobler, W. R. (1970). A Computer Movie Simulating Urban Growth in the Detroit Region. Economic Geography, 46, 234. https://doi.org/10.2307/143141