R Geospatial Training, 4 Nov 2020

Intro and "sf"

sf (simple features) objects

sf is a package providing a class system for geographic vector data. - sf can represent all common vector geometry types: points, lines, polygons and their respective ‘multi’ versions.

load packages we will use:

library(sf)

lets set some variables:

# coordinate system for WGS84
WGS84 <- 4326

# Reach colour vars
reach_red <- "#EE5859"
reach_grey <- "#58585A"
reach_lt_grey <- "#D1D3D4"
reach_color_palette <- c('#FEF2F2','#EE5859', '#772B2C')

st_read: loading vector data into R
st_as_sf: casting as sf object in R

st_read can be used for most spatial data types (shp, geojson, kml etc)

# read in shp and define coordinate ref system (4326 = WGS84)
irq_dist <- st_read("data/shp/boundaries/irq_admbnda_adm2_cso_20190603.shp", crs = 4326)
# examine irq_dist
irq_dist

# read in roads shp and define coordinate ref system (4326 = WGS84)
irq_prim_roads <- st_read("data/shp/roads/primary_routes.shp", crs = 4326)

plot sf object

we can use the plot function to display spatial data with sf:

plot(irq_dist)

we can see that the result here renders a new plot for each attribute. None of these other attributes make sense to plot in this manner, as we simply want to plot the geometruy

In R we access individual attributes with the '$' operator. This is also the case when we plot a layer geometry with sf.

# to plot the geometry attribute
plot(irq_dist$geometry)

We pass additional options for a specific function as arguments. to add a layer to the map: use the add = T argument colour (col) and line width (lwd) are 2 other arguments we can add

# plot district layer
plot(irq_dist$geometry, col= reach_lt_grey)
# add roads layer in a brown colour with a width of 1
plot(irq_prim_roads$geometry, add = T, col= "brown", lwd = 1)

We can use conditionals like ifelse to highlight one governorate

# we can highlight one governorate with an ifelse statement for the color. Lets try highlighting Erbil
plot(irq_dist$geometry, col = ifelse(irq_dist$ADM1_EN == "Erbil", reach_red, reach_lt_grey))

Interactive maps with Leaflet

load packages we will use

library(sf)
library(raster)
library(leaflet)
library(dplyr)
library(rmapshaper)

again loading vector data into R:

# read in shp and define coordinate ref system (4326 = WGS84)
irq_dist <- st_read("data/shp/boundaries/irq_admbnda_adm2_cso_20190603.shp")
# define coordinate ref system (crs) and cast as sf object 
irq_dist <- st_as_sf(irq_dist, coords = c("long", "lat"), crs = 4326)

# read in roads shp and define coordinate ref system (4326 = WGS84)
irq_prim_roads <- st_read("data/shp/roads/primary_routes.shp")
# define crs cast as sf object 
irq_prim_roads <- st_as_sf(irq_prim_roads, coords = c("long", "lat"), crs = 4326)

call leaflet function and add pipe (%>%) for each new layer

style property options are similar to web styling and can be seen by entering ?addPolygons to the console. The order that the layers are called in the function is the order they appear on the map. Therefore the layer called last will appear on top.

irq_leaflet_map <-  leaflet() %>% 
  addPolygons(
              data = irq_dist, # select sf dist polygon layer
              color = reach_grey, # stroke color
              fillColor =  reach_lt_grey, 
              weight = 0.3,
              smoothFactor = 0.5,
              opacity = 1.0,
              fillOpacity = 0.8
              ) 

# load map object
irq_leaflet_map

now we want to add governorate boundaries and labels.

# read in gov shp and define coordinate ref system (4326 = WGS84)
irq_gov <- st_read("data/shp/boundaries/irq_admbnda_adm1_cso_20190603.shp")

# convert gov boundaries to lines to display on map so they do not clash with district interactivity
irq_gov <- st_cast(irq_gov,"MULTILINESTRING")
irq_gov <- st_cast(irq_gov,"LINESTRING")

plot map again with governorate boundaries

irq_leaflet_map <-  leaflet() %>% 
  addPolygons(
              data = irq_dist, # select sf dist polygon layer
              color = reach_grey, # stroke color
              fillColor = reach_lt_grey, 
              weight = 0.3,
              smoothFactor = 0.5,
              opacity = 1.0,
              fillOpacity = 0.8
              ) %>%
    addPolylines(
               data = irq_gov,
               color = "black",
               weight= 0.8,
               opacity = 0.8) 

