Intro2R-spatial

.title[
# Intro2R-spatial
]
.subtitle[
## Analysing spatial patterns of the landscape
]
.author[
### <a href="https://xp-song.github.io">Song, Xiao Ping</a>
]
.institute[
### <a href="mailto:xp.song@u.nus.edu" class="email">xp.song@u.nus.edu</a> Course materials: <a href="https://github.com/xp-song/Intro2R-spatial" class="uri">https://github.com/xp-song/Intro2R-spatial</a>
]
.date[
### updated 2020-10-09
]

---

# Intro2R-spatial

1. Navigate to course webpage for instructions and background information  
*https://github.com/xp-song/Intro2R-spatial*

2. [Download](https://github.com/xp-song/Intro2R-spatial/archive/master.zip) workshop materials  
*(green button on webpage)*

---

# Outline  
 
__[Why analyse spatial patterns?](#why)__

[Landscape ecology: Conceptual models](#concepts)

[Recap: R Environment](#renviron)

[Land cover classification](#classify)

[Landscape metrics](#lsm)

[Useful resources](#resources)

???
- First off, I will explain what I mean by spatial patterns
- Recap on the basics of R
- Get into the analysis

---
class: left
background-image: url("images/sg_urech2017.jpg")
background-position: center
background-size: cover

# Spatial pattern of a landscape

### A function of the *__composition__* and *configuration* of the physical landscape

???
- What is this landscape composed (consist) of?

---
class: left
background-image: url("images/sg_urech2017.jpg")
background-position: center
background-size: cover

# Spatial pattern of a landscape

### **Spatial composition**
- Variety - What types/classes?
- Abundance - How much?

???
- E.g. Basic type of classification: as green, blue, grey
- Abundance: Absolute (i.e. total area) or relative (%)

---

class: left
background-image: url("images/sg_urech2017.jpg")
background-position: center
background-size: cover

# Spatial pattern of a landscape

### A function of the *composition* and *__configuration__* of the physical landscape

???
- Quantifying this is more complex
- What are some examples of configuration?

---
class: left
background-image: url("images/sg_urech2017.jpg")
background-position: center
background-size: cover

# Spatial pattern of a landscape

### **Spatial configuration**
- Position
- Orientation
- Shape/arrangement

???
What if...
1. All greenery was on the left... buildings/roads on the right?
2. Greenery finely interspersed/intertwined with the grey cover?

- Focusing on configuration is where the idea of "connectivity" becomes relevant
  - Connectivity for... Wind, biodiverisity, recreation (e.g. cycling networks)

---

# Why is spatial pattern important?

### The study of landscape ecology

???
- One area of study tt focuses very much on spatial pattern...  Landscape Ecology

- Interested in how spatial pattern affects ecological processes & vice versa (reciprocal)

- __What are some egs of this pattern-process relationship?__
  - Pattern of topography/altitude -> _downstream hydrological processes_
  - Pattern/Arrangement of greenery -> _temperature? distribution of noise/pollutants? visual appeal?_

---

# Why is spatial pattern important?

### The study of landscape ecology

__Pattern-process relationships__  
  
- How do we optimise spatial *patterns* to achieve desired ecological *processes*?  
- How can we incorporate such thinking into policy & design?

???
- Some qns...

- Ultimately, we want to link these relationships to larger goals.
- As designers, your job is to produce sth that your stakeholders will like.
- Key qn: How do we go beyond the looks, and also improve the __performance__ of greenery?

---

# Outline  
 
[Why analyse spatial patterns?](#why)

__[Landscape ecology: Conceptual models](#concepts)__

[Recap: R Environment](#renviron)

[Land cover classification](#classify)

[Landscape metrics](#lsm)

[Useful resources](#resources)

---

class: center
background-image: url("images/sg_urech2017.jpg")
background-position: center
background-size: cover

# Landscape ecology: Conceptual models

__Mosaic model vs. Gradient model__

???
- Previously, we talked about green/grey cover
- Classified the landscape into discrete classes > __Mosaic model__  
  
- But... Canopy above road.... Are we oversimplifying...?
- __Gradient model__ > ... Soft edges between classes

---

# Landscape ecology: Conceptual models

???
- PCM: Patch Corridor Matrix
  - Patches & corridors are places tt animals can live/survive.

- Note: Both images can be in raster format!
  
---

class: center
background-image: url("images/sg_urech2017.jpg")
background-position: center
background-size: cover

# Landscape ecology: Conceptual models

__Do you see the patches & corridors in this landscape?__

???

- Which is the patch/corridor for
  - Humans? A wild boar?
  - An insect - The same? (corridor - canopy)
  - A bird - larger or smaller? (much larger - more mobile!)
  
- PCM varies depending on species
- For today, we simply classify our landscape based on land cover

---

# Outline  
 
[Why analyse spatial patterns?](#why)

[Landscape ecology: Conceptual models](#concepts)

__[Recap: R Environment](#renviron)__

[Land cover classification](#classify)

[Landscape metrics](#lsm)

[Useful resources](#resources)

???
- Time to get our hands dirty
- Before we get into the analysis... Quick recap to make sure we are all on the same page

---

# Course materials

On your computer, navigate to folder downloaded from  
*https://github.com/xp-song/Intro2R-spatial*

_notes.Rmd_ 
_notes.html_ 1 
_/slides_ 1 
_/data_ 
_/clean_data_ 
_Intro2R-spatial.Rproj_

???

- clean_data: Processed data. We will specify in our code to export our results here
- notes.html is the output file of notes.Rmd
- Open the RStudio Project file. This boots up RStudio (demo)

---

# R Studio Client

- __Console:__ Command line input/output  
- __Script editor:__ View/edit files that contain code
- __Environment/History__  
- __Files/Plots/Packages/Help/Viewer__

???
- Console is mainly for testing/real-time interaction (input/output)

- Impt to work primarily within scripts (File > New File > R script)

- R Notebook (.Rmd) - provide more functions (Word vs. Texteditor)
  - (File > New File > R Notebook) - Follow instructions

---

# RStudio Projects

- RStudio Projects help organise your work into separate 'R sessions'.

- Each project has it's own workspace a.k.a. 'working directory' (separate configuration, history, etc.)

???
- Workspace = Table in my office 
- How do I know which is my table?

- The location of the '_.RProj_' file defines the 'working directory'

???
- E.g. tissue paper
- Likewise, .RProj file tells R that this is my workspace 
  - Address: see window header, type getwd()  
  
- Each project will have a diff wd

## 🌟 Best Practice

- Use _relative_ paths in your script, based on _.RProj_ file location   
  
- Keep all project items in the working directory

???

- By using RStudio Projects... don't need to write out initial path
- Demo: list.files("data") vs. absolute path
- Scripts work across different computers

---

# Set up

Use a new R notebook for our analysis 1

Open the R notebook _notes.Rmd_ file for reference

.footnote[[1] Ensure that it is in the working directory (i.e. same folder as the _Intro2R-spatial.Rproj_ file)
]

???
- We already created one earlier, can save it as a new file.

---

# Set up

### Install packages

- [tidyverse](https://www.tidyverse.org/packages/): A collection of packages commonly used for data analyses  
- [raster](https://www.rdocumentation.org/packages/raster/versions/3.0-7): Analyse raster data
- [sf](https://r-spatial.github.io/sf/articles/sf2.html): Analyse spatial data as "simple features" (i.e. tabular 2D data)   
- [landscapetools](https://cran.r-project.org/web/packages/landscapetools/vignettes/overview.html): Collection of functions for landscape analysis  
- [landscapemetrics](https://r-spatialecology.github.io/landscapemetrics/index.html): Analyse spatial patterns of discrete landscape classes

???

- Packages = apps (each w diff set of functions)
- Packages related to the tidyverse are designed to follow a specific workflow
  - E.g. how code is written... how data is structured... provides an intuitive way to work

```r
install.packages("tidyverse", dependencies = TRUE)
install.packages("raster", dependencies = TRUE)
install.packages("sf", dependencies = TRUE)
install.packages("landscapetools", dependencies = TRUE)
install.packages("landscapemetrics", dependencies = TRUE)
```

.small[
Note:  
- Enter `n` if you get the following prompt: `Do you want to install from sources the package which needs compilation?`  
- Click 'Yes' if you are asked to restart R  
]
???
- If stuck in restart loop, click no
- If lib not writable, use personal library (admin permissions?)

---

# Outline  
 
[Why analyse spatial patterns?](#why)

[Landscape ecology: Conceptual models](#concepts)

[Recap: R Environment](#renviron)

__[Land cover classification](#classify)__

[Landscape metrics](#lsm)

[Useful resources](#resources)

???
- In our workshop: Classify landscape into discrete land cover classes (mosaic model)

---

# Land cover classification

- Remotely sensed data (i.e. satellites, drones, planes)
- Historical maps
- Field surveys
]

???
- Google maps + photoshop
- We want to use a more robust method

.right-column[ 
We will be using freely available [Landsat](https://landsat.gsfc.nasa.gov/a-landsat-timeline/) satellite images from the public database https://earthexplorer.usgs.gov/
.small[
- Sign up for an account
- Specify search criteria (i.e. location & date range)
- Specify dataset (Landsat > Landsat Collection 1 Level-1 > Landsat 8 OLI/TIRS C1 Level-1)
- Specify additional criteria (e.g. < 10% cloud cover; Day)
- Order scene and checkout
- Submit Order to download scenes when status is "Processing Complete"
- Install Bulk Download Application (BDA) on your computer
- Login & download ordered scenes
]
]

???
- Landsat: Specific type of satellite... offers free imagery

- OLI: Operational Land Imager
  - Spectral bands to detect clouds & NIR
- TIRS: Thermal Infrared Scanner

---

# Land cover classification

.right-column[ 
__Previewing images on USGS Earth Explorer__
 
<img src="images/usgs.png" width="600" style="display: block; margin: auto;" />
]

???
- Looking for the right image
- Try to find one that covers the entire city boundary
- LOW cloud cover, 
- Taken during SUMMER (winter no vegetation)!

---

# Land cover classification

- Satellites have multiple sensors that measure surface _reflectance_

---

# Land cover classification

- Satellites have multiple sensors that measure surface _reflectance_ 
- Specific bands of wavelengths measured along the electromagnetic spectrum (spatial resolution varies)

???
- What we see on Google Earth: RGB bands

- Satellites can detect many more bands
  - Multispectral images: IR/NIR bands
  - Hyperspectral: 100-1000 bands
- Having more bands = classify land cover more accurately

---

# Land cover classification

- Satellites have multiple sensors that measure surface _reflectance_ 
- Specific bands of wavelengths measured along the electromagnetic spectrum (spatial resolution varies)  
- Data stored as multiple layers of raster images, a.k.a. "scenes" or "bands"
]

---

# Land cover classification

.left-column[ 
### Images
]
.right-column[ 
__We're using the Landsat-8 scene of Singapore on 2018/05/24__  
(30m resolution for RGB, NIR & SWIR)

- Navigate to the *data* folder

???
- SWIR: Shortwave IR

- Contains 3 preview images: 
  - Natural color
  - Thermal infrared
  - Band quality (shows cloud cover)

---

# Land cover classification

.left-column[ 
### Images
]
.right-column[ 
__We're using the Landsat-8 scene of Singapore on 2018/05/24__  
(30m resolution for RGB, NIR & SWIR)

- Navigate to *data/LC08_L1TP_125059_20180524_20180605_01_T1* (named according to [product identifier](https://www.usgs.gov/faqs/what-naming-convention-landsat-collections-level-1-scenes?qt-news_science_products=0#qt-news_science_products))
- Unzip each file into individual bands (*.TIF* rasters)
- Unzipped folder will contain 5 multispectral bands of the same region 
]

---

# Land cover classification

.left-column[ 
### Images
]
.right-column[ 
__Bands:__ 
- __B1:__ Ultra Blue
- __B2:__ Blue
- __B3:__ Green
- __B4:__ Red
- __B5:__ Near Infrared (NIR)
- __B6:__ Shortwave Infrared (SWIR) 1
- __B7:__ SWIR 2
- __B8:__ Panchromatic
- __B9:__ Cirrus
- __B10:__ Thermal Infrared (TIRS) 1
- __B11:__ TIRS 2
- __BQA:__ Band quality assessment (based on cloud cover)
] 
???
- Because file sizes are too large, I've only uploaded bands 1-5

---

# Land cover classification

.left-column[ 
### Images
### Import
]
.right-column[ 
__Let's import the visible light bands into R!__

```r
library(raster)

#blue
b2 <- raster('data/Landsat 8 OLI_TIRS C1 Level-1/LC08_L1TP_125059_20180524_20180605_01_T1/LC08_L1TP_125059_20180524_20180605_01_T1_B2.tif')

#green
b3 <- raster('data/Landsat 8 OLI_TIRS C1 Level-1/LC08_L1TP_125059_20180524_20180605_01_T1/LC08_L1TP_125059_20180524_20180605_01_T1_B3.tif')

#red
b4 <- raster('data/Landsat 8 OLI_TIRS C1 Level-1/LC08_L1TP_125059_20180524_20180605_01_T1/LC08_L1TP_125059_20180524_20180605_01_T1_B4.tif')
```

]

???
- Let's do a quick visualisation
- Print raster for details, see in environment

---

# Land cover classification

```r
landsatRGB <- stack(b4, b3, b2) #order is impt

plotRGB(landsatRGB, 
        stretch = "lin")
```

<img src="slides_files/figure-html/unnamed-chunk-13-1.png" style="display: block; margin: auto;" />
.center[.small[Landsat-8 true color composite (RGB). Source: U.S. Geological Survey.]]
]

???
- `?plotRGB`
- Try changing `stretch = ` argument to `"hist"`

- Notice that the image we have covers a very large area
- Let's crop it according to the boundaries of SG

---

# Land cover classification

.left-column[ 
### Images
### Import
]
.right-column[ 
__Import bands 1-5 as a RasterStack named `landsat`__

```r
filenames <- paste0('data/Landsat 8 OLI_TIRS C1 Level-1/LC08_L1TP_125059_20180524_20180605_01_T1/LC08_L1TP_125059_20180524_20180605_01_T1_B', 1:5, ".tif")
landsat <- stack(filenames) 
```

] 
???
- First, lets import more bands as a stack

- `1:5` creates a sequence

```r
names(landsat) <- c('ultra-blue', 'blue', 'green', 'red', 'NIR')
```
]

```r
crs(landsat)
```

```
## Coordinate Reference System:
## Deprecated Proj.4 representation:
##  +proj=utm +zone=48 +datum=WGS84 +units=m +no_defs 
## WKT2 2019 representation:
## PROJCRS["unknown",
##     BASEGEOGCRS["unknown",
##         DATUM["World Geodetic System 1984",
##             ELLIPSOID["WGS 84",6378137,298.257223563,
##                 LENGTHUNIT["metre",1]],
##             ID["EPSG",6326]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8901]]],
##     CONVERSION["UTM zone 48N",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",105,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9996,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",500000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]],
##         ID["EPSG",16048]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1,
##                 ID["EPSG",9001]]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1,
##                 ID["EPSG",9001]]]]
```
]
???

- Projection: Universal Transverse Mercator
- Datum: World Geodetic System 1984 (EPSG code 4326) - used by Google Earth

---

# Land cover classification

.left-column[ 
### Images
### Import
### Crop
]
.right-column[ 
__Use city boundaries to crop out the landscape of interest__

Polygon file in *data* folder, downloaded from [data.gov.sg](https://data.gov.sg/dataset/master-plan-2014-region-boundary-web)

<img src="slides_files/figure-html/unnamed-chunk-20-1.png" style="display: block; margin: auto;" />
.center[.small[Singapore city boundaries based on Regional Master Plan 2014. Source: Urban Redevelopment Authority (URA)]]

]

---

# Land cover classification

.left-column[ 
### Images
### Import
### Crop
]
.right-column[ 
__Import polygon of Singapore boundaries as `sgshp`__

```r
library(sf)
sgshp <- st_read("data/master-plan-2014-region-boundary-web-shp/MP14_REGION_WEB_PL.shp")
```
]
???
- sgshp spatial info stored in "geometry" col

--
.right-column[ 
__Check CRS of `sgshp`__

```r
crs(sgshp)
```

```
## [1] "PROJCRS[\"SVY21\",\n    BASEGEOGCRS[\"SVY21[WGS84]\",\n        DATUM[\"World Geodetic System 1984\",\n            ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n                LENGTHUNIT[\"metre\",1]],\n            ID[\"EPSG\",6326]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"Degree\",0.0174532925199433]]],\n    CONVERSION[\"unnamed\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",1.36666666666667,\n            ANGLEUNIT[\"Degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",103.833333333333,\n            ANGLEUNIT[\"Degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",1,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",28001.642,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",38744.572,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1,\n                ID[\"EPSG\",9001]]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1,\n                ID[\"EPSG\",9001]]]]"
```
]

???
- CRS is different!
--
.right-column[ 
__Transform CRS of `sgshp` to match `landsat`__

```r
sgshp <- st_transform(sgshp, crs = crs(landsat))
```
]

???
- `?st_transform`

---

# Land cover classification

.left-column[ 
### Images
### Import
### Crop
]
.right-column[ 
__Crop RasterStack `landsat` to the extent of `sgshp`__ 1

```r
landsat <- crop(landsat, sgshp) # crop to rectangle
```
.footnote[[1] The polygon should be in 2D format. Use `st_zm()` to drop the Z-dimension if it is a 3D object.]
]
--
.right-column[ 
__Mask values in `landsat` based on shape of `sgshp`__

```r
landsat <- mask(landsat, sgshp)
```
]

---

# Land cover classification

```r
landsatRGB <- subset(landsat, c(4,3,2))
plotRGB(landsatRGB,
 stretch = "lin")
```

<img src="slides_files/figure-html/landsat rgb-1.png" style="display: block; margin: auto;" />
.center[.small[Landsat-8 true color composite (USGS, 2018) cropped to city boundaries (URA, 2014)]]
]

???
- Extract the RBG bands again (order is impt)

---

# Land cover classification

- Use [vegetation indices](http://www.un-spider.org/links-and-resources/data-sources/daotm/daotm-vegetation) such as the [Normalized Difference Vegetation Index (NDVI)](https://gisgeography.com/ndvi-normalized-difference-vegetation-index/) 
- Ranges from -1 to 1 (higher = greener)
]
???
- Image classification is a huge field: Supervised/unsupervised...
- Our workflow is quite rudimentary, low accuracy

- __What wavelengths (colors) of light do plants absorb?__
--
.right-column[
- Vegetation reflects NIR & absorbs red light, captured in the formula:

$$
NDVI = \dfrac{NIR - Red}{NIR + Red}
$$
]

---

# Land cover classification

$$
NDVI = \dfrac{NIR - Red}{NIR + Red}
$$
]
 
--
.right-column[

```r
ndvi <- function(x, y) {
 (x - y) / (x + y)
}
```
]
???
- Revise: general structure of function

---

# Land cover classification

```r
landsatNDVI <- overlay(landsat[[5]], landsat[[4]], fun = ndvi)
```
]

???
- `?overlay`
- `range(landsatNDVI)`, `hist(landsatNDVI)`
- There can be abnormal values from atmospheric correction... must read product guide
--
.right-column[
__Limit range of values to be from -1 to 1__

```r
landsatNDVI <- reclassify(landsatNDVI, c(-Inf, -1, -1)) # <-1 becomes -1
landsatNDVI <- reclassify(landsatNDVI, c(1, Inf, 1)) # >1 becomes 1
```
]
???
- `?reclassify` (need vector of c("from", "to", "value"))
- `hist(landsatNDVI)` again

---

# Land cover classification

```r
plot(landsatNDVI, 
     col = rev(terrain.colors(10)), 
     main = "Landsat 8 NDVI",
     axes = FALSE, box = FALSE)
```

<img src="slides_files/figure-html/landsat ndvi-1.png" style="display: block; margin: auto;" />
]
???
- We're plotting a single raster
- `plot(landsatNDVI)`

---

# Land cover classification

.left-column[ 
### Images
### Import
### Crop
### Classify
]
.right-column[
__NDVI ranges for land cover classes derived from Landsat-8__1 
<table class="table" style="">
 <thead>
 <tr>
 <th style="text-align:left;"> Class </th>
 <th style="text-align:left;"> NDVI </th>
 </tr>
 </thead>
<tbody>
 <tr>
 <td style="text-align:left;"> Water </td>
 <td style="text-align:left;"> −0.28–0.015 </td>
 </tr>
 <tr>
 <td style="text-align:left;"> Built-up </td>
 <td style="text-align:left;"> 0.015–0.14 </td>
 </tr>
 <tr>
 <td style="text-align:left;"> Barren Land </td>
 <td style="text-align:left;"> 0.14–0.18 </td>
 </tr>
 <tr>
 <td style="text-align:left;"> Shrub and Grassland </td>
 <td style="text-align:left;"> 0.18–0.27 </td>
 </tr>
 <tr>
 <td style="text-align:left;"> Sparse Vegetation </td>
 <td style="text-align:left;"> 0.27–0.36 </td>
 </tr>
 <tr>
 <td style="text-align:left;"> Dense Vegetation </td>
 <td style="text-align:left;"> 0.36–0.74 </td>
 </tr>
</tbody>
</table>

???
- Water: negative?
- Veg: >0.2?

---

# Land cover classification

```r
hist(landsatNDVI,
     main = "Distribution of NDVI values", xlab = "NDVI", 
     xlim = c(-1, 1), breaks = 100, yaxt = 'n')
abline(v=0.2, col="red", lwd=2)
abline(v=0, col="red", lwd=2)
```

<img src="slides_files/figure-html/ndvi histogram-1.png" style="display: block; margin: auto;" />
]
???
- `hist(landsatNDVI, breaks = 100)` - see in greater detail
- Distinct peaks

--
.right-column[
__Let's define vegetation cover by threshold value 0.2__
]
???
- Values >0.2 are typically vegetation

---

# Land cover classification

```r
landsatGreen <- reclassify(landsatNDVI, c(-1, 0.2, NA)) #-1 to 0.2 becomes NA
plot(landsatGreen, col = "darkgreen", 
 axes = FALSE, box = FALSE, legend = FALSE)
```
]
--
.right-column[
<div class="figure" style="text-align: center">
<img src="slides_files/figure-html/landsat ndvi green-1.png" alt="Looks like vegetation cover" />
Looks like vegetation cover
</div>
]

---

# Land cover classification

```r
landsatBlue <- reclassify(landsatNDVI, c(0, 1, NA)) # 0 to 1 becomes NA
plot(landsatBlue, col = "blue", 
 axes = FALSE, box = FALSE, legend = FALSE,)
```

<div class="figure" style="text-align: center">
<img src="slides_files/figure-html/landsat ndvi blue-1.png" alt="Looks like water cover" />
Looks like water cover
</div>
]

???
- Visualise only the NDVI with negative values

---

# Land cover classification

<img src="slides_files/figure-html/unnamed-chunk-35-1.png" style="display: block; margin: auto;" />
]
???
- A Mosaic model with discrete classes

---

# Land cover classification

```r
reclass_m <- matrix(c(-Inf, 0, 1, #water
 0, 0.2, 2, #urban
 0.2, Inf, 3), #veg
 ncol = 3, byrow = TRUE)
reclass_m
```

```
##      [,1] [,2] [,3]
## [1,] -Inf  0.0    1
## [2,]  0.0  0.2    2
## [3,]  0.2  Inf    3
```
]
???
- A bit more complicated
- This time, we need a matrix (a table)
- The matrix should have 3 columns (`from=, to=, becomes=`)
- `?reclassify`

--
.right-column[
__Create a new raster `landsatCover` using the threshold values defined in `reclass_m`__

```r
landsatCover <- reclassify(landsatNDVI, rcl = reclass_m)
```

]

---

# Land cover classification

__Plot `landsatCover` with custom colours__

```r
plot(landsatCover,
     col = c("blue", "grey", "darkgreen"),
     legend = FALSE,
     axes = FALSE,
     box = FALSE,
     main = "Land cover (mosaic) in Singapore")
legend("bottomright",
       legend = c("Water", "Urban", "Vegetation"),
       fill = c("blue", "grey", "darkgreen"),
       border = FALSE,
       bty = "n")
```
]

???
- Remember how to plot? `plot(landsatCover)`
  - `?plot` for customisation (i.e. `col = c("blue", "grey", "darkgreen")`)  
  
- Proper analysis: need validation with a reference map

---

# Land cover classification

---

# Land cover classification

.left-column[ 
### Images
### Import
### Crop
### Classify
### Save
]
.right-column[
__Save raster `landsatCover` to disk__

```r
writeRaster(landsatCover, 
            filename = "clean_data/landsat_landcover.tif", 
            overwrite = TRUE)

#save as GeoTiff format
```
]
???
- `overwrite = TRUE` overwrites existing file (do so with caution)

- If not, you have to download the course materials again
  - Can save using another file name

---

# Land cover classification

##💭 What if the satellite image does not cover the whole city?

1) Calculate the NDVI for multiple images (each a *single* raster layer)

2) Combine the rasters by using the mean NDVI between images

3) Use combined image to derive land cover classes

???
- FAQ
- Instead of having to deal with raster stacks

---

# Land cover classification

##💭 What if the satellite image does not cover the whole city?

__If CRS of images are different, re-project one to match the other__

```r
NDVIraster2 <- projectRaster(NDVIraster2, NDVIraster1) #overwrite NDVIraster2
```

__Calculate the mean NDVI for overlapping sections between images__

```r
NDVIcombined <- mosaic(NDVIraster1, NDVIraster2, fun = mean) #mean() function applied
```

???
- Remember to assign different names to the data!
- `?projectRaster` for more details

---

# Outline  
 
[Why analyse spatial patterns?](#why)

[Landscape ecology: Conceptual models](#concepts)

[Recap: R Environment](#renviron)

[Land cover classification](#classify)

__[Landscape metrics](#lsm)__

[Useful resources](#resources)

???
- Use numbers to analyse and tell the story about landscapes

---

# Landscape metrics

- Quantifies spatial patterns in a patch-corridor-matrix (mosaic) model

???
- Discrete landscape classes

- We will be using the package `landscapemetrics`, which is based on the software [FRAGSTATS](http://www.umass.edu/landeco/research/fragstats/fragstats.html)

???
- There are tons of metrics we can use.
- But, they are meaningless numbers... unless you can link it to an explicit ecological process
- This is the question u should be asking yourself...

---

# Landscape metrics

__3 levels of analyses:__ 
- Patch `lsm_p_<metric>()` 
- Class `lsm_c_<metric>()` 
- Landscape `lsm_l_<metric>()`

???
- Patch: 
  - Every single patch will have a result (regardless of class)
  - Most "zoomed-in"
- Class: 
  - Every class will have a result (summarises all patches per class)  
- Landscape: 
  - Whole landscape has one result/value

- Functions follow a consistent naming convention

---

# Landscape metrics

<table class="table" style="">
 <thead>
 <tr>
 <th style="text-align:left;"> Composition </th>
 <th style="text-align:left;"> Configuration </th>
 </tr>
 </thead>
<tbody>
 <tr>
 <td style="text-align:left;"> Proportional abundance </td>
 <td style="text-align:left;"> Patch area &amp; edge </td>
 </tr>
 <tr>
 <td style="text-align:left;"> Richness </td>
 <td style="text-align:left;"> Patch shape complexity </td>
 </tr>
 <tr>
 <td style="text-align:left;"> Evenness </td>
 <td style="text-align:left;"> Core area </td>
 </tr>
 <tr>
 <td style="text-align:left;"> Diversity </td>
 <td style="text-align:left;"> Isolation/proximity </td>
 </tr>
 <tr>
 <td style="text-align:left;"> </td>
 <td style="text-align:left;"> Contrast </td>
 </tr>
 <tr>
 <td style="text-align:left;"> </td>
 <td style="text-align:left;"> Dispersion </td>
 </tr>
 <tr>
 <td style="text-align:left;"> </td>
 <td style="text-align:left;"> Contagion &amp; interspersion </td>
 </tr>
 <tr>
 <td style="text-align:left;"> </td>
 <td style="text-align:left;"> Subdivision </td>
 </tr>
 <tr>
 <td style="text-align:left;"> </td>
 <td style="text-align:left;"> Connectivity </td>
 </tr>
</tbody>
</table>

.center[.small[More info available in [documentation](http://www.umass.edu/landeco/research/fragstats/documents/Conceptual%20Background/Landscape%20Metrics/Landscape%20Metrics.htm) by
UMass Landscape Ecology Lab
]]

???
- An overview of the subject areas/metrics tt can be quantified

- Grouped according to composition/configuration
  - ...each of which can be patch/class/landscape-lvl (many diff combinations)

---

# Landscape metrics

<img src="slides_files/figure-html/unnamed-chunk-49-1.png" style="display: block; margin: auto;" />
__Patch__

.small[
__Area & edge:__ Area `area`, Perimeter `perim`, Perimeter-area-ratio `para`  
__Shape:__ Shape `shape`   
]

.small[
*Note:* See [R documentation](https://cran.r-project.org/web/packages/landscapemetrics/landscapemetrics.pdf) for all functions. More egs grouped by _subject area_ & _level of analysis_: [Area/density/edge](http://www.umass.edu/landeco/research/fragstats/documents/Metrics/Area%20-%20Density%20-%20Edge%20Metrics/AREA%20-%20DENSITY%20-%20EDGE%20METRICS.htm), [Shape](http://www.umass.edu/landeco/research/fragstats/documents/Metrics/Shape%20Metrics/SHAPE%20METRICS.htm), [Diversity](http://www.umass.edu/landeco/research/fragstats/documents/Metrics/Diversity%20Metrics/DIVERSITY%20METRICS.htm), [Contagion](http://www.umass.edu/landeco/research/fragstats/documents/Metrics/Contagion%20-%20Interspersion%20Metrics/CONTAGION%20-%20INTERSPERSION%20METRICS.htm).  
]

???
- These are only some basic examples!

---

# Landscape metrics

<img src="slides_files/figure-html/unnamed-chunk-50-1.png" style="display: block; margin: auto;" />
__Class__

.small[
__Area & edge:__ Mean patch area `area_mn`,  Edge density `ed`, Percentage of landscape  `pland`, Largest patch index `lpi`  
__Shape:__ Mean perimeter-area ratio `para_mn`, Mean shape index `shape_mn`   
__Aggregation:__ No. of patches `np`, Patch density `pd`, Mean of euclidean nearest-neighbor distance `enn_mn`, Landscape shape index `lsi`   
]

???
- Edge density: edge length/area

---

# Landscape metrics

<img src="slides_files/figure-html/unnamed-chunk-51-1.png" style="display: block; margin: auto;" />
__Landscape__

.small[
__Area & edge:__ Total area `ta`, Mean patch area `area_mn`, Total edge `te`, Edge density `ed`  
__Shape:__ Mean perimeter-area ratio `para_mn`, Mean shape index `shape_mn`  
__Aggregation:__ No. of patches `np`, Patch density `pd`, Mean of euclidean   nearest-neighbor distance `enn_mn`, Landscape shape index `lsi`  
__Diversity:__ Patch richness `pr`, Shannon's diversity index `shdi`, Simpson's diversity index `sidi`  
]

.small[
*Note:* See [R documentation](https://cran.r-project.org/web/packages/landscapemetrics/landscapemetrics.pdf) for all functions. More egs grouped by _subject area_ & _level of analysis_: [Area/density/edge](http://www.umass.edu/landeco/research/fragstats/documents/Metrics/Area%20-%20Density%20-%20Edge%20Metrics/AREA%20-%20DENSITY%20-%20EDGE%20METRICS.htm), [Shape](http://www.umass.edu/landeco/research/fragstats/documents/Metrics/Shape%20Metrics/SHAPE%20METRICS.htm), [Diversity](http://www.umass.edu/landeco/research/fragstats/documents/Metrics/Diversity%20Metrics/DIVERSITY%20METRICS.htm), [Contagion](http://www.umass.edu/landeco/research/fragstats/documents/Metrics/Contagion%20-%20Interspersion%20Metrics/CONTAGION%20-%20INTERSPERSION%20METRICS.htm). 
]

???
- Patch richness: no. of classes

---

# Landscape metrics

```r
landsatCover <- raster('clean_data/landsat_landcover.tif')
```

]
???
- If it's not in your workspace
--
.right-column[
__Load required packages__

```r
library(landscapemetrics)
library(landscapetools)
```
]

---

# Landscape metrics

```r
show_landscape(landsatCover, discrete = TRUE) 
```

<img src="slides_files/figure-html/unnamed-chunk-55-1.png" style="display: block; margin: auto auto auto 0;" />
]
???
- from `library(landscapetools)`
- Remind about `landscapetools::`
- Default colors shown - not a big deal
- Focus on 3 (vegetation) & 2 (urban)

---

# Landscape metrics

__Check if `landsatCover` is in the right format__

```r
check_landscape(landsatCover)
```

```
##   layer       crs units   class n_classes OK
## 1     1 projected     m integer         3  ✔
```
]

???
- `landscapemetrics::`
- Units: m
- Landscape classes need to be integers
- Will throw warning if too many classes

---

# Landscape metrics

```r
list_lsm()
```

```
## # A tibble: 133 × 5
## metric name type level function_name
## <chr> <chr> <chr> <chr> <chr> 
## 1 area patch area area and edge… patch lsm_p_area 
## 2 cai core area index core area met… patch lsm_p_cai 
## 3 circle related circumscribing circle shape metric patch lsm_p_circle 
## 4 contig contiguity index shape metric patch lsm_p_contig 
## 5 core core area core area met… patch lsm_p_core 
## 6 enn euclidean nearest neighbor distance aggregation m… patch lsm_p_enn 
## 7 frac fractal dimension index shape metric patch lsm_p_frac 
## 8 gyrate radius of gyration area and edge… patch lsm_p_gyrate 
## 9 ncore number of core areas core area met… patch lsm_p_ncore 
## 10 para perimeter-area ratio shape metric patch lsm_p_para 
## # ℹ 123 more rows
```
]
???
- Tibble is basically a df (more stringent rules)
- Save as variable `list <- list_lsm()` to view separately

---

# Landscape metrics

```r
list_lsm(level = "patch", type = "area and edge metric")
```

```
## # A tibble: 3 × 5
## metric name type level function_name
## <chr> <chr> <chr> <chr> <chr> 
## 1 area patch area area and edge metric patch lsm_p_area 
## 2 gyrate radius of gyration area and edge metric patch lsm_p_gyrate 
## 3 perim patch perimeter area and edge metric patch lsm_p_perim
```
]
???
- `?list_lsm` for possible arguments
- Most basic measure of composition: area

---

# Landscape metrics

```r
lsm_p_area(landsatCover)
```

```
## # A tibble: 12,580 × 6
## layer level class id metric value
## <int> <chr> <int> <int> <chr> <dbl>
## 1 1 patch 1 1 area 14.6 
## 2 1 patch 1 2 area 0.09
## 3 1 patch 1 3 area 27.3 
## 4 1 patch 1 4 area 0.09
## 5 1 patch 1 5 area 0.18
## 6 1 patch 1 6 area 0.36
## 7 1 patch 1 7 area 0.09
## 8 1 patch 1 8 area 0.09
## 9 1 patch 1 9 area 0.09
## 10 1 patch 1 10 area 0.09
## # ℹ 12,570 more rows
```

.small[*Note:* Use `?lsm_p_area` for more info on the function]  
]
???
- Example of patch-level: area
- class 1: water
- What unit is the value in? Check with `?lsm_p_area`
  - Hectares

---

# Landscape metrics

```r
lsm_c_ca(landsatCover) #total class area
```

```
## # A tibble: 3 × 6
## layer level class id metric value
## <int> <chr> <int> <int> <chr> <dbl>
## 1 1 class 1 NA ca 6742.
## 2 1 class 2 NA ca 28496.
## 3 1 class 3 NA ca 42933.
```

```r
lsm_c_area_mn(landsatCover) #mean area
```

```
## # A tibble: 3 × 6
## layer level class id metric value
## <int> <chr> <int> <int> <chr> <dbl>
## 1 1 class 1 NA area_mn 3.22
## 2 1 class 2 NA area_mn 4.52
## 3 1 class 3 NA area_mn 10.3
```
]

---

# Landscape metrics

```r
lsm_c_pland(landsatCover) #percentage of landscape
```

```
## # A tibble: 3 × 6
## layer level class id metric value
## <int> <chr> <int> <int> <chr> <dbl>
## 1 1 class 1 NA pland 8.62
## 2 1 class 2 NA pland 36.5 
## 3 1 class 3 NA pland 54.9
```
]
--
.right-column[
Vegetation cover of ~50% is relatively similar to results reported in Gaw et al. (2019), Data, 4(3), 116.
]

???
- Our image was taken in 2018
- `pland` of class 3 (vegetation)

--
.right-column[
__What are some reasons why our results may differ from actual land cover?__
]
???
- Landset resolution (small patches not captured)
- NDVI threshold value varies between images/satellite/regions
- Cloud cover
- City boundaries vary/not standardised
- Time of year (autumn/winter)

---

# Landscape metrics

```r
lsm_l_ta(landsatCover) #total area of landscape
```

```
## # A tibble: 1 × 6
## layer level class id metric value
## <int> <chr> <int> <int> <chr> <dbl>
## 1 1 landscape NA NA ta 78170.
```
]
???
- 780 sqkm

---

# Landscape metrics

```r
calculate_lsm(landsatCover, 
              what = c("lsm_l_ta", #total area of landscape
                       "lsm_c_ca", #total area of each class
                       "lsm_c_pland"), #proportional area of each class
              full_name = TRUE) #include info about metrics in results
```

```
## # A tibble: 7 × 9
## layer level class id metric value name type function_name
## <int> <chr> <int> <int> <chr> <dbl> <chr> <chr> <chr> 
## 1 1 class 1 NA ca 6742. total (class)… area… lsm_c_ca 
## 2 1 class 2 NA ca 28496. total (class)… area… lsm_c_ca 
## 3 1 class 3 NA ca 42933. total (class)… area… lsm_c_ca 
## 4 1 class 1 NA pland 8.62 percentage of… area… lsm_c_pland 
## 5 1 class 2 NA pland 36.5 percentage of… area… lsm_c_pland 
## 6 1 class 3 NA pland 54.9 percentage of… area… lsm_c_pland 
## 7 1 landscape NA NA ta 78170. total area area… lsm_l_ta
```
.small[
__Note:__ use `options(scipen=999)` to disable scientific notation for reported values
]
]
???
- `?calculate_lsm`

---

# Landscape metrics

.left-column[ 
### Load
### Analyse
]
.right-column[
__Calculate multiple metrics by level and type__

```r
calculate_lsm(landsatCover, 
              level = "class",
              type = "aggregation metric",
              full_name = TRUE) #include info about metrics in results
```

```
## # A tibble: 45 × 9
## layer level class id metric value name type function_name
## <int> <chr> <int> <int> <chr> <dbl> <chr> <chr> <chr> 
## 1 1 class 1 NA ai 88.4 aggregation index aggr… lsm_c_ai 
## 2 1 class 2 NA ai 79.5 aggregation index aggr… lsm_c_ai 
## 3 1 class 3 NA ai 87.8 aggregation index aggr… lsm_c_ai 
## 4 1 class 1 NA clumpy 0.889 clumpiness index aggr… lsm_c_clumpy 
## 5 1 class 2 NA clumpy 0.687 clumpiness index aggr… lsm_c_clumpy 
## 6 1 class 3 NA clumpy 0.735 clumpiness index aggr… lsm_c_clumpy 
## 7 1 class 1 NA cohesion 95.5 patch cohesion i… aggr… lsm_c_cohesi…
## 8 1 class 2 NA cohesion 99.2 patch cohesion i… aggr… lsm_c_cohesi…
## 9 1 class 3 NA cohesion 99.8 patch cohesion i… aggr… lsm_c_cohesi…
## 10 1 class 1 NA division 1.00 division index aggr… lsm_c_divisi…
## # ℹ 35 more rows
```
]

???
- These tell you sth about connectivity between patches

---

# Landscape metrics

.left-column[ 
### Load
### Analyse
]
.right-column[
__Calculate multiple metrics by level and type__

```r
calculate_lsm(landsatCover, 
              level = "landscape",
              type = "diversity metric",
              full_name = TRUE)
```

```
## # A tibble: 6 × 9
## layer level class id metric value name type function_name
## <int> <chr> <int> <int> <chr> <dbl> <chr> <chr> <chr> 
## 1 1 landscape NA NA msidi 0.817 modified simp… dive… lsm_l_msidi 
## 2 1 landscape NA NA msiei 0.743 modified simp… dive… lsm_l_msiei 
## 3 1 landscape NA NA pr 3 patch richness dive… lsm_l_pr 
## 4 1 landscape NA NA prd 0.00384 patch richnes… dive… lsm_l_prd 
## 5 1 landscape NA NA rpr NA relative patc… dive… lsm_l_rpr 
## 6 1 landscape NA NA shdi 0.908 shannon's div… dive… lsm_l_shdi
```
]
--
.right-column[
__What do these "diversity metrics" actually mean?__
]

???
- Are they referring to biodiversity?

---

# Landscape metrics

.left-column[ 
### Load
### Analyse
### Discuss
]
.right-column[
__How can we use our understanding of spatial patterns to improve landscape planning and design?__ 
]

---

# Landscape metrics

###Ecosystem services framework

__Discuss how landscape metrics link to ecological processes that provide benefits or 'services' to humans__

---

# Landscape metrics

###Ecosystem services framework

__Discuss how landscape metrics link to ecological processes that provide benefits or 'services' to humans__

???
- Range of ecosystem services that are associated with city landscapes

---

# Landscape metrics

### E.g. Vegetation cover and temperature regulation

### Correlation between metric & land surface temp (LST) 1

__Percentage of landscape `pland`__: -0.7

- Larger proportional area of greenery associated with a lower LST

- Beijing should maintain existing amounts of greenery and adopt appropriate urban greening measures

???
- Most basic type: correlation analysis
- This study was done in Beijing, China

---

# Landscape metrics

### E.g. Vegetation cover and temperature regulation

### Correlation between metric & land surface temp (LST) 1

__Mean patch area `area_mn`__: -0.35

- Larger patches of greenery (less disturbed) could help lower temperatures

- Greening programs should reduce the fragmentation of green patches

---

# Landscape metrics

### E.g. Vegetation cover and temperature regulation

### Correlation between metric & land surface temp (LST) 1

__Patch density `pd`__: 0.39  
__Edge density `ed`__: 0.31

- Increased fragmentation/shape complexity of greenery is related to higher temperatures

- Policy measures should aim to decrease fragmentation

???
- PD: total no. of patches/area
- ED: total perimeter/area

---

# Landscape metrics

###💭 How can we compare metrics between different cities?

- Urban development: Sprawling vs. compact cities – Big difference in total area?
- City boundaries: Rural hinterlands included?

???
- Food for thought

--

__How can we account for such differences?__

???
- Analyse city core only?
- Does the metric account for total area (if not, / total area)?
- / population size?

---

# Questions?

[Why analyse spatial patterns?](#why)

[Landscape ecology: Conceptual models](#concepts)

[Recap: R Environment](#renviron)

[Land cover classification](#classify)

[Landscape metrics](#lsm)

[Useful resources](#resources)

---

# Useful Resources

__Free satellite imagery__   
- [Overview](https://gisgeography.com/free-satellite-imagery-data-list/)  
- http://earthexplorer.usgs.gov/
- https://lpdaacsvc.cr.usgs.gov/appeears/
- https://scihub.copernicus.eu/
- https://search.earthdata.nasa.gov/search

__R Resources__
- [Spatial data science with R](https://www.rspatial.org/raster/index.html)
- [Remote Sensing Image Analysis in R](https://www.rspatial.org/raster/rs/index.html)
- [R packages for spatial data analyses](https://cran.r-project.org/web/views/Spatial.html)
- [Interactive visualisations using `mapview`](https://r-spatial.github.io/mapview/index.html)