---
title: "Technical Guideline Series"
output:
  html_document:
    toc: yes
    toc_float: yes
    toc_depth: 4
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r, echo=FALSE}
htmltools::img(src = knitr::image_uri("C:/Users/jkumagai/Documents/IPBES/R/Geoinformatics/Technical Guidelines Series/Mapping Guidelines v2/logo2.png"), 
               alt = 'logo', 
               style = 'position:absolute; top:0; right:0; padding:10px; width: 600px; height: 135px;') 
```


## Part 2 - Preparing and mapping data
##          to IPBES Regions and Sub-regions  

**Prepared by Joy Kumagai - Technical Support Unit of Knowledge and Data**  
**Reviewed by Aidin Niamir - Head of the Technical Support Unit of Knowledge and Data**  
*For any inquires please contact tsu.data@ipbes.net*  

Version: 2.3  
Last Updated: 15 August 2022

DOI: [10.5281/zenodo.6992546](https://doi.org/10.5281/zenodo.6992546)
  
The guide will show how to aggregate and map FAO data according to the IPBES Regions and Sub-Regions polygons using R. For this exercise, we chose the FAO population data but any FAOSTAT dataset can be used.

Let's begin by loading the following packages. 

```{r message = FALSE, warning=FALSE}
library(sf) 
library(dplyr)
library(magrittr)
library(FAOSTAT)
library(httr) # to download data off of Zenodo
library(rnaturalearth) # download ocean data from natural earth 
library(graticule) # for mapping 
```

### I. Downloading Necessary Data  

#### *A. Downloading FAO data*  

The first step is to download the FAO data using the FAOSTAT package. The url can be found in FAO STAT's data description file [here.](http://fenixservices.fao.org/faostat/static/bulkdownloads/datasets_E.xml)  

```{r}
FAOSTAT::download_faostat_bulk("http://fenixservices.fao.org/faostat/static/bulkdownloads/Population_E_All_Data_(Normalized).zip", getwd())
```

#### *B. Downloading IPBES regions and subregions*   

Now we will download the shapefile of the IPBES Regions and Sub-regions off of Zenodo. This can be accomplished manually or through a few lines of code.   

To download the shapefile manually, please go to the [IPBES Regions and Sub-Regions Zenodo entry](https://doi.org/10.5281/zenodo.3923633).  

To do this through a script, first identify the record ID of the Zenodo entry, which is the numbers following "*zenodo.*" at the end of the URL. We then create a URL with the record ID and query the API for information about the record.  

```{r collapse=TRUE, message='hide'}
recordID <- "3923633"
url_record <- paste0("https://zenodo.org/api/records/", recordID)
record <- httr::GET(url_record)
record # Status 200 indicates a successful download
```

Now, we can inspect the contents downloaded with the function content()  
```{r}
View(content(record)) # view displays the output in a human readable form within R Studio
```

![The picture above shows the resulting R Studio window which displays what was downloaded in a human readable form.](C:/Users/jkumagai/Documents/IPBES/R/Geoinformatics/Technical Guidelines Series/Mapping Guidelines v2/view_content_zendodo_record.png)
  
This information we received contains metadata for the record, and within this we can find the specific URL to download the IPBES regions and sub-regions shapefile. We then use this URL and the function GET() to download the shapefile.  
```{r collapse=TRUE, message='hide'}
# Contains the url to download the shapefile
url_shape <- content(record)$files[[5]]$links$download 

httr::GET(url_shape, write_disk("ipbes_regions_subregions.zip", overwrite = T)) # Downloads shapefile
unzip("ipbes_regions_subregions.zip") # unzips shapefile
```
  
### II. Uploading data into R Studio 
Now that our data is on our computer, we need to upload the data into R studio and project the spatial data.  
  
```{r}
pop_raw <- FAOSTAT::read_faostat_bulk("Population_E_All_Data_(Normalized).zip") # load the population data using FAOSTAT's built in function 
shape <- sf::st_read("IPBES_Regions_Subregions2.shp") # shapefile
```

We chose to project the data into the Robinson projection as it minimizes distortions in both area and distance. To find the proj4 notation please visit [this link](https://epsg.io/54030).    
  
```{r}
crs_robin <-  "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
shape <- sf::st_transform(shape, crs_robin)
```
   
To plot the ocean in our maps, we will also download ocean data from the rnaturalearth package and project it  
```{r}
ocean <- rnaturalearth::ne_download(scale = 10, type = 'ocean', category = 'physical', returnclass = "sf")
ocean <- sf::st_transform(ocean, crs = crs_robin) # changes the projection
ocean <- ocean[,1]
```
  
### III. Cleaning the data  
The next important step is to clean the data to ensure it can be joined and mapped easily. For this example, we will filter to only include the total population for each country in 2018.    
```{r}
pop_2018 <- pop_raw %>% 
  dplyr::filter(element == "Total Population - Both sexes" &
           year == 2018) %>% 
  dplyr::select(area_code, # these columns are selected from the original data
         area,
         element,
         year,
         unit,
         value)
