In the last lesson we learned about some basic concepts related to spatial data (e.g. projections and data types) and some of the basics of handling and displaying spatial data in R. In that lesson we saw how to generate maps using data from the maps package and the plotting functions from the ggplot2 package. We also learned how to bring layer data on to a map – we created a point map showing recent earthquakes and a choropleth map showing temperature data for each US state.

As we saw the ggplot2 package allows you to create high-quality maps (there is much more functionality than we were able to cover in one class) and you could easily do much of your map making work in ggplot2. However, there are many other packages/tools available for creating maps in R and it is worthwhile to be introduced to two highly useful and widely used packages – leaflet and tmap. The leaflet and tmap packages, allow you to create both interactive (leaflet and tmap) and static maps (tmap) and will expand the map making tools at your disposal.

Again, I strongly recommend those of you who are interested to map making (and spatial analysis more broadly) to consult the excellent and freely available textbook Geocomputation with R.


Interactive maps with leaflet

Leaflet is a very popular, open-source library for creating interactive maps. You’ve almost certainly interacted Leaflet maps before as they are widely used by The US National Parks Service, The New York Times and The Washington Post, along with a host of other popular sites. You can easily generate Leaflet maps in R using the leaflet package. These interactive maps will be nicely embedded in your knit html documents, thus allowing you to easily share them.

Let’s first load in the leaflet package, which we will use to make interactive maps. We’ll also load in the tidyverse package, since we will use some of its data wrangling/handling packages.

library(leaflet)
library(tidyverse)

With Leaflet you can quickly load in a basemap image (tile) into the interactive map widget.

To do this you call the leaflet() function, which create the map widget and then call the addTiles() function, which loads in a basemap tile (i.e. a basemap image).

Let’s give that a try. The code below will bring up an interactive map in your Viewer pane. After running the code take a minute to interact with the map widget to get familiar with its functionality.

leaflet() %>% 
  addTiles()

Notice that unlike in ggplot2, we add additional functions to the leaflet call using the pipe %>% operator and NOT a +.


The default basemap tile (OpenStreetMap) is suitable for some applications, though in many cases you will want to use a different map tile. Leaflet provides access to a range of different basemap tiles – providing basemaps that suit a wide-range of applications.

The basemap tiles available in Leaflet are listed here. Take a look at the link to see all of the options.

The addProviderTiles() function allows you to specify the basemap tile you want to use. Let’s illustrate this with an example. We’ll us the National Geographic basemap.

leaflet() %>% 
  addProviderTiles(providers$Esri.NatGeoWorldMap)

Try out a few other map tiles of interest. Note: some provider tiles do not knit properly or show up in the Viewer pane. It seems to be an RStudio issue, as the map will generally appear properly when you show it in a web browser window.


Adding layers/geometries to Leaflet maps

So far we’ve just generated a basemap, however we almost always want to display additional layers of data on top of the basemap. We can easily use data from a range of data formats (e.g. data frames, sp and sf objects, GeoJSON) display a range of geometries (e.g. point markers, polygons, lines, …).

Markers

Let’s load in the locations of the US state capitals and save the data to a data frame.

state_cap_locs <- read_csv("https://stahlm.github.io/ENS_215/Data/us_capitals_locs.csv", skip = 2)

state_cap_locs <- state_cap_locs %>% 
  filter(state_cd %in% state.abb) # remove locations that aren't US state list (e.g. Puerto Rico, Guam, ...)

We can now generate a leaflet map where we show the location of each capital as a marker.

map_cap <- leaflet(data = state_cap_locs) %>% 
  addProviderTiles(providers$Esri.NatGeoWorldMap) %>% 
  addMarkers(lng = ~long_deg, lat = ~lat_deg, label = ~state_cd)

Note the syntax in the addMarkers() function. We assign the latitude and longitude aesthetics using lng = ~Name_of_dataframe_col. Where Name_of_dataframe_col is the name of the column in our state_cap_locs data frame that has the variable interest. We also assigned the variable state_cd from the state_cap_locs data frame to the label aesthetic. You will see the state abbreviation will appear when you hover over a given state capital’s marker .

Note that we assign the variable to the aesthetic using = ~. In leaflet variables are mapped to aesthetics using that syntax.