# load map object
irq_leaflet_map

Labelling

create polygon centroids for the governorates in order to create labels

plot map again with governorate labels using the 'addLabelOnlyMarkers' function Lets also try highlighting Erbil by giving it a different color than the rest

irq_leaflet_map <-  leaflet() %>% 
  addPolygons(
              data = irq_dist, # select sf dist polygon layer
              color = ifelse(irq_dist$ADM1_EN == "Erbil", "#white",reach_grey), # stroke color
              fillColor = ifelse(irq_dist$ADM1_EN == "Erbil", reach_red, reach_lt_grey), 
              weight = 0.3,
              smoothFactor = 0.5,
              opacity = 1.0,
              fillOpacity = 0.8
              ) %>%
    addPolylines(
               data = irq_gov,
               color = "black",
               weight= 0.8,
               opacity = 0.8) %>%
    addLabelOnlyMarkers(
               data = irq_gov_cent,
               label= irq_gov_cent$ADM1_EN, # this is the attribute used to determine the labels
               labelOptions = labelOptions(noHide = T,
                            direction = 'center',
                            textOnly = T,
                            style = list(
                              "font-family" = "Arial"
                            ))) 

# load map object
irq_leaflet_map

Point data in R

read in excel data and convert to sf object

library(readxl)
library(dplyr)
library(sf)
library(leaflet)

# read in DTM returnee excel dataset 
dtm_returnee<- read_xlsx("data/table/202092048574_Round117_Master_List_Returnee_2020_08_31_IOM_DTM.xlsx", sheet = "RETURNEE DATASET") %>% 
  filter(!is.na(Latitude), !is.na(Longitude))

# convert to sf object and define coordinate columns in order to plot in leaflet
dtm_returnee_sf <- dtm_returnee %>%
  st_as_sf(
    coords = c("Longitude", "Latitude"), # choose lat and long columns in table 
    crs = WGS84) 

add the points to the leaflet map:

irq_leaflet_map <-  leaflet() %>% 
  addPolygons(
              data = irq_dist, # select sf dist polygon layer
              color = reach_grey, # stroke color
              fillColor = reach_lt_grey, 
              weight = 0.3,
              smoothFactor = 0.5,
              opacity = 1.0,
              fillOpacity = 0.8
              ) %>%
    addPolylines(
               data = irq_gov,
               color = "black",
               weight= 0.8,
               opacity = 0.8) %>%
  # add points to map
  addCircleMarkers( 
               data = dtm_returnee_sf, 
               color = "#301212",
               fillColor = "#EE5859",
               radius=3,
               weight=0.2,
               fillOpacity = 0.8,
               popup = dtm_returnee_sf$`Place id`) %>% 
  # add governorate labels
    addLabelOnlyMarkers(
               data = irq_gov_cent,
               label= irq_gov_cent$ADM1_EN, # this is the attribute used to determine the labels
               labelOptions = labelOptions(noHide = T,
                            direction = 'center',
                            textOnly = T,
                            style = list(
                              "font-family" = "Arial"
                            ))) 

# load map object
irq_leaflet_map

Spatial Joins and choropleth maps

packages

library(dplyr)
library(sf)
library(leaflet)

Spatial joins in R

We will work with the DTM returnee dataset here The district and governorate names are spelled different in the DTM dataset than in the OCHA admin layers. Lets do a spatial join of the points to the district layer so that the attributes from the districts can be carried over to the points.

# spatial join points to get true districts of GPS locations
dtm_returnee_dist_join <- st_join(dtm_returnee_sf, irq_dist )

# now we want to take the join, convert to a table (data frame) and select the columns we want
dtm_returnee_dist_join <- as.data.frame(dtm_returnee_dist_join) %>% 
  dplyr::select(`Place id`, ADM2_EN, Individuals)

# now we need to sum the population by district so we use dplyr functions group by and summarise
dtm_returnee_dist_join <- dtm_returnee_dist_join %>% 
  dplyr::group_by(ADM2_EN) %>% 
  dplyr::summarise(returnee_pop = sum(Individuals))

# now we can join this table to the district layer to have the population in the spatial layer
irq_dist <- left_join(irq_dist, dtm_returnee_dist_join, by = "ADM2_EN")

choropleth maps in R

We will now use the returnee population stats to make s choropleth map in leaflet