```
  
By examining the Area names within the dataset, one will notice that every name after Zimbabwe refers to aggregated data, therefore we will remove these from our analysis.
```{r}
tail(pop_2018$area, 34)
pop_2018 <- pop_2018[1:237, ] # Selects the first 237 records, thus removing the last 34 which are aggregated data
```
    
Finally, we need to add the ISO3 codes onto the dataframe, so we can easily join it to the IPEBS Regions and Sub-Regions data. The translateCountryCode() function provided by the FAOSTAT package allows us to easily do this.   

```{r collapse=TRUE}
pop_2018 <- FAOSTAT::translateCountryCode(data = pop_2018, from = "FAOST_CODE", to = "ISO3_CODE", "area_code") # Add's the ISO Code to the data
```
    
There are two records where no ISO3 Code was assigned: China (including mainland, Hong Kong SAR, Macao SAR, and Taiwan) and South Sudan. "China mainland" refers to the same area as "China" in our dataset, so we are safe to exclude the China (including mainland, Hong Kong SAR, Macao SAR, and Taiwan) from our analysis.  
  
For South Sudan, we will add the same ISO-3 Code we have in the IPBES Regions and Sub-regions dataset.   
```{r}
pop_2018[230,2] <- "SSD" # South Sudan
pop_2018 <- na.omit(pop_2018) # removes China (including other areas)
```
    
### IV: Joining and Aggregating 
We have all of our data downloaded locally, uploaded into R, and formatted properly. The last step is to join and aggregate the data to the IPBES regions and sub-regions shapefile.  
  
First, we join the IPBES regions and sub-regions attributes to our data table. I drop the spatial attributes of the IPBES regions and sub-regions dataset to speed up the process.   
```{r}
colnames(shape)[2] <- "ISO3_CODE" 
regions <- shape %>%   
  as.data.frame() %>% # drops the spatial attributes
  dplyr::select(ISO3_CODE, Region, Sub_Region)  # filters the columns

pop_2018 <- dplyr::left_join(x = pop_2018, y = regions, by = "ISO3_CODE") %>% # Joins data
 tidyr::drop_na() # conveenient function from tidyr package 
```
    
Secondly, we aggregate the data per IPBES regions and sub-regions. In our example, I calculate the total population per region and per sub-region using the group_by() function.   
```{r}
pop_2018 <- pop_2018 %>% 
  dplyr::group_by(Region) %>% # Grouping by regions
  dplyr::mutate(region_pop = sum(value)/1000) %>% # calculates total population (millions) per region
  dplyr::ungroup() %>% 
  dplyr::group_by(Sub_Region) %>% # Grouping by sub-region 
  dplyr::mutate(sub_region_pop = sum(value)/1000) %>% # calculates total population (millions) per sub-region
  dplyr::ungroup() 
pop_2018
```
    
Finally, we join the formatted FAO data to the spatial data we originally had so we can create maps. 
```{r}
data <- dplyr::full_join(x = shape, y = pop_2018, by = "ISO3_CODE")
```
  
### V. Mapping 
  
All that is left to do is to map the data per region and sub-region. We begin by dissolving the spatial data per region and subregion so country borders are not included.    
```{r}
data_region <- data %>% # this dissolves the data by region
  dplyr::group_by(Region.x) %>% 
  dplyr::summarise(region_pop2 = sum(value, na.rm = T)/1000) %>% 
  sf::st_cast() 

data_subregion <- data %>% # this dissolves the data by subregion
  dplyr::group_by(Sub_Region.x) %>% 
  dplyr::summarise(sub_region_pop2 = sum(value, na.rm=T)/1000) %>% 
  sf::st_cast()
```
  
Now, we choose the palette and plot by region.   
```{r}
data_region$region_pop2 <- as.character(round(data_region$region_pop2 )) # Treats the values as groups so the legend displays correctly 
data_region$region_pop2[5] <- "0927" # Ensures the legend displays correctly

# Plotting by regions
palette <- c("grey","aliceblue", "lightskyblue", "dodgerblue", "dodgerblue4") # colors


plot(data_region[,2], pal = palette, main = "Total population (millions) in 2018 per region")
```
   
If you would like to add graticules and an ocean background, follow this example. First, we will set up the graticules we will plot   
```{r message = FALSE, warning = FALSE}
# Creates latitude and longitude labels and graticules
lat <- c(-90, -60, -30, 0, 30, 60, 90)
long <- c(-180, -120, -60, 0, 60, 120, 180)
labs <- graticule::graticule_labels(lons = long, lats = lat, xline = -180, yline = 90, proj = crs_robin) # labels for the graticules 
lines <- graticule::graticule(lons = long, lats = lat, proj = crs_robin) # graticules 
```
   
Then, we will plot the graticules, ocean data, then region data, and finally the text for latitude and longitude lines, legend, and surrounding box.   
```{r}
par(mar = c(2,3,1,2)) # Adjusts the edges of the frame 
plot(lines, lty = 5, col = "lightgrey",  main = "Total population (millions) in 2018 per region") # plots graticules 
plot(ocean, col = ggplot2::alpha("slategray1", 0.3), add = TRUE) 
plot(data_region[,2], pal = palette, add = TRUE)
text(subset(labs, labs$islon), lab = parse(text = labs$lab[labs$islon]), pos = 3, xpd = NA) # plots longitude labels
text(subset(labs, !labs$islon), lab = parse(text = labs$lab[!labs$islon]), pos = 2, xpd = NA) # plots latitude labels
legend("bottom", # adding the legend last 
       legend = (data_region %>% pull(region_pop2) %>% sort()),
       fill = palette, 
       horiz = TRUE, bty = "n")
box(which = "plot", lty = "solid") # Map frame 
```
   
We can also plot by subregion.  
```{r}
plot(data_subregion[,2], main = "Total population (millions) in 2018 per sub-region", breaks = "quantile")
```
   
Your feedback on this content is welcome. Let us know what other useful material would you like to see here by emailing tsu.data@ipbes.net  