Also note that just like graphics generated with ggplot2, we can similarly assign the graphic to an R object, in this case map_cap. We’ll now display map_cap

map_cap


There are a number of other marker features that you can adjust, including the marker icon used and marker size. For additional information and examples, check out the Leaflet for R website.

  • Create another version of the above map where you adjust some of the marker features.


Circles

Another useful geometry to add to leaflet maps are circles, which are added using the addCircles() function. Circles can be particularly useful in certain situations, given that the radii of the circles are specified in meters and thus the circle will scale accordingly as you zoom in and out. This differs from markers, which have a fixed pixel size and thus the size of the marker area to a given geographical area changes as you zoom in and out.

One particularly useful application of the circles geometry in Leaflet is for creating circular buffers around regions. For instance you might want to map a 10 km radius around a nuclear power plant to help visualize the areas most likely to be affected during an emergency.

Let’s illustrate the use of circles with an example. Around each US state capital we will create circles with at 30 km radius. This could be useful if we interested in seeing the regions with relatively easy travel distance to each US state capital. Note that the radius is set in units of meters.

leaflet(state_cap_locs) %>% 
  addTiles() %>% 
  addCircles(lng = ~long_deg, lat = ~lat_deg, radius = 30*10^3, fillColor = "transparent")


Polygons

In addition to markers and circles, displaying polygons is often required. A common situation where this is needed is when displaying border/boundaries of regions (e.g. states, countries, national parks,…).

Let’s load in a dataset containing polygons outlining the border of each US state. The data we will load in is in GeoJSON format. GeoJSON is a standardized format for representing simple geographic features (_e.g. points, lines, polygons) along with their non-spatial attributes (_e.g. location names, populations,…). When you download spatial data in your work you will find that it is often in GeoJSON or shapefile format (another very common file format for spatial data).

In order to load GeoJSON data into R we will need to use the geojsonio package.

library(geojsonio)

Now we will use the geojson_read() function to load in the .geojson file that has the state boundaries. You can see that we’ve set the what argument equal to "sp". This tells the geojson_read() function to convert the GeoJSON file into a Spatial data-structure that can be easily computed on, manipulate, and plotted in many R packages including Leaflet.

state_borders <- geojson_read("https://stahlm.github.io/ENS_215/Data/gz_2010_us_040_00_5m.geojson", what = "sp")

# geojson state boundaries from http://eric.clst.org/tech/usgeojson/

Let’s now add the state polygons to our map with the state capitals. We can add the polygons using the addPolygons() function in leaflet. Notice that we need to specify within the addPolygons() function call the dataset where the polygon data is stored.

map_cap %>% 
  addPolygons(data = state_borders)


Controlling map layers

Oftentimes you will want to give yourself (or other users) the ability to interactively toggle map layers on and off and to change the basemap between different map tiles. You can easily do this with the addLayersControl() function in leaflet.

Let’s create our map with the state border polygons and state capital markers. This time we will also load in three different basemap tiles. We will use the addLayersControl() function to allow us to toggle between the basemap and layers in the map we create. Note that each of the geometry layers (markers and polygons) and each of the tile layers are assigned to a group with a name of our choosing. The group provides a unique identifier for each layer which is used in the addLayersControl() to specify the layers that we will allow to be toggled between.

leaflet(state_cap_locs) %>% 
  addTiles(group = "OSM") %>%
  addProviderTiles(providers$Stamen.Toner, group = "Toner") %>% 
  addProviderTiles(providers$Stamen.Watercolor, group = "Watercolor") %>% 
  addMarkers(lng = ~long_deg, lat = ~lat_deg, label = ~state_cd, group = "Capitals") %>% 
  addPolygons(data = state_borders, group = "Borders") %>% 
  
  addLayersControl(
    baseGroups = c("OSM", "Toner","Watercolor"),   # specify the base map groups that we can toggle between
    overlayGroups = c("Capitals", "Borders"),   # specify the geometry layer groups that we can toggle between
    options = layersControlOptions(collapsed = FALSE) 
  )
  • Try creating one or two maps of your own using the techniques learned in the section above. Many of you will likely want to create an interactive map for your final project (now’s your chance to give it a try!).


Making maps with the tmap package

