---
title: "Maps in R: Additional topics"
author: "ENS-215"
date: "12-Feb-2025"
output:
  html_document:
    df_print: paged
    theme: journal
    toc: TRUE
    toc_float: TRUE
---


<br/>

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](https://r.geocompx.org/).  

<br/> 

## 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.  

```{r warning=FALSE, message=FALSE}
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.  

```{r}
leaflet() %>% 
  addTiles()
```

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

<br/>

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](http://leaflet-extras.github.io/leaflet-providers/preview/).  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.  

```{r}
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.  


<br/>


### 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.
```{r message=FALSE}
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. 
```{r eval=T}
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`
```{r}
map_cap
```

<br/> 

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](http://rstudio.github.io/leaflet/index.html) website.

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

<br/>

#### 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.

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


<br/>

#### 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.  
```{r message=FALSE, warning=FALSE}
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.  

```{r}
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.  

```{r}
map_cap %>% 
  addPolygons(data = state_borders)
```

<br/> 

### 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.  
```{r}
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!).

<br/>


## 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 features 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](https://r-spatial.github.io/sf/).  

```{r message=FALSE, warning=FALSE}
library(tmap)
library(spData)
library(sf)
```

<br/> 

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

```{r}
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. 
```{r}
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.
```{r}
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"`).  

<br/>


Just like in `ggplot2` you can facet map graphics, thus allowing you to quickly and easily create sets of maps.
```{r}
tm_shape(world_map_data) +
  tm_fill("continent") +
  tm_borders() +
  tm_facets(by = "continent", free.coords = TRUE)
```

<br/> 

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).
```{r eval = T}
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](https://en.wikipedia.org/wiki/List_of_map_projections)).  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](https://jjimenezshaw.github.io/crs-explorer/), though a few commonly used ones can be found [here](https://guides.library.duke.edu/r-geospatial/CRS). 
 
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](https://r.geocompx.org/spatial-class#crs-intro) and [chapter 7](https://r.geocompx.org/reproj-geo-data) of **Geocomputation with R**.  


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

<br/>  

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.  

```{r}
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](https://tmap.geocompx.org/) and also chapter 8 of [Geocomputation with R](https://r.geocompx.org/adv-map).


<br/>

#### 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  

<br/>

The `st_read()` function will read in the river shapefile and convert it to a simple feature object.  
```{r}
shape_river <- st_read("./Data/Rivers/ne_50m_rivers_lake_centerlines.shp", stringsAsFactors = FALSE)

```

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.  
```{r}
map_world_rivers <- tm_shape(world) +
  tm_fill() +
  tm_borders() +
  tm_style("natural") +
  tm_shape(shape_river) +
  tm_lines(col = "blue") 
  
map_world_rivers  
```

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

<br/> 

#### 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.    

```{r}
tmap_mode("view") # set viewing mode to interactive
map_world_rivers 

tmap_mode("plot") # set viewing mode to static, so that any maps made later will be in static plot mode
```


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](https://geocompr.robinlovelace.net/adv-map.html#interactive-maps) 


<br/>  

## 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 exact 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](https://github.com/walkerke/tigris).  


```{r warning= F, message= F}
library(tigris)
```


<br/>

Let's get the county borders for New York state.  To do this we will use the `counties()` function.  
```{r message=FALSE, warning=FALSE, results='hide'}
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.  
```{r}
NY_counties <- st_as_sf(NY_counties) # convert to an sf object
```


```{r message=FALSE, warning=FALSE}
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
```

<br/>

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

```{r message=FALSE, warning=FALSE}
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.  

<br/>

```{r message=FALSE, warning=FALSE, results='hide'}
Water_Schenectady <- area_water(state = "NY", county = "Schenectady")  # download the data

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

<br/>

Now let's add the water features to the Schenectady map
```{r}
Schenectady_map <- Schenectady_map +
  tm_shape(Water_Schenectady) +
  tm_polygons(col = "blue", border.alpha = 0)

Schenectady_map
```

<br/>

Let's also add linear water features using the `linear_water()` function
```{r message=FALSE, warning=FALSE, results='hide'}
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
```

```{r}
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)
```


<br/>

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.  