# first we will define our color palette in the following format
returnee_pal <- colorNumeric(
  palette = colorRampPalette(c('#FEF2F2','#EE5859', '#772B2C'))(length(irq_dist$returnee_pop)), 
  domain = irq_dist$returnee_pop, na.color = reach_lt_grey)

irq_leaflet_map <-  leaflet() %>% 
  addPolygons(
              data = irq_dist, # select sf dist polygon layer
              color = reach_grey, # stroke color
              fillColor = ~returnee_pal(returnee_pop), 
              popup = irq_dist$ADM2_EN,
              weight = 0.3,
              smoothFactor = 0.5,
              opacity = 1.0,
              fillOpacity = 0.8
              ) %>%
  # add governorate boundaries
    addPolylines(
               data = irq_gov,
               color = "black",
               weight= 0.8,
               opacity = 0.8) %>%
  # add governorate labels
    addLabelOnlyMarkers(
               data = irq_gov_cent,
               label= irq_gov_cent$ADM1_EN, # this is the attribute used to determine the labels
               labelOptions = labelOptions(noHide = T,
                            direction = 'center',
                            textOnly = T,
                            style = list(
                              "font-family" = "Arial"
                            ))) %>%
  # now we add a legend to the bottom left
      addLegend("bottomleft",
                pal = returnee_pal,
                values = irq_dist$returnee_pop,
                title = "Returnee population",
                opacity = 1)

# load map object
irq_leaflet_map

Advanced Leaflet

packages

library(dplyr)
library(sf)
library(leaflet)

Custom labels

More advanced labels can be applied by applying css style elements

irq_leaflet_map <-  leaflet() %>% 
  addPolygons(
              data = irq_dist, # select sf dist polygon layer
              color = ifelse(irq_dist$ADM1_EN == "Erbil", "#white", reach_grey), # stroke color
              fillColor = ifelse(irq_dist$ADM1_EN == "Erbil", reach_red, reach_lt_grey), 
              weight = 0.3,
              smoothFactor = 0.5,
              opacity = 1.0,
              fillOpacity = 0.8
              ) %>%
    addPolylines(
               data = irq_gov,
               color = "black",
               weight= 0.8,
               opacity = 0.8) %>%
    addLabelOnlyMarkers(
               data = irq_gov_cent,
               label= irq_gov_cent$ADM1_EN, 
               labelOptions = labelOptions(noHide = T,
                                            direction = 'center',
                                            textOnly = T,
                                            style = list(
                                              "color" = "#FFF",
                                              "font-weight"= "bold",
                                              "padding" = "3px 8px",
                                              "font-family" = "Arial",
                                              "text-shadow"="2px 2px 5px #58585A"
                                            ))) # here we have changed the options for a nicer label with drop shadow

# load map object
irq_leaflet_map

Adding tooltips and colors on mouse hover

Lets add a tooltip that appears when we hover over the districts. For this we use an HTML

# district tooltip
district_tooltip <- sprintf(
  '<strong><span style="font-size: 8px; color: #EE5859;font-family:Arial;margin-bottom: 1px;">%s District</span><br><span style="font-size: 7px; color: #58585A;"> %s Governorate</strong></span>',
  irq_dist$ADM2_EN, irq_dist$ADM1_EN)%>% 
  lapply(htmltools::HTML)

# now our map again with the label for the tool 
irq_leaflet_map <-  leaflet() %>% 
  addPolygons(
              data = irq_dist, 
              color = ifelse(irq_dist$ADM1_EN == "Erbil", "white", reach_grey), 
              fillColor = ifelse(irq_dist$ADM1_EN == "Erbil", reach_red, reach_lt_grey), 
              weight = 0.3,
              smoothFactor = 0.5,
              opacity = 1.0,
              fillOpacity = 0.8,
              label = district_tooltip, # here we link to the tooltip mark up described above
              highlightOptions = highlightOptions(fillColor="white",
                                   fillOpacity = 0.5,
                                   color="white",
                                   weight = 1.4,
                                   bringToFront = F) # here we add a white color for the distr hover
              ) %>%
    addPolylines(
               data = irq_gov,
               color = "black",
               weight= 0.8,
               opacity = 0.8) %>%
    addLabelOnlyMarkers(
               data = irq_gov_cent,
               label= irq_gov_cent$ADM1_EN, 
               labelOptions = labelOptions(noHide = T,
                                            direction = 'center',
                                            textOnly = T,
                                            style = list(
                                              "color" = "#FFF",
                                              "font-weight"= "bold",
                                              "padding" = "3px 8px",
                                              "font-family" = "Arial",
                                              "text-shadow"="2px 2px 5px #58585A"
                                            )))