The tmap package is a widely used package for creating both static and interactive maps in R. The tmap package have a wide range of functionality and flexibility and allows for the display of more spatial classes than ggplot2 – thus allowing for certain types of maps that are difficult to create with ggplot2. A nice feature of the tmap package is that it is based on the “grammar of graphics” and thus the syntax is very similar to ggplot2. This makes learning tmap relatively easy for those who already have experience with ggplot2 (i.e. you).

Let’s load in the tmap package, along with the spData package which has provides some nice spatial datasets that we can use to create example maps. We will also load in the sf package, which provides functionality in R for simple features (sf) data types. Simple features are a standard for storing storing spatial data. Simple featurs are widely used in the field of spatial analysis and given that they are spatial data frames, they are well suited to use in R. You can learn more about simple features and R here.

library(tmap)
library(spData)
library(sf)


We will use the tmap package to create a world map based on the world country polygons simple features dataset from the spData package.

world_map_data <- world # save the world data from the spData package to our own R object

Type View(world_map_data) in your console. You will see that the world_map_data simple features object looks like a typical data frame, except that it has a geom column which stores the spatial data associated with the other variables.

Now let’s create a choropleth map of life expectancy using tmap and the world_map_data dataset.

tm_shape(world_map_data) +
  tm_polygons("lifeExp")

With the tm_shape() function we define the spatial dataset to be used. Then we plot the spatial data by adding geometry layers (just like we do in ggplot2). Here we’ve added polygons using the tm_polygon() function and we specified the variable (“lifeExp”) from the spatial dataset that will be used to color the polygons.

There are other geometries (e.g. border, fill, lines, dots,…) that you can layer on top of your map.

Let’s step back for a moment and demonstrate how we can build up a map from it’s basic components.

map_1 <- tm_shape(world_map_data) +
  tm_fill()

map_2 <- tm_shape(world_map_data) +
  tm_fill("lifeExp") 

map_3 <- tm_shape(world_map_data) +
  tm_fill("lifeExp") +
  tm_borders("black") 

tmap_arrange(map_1, map_2, map_3) # arrange the three maps

You can see from the above code, that we have control of each of the individual geometries – in this example fill tm_fill() and borders tm_borders(). We can also set the aesthetics of each of the geometries. Note that we set the aesthetic with a character string specifying the variable to use (e.g. "lifeExp").


Just like in ggplot2 you can facet map graphics, thus allowing you to quickly and easily create sets of maps.

tm_shape(world_map_data) +
  tm_fill("continent") +
  tm_borders() +
  tm_facets(by = "continent", free.coords = TRUE)


Since the world_map_data is a spatial data frame, we can use dplyr functionality to filter(), select(), pipe, etc it.

Try creating a world map that has every country’s borders shown where only countries with a per capita GDP > 10,000 are colored lightblue. Hint you can add multiple layers to a map by using a tm_shape() function call for each layer.

My map looks like this (I’ve added a few extra features, such as graticules (lat/long lines), but you don’t need to add these).

world_map_data %>% 
  filter(gdpPercap > 10000) %>% 
  tm_shape(projection= "+proj=moll") +
  tm_borders() +
  tm_graticules(labels.size = 0.5, n.x = 20, alpha = 0.5) +

  tm_fill(col = "lightblue") +
  tm_shape(world_map_data) +
  tm_borders()

You’ll notice that in the tm_shape() function call I set projection = "+proj=moll". This set the map’s projection (i.e., how the geographic data from the Earth’s curved surface is transformed to be displayed on a flat surface). There are many, many different projections that can be used when making maps (see some common ones here). Different projections have different strengths and weaknesses (e.g., how well they preserve areas or distances and how well they preserve shapes). Some projections are only useful when mapping a particular region (e.g., polar region vs. the whole world).

In many cases where you obtain data that is distributed as spatial data (e.g., shapefile, sf) it may already have a projection assigned to it. In these cases, tmap will use the assigned projection.

If you want to reproject the data to better meet the goals of your map you can specify what projection you want to use (as we did in the code above). You can specify the projection a number of ways. One way is to use the EPSG code, which is a numeric code that refers to particular projection. There are thousands of EPSG codes, though a few commonly used ones can be found here.

You can also specify the projection using what is called a proj4string. That is what we did above when we set projection = "+proj=moll".

