edwebinar_mar19_ornldaac_supplemental.Rmd

---
title: "Introduction to Geospatial Analysis in R"
subtitle: "Supplemental Tutorial: Neighborhood analysis"
author: 'Presented by the ORNL DAAC  https://daac.ornl.gov'
date: "January 9, 2024"
output:
  html_document:
    keep_md: yes
    number_sections: yes
    toc: yes
  html_notebook:
    number_sections: yes
    toc: yes
editor_options:
  chunk_output_type: inline
---
***

<!--------------------------- LOAD THIS FILE INTO RSTUDIO --------------------------->

This supplemental tutorial will demonstrate how perform a **neighborhood analysis** to identify locations disturbed by insects are adjacent to fire disturbance. The tutorial will use two raster maps created during the main tutorial.

# Install and Load Packages  

As with the main portion of the tutorial, the *terra* and *sf* must be installed along with their dependencies.

```{r install_packages, message=FALSE, warning=FALSE, eval=FALSE}
install.packages("terra", dependencies = TRUE)  
install.packages("sf", dependencies = TRUE)  
```

Load the *terra* and *sf* packages.

```{r load_packages, message=FALSE, warning=FALSE}
library(terra)       # provides most functions
library(sf)           # provides function for saving files
```

# Load Data  

> *Functions featured in this section:*  
> **sf::st_read()**  
> reads a shapefile and loads it as a simple feature object  

We must load two objects, *threeStates* and *prjFireInsect*, that were created from manipulations described in the main part of the tutorial and saved in the "./data" folder. To load a shapefile, use the `st_read()` function from the *sf* package. We will create a third object, *focalFireInsect*, later in this supplemental tutorial.

```{r read_maps, message=FALSE, warning=FALSE, results="hide"}
threeStates <- st_read(dsn="./data/threeStates.shp", quiet=TRUE)
prjFireInsect <- rast("./data/prjFireInsect.tif")
```

# Perform a Focal Analysis: a simple demonstration  

> *Functions featured in this section:*  
> **terra::focal()**  
> calculates values from the neighborhood of cells ("moving window") around a target cell using a matrix  

A neighborhood analysis computes the value of a single focal cell using information from nearby cells.  Imagine a donut (i.e., an annulus) where the focal cell is in the donut hole and values from the ring of surrounding cells  are used to derive a value for the focal cell. 

Users can specify neighborhoods that best fit the spatial processes under analysis. The neighborhood may or may not include the focal cell, and the width of the band of surrounding cells can be specified to alternatively enlarge or shrink the neighborhood. Moreover, the neighborhood does not have to take the form of ring. It can be rectangular and/or expanded on one side but not the others.  For example, imagine estimating air pollution values for a cells in a landscape with westerly winds. The analysis could use a triangular neighborhood that fans out to the west of the focal cell to use values from upwind cells while ignoring downwind values to the east. All cells within the neighborhood can be treated equally. Alternatively, the relative influence of cells can be weighted by their distance from the focal cell.

The `focal()` function is used in the *terra* package for neighborhood analysis.  We begin with a demonstration of how the `focal()` behaves using a small SpatRaster object.

First, we create a matrix object similar to the raster grid of a SpatRaster object. Like *prjFireInsect*, the matrix object is made up of cell values (1's code for insect damage and 2's code for fire damage) and NA's (non-values). Unlike *prjFireInsect*, it is much smaller with only eight rows and eight columns.

```{r demo_1, message=FALSE, warning=FALSE}
demo_mat <- matrix(data = c(1, NA, NA, NA, NA, NA, NA, 1,
                            1, NA, NA, NA, 1, NA, 2, 1,
                            NA, 2, 1, NA, 1, NA, 1, NA,
                            NA, 1, NA, NA, NA, 2, NA, 1,
                            2, 2, NA, 1, 2, NA, NA, 1,
                            NA, NA, NA, 1, NA, 2, NA, NA,
                            1, NA, 1, 1, 2, NA, NA, NA,
                            1, NA, NA, NA, NA, NA, NA, 1), 
                   ncol = 8, byrow = TRUE)
print(demo_mat)
```

Printing the matrix object shows how it resembles an 8 by 8 grid.

Next, convert the matrix object to a raster object using *terra*'s `rast()` function. Calling `print()` shows the properties of a SpatRaster object with no CRS.

```{r demo_2, message=FALSE, warning=FALSE}
demo_rast <- rast(demo_mat)
print(demo_rast)
```