# load map object
irq_leaflet_map

Adding toggle between arabic and english names

For labels and tooltips in both languages we use layer groups

# district tooltip - English - same as previously
district_tooltip <- sprintf(
  '<strong><span style="font-size: 8px; color: #EE5859;font-family:Arial;margin-bottom: 1px;">%s District</span><br><span style="font-size: 7px; color: #58585A;"> %s Governorate</strong></span>',
  irq_dist$ADM2_EN, irq_dist$ADM1_EN)%>% 
  lapply(htmltools::HTML)

# district tooltip - Arabic 
district_tooltip_arabic <- sprintf(
  '<strong><span style="font-size: 10px; color: #EE5859;font-family:Arial;margin-bottom: 1px;">%s District</span><br><span style="font-size: 9px; color: #58585A;"> %s Governorate</strong></span>',
  irq_dist$ADM2_AR, irq_dist$ADM1_AR)%>%  # here we change the attributes to the arabic columns
  lapply(htmltools::HTML)

# now our map again with the label for the tool 
irq_leaflet_map <-  leaflet() %>% 
  addPolygons( # english
              data = irq_dist, 
              group="English", # here we specify the group
              color = ifelse(irq_dist$ADM1_EN == "Erbil", "white", reach_grey), 
              fillColor = ifelse(irq_dist$ADM1_EN == "Erbil", reach_red, reach_lt_grey), 
              weight = 0.3,
              smoothFactor = 0.5,
              opacity = 1.0,
              fillOpacity = 0.8,
              label = district_tooltip, # english tooltip - same as before
              highlightOptions = highlightOptions(fillColor="white",
                                   fillOpacity = 0.5,
                                   color="white",
                                   weight = 1.4,
                                   bringToFront = F) 
              ) %>%
    addPolygons( # arabic
              data = irq_dist, 
              group="Arabic",
              color = ifelse(irq_dist$ADM1_EN == "Erbil", "white", reach_grey), 
              fillColor = ifelse(irq_dist$ADM1_EN == "Erbil", reach_red, reach_lt_grey), 
              weight = 0.3,
              smoothFactor = 0.5,
              opacity = 1.0,
              fillOpacity = 0.8,
              label = district_tooltip_arabic, # here we link to the arabic tooltip
              highlightOptions = highlightOptions(fillColor="white",
                                   fillOpacity = 0.5,
                                   color="white",
                                   weight = 1.4,
                                   bringToFront = F) # here we add a white color for the distr hover
              ) %>%
    addPolylines(
               data = irq_gov,
               group = "English",
               color = "black",
               weight= 0.8,
               opacity = 0.8) %>%
    addPolylines(
               data = irq_gov,
               group = "Arabic",
               color = "black",
               weight= 0.8,
               opacity = 0.8) %>%
    addLabelOnlyMarkers(
               data = irq_gov_cent,
               group = "English",
               label= irq_gov_cent$ADM1_EN, # english labels
               labelOptions = labelOptions(noHide = T,
                                            direction = 'center',
                                            textOnly = T,
                                            style = list(
                                              "color" = "#FFF",
                                              "font-weight"= "bold",
                                              "padding" = "3px 8px",
                                              "font-family" = "Arial",
                                              "text-shadow"="2px 2px 5px #58585A"
                                            ))) %>% 
      addLabelOnlyMarkers(
               data = irq_gov_cent,
               group = "Arabic",
               label= irq_gov_cent$ADM1_AR, # arabic labels
               labelOptions = labelOptions(noHide = T,
                                            direction = 'center',
                                            textOnly = T,
                                            style = list(
                                              "color" = "#FFF",
                                              "font-weight"= "bold",
                                              "padding" = "3px 8px",
                                              "font-family" = "Arial",
                                              "text-shadow"="2px 2px 5px #58585A"
                                            ))) %>% 
  # we add a layers control window to allow for the toggling
            addLayersControl(
            baseGroups = c("English","Arabic"),
            options = layersControlOptions(collapsed = F)
          )

# load map object
irq_leaflet_map