#to avoid messages of not recognising a file by R
load('E:/Documents/R Studio is here/R Studio files projects/gis_coding/.RData')

All along we’ve been dealing with creating print maps showing all manner of acrobatics to prove we are a cool cartographer or a special spatial analyst. However, being cool or special is not worth the salt if we cannot create a web map to represent the result of our operations. Therefore, to prove ourselves, we will create a webmap to present a few findings. Of course your personal worth is not pegged on your successes.

#load leaflet. Leaflet is a language that enables plotting of interactive maps
library(leaflet)
## Warning: package 'leaflet' was built under R version 4.1.3
#load this other packages necessary for R spatial analysis
library(rgdal)
## Loading required package: sp
## Please note that rgdal will be retired by the end of 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## 
## rgdal: version: 1.5-28, (SVN revision 1158)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.2.1, released 2020/12/29
## Path to GDAL shared files: C:/Users/User/Documents/R/win-library/4.1/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: C:/Users/User/Documents/R/win-library/4.1/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.4-6
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
## Overwritten PROJ_LIB was C:/Users/User/Documents/R/win-library/4.1/rgdal/proj
library(rgeos)
## rgeos version: 0.5-9, (SVN revision 684)
##  GEOS runtime version: 3.9.1-CAPI-1.14.2 
##  Please note that rgeos will be retired by the end of 2023,
## plan transition to sf functions using GEOS at your earliest convenience.
##  GEOS using OverlayNG
##  Linking to sp version: 1.4-6 
##  Polygon checking: TRUE
library(sf)
## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1; sf_use_s2() is TRUE
library(sp)

