01_z07_nunn_2012.rmd

### Project Overview

---

**Objective:**

@nunn2012ruggedness showed empirically that the ruggedness of the terrain has had a positive impacts on economic developement in African countries. In this demonstration, we calculate Terrain Ruggedness Index (TRI) for African countries from the world elevation data.


---

**Datasets**

* World elevation data 
* World country borders

---

**GIS tasks**

* read a raster file 
  - use `terra::rast()`
* import world country border data
  - use `rnaturalearth::ne_countries()`
* crop a raster data to a particular region
  - use `terra::crop()`
* replace cell values
  - use `terra::subst()`
* calculate TRI 
  - use `terra::focal()`
+ create maps 
  * use the `tmap` package 

---

**Preparation for replication**

Run the following code to install or load (if already installed) the `pacman` package, and then install or load (if already installed) the listed package inside the `pacman::p_load()` function.

```{r demo_nunn_2012_packages}
if (!require("pacman")) install.packages("pacman")
pacman::p_load(
  sf, # vector data operations
  tidyverse, # data wrangling
  stars,
  tmap,
  raster,
  rnaturalearth,
  skimr
)
```

### Project Demonstration

We first read the world elevation data using `terra::rast()`.

```{r read-dem}
(
  dem <- terra::rast("Data/nunn_2012/GDEM-10km-BW.tif")
)
```

```{r map-world-elev, cache = TRUE}
tm_shape(dem) +
  tm_raster()
```

In this dataset, the elevation of the ocean floor is recorded as 0, so let's replace an elevation of 0 with `NA`s. This avoids a problem associated with calculating TRI later. 

```{r replace-zero}
dem <- terra::subst(dem, 0, NA)
```

```{r map-world-elev-no-zero, cache = TRUE}
tm_shape(dem) +
  tm_raster()
```

Now, since our interest is in Africa, let's just crop the raster data to its African portion. To do so, we first get the world map using `rnaturalearth::ne_countries()` and filter out the non-African countries.

```{r get-africa-sf}
africa_sf <-
  #--- get an sf of all the countries in the world ---#
  rnaturalearth::ne_countries(scale = "medium", returnclass = "sf") %>%
  #--- filter our non-African countries ---#
  filter(continent == "Africa")
```

We can now apply `terra::crop()` to `dem` based on the bounding box of `africa_sf`.

```{r crop-to-africa}
africa_dem <- terra::crop(dem, africa_sf)
```

Here is a map of the crop data.

```{r map-africa-elev, cache = TRUE}
tm_shape(africa_dem) +
  tm_raster()
```

Now, we are ready to calculte TRI, which is defined as

$$
  Ruggedness_{r,c} = \sqrt{\sum_{i=r-1}^{r+1} \sum_{j= r-1}^{r+1} (e_{i,j} - e_{r,c})^2}
$$

We are going to loop through the raster cells and calculate TRI. To do so, we make use of `terra::focal()` It allows you to apply a function to every cell of a raster and bring a matrix of the surrounding values as an input to the function. For example `terra::focal(raster, w = 3, fun = any_function)` will pass to your `any_function` a 3 by 3 matrix of the raster values centered at the point.

Let's define a function that calculates TRI for a given matrix.

```{r define-calc-tri}
calc_tri <- function(matr) {
  # matr is a length 9 matrix
  center <- matr[5]
  sum_squares <- sum((matr - center)^2, na.rm = TRUE)
  return(sqrt(sum_squares))
}
```

Now, let's calculate TRI.

```{r calc-tri}
tri_africa <-
  terra::focal(
    africa_dem,
    w = 3,
    fun = calc_tri
  )
```

Here is the map of the calculated TRI.

```{r map-tri}
tm_shape(tri_africa) +
  tm_raster()
```