Review
Today’s lab is relatively short, before we begin we will do some interactive review of the lab material we have covered to date. Open the following link on one of the lab computers to follow along https://poll.gis.unbc.ca/p/31359493
Moran’s I
In today’s lab we will be looking at calculating Moran’s I, and Local Moran’s I using R; while in R we will also take a look at how we can produce reproducible reports using R-Markdown documents. For those unfamiliar Markdown is a tool used by various programs that allow for documents to be written in plain text, and rendered with various style options that were defined in said text.
In RStudio you can create your first Markdown document by going into File, New and choosing R Notebook.
![](https://gis.unbc.ca/wp-content/uploads/2022/11/image.png)
The main parts of the new new document are
- The heading starts and ends with —
- Markdown (Think of this as plain text
- R Code surrounded by “`
{r}
“`
It is generally considered proper form to include your name and date in the file header.
--- title: "Morans-I" author: "Matt McLean" date: "07/11/2022" output: html_notebook: default pdf_document: default ---
With that in place you can preview your document with the preview button
![](https://gis.unbc.ca/wp-content/uploads/2022/11/image-1.png)
Finally the Table of Contents is another very helpful tool allowing the document to get some structure, outside of r code this is done with the # symbol, for example,
# Heading 1 ## Heading 2 # H1 Two ## H2 Two ### Heading 3
Will produce the following headings, these will both appear in the final document, and can also be used for navigating within the document inside of R Studio, by clicking on the desired heading.
![](https://gis.unbc.ca/wp-content/uploads/2022/11/image-2.png)
Preparing for Moran’s-I
For today’s lab we will be calculating Moran’s-I on some trials conducted by the Ministry of Forests, and will be using a package not normally available in R, this will be installed in the R command line with the following command.
install.packages('bcmapsdata', repos='https://bcgov.github.io/drat/')
After this is installed you can import the libraries that we will be using today. Note there is an argument include=FALSE, this will prevent all these import statements from appearing in our final report.
```{r, include=FALSE} library(bcdata) library(bcmaps) library(bcmapsdata) library(tidyverse) library(sf) library(ggplot2) library(ape) library(units) library(lubridate) library(spdep) library(sp) library(tmap) library(rgdal) ```
Download the trial data from GitHub https://github.com/FLNRO-Smithers-Research/OffSite_Trial_Reporting/blob/master/trial_data.Rdata
After you have downloaded into your working directory you can open with the following. It is relatively large code, as it expands the data from a single data-frame to multiple data-frames as well as calculate the extent of data (ie how many years does the data cover).
```{r} load('trial_data.Rdata') seedlots = all_data$seedlots trees = all_data$trees climate = all_data$climate plots = all_data$plots meta = all_data$meta sites = distinct(select(merge(meta, plots, by="FID"), Location, FID, Lat, Lon), Location, FID, .keep_all = TRUE) #Location based on a plot not center of entire site climate_sites = distinct(climate, ClimateStationID) #List of years for which tree data is present, don't try to plot dates outside this range years = year(as.character(min(as.Date(as.character(trees$Date))))):year(as.character(max(as.Date(as.character(trees$Date))))) #List of years for climate data can differ from plot data. climate_years = year(as.character(min(as.Date(as.character(climate$Date))))):year(as.character(max(as.Date(as.character(climate$Date))))) #This function was used by many of the plots, moved here to not repeat tree_summary=list() for(year in years){ tree_summary[[year]] <- trees %>% filter(year(as.character(Date)) == year) %>% # need to convert seedlot to factor mutate(Seedlot=factor(Seedlot)) %>% # Relevel condition factors mutate(Condition=fct_relevel(Condition,"Missing","Dead", "Moribund","Poor","Fair","Good")) %>% # reorder levels # Create plot summaries group_by(FID,Plot,Species,Seedlot,Condition) %>% summarise(n=n()) %>% # summarize number of trees at FID/Plot/Species/Condition ungroup() } ```
To prepare for calculating the statistics the data must be opened, this sample data set contains 3 attributes that we can use to check for spatial auto correlation with Species, Seedlot (which tree farm are the seedlings from), and condition. You may note that there are also height and diameter columns, while potentially interested this data has been scrubbed from the open source side of the project.
attribute_list <- c('Species', 'Seedlot', 'Condition')
Next the data needs to be sorted by trial site, (Identified by FID)
# Get list of unique trial sites fids = as.data.frame(trees %>% distinct(FID) ) # Fore each trial site for(id in 1:nrow(fids)) { # Convert iterator back to factor feature <- as.vector(fids$FID[[id]]) # Select trees at trial site and that have a location tree_list = as.data.frame(trees %>% filter(FID == feature & (Lon != 0 | Lat != 0))) if (nrow(tree_list) > 0){ # Print trial site name cat(" \n####", as.character(feature), " \n") } else{ # No trees were found in filter cat(" \n", as.character(feature)", tree positions unknown") } }
To compute Moran’s I we need to convert the Trees to SF features and project into UTM
tree_sf <- st_as_sf(tree_list, coords = c('Lat', 'Lon')) st_crs(tree_sf) <- 4326 tree_utm <- st_transform(tree_sf, 32610)
Once the data is sf the next stage is to calculate the inverse distance of all points
# Compute distance between points knn <- knn2nb(knearneigh(tree_utm)) # Determine maximum Neighbor Distances all.linked<- max(unlist(nbdists(knn, tree_utm))) # Make Nerest Neighbor List tree_dist <- dnearneigh(tree_utm, 0, all.linked) # Style W is row-wise normalization lw <- nb2listw(tree_dist, style="W", zero.policy = T)
When working with categorical data it will be converted to factor in R, for proper calculation we would like this numbers to have a logical pattern, as such we will re-level the condition attribute
tree_list %>% mutate(Condition=fct_relevel(Condition,"Missing","Dead", "Moribund","Poor","Fair","Good"))
Finally we want to plot the attribute and calculate Moran’s I
for (att in attribute_list){ p <- ggplot(tree_utm) + geom_sf(aes(col = .data[[att]])) + coord_sf(datum = st_crs(tree_utm)) + theme(axis.text.x = element_text(angle = 45, vjust = 0.5)) + ggtitle(paste(feature, "-", att)) print(p) # Test Morans I test <- moran.test(as.numeric(tree_utm[[att]]), lw) print(test) # Monty Carlo method MI <- moran.mc(as.numeric(tree_utm[[att]]), lw, nsim = 599, zero.policy = T) print(MI) }
This is put all together as the following function
attribute_list <- c('Species', 'Seedlot', 'Condition') fids = as.data.frame(trees %>% distinct(FID) ) for(id in 1:nrow(fids)) { feature <- as.vector(fids$FID[[id]]) tree_list = as.data.frame(trees %>% filter(FID == feature & (Lon != 0 | Lat != 0))) if (nrow(tree_list) > 0){ cat(" \n####", as.character(feature), " \n") tree_sf <- st_as_sf(tree_list, coords = c('Lat', 'Lon')) st_crs(tree_sf) <- 4326 tree_utm <- st_transform(tree_sf, 32610) # Compute distance between points knn <- knn2nb(knearneigh(tree_utm)) # Determine maximum Neighbor Distances all.linked<- max(unlist(nbdists(knn, tree_utm))) # Make Nerest Neighbor List tree_dist <- dnearneigh(tree_utm, 0, all.linked) # Style W is row-wise normalization lw <- nb2listw(tree_dist, style="W", zero.policy = T) tree_list %>% mutate(Condition=fct_relevel(Condition,"Missing","Dead", "Moribund","Poor","Fair","Good")) for (att in attribute_list){ p <- ggplot(tree_utm) + geom_sf(aes(col = .data[[att]])) + coord_sf(datum = st_crs(tree_utm)) + theme(axis.text.x = element_text(angle = 45, vjust = 0.5)) + ggtitle(paste(feature, "-", att)) print(p) # Test Morans I test <- moran.test(as.numeric(tree_utm[[att]]), lw) print(test) # Monty Carlo method MI <- moran.mc(as.numeric(tree_utm[[att]]), lw, nsim = 599, zero.policy = T) print(MI) } } else{ cat(" \nTree positions unknown") } cat(" \n") }
Assignment
Create a Report on the Spatial Auto-correlation between Annual Precipitation and location of weather station using the stations data found in L:\GEOG413\lab08
You should be able to
- Import the data, and confirm the CRS / Translate if Needed
- Plot the data on a basemap (note that as you are making a PDF report interactive leaflet will not be an option)
- Report on if there is a correlation, and if it is significant.