The `focal()` function is useful for raster analyses because it allows the user to compute the value of a focal (target) cell based upon the values of nearby cells. This utility has two components: (a) a definition of the window that identifies which cells to include in the neighborhood and (b) the function that evaluates the values of the nearby cells.

First , we define the "neighborhood", that is, the locations and number of surrounding cells considered as a "neighbor" of the focal cell. We will use a rectangular neighborhood of eight cells so that every cell adjacent to the focal cell is considered a neighbor. This neighborhood is represented by a 3 x 3 matrix; it includes the focal cell in the center and 8 adjacent cells. 

```{r demo_3, message=FALSE, warning=FALSE}
demo_neigh <- matrix(rep(1,9), ncol = 3)   # a 3 x 3 matrix of 1's
print(demo_neigh)
```

We use 1's to represent all the neighboring cells because we want all cells weighted equally, and we also include the focal cell in our analysis. The *demo_neigh* matrix is passed to the `focal()` command below following "w = " argument below; the "w" is a reference to "moving window", a synonym for neighborhood.  

The options "expand = TRUE" and "fillValue = 0" are used because we want to consider **all** cells in *demo_rast*, even the "edge" cells that do not have a complete neighborhood of cells. In the case of edges, the "missing" neighbors will be considered zeros. 

The function used to evaluate cells in the neighborhood is set by "fun = 'max'", which tells R to find the maximum value for all cells in the window. The "na.rm = TRUE" specifies that NA cells in the neighborhood will be skipped. The "na.policy='omit'" tells `focal()` to not to compute values for focal cells that are NA.

```{r demo_4, message=FALSE, warning=FALSE}
demo_foc <- focal(x=demo_rast, w = demo_neigh, fun = 'max', na.rm=TRUE, 
                  na.policy='omit', expand = TRUE, fillValue = 0)
print(matrix(demo_foc[], ncol=8))
```

Compare the *demo_foc* output with *demo_mat*. If a focal cell was 1, but adjacent to a 2, the focal cell became a 2. A focal cell that was a 2 remains a 2 because 2 is the maximum value of neighboring cells in *demo_rast*.

Suppose we want to identify the cells in *demo_rast* that are 1's having at least one 2 in its neighborhood. We can draw on the information held in both *demo_rast* and *demo_foc*.  First, combine the two rasters using raster algebra (i.e., adding the two SpatRaster objects).

```{r demo_5, message=FALSE, warning=FALSE}
demo_temp <- demo_rast + demo_foc   # adds values in the two rasters, cell by cell
print(matrix(demo_temp[], ncol=8))
```

Every cell that was a 1 in *demo_rast* but changed to a 2 in *demo_foc* is now a 3 (e.g., *demo_rast* = 1 and *demo_foc* = 2; therefore, *demo_rast* + *demo_foc* = 1 + 2 = 3). A cell that was a 2 in *demo_rast* and stayed a 2 in *demo_foc* is now a 4. So, the cells with a 3 in *demo_tmp* are the ones in the *demo_rast* cells with value 1 and at least one 2 in their neighborhood.

We will use *terra*'s `app()`function (apply) to reclassify the values of *demo_temp* to get the final values we'd like to see represented. In the code below, we tell R that if a *demo_temp* value is 3, change it to a 1, and save  every other value as NA.

```{r demo_6, message=FALSE, warning=FALSE}
demo_tar <- app(demo_temp, 
                 fun = function(x) {ifelse( x == 3, 1, NA) } )
print(matrix(demo_tar[], ncol=8))
```

This raster is our final product for the demonstration. The cells with value of 1 indicate 1-cells in *demo_rast* that were adjacent to a 2. We will apply this same logic to the "real" data.

# Perform a Focal Analysis: identify locations adjacent fire and insect disturbances 

Using *prjFireInsect*, we will determine which forest locations were damaged by insects and were adjacent to forest locations damaged by fire. In other words, we will locate cells with the value 1 that are next to cells with the value 2. As explained in the main part of the tutorial, there are no cells that had both *fire* and *insect* disturbance, which is why we will use `focal()` here.

The code below is the same as demonstrated above, but it uses a much larger SpatRaster object, *prjFireInsect*.  
The code runs the `focal()` analysis to create *focalFireInsect*.  Then, that raster is reclassified to create *tarFireInsect*. 

