Welcome to the first script of Module 6! Here, we will visualize one of the five climate data sets in this folder with the package leaflet. This is the basis for for the app in this second part of the module.

Please read the following code instructions carefully, try to understand the code that follows each instruction, execute it and see what happens. Do not hesitate to change the code or write your own code, and execute it as well!

1. Preparations & data inspection

As usual, we start by loading the packages and data we need, and by defining our working directory.

library(randomForest)
library(dplyr)
library(ggplot2)
library(tigris)
library(rgdal)
library(viridis)

dir <- "C:/Users/Max/Desktop/IAMO_neu/eLearning/Module_6_Shiny_Apps/M06_Part02_Shiny_with_leaflet/"

marzes <- readOGR(paste0(dir, "shapefiles/ARM_marzes.shp"), verbose=F)
futyield <- read.csv(paste0(dir, "ARM_futyield.csv"))
## convert values to percent
futyield$value <- futyield$value*100

The object “futyield” contains predicted future yield changes in % for each marz and six different crops, for two different RCPs and two different future time periods:

head(futyield)
region crop value RCP period
Aragatsotn Cucumber -8.4752644 RCP 4.5 2040 to 2060
Ararat Cucumber 0.3734708 RCP 4.5 2040 to 2060
Armavir Cucumber 1.3337670 RCP 4.5 2040 to 2060
Gegharkunik Cucumber -10.8317694 RCP 4.5 2040 to 2060
Kotayk Cucumber -10.1736820 RCP 4.5 2040 to 2060
Lori Cucumber -24.3081504 RCP 4.5 2040 to 2060 
unique(futyield$crop)
## [1] "Cucumber"      "Silage Maize"  "Potato"        "Spring Barley"
## [5] "Tomato"        "Winter Wheat"
unique(futyield$RCP)
## [1] "RCP 4.5" "RCP 8.5"
unique(futyield$period)
## [1] "2040 to 2060" "2080 to 2100"

The file “app.R” of this second part of module 6 contains a shiny app in which which the user can select a certain crop, RCP and future time period, and in which an interactive map shows the selected data in a map. In this scrip, we focus on laying the basis for this app and explore how we can create the interactive map outside a shiny app.

Let us first define some variables that we use to make a selection from “futyield”:

crop_sel   <- "Cucumber"
rcp_sel    <- "RCP 4.5"
period_sel <- "2040 to 2060"

futyield_sel <- filter(futyield, crop == crop_sel, RCP == rcp_sel, period == period_sel)

futyield_sel
region crop value RCP period
Aragatsotn Cucumber -8.4752644 RCP 4.5 2040 to 2060
Ararat Cucumber 0.3734708 RCP 4.5 2040 to 2060
Armavir Cucumber 1.3337670 RCP 4.5 2040 to 2060
Gegharkunik Cucumber -10.8317694 RCP 4.5 2040 to 2060
Kotayk Cucumber -10.1736820 RCP 4.5 2040 to 2060
Lori Cucumber -24.3081504 RCP 4.5 2040 to 2060
Shirak Cucumber -18.6350460 RCP 4.5 2040 to 2060
Syunik Cucumber -6.0147475 RCP 4.5 2040 to 2060
Tavush Cucumber -8.1875939 RCP 4.5 2040 to 2060
Vayotz Dzor Cucumber -4.6316975 RCP 4.5 2040 to 2060
Yerevan Cucumber NA RCP 4.5 2040 to 2060

Join the selected data with the shapefile using the function geo_join() from the tigris package:

marzes_sel <- geo_join(marzes, futyield_sel, "region", "region", how="left")

We could now plot the data in the shapefile in a static manner as we have done several times before:

spplot(marzes_sel, "value", main = "Cucumber yield change under RCP 4.5, 2040-2060")

2. Leaflet

However, we will use a different map plotting package now, which allows zooming, highlighting polygons, etc. - it is called leaflet. There is a very good overview about the leaflet package here: https://rstudio.github.io/leaflet/

Similarly to ggplot(), we sequentially build a leaflet() map by adding layers to each other. However, instead of using the sign “+”, we here have to connect the layers with the operator “%>%”. A simple leaflet map would just consist of some Open Street Map tiles that you can add with the function addTiles(). For other tiles, you have to use the function addProviderTiles(). You can find previews of many available tiles here: https://leaflet-extras.github.io/leaflet-providers/preview/

library(leaflet)

leaflet() %>%
      addTiles()
leaflet() %>%
      addProviderTiles(providers$CartoDB.Positron)

You have seen that by default, we see the entire world. We can use the function setView() to zoom in and set the center coordinate of our window to Armenia - in that way, we decide which part of the world should be seen when the map is loaded:

leaflet() %>%
      setView(45, 40, zoom = 7) %>%
      addTiles()

We can also add different tile sets and use addLayersControl() to create a little window in the upper right corner to select one of them, and add a mini world map in the lower right corner:

leaflet() %>%
      setView(45, 40, zoom = 7) %>%
      addTiles(group = "OSM") %>%
      addProviderTiles(providers$CartoDB.Positron, group = "Positron") %>%
      addLayersControl(baseGroups = c("OSM", "Positron"),
                       options = layersControlOptions(collapsed = FALSE)) %>%
      addMiniMap()

Let us now add the shapefile “marzes_sel” to the map. We use addPolygons() for this. Whilst the argument “fillColor” determines the color of the area of the polygons, the argument “color” determines the color of the polygon borders. “Weight” determines the boldness of these borders. The argument “fillOpacity” determines the level of transparency of the polygon color:

leaflet() %>%
      setView(45, 40, zoom = 7) %>%
      addTiles() %>%
      addPolygons(data = marzes_sel,
                  fillColor = "red", 
                  weight = 1, 
                  color = "black",
                  fillOpacity = 0.6)

Since leaflet maps are interactive, we can define actions that happen when the mouse cursor hovers over a feature in the map. Let us first define some parameters to highlight polygons. The function highlightOptions() is inside addPolygons() but redefines the arguments defined therein for when the mouse cursor passes over a polygon:

leaflet() %>%
      setView(45, 40, zoom = 7) %>%
      addTiles() %>%
      addPolygons(data = marzes_sel,
                  fillColor = "red", 
                  weight = 1, 
                  color = "black",
                  fillOpacity = 0.6,
                  highlight = highlightOptions(weight = 3, color = "white", fillOpacity = 0.9, bringToFront = FALSE))

Let us now add labels to the polygons. We would like to show the name of each marz when we pass over it - for this, we can use the column “region” in “marzes_sel”. The argument “label” would be sufficient to achieve that, but we can further format the labels with “labelOptions”

leaflet() %>%
      setView(45, 40, zoom = 7) %>%
      addTiles() %>%
      addPolygons(data = marzes_sel,
                  fillColor = "red", 
                  weight = 1, 
                  color = "black",
                  fillOpacity = 0.6,
                  highlight = highlightOptions(weight = 3, color = "white", fillOpacity = 0.9, bringToFront = FALSE),
                  label = marzes_sel$region,
                  labelOptions = labelOptions(style = list("font-weight" = "bold", 
                                                              padding = "3px 8px"),
                                                              textsize = "15px", 
                                                              direction = "auto"))

Coloring the polygons according to the values in the column “value” in the attribute table is a bit more complicated in leaflet compared to spplot. To do so, we first have to define the bins, i.e. the intervals that define which values should be displayed with the same color. Let us define 5 equal intervals from the minimum to the maximum value of the data range in the column “value”. The resulting object “bins” is a vector of six elements that determine the interval limits. We also then have to define a color palette function “pal” that is defined using “bins”.

bins <- seq(from = min(na.omit(marzes_sel$value)), 
            to = max(na.omit(marzes_sel$value))+0.001, 
            by = (max(na.omit(marzes_sel$value))+0.001 - min(na.omit(marzes_sel$value)))/5)

bins 
## [1] -24.308150 -19.179567 -14.050983  -8.922400  -3.793816   1.334767
pal <- colorBin(magma(5), domain = na.omit(marzes_sel$value), bins = bins)

We now just have to supply the data column of the shapefile to “pal”, and supply this to the argument “fillColor” in addPolygons(). For each polygon, the palette function will allocate the data value of the polygon to one of the 5 intervals defined by “bins”, and then color the polygon according to color ramp defined in “pal”:

leaflet() %>%
      setView(45, 40, zoom = 7) %>%
      addTiles() %>%
      addPolygons(data = marzes_sel,
                  fillColor = ~pal(marzes_sel$value), 
                  weight = 1, 
                  color = "black",
                  fillOpacity = 0.6,
                  highlight = highlightOptions(weight = 3, color = "white", fillOpacity = 0.9, bringToFront = FALSE),
                  label = marzes_sel$region,
                  labelOptions = labelOptions(style = list("font-weight" = "bold", 
                                                              padding = "3px 8px"),
                                                              textsize = "15px", 
                                                              direction = "auto"))

We can also expand the labels to include the respective value in a second row, in the following way:

labels <- sprintf(paste0("<strong>", marzes_sel$region, "</strong><br/>Change: %s ", "%s"), round(marzes_sel$value, digits=2), "percent") %>% lapply(htmltools::HTML)

leaflet() %>%
      setView(45, 40, zoom = 7) %>%
      addTiles() %>%
      addPolygons(data = marzes_sel,
                  fillColor = ~pal(marzes_sel$value), 
                  weight = 1, 
                  color = "black",
                  fillOpacity = 0.6,
                  highlight = highlightOptions(weight = 3, color = "white", fillOpacity = 0.9, bringToFront = FALSE),
                  label = labels,
                  labelOptions = labelOptions(style = list("font-weight" = "normal", 
                                                              padding = "3px 8px"),
                                                              textsize = "15px", 
                                                              direction = "auto"))

Finally, let us add a legend to the map:

leaflet() %>%
      setView(45, 40, zoom = 7) %>%
      addTiles() %>%
      addPolygons(data = marzes_sel,
                  fillColor = ~pal(marzes_sel$value), 
                  weight = 1, 
                  color = "black",
                  fillOpacity = 0.6,
                  highlight = highlightOptions(weight = 3, color = "white", fillOpacity = 0.9, bringToFront = FALSE),
                  label = labels,
                  labelOptions = labelOptions(style = list("font-weight" = "normal", 
                                                              padding = "3px 8px"),
                                                              textsize = "15px", 
                                                              direction = "auto")) %>%
      addLegend(data = marzes_sel, 
                pal = pal, 
                values = marzes_sel$value, 
                opacity = 0.9,
                title = "Change in %",
                position = "bottomleft")

Congratulations, you made it to the end of the script, and almost to the end of the course! You can now go to the file “app.R” in this folder, and then go on with solving the exercises!