For this exercise, we will use a shapefile containing Kenyan towns. To save on time, this particular link contains a list of all major Kenyan towns georeferenced (http://datasets.wri.org/dataset/27f2f740-4fdb-474d-9c77-b78072beb807/resource/3c780407-0046-4fc8-988a-247317ae9260/download/kemajor-towns.zip) However, these towns in the shapefile lack the lat-long coordinates that are needed by leaflet for a couple of operations, say when adding Markers. You will see what I mean later. Therefore, using Qgis, the towns outside the top five populous counties were clipped out, and the lat-long coordinates for each town were auto-assigned using the Qgis ‘Add X/Y fields to layer’. This modified file to be used for the rest of the analysis is available from Dropbox here: https://www.dropbox.com/s/uayzj6ybi8cisgp/kenya_alltowns.zip?dl=0

#call the shapefile into R
towns <- readOGR('E:/Documents/R Studio is here/R Studio files projects/gis_coding/kenya_alltowns.shp')
## OGR data source with driver: ESRI Shapefile 
## Source: "E:\Documents\R Studio is here\R Studio files projects\gis_coding\kenya_alltowns.shp", layer: "kenya_alltowns"
## with 13 features
## It has 7 fields
## Integer64 fields read as strings:  TOWN_ID
#check the attributes of the towns shapefile
head(towns@data)
##   AREA PERIMETER TOWN_NAME TOWN_ID    TOWN_TYPE        x           y
## 1    0         0    Kitale       1 Municipality 35.00683  1.02031000
## 2    0         0    Kisumu       1 Municipality 34.76333 -0.09653389
## 3    0         0   Kericho       1 Municipality 35.28588 -0.36461039
## 4    0         0  Kakamega       1 Municipality 34.75269  0.29018246
## 5    0         0   Eldoret       1 Municipality 35.26689  0.52684920
## 6    0         0   Bungoma       1 Municipality 34.63162  0.54046363
#plot the towns overlying the default OSM m
leaflet() %>%
  addTiles() %>% #this adds the OSM layer
  fitBounds(towns, lng1 = 34.6316215090692410, lat1 = -4.0564121656444732,
                                                         lng2 =  40.1165678143107414, lat2 = 1.0203100000000001) %>% #we fit the OSM to the bounds of our point shapefile
  addCircles(data = towns, ~x,  ~y) # the ~x and ~y stand for longitude and latitude arguments for the addCircles function.

Your webmaps appear on the Viewer tab as shown in the image below.

viewer tab

However, we can make our map more interactive. How about we make the towns clickable, and a small textbox of the name of the town shows up? Therefore, we can’t be lost, literally.

#make the towns info ie. name appear when mouse clicks on them
leaflet(towns) %>%
  addTiles() %>% #this adds the OSM layer
  fitBounds(towns, lng1 = 34.6316215090692410, lat1 = -4.0564121656444732,
                                                         lng2 =  40.1165678143107414, lat2 = 1.0203100000000001) %>% #we fit the OSM to the bounds of our point shapefile
  addCircles( ~x,  ~y, radius = 15, popup = towns@data$TOWN_NAME, options = popupOptions(maxWidth = 15,  textsize = 3, closeButton = F, closeOnClick = F)) # the ~x and ~y stand for longitude and latitude arguments for the addCircles function.

Now what just happened here? Let’s start from the very beginning. The leaflet function creates a map widget in which your map will be drawn. It is what activates the Viewer tab on the bottom right of your R window. addTiles adds the OSM layer, fitBounds specifies the maximum spatial extent that should be displayed, though its not always necessary. addCircles adds the graphics we would like to see on the map. Because we had already specified the towns data (which is the shapefile containing all our data), we need not include it in the addCircles(data = towns). This is because it is already specified in the leaflet function, our parent who enable us to create this webpages in the first place. You can think of leaflet as the magic word to open our portal to webmapping, figuratively speaking. If we had simply put in leaflet()%>% without anything in the brackets, addCircles would need the data = specified.

In the second map which contains our clickable towns, we specified our popup parameters. popup specifies the field that should be displayed on click, while popupOptions goes a step further by specifying the textbox parameters such as the size of our popup rectangle and textsize.

Looks beautiful, but it was like chiseling out diamond from rock with some tender…loving…care.

Why not use the addMarkers? This would have needed one to specify the icon. You would have to call the icon from your directory, or from the web, specify the width, height and other parameters, all too much of a hustle when there are other less tedious options available. That’s why we settled for addCircles. If interested in addMarkers, check out from the creator of this leaflet tool here: https://rstudio.github.io/leaflet/markers.html

In the following map, we will try to differentiate the circles by a variable. We will also use the label argument this time round because it gives us the option to find the name of a town simply by hovering over it.

#this is to create a range of colors differentiated by the TOWN_NAME. Any other variable can be used but TOWN_NAMES has been used because it is more varied than TOWN_TYPE which only has two. 
mapcolors <- colorFactor(palette = topo.colors(13, 0.7, rev = F), domain = towns@data$TOWN_NAME) 
##other better parameters to differentiate color would have been like year of establishment, population etc. but are not available in this dataset

#create the map
leaflet(towns) %>%
  addTiles() %>% #this adds the OSM layer
  fitBounds(towns, lng1 = 34.6316215090692410, lat1 = -4.0564121656444732,
                                                         lng2 =  40.1165678143107414, lat2 = 1.0203100000000001) %>% #we fit the OSM to the bounds of our point shapefile
  addCircles( ~x,  ~y, radius = 15, color = ~mapcolors(TOWN_NAME),  label = towns@data$TOWN_NAME, options = labelOptions(interactive = T, clickable = T, noHide = T, textsize = 9))

In the third web map, try to settle on a town without clicking. I know it’s tricky, but when you see the mouse pointer shift to a hand you’ll know leaflet has registered that as the point of interest. Let the mouse hand rest there briefly, without clicking on the mouse. The name of that town will appear all by itself.

Next, we will proceed to adding the polygons of our county_hs_assets shapefile. For one, we just want to experiment adding polygons to a webmap. Secondly having something other than plain name towns gives the webmap some educative edge. The county_hs_assetts polygon data will underlay the town points so as not to obstruct them.

#check attributes of county_hs_assets
summary(county_hs_assets@data)
##    ADM1_EN            Shape_Leng       Shape_Area      ADM1_PCODE       
##  Length:6           Min.   : 1.658   Min.   :0.0575   Length:6          
##  Class :character   1st Qu.: 3.095   1st Qu.:0.2066   Class :character  
##  Mode  :character   Median : 3.193   Median :0.2255   Mode  :character  
##                     Mean   : 4.305   Mean   :0.3740                     
##                     3rd Qu.: 3.892   3rd Qu.:0.2449                     
##                     Max.   :10.599   Max.   :1.2836                     
##    ADM1_REF          ADM1ALT1EN         ADM1ALT2EN          ADM0_EN         
##  Length:6           Length:6           Length:6           Length:6          
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##   ADM0_PCODE            date             validOn            validTo         
##  Length:6           Length:6           Length:6           Length:6          
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##  County...Sub.County Conventional.Households Stand...alone.Radio
##  Length:6            Length:6                Min.   :53.40      
##  Class :character    Class :character        1st Qu.:61.12      
##  Mode  :character    Mode  :character        Median :61.30      
##                                              Mean   :60.20      
##                                              3rd Qu.:61.85      
##                                              Max.   :62.10      
##  Desk.Top.Computer..Laptop..Tablet Functional.Television Analogue.Television
##  Min.   : 4.100                    Min.   :25.40         Min.   :3.600      
##  1st Qu.: 5.625                    1st Qu.:35.02         1st Qu.:4.650      
##  Median :14.600                    Median :60.60         Median :6.000      
##  Mean   :14.783                    Mean   :53.27         Mean   :5.733      
##  3rd Qu.:21.925                    3rd Qu.:69.67         3rd Qu.:6.825      
##  Max.   :28.400                    Max.   :73.80         Max.   :7.500      
##     Internet        Bicycle       Motor.Cycle      Refrigerator  
##  Min.   : 7.20   Min.   :12.50   Min.   : 4.300   Min.   : 2.80  
##  1st Qu.:11.68   1st Qu.:17.00   1st Qu.: 5.475   1st Qu.: 4.70  
##  Median :30.00   Median :19.15   Median : 8.900   Median :14.25  
##  Mean   :27.93   Mean   :19.62   Mean   : 8.417   Mean   :15.05  
##  3rd Qu.:41.35   3rd Qu.:23.10   3rd Qu.:11.350   3rd Qu.:22.52  
##  Max.   :49.70   Max.   :26.20   Max.   :11.900   Max.   :32.20  
##       Car         Truck..Lorry..Bus..Three.Wheeler.truck    Tuk.Tuk      
##  Min.   : 3.100   Min.   :0.600                          Min.   :0.4000  
##  1st Qu.: 4.250   1st Qu.:0.700                          1st Qu.:0.4000  
##  Median : 9.750   Median :1.100                          Median :0.4000  
##  Mean   : 9.983   Mean   :1.200                          Mean   :0.4333  
##  3rd Qu.:12.775   3rd Qu.:1.425                          3rd Qu.:0.4750  
##  Max.   :21.100   Max.   :2.300                          Max.   :0.5000  
##     names          
##  Length:6          
##  Class :character  
##  Mode  :character  
##                    
##                    
## 
#remember color differentiation for towns was - mapcolors <- colorFactor(palette = topo.colors(13, 0.7, rev = F), domain = towns@data$TOWN_NAME)
county_color <- colorFactor(hcl.colors(6, palette = 'viridis', 0.9, rev = F, fixup = T), 
                            domain = county_hs_assets@data$names) #create a palette to differentiate the values of a particular shapefile field.

leaflet() %>%
  addTiles() %>% #this adds the OSM layer
  fitBounds(towns, lng1 = 34.6316215090692410, lat1 = -4.0564121656444732,
                                                         lng2 =  40.1165678143107414, lat2 = 1.0203100000000001) %>% #we fit the OSM to the bounds of our point shapefile
  addPolygons(data = county_hs_assets, color = ~county_color(names), popup = paste('County: ', county_hs_assets@data$County...Sub.County, '<br>', 'No. Households: ', county_hs_assets@data$Conventional.Households),   highlightOptions = highlightOptions(stroke = T, color = 'white', bringToFront = T)) %>% #popup indicates that when we click on the polygons, county names and no. of households are displayed
  addCircleMarkers(data = towns, ~x,  ~y, radius = 10, color =  ~mapcolors(TOWN_NAME), fillColor = ~mapcolors(TOWN_NAME), fill = T, fillOpacity = 1, stroke = F, label = towns@data$TOWN_NAME, options = labelOptions(interactive = T, clickable = T, noHide = T, textsize = 9))

You may wonder where somethings like paste came from. The idea here was to display multiple values on the popup for the county polygons layer i.e ‘county: name xxx’ a break and ‘population: number yyy’. Initially I had tried popup = c(...$names, ...$conventionalhouseholds) but R only displayed the values of the enumerated households as though the call function was persona nongrata (aka unrecognisable) this time round. This link here saved me a ton - https://stackoverflow.com/questions/31562383/using-leaflet-library-to-output-multiple-popup-values.

The highlightOptions was used to specify how the polygons should appear when hovered. This makes them distinct from the rest of the map, including the town points-when hovered over. This helps one know which layer is active especially where there are two or more at the same point.

Our map cannot make it to the hall of fame without a legend. The addLegend function is the handy here. In the addLegend function, we use pal to call the color palette object used in displaying a particular layer, say ‘mapcolors’ for towns. Thereafter, for addLegend to associate which color palette with which layer, the pal argument is used to specify the color palette. You don’t want your legend to dispaly a layer matched with the wrong colors. Now that pal understands which color palette is hooked to which layer, the values argument is used to specify which layer variable was used to create the color palette. Does it sound like gibberish? (yes) In very simple terms, place your customized color palette to pal and your domain to the values argument. This sounds better. (hope it did)

#remember color differentiation for towns was - mapcolors <- colorFactor(palette = topo.colors(13, 0.7, rev = F), domain = towns@data$TOWN_NAME)
#remember for counties - county_color <- colorFactor(hcl.colors(6, palette = 'viridis', 0.9, rev = F, fixup = T),  domain = county_hs_assets@data$names)
#####

leaflet() %>%
  addTiles() %>% #this adds the OSM layer
  fitBounds(towns, lng1 = 34.6316215090692410, lat1 = -4.0564121656444732,
                                                         lng2 =  40.1165678143107414, lat2 = 1.0203100000000001) %>% #we fit the OSM to the bounds of our point shapefile
  addPolygons(data = county_hs_assets, color = ~county_color(names), popup = paste('County: ', county_hs_assets@data$County...Sub.County, '<br>', 'No. Households: ', county_hs_assets@data$Conventional.Households),   highlightOptions = highlightOptions(stroke = T, color = 'white', bringToFront = T)) %>% #popup indicates that when we click on the polygons, countynames are displayed
  addCircleMarkers(data = towns, ~x,  ~y, radius = 10, color =  ~mapcolors(TOWN_NAME), fillColor = ~mapcolors(TOWN_NAME), fill = T, fillOpacity = 1, stroke = F, label = towns@data$TOWN_NAME, options = labelOptions(interactive = T, clickable = T, noHide = T, textsize = 9)) %>%
  addLegend(position = 'bottomright', pal = mapcolors, values = towns@data$TOWN_NAME, title = 'Major towns') %>%
  addLegend(position = 'bottomleft', pal = county_color, values = county_hs_assets@data$names, title = 'Counties')

There is one more decoration to make this work legendary - the toggle layers on/off widget.

#lets shorten our long codes into a single variable called 'mymap'. Note that for the addPolygons and addCircles, group argument has been added. 
mymap <- leaflet() %>%
  addTiles(group = 'OSM') %>% #this adds the OSM layer
  fitBounds(towns, lng1 = 34.6316215090692410, lat1 = -4.0564121656444732,
                                                         lng2 =  40.1165678143107414, lat2 = 1.0203100000000001) %>% #we fit the OSM to the bounds of our point shapefile
  addPolygons(data = county_hs_assets, color = ~county_color(names), popup = paste('County: ', county_hs_assets@data$County...Sub.County, '<br>', 'No. Households: ', county_hs_assets@data$Conventional.Households),   highlightOptions = highlightOptions(stroke = T, color = 'white', bringToFront = T), group = 'Counties') %>% #popup indicates that when we click on the polygons, countynames are displayed
  #group name has been specified to indicate which group the created layer belongs to
  addCircleMarkers(data = towns, ~x,  ~y, radius = 10, color =  ~mapcolors(TOWN_NAME), fillColor = ~mapcolors(TOWN_NAME), fill = T, fillOpacity = 1, stroke = F, label = towns@data$TOWN_NAME, options = labelOptions(interactive = T, clickable = T, noHide = T, textsize = 9), group = 'Towns') %>%
  addLegend(position = 'bottomright', pal = mapcolors, values = towns@data$TOWN_NAME, title = 'Major towns') %>%
  addLegend(position = 'bottomleft', pal = county_color, values = county_hs_assets@data$names, title = 'Counties')

mymap

Although when saving the webmap code to an object made the entire work shorted to a single word, adding more decorations was becoming problematic. For example, adding the addLayerOptions to mymap variable was simply giving me a static map with interactive control widgets. Seems you have to add it to the long code. But you get the drift, you can shorten long codes into a word and add more functionalities to the object ie. `mymap + addPolygons()…’. Try this:

mymap + addLayersControl(baseGroups = 'OSM', overlayGroups = c('Counties', 'Towns'), position = 'topright', options = layersControlOptions(collapsed = T))

In my case I got the following error:

Error in getMapData(map) : argument "map" is missing, with no default
#draw my map
mymap

Lets use the long route anywhere. At least the end is in sight.

#use addLayerOptions to create a layer widget. This is where the groups come in handy. Group names will be used to store the info of a particular element in the layers widget
leaflet() %>%
  addTiles(group = 'OSM') %>% #this adds the OSM layer
  fitBounds(towns, lng1 = 34.6316215090692410, lat1 = -4.0564121656444732,
                                                         lng2 =  40.1165678143107414, lat2 = 1.0203100000000001) %>% #we fit the OSM to the bounds of our point shapefile
  addPolygons(data = county_hs_assets, color = ~county_color(names), popup = paste('County: ', county_hs_assets@data$County...Sub.County, '<br>', 'No. Households: ', county_hs_assets@data$Conventional.Households),   highlightOptions = highlightOptions(stroke = T, color = 'white', bringToFront = T), group = 'Counties') %>% #popup indicates that when we click on the polygons, countynames are displayed
  addCircleMarkers(data = towns, ~x,  ~y, radius = 10, color =  ~mapcolors(TOWN_NAME), fillColor = ~mapcolors(TOWN_NAME), fill = T, fillOpacity = 1, stroke = F, label = towns@data$TOWN_NAME, options = labelOptions(interactive = T, clickable = T, noHide = T, textsize = 9)) %>%
  addLegend(position = 'bottomright', pal = mapcolors, values = towns@data$TOWN_NAME, title = 'Major towns', group = 'Towns') %>%
  addLegend(position = 'bottomleft', pal = county_color, values = county_hs_assets@data$names, title = 'Counties') %>%
  addLayersControl(baseGroups = 'OSM', overlayGroups = c('Counties', 'Towns'), position = 'topright', options = layersControlOptions(collapsed = T))

Open StreetMap (OSM) is not the only tile layer callable in leaflet. There are others from ESRI, USGS et cetera and the leaflet creator was very generous in putting them all in one basket. Check them out: http://leaflet-extras.github.io/leaflet-providers/preview/index.html.

#new tiles will be added, passed to their respective groups and update the layers widget
leaflet() %>%
  addTiles(group = 'OSM') %>% #this adds the OSM layer
  fitBounds(towns, lng1 = 34.6316215090692410, lat1 = -4.0564121656444732,
                                                         lng2 =  40.1165678143107414, lat2 = 1.0203100000000001) %>% #we fit the OSM to the bounds of our point shapefile
  addTiles(urlTemplate = 'https://server.arcgisonline.com/ArcGIS/rest/services/World_Street_Map/MapServer/tile/{z}/{y}/{x}', group = 'ESRI')%>% #add the ESRI map tile 
  addTiles(urlTemplate = 'https://basemap.nationalmap.gov/arcgis/rest/services/USGSImageryOnly/MapServer/tile/{z}/{y}/{x}', group = 'USGS') %>% #add the USGS map tile
  addPolygons(data = county_hs_assets, color = ~county_color(names), popup = paste('County: ', county_hs_assets@data$County...Sub.County, '<br>', 'No. Households: ', county_hs_assets@data$Conventional.Households),   highlightOptions = highlightOptions(stroke = T, color = 'white', bringToFront = T), group = 'Counties') %>% #popup indicates that when we click on the polygons, countynames are displayed
  addCircleMarkers(data = towns, ~x,  ~y, radius = 10, color =  ~mapcolors(TOWN_NAME), fillColor = ~mapcolors(TOWN_NAME), fill = T, fillOpacity = 1, stroke = F, label = towns@data$TOWN_NAME, options = labelOptions(interactive = T, clickable = T, noHide = T, textsize = 9)) %>%
  addLegend(position = 'bottomright', pal = mapcolors, values = towns@data$TOWN_NAME, title = 'Major towns', group = 'Towns') %>%
  addLegend(position = 'bottomleft', pal = county_color, values = county_hs_assets@data$names, title = 'Counties') %>%
  addLayersControl(baseGroups = c('OSM', 'USGS', 'ESRI'),  overlayGroups = c('Counties', 'Towns'), position = 'topright', options = layersControlOptions(collapsed = T))

This tutorial has been made possible by standing on the shoulders of giants. As an extra exercise, what if I want the circle markers differentiated by size? To do this, one has to choose a reasonable numericable variable. In the case our towns file, no such reasonable variable exists. However, we can make honey from a beehive by getting the length of our town names multiplied by 1.2. Since the names of our towns under the TOWN_NAME variable all vary slightly, multiplying by 1.2 will exaggerate the total value of their characters so as to get distinct sizes for each town.

#make the sizes of the town points dependent on length of town name x1.2
leaflet() %>%
  addTiles(group = 'OSM') %>% #this adds the OSM layer
  fitBounds(towns, lng1 = 34.6316215090692410, lat1 = -4.0564121656444732,
                                                         lng2 =  40.1165678143107414, lat2 = 1.0203100000000001) %>% #we fit the OSM to the bounds of our point shapefile
  addTiles(urlTemplate = 'https://server.arcgisonline.com/ArcGIS/rest/services/World_Street_Map/MapServer/tile/{z}/{y}/{x}', group = 'ESRI')%>% #add the ESRI map tile 
  addTiles(urlTemplate = 'https://basemap.nationalmap.gov/arcgis/rest/services/USGSImageryOnly/MapServer/tile/{z}/{y}/{x}', group = 'USGS') %>% #add the USGS map tile
  addPolygons(data = county_hs_assets, color = ~county_color(names), popup = paste('County: ', county_hs_assets@data$County...Sub.County, '<br>', 'No. Households: ', county_hs_assets@data$Conventional.Households),   highlightOptions = highlightOptions(stroke = T, color = 'white', bringToFront = T), group = 'Counties') %>% #popup indicates that when we click on the polygons, countynames are displayed
  addCircleMarkers(data = towns, ~x,  ~y, radius = ~nchar(towns@data$TOWN_NAME)*1.2, #nchar() is a function to check for length of character values in a string
                   color =  ~mapcolors(TOWN_NAME), fillColor = ~mapcolors(TOWN_NAME), fill = T, fillOpacity = 1, stroke = F, label = towns@data$TOWN_NAME, options = labelOptions(interactive = T, clickable = T, noHide = T, textsize = 9), group = 'Towns') %>%
  addLegend(position = 'bottomright', pal = mapcolors, values = towns@data$TOWN_NAME, title = 'Major towns') %>%
  addLegend(position = 'bottomleft', pal = county_color, values = county_hs_assets@data$names, title = 'Counties') %>%
  addLayersControl(baseGroups = c('OSM', 'USGS', 'ESRI'),  overlayGroups = c('Counties', 'Towns'), position = 'topright', options = layersControlOptions(collapsed = T))

A clean looking map.