```{r focal_FI, message=FALSE, warning=FALSE}
# this could take several minutes to run; use saved file if available.
if(file.exists("./data/focalFireInsect.Rds")) { 
  focalFireInsect <- readRDS("./data/focalFireInsect.rds") 
  }else{ 
    fireinsect_neigh <- matrix(rep(1,9), ncol = 3)
    focalFireInsect <- focal(prjFireInsect, w = fireinsect_neigh, 
                             fun = 'max', na.rm=TRUE, na.policy='omit', 
                             expand = TRUE, fillValue = 0)
      }
print(focalFireInsect)
```

```{r tarFI_map, message=FALSE, warning=FALSE, results="hide"}
# this could take a few minutes to run
temp <- focalFireInsect + prjFireInsect
tarFireInsect <- app(temp, fun = function(x) {ifelse( x == 3, 1, NA) } )
print(tarFireInsect)
cat("\nSummary\n", summary(tarFireInsect[]), "\n")
rm(temp)    # clean up
```

Notice that *tarFireInsect* has only the value 1 and NA's.  The 1's identify the cells that meet the criterion: insect damage adjacent to at least one cell with fire damage.


# Get Cell Coordinates  

> *Functions featured in this section:*  
> **terra::xyFromCell**  
> gets coordinates of the center of raster cells for a SpatRaster object  

In this section, we will get the geographic coordinates of the cells that were identified in the last section. First, we create a vector that stores the index of *tarFireInsect* cells that are 1. That is, we store the matrix  or grid "locations" of the values across the extent of the SpatRaster object, as if counting from one to the total number of cells across the raster grid. The cell index number starts with 1 in upper left corner of raster grid. Index numbers increase from left to right across the row, and then from top to bottom. The cell index at the bottom right corner is equal to `ncell(tarFireInsect)`, the total number of cells in the raster layer.

```{r get_coord_1, message=FALSE, warning=FALSE}
val_tar <- which(tarFireInsect[]==1)                  # returns index of cells with value of 1
cat("Number of 1-cells: ", length(val_tar), "\n\n")   # print number of cells with value of 1
cat("First six index values: ", head(val_tar) )
```

There are 11,487 cells that have value of 1.  We use the function `head()` to view the first six values of *val_tar*. For example, the output tells us that the cell at position 18,468,526 of the "grid" has the value 1

Using the cell "locations" provided by *val_tar* and the `xyFromCell()` function, we can extract the coordinates of those cells. The function outputs a matrix object, *loc_tarFireInsect*, and we name its columns. Because the CRS of this SpatRaster uses geographic coordinates, the x,y coordinates will be in longitude, latitude.

```{r get_coord_2, message=FALSE, warning=FALSE}
loc_tarFireInsect <- xyFromCell(tarFireInsect, val_tar)
colnames(loc_tarFireInsect) <- c("lon","lat")
head(loc_tarFireInsect)
```

The first six rows of *loc_traFireInsect* shows the longitude and latitude of cells that have the value 1. In other words, we now know the exact locations of forested areas throughout Idaho, Montana, and Wyoming that had insect damage **and** were adjacent to a location that had fire damage.

We save these coordinates in a comma-separated values (CSV) format for future use. If we wanted to, we could visit these locations in person! Try looking up one of these locations using a web-based map, like Google Maps.

```{r save_coord, message=FALSE, warning=FALSE, eval=FALSE}
write.csv(loc_tarFireInsect, file = "loc_tarFireInsect.csv", row.names = FALSE)
```


# Plot Fire-Insect locations 

> *Functions featured in this section:*  
> **sf::st_as_sf**  
> gets coordinates of the center of raster cells for a SpatRaster object

Using the coordinate, we can plot these locations on a dot map. The `st_as_sf()` function, as used below, creates a simple feature point object from a data frame holding point coordinates.  The "coords=c(1,2)" and "crs =" arguments specify which columns of the data frame hold the x,y coordinates as well as their CRS. 

```{r plot_FI_loc, message=FALSE, warning=FALSE}
# create simple feature point object
fi_pt <- st_as_sf(as.data.frame(loc_tarFireInsect), 
                  coords=c(1,2), crs = "+proj=longlat +datum=WGS84")
plot(threeStates$geometry, main="Locations with fire-insect damage")
plot(fi_pt, add=T, col="blue", cex=0.7)

```



<!--------------------------------- END OF TUTORIAL --------------------------------->