This tutorial covers how to create interactive maps in R, focusing on two types of visualizations:
There are many ways to create maps using R, but we will focus on leaflet, a popular open-source JavaScript library that is used by many professional organizations to create interactive maps.
In this tutorial we will require the following packages:
library(leaflet) # The map-making package
library(geojsonio) # A package for geographic and spatial data, requires the latest version of dplyr
library(dplyr) # Used for data manipulation and merging
library(htmltools) # Used for constructing map labels using HTMLIn this tutorial we will use two datasets:
data(quakes)
head(quakes)## lat long depth mag stations
## 1 -20.42 181.62 562 4.8 41
## 2 -20.62 181.03 650 4.2 15
## 3 -26.00 184.10 42 5.4 43
## 4 -17.97 181.66 626 4.1 19
## 5 -20.42 181.96 649 4.0 11
## 6 -19.68 184.31 195 4.0 12happiness <- read.csv("https://raw.githubusercontent.com/skuiper/datascience/master/datasets/HappyPlanet.csv", stringsAsFactors = FALSE)
head(happiness)## Country Region Happiness LifeExpectancy Footprint HLY HPI HPIRank
## 1 Albania 7 5.5 76.2 2.2 41.7 47.91 54
## 2 Algeria 3 5.6 71.7 1.7 40.1 51.23 40
## 3 Angola 4 4.3 41.7 0.9 17.8 26.78 130
## 4 Argentina 1 7.1 74.8 2.5 53.4 58.95 15
## 5 Armenia 7 5.0 71.7 1.4 36.1 48.28 48
## 6 Australia 2 7.9 80.9 7.8 63.7 36.64 102
## GDPperCapita HDI Population
## 1 5316 0.801 3.15
## 2 7062 0.733 32.85
## 3 2335 0.446 16.10
## 4 14280 0.869 38.75
## 5 4945 0.775 3.02
## 6 31794 0.962 20.40Generally speaking, there are three steps to creating a map using leaflet:
leaflet functionaddTiles, addMarkers, and addPolygonsLayers are usually added by “piping” via the %>% operator. Put simply, this means passing the resulting object on the left of the %>% symbol as the first argument to the function on the right of the %>%, thereby chaining together a series of functions that each operate on the results of the prior step.
## Step 1 call leaflet
map1 <- leaflet() %>%
addTiles() ## Step 2 add a layer using %>%
## Step 3 print the map
map1Leaflet constructs a base map using map tiles, an approach popularized by google maps that is now used by nearly every interactive map online.
Calling addTiles() with no arguments will generate a base map using OpenStreetMap tiles. For most purposes OpenStreetMap tiles are all you’ll need; however, the leaflet package also comes with dozens of different tile options created by third party providers and stored in an object named providers. Provider tile sets can be used to create a base map via the addProviderTiles function:
names(providers) ## Providers contains a list of different tiles## [1] "OpenStreetMap"
## [2] "OpenStreetMap.Mapnik"
## [3] "OpenStreetMap.DE"
## [4] "OpenStreetMap.CH"
## [5] "OpenStreetMap.France"
## [6] "OpenStreetMap.HOT"
## [7] "OpenStreetMap.BZH"
## [8] "OpenSeaMap"
## [9] "OpenPtMap"
## [10] "OpenTopoMap"
## [11] "OpenRailwayMap"
## [12] "OpenFireMap"
## [13] "SafeCast"
## [14] "Thunderforest"
## [15] "Thunderforest.OpenCycleMap"
## [16] "Thunderforest.Transport"
## [17] "Thunderforest.TransportDark"
## [18] "Thunderforest.SpinalMap"
## [19] "Thunderforest.Landscape"
## [20] "Thunderforest.Outdoors"
## [21] "Thunderforest.Pioneer"
## [22] "Thunderforest.MobileAtlas"
## [23] "Thunderforest.Neighbourhood"
## [24] "OpenMapSurfer"
## [25] "OpenMapSurfer.Roads"
## [26] "OpenMapSurfer.Hybrid"
## [27] "OpenMapSurfer.AdminBounds"
## [28] "OpenMapSurfer.ContourLines"
## [29] "OpenMapSurfer.Hillshade"
## [30] "OpenMapSurfer.ElementsAtRisk"
## [31] "Hydda"
## [32] "Hydda.Full"
## [33] "Hydda.Base"
## [34] "Hydda.RoadsAndLabels"
## [35] "MapBox"
## [36] "Stamen"
## [37] "Stamen.Toner"
## [38] "Stamen.TonerBackground"
## [39] "Stamen.TonerHybrid"
## [40] "Stamen.TonerLines"
## [41] "Stamen.TonerLabels"
## [42] "Stamen.TonerLite"
## [43] "Stamen.Watercolor"
## [44] "Stamen.Terrain"
## [45] "Stamen.TerrainBackground"
## [46] "Stamen.TerrainLabels"
## [47] "Stamen.TopOSMRelief"
## [48] "Stamen.TopOSMFeatures"
## [49] "TomTom"
## [50] "TomTom.Basic"
## [51] "TomTom.Hybrid"
## [52] "TomTom.Labels"
## [53] "Esri"
## [54] "Esri.WorldStreetMap"
## [55] "Esri.DeLorme"
## [56] "Esri.WorldTopoMap"
## [57] "Esri.WorldImagery"
## [58] "Esri.WorldTerrain"
## [59] "Esri.WorldShadedRelief"
## [60] "Esri.WorldPhysical"
## [61] "Esri.OceanBasemap"
## [62] "Esri.NatGeoWorldMap"
## [63] "Esri.WorldGrayCanvas"
## [64] "OpenWeatherMap"
## [65] "OpenWeatherMap.Clouds"
## [66] "OpenWeatherMap.CloudsClassic"
## [67] "OpenWeatherMap.Precipitation"
## [68] "OpenWeatherMap.PrecipitationClassic"
## [69] "OpenWeatherMap.Rain"
## [70] "OpenWeatherMap.RainClassic"
## [71] "OpenWeatherMap.Pressure"
## [72] "OpenWeatherMap.PressureContour"
## [73] "OpenWeatherMap.Wind"
## [74] "OpenWeatherMap.Temperature"
## [75] "OpenWeatherMap.Snow"
## [76] "HERE"
## [77] "HERE.normalDay"
## [78] "HERE.normalDayCustom"
## [79] "HERE.normalDayGrey"
## [80] "HERE.normalDayMobile"
## [81] "HERE.normalDayGreyMobile"
## [82] "HERE.normalDayTransit"
## [83] "HERE.normalDayTransitMobile"
## [84] "HERE.normalDayTraffic"
## [85] "HERE.normalNight"
## [86] "HERE.normalNightMobile"
## [87] "HERE.normalNightGrey"
## [88] "HERE.normalNightGreyMobile"
## [89] "HERE.normalNightTransit"
## [90] "HERE.normalNightTransitMobile"
## [91] "HERE.reducedDay"
## [92] "HERE.reducedNight"
## [93] "HERE.basicMap"
## [94] "HERE.mapLabels"
## [95] "HERE.trafficFlow"
## [96] "HERE.carnavDayGrey"
## [97] "HERE.hybridDay"
## [98] "HERE.hybridDayMobile"
## [99] "HERE.hybridDayTransit"
## [100] "HERE.hybridDayGrey"
## [101] "HERE.hybridDayTraffic"
## [102] "HERE.pedestrianDay"
## [103] "HERE.pedestrianNight"
## [104] "HERE.satelliteDay"
## [105] "HERE.terrainDay"
## [106] "HERE.terrainDayMobile"
## [107] "FreeMapSK"
## [108] "MtbMap"
## [109] "CartoDB"
## [110] "CartoDB.Positron"
## [111] "CartoDB.PositronNoLabels"
## [112] "CartoDB.PositronOnlyLabels"
## [113] "CartoDB.DarkMatter"
## [114] "CartoDB.DarkMatterNoLabels"
## [115] "CartoDB.DarkMatterOnlyLabels"
## [116] "CartoDB.Voyager"
## [117] "CartoDB.VoyagerNoLabels"
## [118] "CartoDB.VoyagerOnlyLabels"
## [119] "CartoDB.VoyagerLabelsUnder"
## [120] "HikeBike"
## [121] "HikeBike.HikeBike"
## [122] "HikeBike.HillShading"
## [123] "BasemapAT"
## [124] "BasemapAT.basemap"
## [125] "BasemapAT.grau"
## [126] "BasemapAT.overlay"
## [127] "BasemapAT.highdpi"
## [128] "BasemapAT.orthofoto"
## [129] "nlmaps"
## [130] "nlmaps.standaard"
## [131] "nlmaps.pastel"
## [132] "nlmaps.grijs"
## [133] "nlmaps.luchtfoto"
## [134] "NASAGIBS"
## [135] "NASAGIBS.ModisTerraTrueColorCR"
## [136] "NASAGIBS.ModisTerraBands367CR"
## [137] "NASAGIBS.ViirsEarthAtNight2012"
## [138] "NASAGIBS.ModisTerraLSTDay"
## [139] "NASAGIBS.ModisTerraSnowCover"
## [140] "NASAGIBS.ModisTerraAOD"
## [141] "NASAGIBS.ModisTerraChlorophyll"
## [142] "NLS"
## [143] "JusticeMap"
## [144] "JusticeMap.income"
## [145] "JusticeMap.americanIndian"
## [146] "JusticeMap.asian"
## [147] "JusticeMap.black"
## [148] "JusticeMap.hispanic"
## [149] "JusticeMap.multi"
## [150] "JusticeMap.nonWhite"
## [151] "JusticeMap.white"
## [152] "JusticeMap.plurality"
## [153] "Wikimedia"
## [154] "GeoportailFrance"
## [155] "GeoportailFrance.parcels"
## [156] "GeoportailFrance.ignMaps"
## [157] "GeoportailFrance.maps"
## [158] "GeoportailFrance.orthos"
## [159] "OneMapSG"
## [160] "OneMapSG.Default"
## [161] "OneMapSG.Night"
## [162] "OneMapSG.Original"
## [163] "OneMapSG.Grey"
## [164] "OneMapSG.LandLot"map1 <- leaflet() %>%
addProviderTiles(providers$Esri.WorldTopoMap)
map1The simplest addition to base map is a marker, or identifiable point. Markers can be added to the base map using addMarkers function. At minimum a marker requires a longitude and latitude, but it is useful to also include labels, which display text when you hover over the marker, or popups, which display a dialog box when you click on the marker.
fiji_reduced <- quakes[1:10,]
map1 <- leaflet(data = fiji_reduced) %>%
addTiles() %>%
addMarkers(lng = ~long, lat = ~lat,
label = ~paste(depth, "meters"),
popup = ~paste(mag, "Richter Magnitude", "<br/>", depth, "Depth"))
map1What is displayed by a popup or label is an interpretted HTML strings, meaning commands like “<br/>” can be used to separate the text onto a new line, and other HTML commands can be used to modify the text’s appearance.
Additionally, in the example above, you should notice the \(\sim\) character that is used in front of the inputs to lng, lat, label, and popup. This tells addMarkers to look for those variable within the data object specified in the call to leaflet. So, in this example, leaflet looks for the variables “lat” and “lat” in fiji_reduced. Later in this tutorial we’ll see examples that don’t do this.
Question #1
Construct a base map using “Esri.OceanBasemap” tiles. Then, filter the quakes dataset to include only observations reported by more than 40 stations and plot the location of these earthquakes on your map with hoverable labels that display the number of stations reporting (ie: text like “47 Stations Reporting”).
Adding a large number of markers to a map is computationally demanding, and it also makes the map difficult to read. To address these and other concerns, leaflet includes a suite of built in functions includingmarkerClusterOptions, which aggregates nearby markers into a group that will split when you click on it to zoom in.
map2 <- leaflet(data = quakes) %>%
addTiles() %>%
addMarkers(lng = ~long, lat = ~lat,
label = ~paste(depth, "meters"),
clusterOptions = markerClusterOptions())
map2The appearance of markers can be changed in two ways. The most straightfoward is to use addCircleMarkers, which plots circular points at the marker locations. These markers have attributes like color, opacity, and fill that can be used to add information to the map. The addCircleMarkers function uses syntax similar to addMarkers, a few of its options are demonstrated in the example below:
quakes$col <- ifelse(quakes$mag > 5, "red", "blue")
map2 <- leaflet(data = quakes) %>%
addTiles() %>%
addCircleMarkers(lng = ~long, lat = ~lat,
label = ~paste(depth, "meters"),
color = ~col,
fillOpacity = 0.25,
stroke = FALSE)
map2Here, we color earthquakes with magnitudes over 5.0 red, while coloring less intense earthquakes blue. We use the argument fillOpacity = 0.25 to make the markers more transparent, allowing viewers to see relative concentrations, and we add the argument stroke = FALSE to remove the solid outline that by default surrounds circular markers.
Custom markers can be added to a map as icons. Icons are created from an image using its web url or local file path. The example below sets up two icons, a ferry icon and a skull with crossbones icon, using images from an online repository.
oceanIcons <- iconList(
ferry = makeIcon(
"http://globetrotterlife.org/blog/wp-content/uploads/leaflet-maps-marker-icons/ferry-18.png"),
danger = makeIcon(
"http://globetrotterlife.org/blog/wp-content/uploads/leaflet-maps-marker-icons/danger-24.png"))
iconid <- ifelse(quakes$mag < 5, "ferry", "danger")
leaflet(data = quakes[1:4,]) %>% addTiles() %>%
addMarkers(~long, ~lat, icon = oceanIcons[iconid])To use an icon, it must first be setup using makeIcon, and then stored in a special list using iconList before being passed to addMarkers.
In this example you should notice that we do not use a \(\sim\) in front of oceanIcons. That is because oceanIcons exists outside of the quakes data frame.
Question #2
Construct a base map using the default map tiles. Then, using your filtered dataset (which contains earthquakes reported by more than 40 stations), construct a map displaying these earthquakes using clustered circle markers which are colored yellow if the earthquake was reported by more than 70 stations and black otherwise. You do not need to add any labels or popups to this map.
A choropleth is a type of map that displays geographic units like countries, states, census areas, or watersheds with each unit colored based upon a particular variable. The first step in creating a choropleth is to define the boundaries of each geographic unit using spatial polygons.
In this example, we load the spatial polygons (an sp object) of country boundaries from a web url using the geojson_read function.
shapeurl <- "https://raw.githubusercontent.com/johan/world.geo.json/master/countries.geo.json"
WorldCountry <- geojson_read(shapeurl, what = "sp")
Map3 <- leaflet(WorldCountry) %>% addTiles() %>% addPolygons()
Map3Learning how to create your own spatial polygons is beyond the scope of this tutorial, but for most applications you can find a shapefile that someone else has created which contains the spatial polygon boundaries you need. Unfortunately I’m not aware of a single repository of different geo.json files, so finding the boundaries you need can sometimes require a bit of searching.
The addPolygons function will add the area boundaries for each geographic unit defined in sp file. To modify the color, or other attributes, of these polygons we can supply additional arguments:
Map3 <- leaflet(WorldCountry) %>% addTiles() %>%
addPolygons(
fillColor = "red",
weight = 2,
opacity = 1,
color = "white",
fillOpacity = 0.7)
Map3You should recognize that the arguments within addPolygons expects you to provide a vector of length 1 (like we did in the example above) or a vector that is the same length as the number of spatial polygons (what is needed to make a choropleth). In the example below, we merge data from the Happy Planet with the countries in our sp file, and then use the variable LifeExpectancy to create a choropleth:
## Join the Happy Planet data to the countries in WorldCountry
CountryHappy <- left_join(data.frame(Name = WorldCountry$name), happiness, by = c("Name" ="Country"))
## Create a color palette (using the "magma" color scheme)
pal <- colorBin("magma", domain = CountryHappy$LifeExpectancy)
## Use the palette and LifeExpectancy to color the spatial polygons
Map3 <- leaflet(WorldCountry) %>% addTiles() %>%
addPolygons(
fillColor = pal(CountryHappy$LifeExpectancy),
weight = 2,
opacity = 1,
color = "white",
fillOpacity = 0.7)
Map3In this example you should notice a few things:
left_join because every spatial polygon (country) needs to have a fillColor.WorldCountry had a matching entry in the Happy Planet data, but addPolygons will use grey to color spatial polygons with an NA value for their fill color.colorBin.To add interactivity to our map we can use the highlight argument to indicate when the mouse is hovering over a particular spatial polygon.
Map3 <- leaflet(WorldCountry) %>% addTiles() %>%
addPolygons(
fillColor = pal(CountryHappy$LifeExpectancy),
weight = 2,
opacity = 1,
color = "white",
fillOpacity = 0.7,
highlight = highlightOptions(
weight = 3,
color = "grey",
fillOpacity = 0.7,
bringToFront = TRUE))
Map3In addition to adding highlights, we can add labels and popups to our spatial polygons.
myLabels <- paste("<strong>", CountryHappy$Name, "</strong>", "<br/>",
"Life Expectancy:", CountryHappy$LifeExpectancy)
myPopups <- paste("Happy Planet Rank", CountryHappy$HPIRank)
Map3 <- leaflet(WorldCountry) %>% addTiles() %>%
addPolygons(
fillColor = pal(CountryHappy$LifeExpectancy),
weight = 2,
opacity = 1,
color = "white",
fillOpacity = 0.7,
highlight = highlightOptions(
weight = 3,
color = "grey",
fillOpacity = 0.7,
bringToFront = TRUE),
label = lapply(myLabels, HTML),
popup = myPopups)
Map3In this example you should notice:
addPolygons to read the HTML commands properly we need to apply (using lapply) the function HTML from the htmltools package to each element of the vector we supply to the label argument.addPolygons that we’re using HTML.As a final step, we might consider adding a legend to our choropleth. The addLegend function will add a legend as a new layer:
myLabels <- paste("<strong>", CountryHappy$Name, "</strong>", "<br/>",
"Life Expectancy:", CountryHappy$LifeExpectancy)
myPopups <- paste("Happy Planet Rank", CountryHappy$HPIRank)
Map3 <- leaflet(WorldCountry) %>% addTiles() %>%
addPolygons(
fillColor = pal(CountryHappy$LifeExpectancy),
weight = 2,
opacity = 1,
color = "white",
fillOpacity = 0.7,
highlight = highlightOptions(weight = 3,
color = "grey",
fillOpacity = 0.7,
bringToFront = TRUE),
label = lapply(myLabels, HTML),
popup = myPopups) %>%
addLegend(pal = pal, values = CountryHappy$LifeExpectancy,
title = "Life Expectancy", position = "bottomright")
Map3Questions:
The us_cities data set in the geojsonio package contains state capitals and all cities in the United States with populations larger than 40,000. The file located at “https://raw.githubusercontent.com/PublicaMundi/MappingAPI/master/data/geojson/us-states.json” contains spatial polygon boundaries for all 50 states (and the District of Columbia and Puerto Rico). The USArrests data set in the datasets package contains violant crime rates by state.
For this question you should create a map with the following features:
#1) Spatial polygons for each state filled according to colors from the “viridis” color scheme and based upon that state’s murder rate in the USArrests data set. #2) Black outlines for each spatial polygon. #3) A legend in the bottom left that depicts the color scheme with an informative title. #4) A highlight that outlines a state in white when a user hovers over it. #5) Popups that show each state’s name in bold, along with its UrbanPop below the name on a new line. #6) Circle markers (with radius = 2) for each city in the us_cities data set where state capitals are colored “red” and all other cities are colored yellow #7) Labels for each city that display the city’s name (it is okay to keep the state abbreviation) and the city’s population
Hints:
bringToFront = FALSE in your highlight to avoid masking the city labels when hovering over a stateshapeurl <- "https://raw.githubusercontent.com/PublicaMundi/MappingAPI/master/data/geojson/us-states.json"
UnitedStates <- geojson_read(shapeurl, what = "sp")
data("USArrests")
data("us_cities")
USArrests$State = row.names(USArrests)
head(us_cities)## name country.etc pop lat long capital
## 1 Abilene TX TX 113888 32.45 -99.74 0
## 2 Akron OH OH 206634 41.08 -81.52 0
## 3 Alameda CA CA 70069 37.77 -122.26 0
## 4 Albany GA GA 75510 31.58 -84.18 0
## 5 Albany NY NY 93576 42.67 -73.80 2
## 6 Albany OR OR 45535 44.62 -123.09 0head(USArrests)## Murder Assault UrbanPop Rape State
## Alabama 13.2 236 58 21.2 Alabama
## Alaska 10.0 263 48 44.5 Alaska
## Arizona 8.1 294 80 31.0 Arizona
## Arkansas 8.8 190 50 19.5 Arkansas
## California 9.0 276 91 40.6 California
## Colorado 7.9 204 78 38.7 Colorado#load us state shapefiles and join with arrest data
shapeurl2 <- "https://raw.githubusercontent.com/PublicaMundi/MappingAPI/master/data/geojson/us-states.json"
Stateshapes <- geojson_read(shapeurl2, what = "sp")
StatesArrests <- left_join(data.frame(name = Stateshapes$name), USArrests, by = c("name" = "State"))
#create color palate for arrest data
pal <- colorBin("viridis", domain = StatesArrests$Murder)
#create state labels and popups
stateLabels <- paste("<strong>", StatesArrests$name, "</strong>", "<br/>", "Urban Population:", StatesArrests$UrbanPop)
statepopups <- paste("<strong>", StatesArrests$name, "</strong>", "<br/>", "Urban Population:", StatesArrests$UrbanPop)
#Sort markers by city category
us_cities$col <- ifelse(us_cities$capital == 0, "yellow", "red")
head(us_cities)## name country.etc pop lat long capital col
## 1 Abilene TX TX 113888 32.45 -99.74 0 yellow
## 2 Akron OH OH 206634 41.08 -81.52 0 yellow
## 3 Alameda CA CA 70069 37.77 -122.26 0 yellow
## 4 Albany GA GA 75510 31.58 -84.18 0 yellow
## 5 Albany NY NY 93576 42.67 -73.80 2 red
## 6 Albany OR OR 45535 44.62 -123.09 0 yellow#Create map
questionmap <- leaflet(Stateshapes) %>%
addProviderTiles(providers$Esri.WorldTopoMap) %>%
addPolygons(
fillColor = pal(StatesArrests$Murder),
weight = 1,
opacity = 1,
color = "black",
fillOpacity = 0.7,
highlight = highlightOptions(
weight = 3,
color = "white",
fillOpacity = 0.7,
bringToFront = FALSE),
label = lapply(stateLabels, HTML),
popup = statepopups) %>%
addCircleMarkers(data = us_cities, lng = ~long, lat = ~lat,
label = ~paste(name, "population:", pop),
color = ~col,
radius = 2,
fillOpacity = 0.25,
stroke = FALSE
) %>%
addLegend(pal = pal, values = StatesArrests$Murder, title = "Murder Rate (per 100,000)", position = "topright")
questionmap