Try out a few other projections by changing the projection = argument in the code above. A few interesting/useful ones to use for world maps are:

  • Winkel tripel projection "+proj=wintri"

  • Lambert azimuthal equal-area projection (centered on lat/long of NYC) "+proj=laea +x_0=0 +y_0=0 +lon_0=-74 +lat_0=40"

You can also try out setting the projection using a EPSG code (e.g. projection = 4326).

To learn more about coordinate reference systems and map projections see chapters 2.4 and chapter 7 of Geocomputation with R.

Ok, now let’s move on to a few more topics with tmap and map making.


There are tons of formatting options for tmap graphics, which allow you to customize your map and to generate research quality graphics.

Let’s make a map of Alaska to demonstrate a few of these options. We’ll use the alaska simple features dataset that is in the spData package.

alaska_data <- alaska # assign alaska data to our own R object

tm_shape(alaska_data) +
  tm_fill("lightgrey") +   # add the fill
  tm_borders("black", lwd = 2) +  # add the borders and set border line width
  
  tm_compass(position = c("left", "top")) +  # add a compass
  tm_scale_bar(breaks = c(0, 250, 500), text.size = 1) +    # add a scale bar
  
  tm_layout(bg.color = "lightblue", title = "Alaska")  # set the background color and title

Above we adjusted a number of different features and stylistic settings. This example highlights just a few of the many options. To learn more about the features/options available, check out the excellent (but still under development) tmap textbook and also chapter 8 of Geocomputation with R.


Loading in and mapping external datasets in tmap

Oftentimes we will have a dataset that we want to load into R and display on our map. We’ve seen some examples already, where we loaded in files and plotted them on maps in ggplot2 and leaflet. When we were loading in point data it was very easy to map them in ggplot2 and leaflet as both of those packages accept tibbles (data frames). It is easy to load in and display data in tmap as well, however we first need to convert it to a compatible format (e.g. simple features).

The st_read() function allows us to read in shapefiles, which are one of the most common spatial data storage formats, and convert them to simple feature objects that be be visualized with tmap.

Paste the link below in your web browser and it will initiate the download of a zipfile containing a shapefile of global rivers. Once the download is complete, you should unzip the data to a folder in your “Lectures” directory. If you unzip the file to a folder called “Rivers” that lives within your “Lectures” folder then the code below should run on your computer. If you named the folder something else, then you will need to update the path below accordingly.

https://github.com/stahlm/stahlm.github.io/raw/master/ENS_215/Data/Rivers.zip


The st_read() function will read in the river shapefile and convert it to a simple feature object.

shape_river <- st_read("./Data/Rivers/ne_50m_rivers_lake_centerlines.shp", stringsAsFactors = FALSE)
## Reading layer `ne_50m_rivers_lake_centerlines' from data source 
##   `C:\Users\stahlm\Documents\Teaching_UnionCollege\Environmental_Data_Analysis\stahlm.github.io\ENS_215\Winter_2025\Lectures\Data\Rivers\ne_50m_rivers_lake_centerlines.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 462 features and 32 fields (with 1 geometry empty)
## Geometry type: MULTILINESTRING
## Dimension:     XY
## Bounding box:  xmin: -165.2439 ymin: -50.24014 xmax: 176.3258 ymax: 73.3349
## Geodetic CRS:  WGS 84

Now that the river data in a format compatible with tmap we can plot this layer. We will first generate our world map with the world polygon data we used earlier and then we will layer on an additional shape shape_river, which we will specify to be represented as lines using the tm_lines() function.

map_world_rivers <- tm_shape(world) +
  tm_fill() +
  tm_borders() +
  tm_style("natural") +
  tm_shape(shape_river) +
  tm_lines(col = "blue") 
  
map_world_rivers  
## Warning: The shape shape_river contains empty units.

Now we have a world map where which has a layer showing the world’s major rivers.


Making interactive maps in tmap

Making interactive maps with tmap is very straightforward and requires just a single additional line of code. The tmap_mode() function allows you to set the viewing mode to “view”, which is the interactive mode, or to “plot” which is the static mode. Note that once you set a particular mode, that mode applies to all subsequently generated maps until you switch to another mode.

tmap_mode("view") # set viewing mode to interactive
## tmap mode set to interactive viewing
map_world_rivers 
## Warning: The shape shape_river contains empty units.
tmap_mode("plot") # set viewing mode to static, so that any maps made later will be in static plot mode
## tmap mode set to plotting

