Technical Guideline Series

Part 2 - Preparing and mapping data

Authors

Affiliation

Rainer M. Krug

Senckenberg

Aidin Niamir

Senckenberg

Joy Kumagai

Senckenberg

Aidin Niamir

Senckenberg

Abstract

TO BE ADDED

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

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.

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.

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.

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.

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()

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.

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.

# 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.

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.

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

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.

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.

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.

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.

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.

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.

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.

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.

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.

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

# 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.

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.

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