These interactive maps will be nicely embedded in your knit html documents, thus allowing you to easily share them.

You can learn more about the interactive features available in tmap here


Obtaining spatial data directly in R

You’ve seen that some packages come pre-loaded with a handful of spatial data frames (e.g. the world map polygons from the spData package that was used last class). You’ve also seen that you can load geospatial data into R that was obtained elsewhere and then saved to your computer (e.g. downloaded to your computer from NASA, USGS, US Forest Service,…).

Now we are going to see that there are some R packages that allow you to directly query sites such as the US Census Bureau and download geospatial data directly in your code! This is amazing in that it allows you to among other thing:

  • Acquire data from the web using R code, which prevents you from having to find the site, click through options, save the data to a folder, and then load it into R
  • Makes your map easily reproducible - anyone can see the exacty query and dataset you obtained
  • Makes your map easily adaptable - with the change of a few bits of code you can select additional/new data
  • Many other benefits related to reproducible science, ease of analysis,…

There are a number of packages that allow you to query different geospatial data repositories, but today we’ll look at the TIGRIS package, which allows you to query US Census Bureau geographic data (e.g. state boundaries, county boundaries, roads, waterbodies,…), many of which are useful across a range of environmental science and policy research questions.

You’ll first need to install the tigris package and then load it in. We are also going to use the sf package below (this package provides functionality for a range of useful spatial operations). An overview of the tigris package is here.

library(tigris)


Let’s get the county borders for New York state. To do this we will use the counties() function.

NY_counties <- counties(state = "NY",
                        cb = TRUE) # cb = TRUE specifies to download a lower resolution file (faster download)

The data is downloaded as a shapefile. To plot in in tmap or ggplot we will need to convert the shapefile to a simple features (sf) spatial data frame. To do this we will use the function st_as_sf() from the sf package.

NY_counties <- st_as_sf(NY_counties) # convert to an sf object
NY_map <- NY_counties %>% 
  tm_shape() +
  tm_polygons(col = "lightgrey", border.col = "black") +
  tm_scale_bar(breaks = c(0, 50, 100), size = 1, position = c("LEFT", "BOTTOM"))     # add a scale bar

NY_map


Let’s also create a map of just Schenectady county.

Schenectady_map <- NY_counties %>% 
  filter(NAME == "Schenectady") %>% 
  tm_shape() +
  tm_polygons(col = "lightgrey", border.col = "black") +
  tm_scale_bar(breaks = c(0, 5, 10), size = 1, position = c("RIGHT", "BOTTOM"))     # add a scale bar

Schenectady_map

That’s pretty cool – we were able to download the county borders for NY and create a map in a few lines of code. This could easily be modified to create a map of other state(s). Let’s check out the function area_water() from tigris that allows us to download waterbodies (e.g. lakes, ponds, wetlands, large rivers). Note that you can download linear water features (streams, canals, ditches,…) using the linear_water() function.


Water_Schenectady <- area_water(state = "NY", county = "Schenectady")  # download the data

Water_Schenectady <- st_as_sf(Water_Schenectady) # convert to an sf object


Now let’s add the water features to the Schenectady map

Schenectady_map <- Schenectady_map +
  tm_shape(Water_Schenectady) +
  tm_polygons(col = "blue", border.alpha = 0)

Schenectady_map


Let’s also add linear water features using the linear_water() function

Linear_water_Schenectady <- linear_water(state = "NY", county = "Schenectady")  # download the data

Linear_water_Schenectady <- st_as_sf(Linear_water_Schenectady) # convert to an sf object
Schenectady_map + 
  tm_shape(Linear_water_Schenectady) +
  tm_lines(col = "blue") +
  tm_compass(position = c("right", "top")) +  # add a compass
  tm_layout(title = "Schenectady hydrography",
            frame = F)


As you can see you can relatively easily make fairly complex maps with your understanding of R and often a few lines of code. Furthermore, you can now integrate spatial analysis directly into your overall analysis workflow! We just scratched the surface today, however these basics alone should provide you with much of what you might need (with respect to maps) for your final projects. We will cover a few additional topics on spatial analysis in the last weeks of class.