03-map-making.Rmd

```{r, echo=FALSE, purl=FALSE, message=FALSE}
knitr::opts_chunk$set(results='hide', comment = "#>", purl = FALSE)

## libraries needed for R code examples
library(sf)
library(classInt)
library(RColorBrewer)
library(ggplot2)
library(leaflet)
library(tmap)
library(dplyr)

```

# Making Maps in R {#mapping}

> Learning Objectives
>
> * plot an `sf` object
> * create a choropleth map with `ggplot`
> * plot a raster map with `ggplot`
> * use `RColorBrewer` to improve legend colors 
> * use `classInt`to improve legend breaks 
> * create a choropleth map with `tmap`
> * plot a raster map with `tmap`
> * create an interactive map with `leaflet`
> * customize a `leaflet` map with popups and layer controls

------------

In the preceding examples we have used the base `plot` command to take a quick look at our spatial objects.

In this section we will explore several alternatives to map spatial data with R. For more packages see the "Visualisation" section of the [CRAN Task View](https://cran.r-project.org/web/views/Spatial.html).


## Choropleth Mapping with `ggplot2`

[`ggplot2`](http://ggplot2.org/) is a widely used and powerful plotting library for R. It is not specifically geared towards mapping, it is possible to create quite nice maps.

For an introduction to `ggplot` check out [this site](http://ggplot2.tidyverse.org/) for more pointers.

`ggplot` can plot `sf` objects directly by using the geom `geom_sf`. So all we have to do is:

```{r ggplot-sf, results='show'}
library(ggplot2)
# if you need to read this in again:
# philly_crimes_sf <- st_read("data/PhillyCrimerate") 
ggplot(philly_crimes_sf) + 
  geom_sf(aes(fill=homic_rate))
```

Homicide rate is a continuous variable and is plotted by `ggplot` as such. If we wanted to plot our map as a 'true' choropleth map we need to convert our continuous variable into a categorical one, according to whichever brackets we want to use. 

This requires two steps: 

- Determine the quantile breaks. 
- Add a categorical variable to the object which assigns each continious vaule to a bracket.

We will use the `classInt` package to explicitly determine the quantile breaks. 

```{r ggplot-sf-getbreaks, results='show'}
library(classInt)
# get quantile breaks. Add .00001 offset to catch the lowest value
breaks_qt <- classIntervals(c(min(philly_crimes_sf$homic_rate) - .00001,
                              philly_crimes_sf$homic_rate), n = 7, style = "quantile")

str(breaks_qt)
```

Ok. We can retrieve the breaks with `breaks$brks`.

We use `cut` to divicde `homic_rate` into intervals and code them according to which interval they are in. 

Lastly, we can use `scale_fill_brewer` and add our color palette. 

```{r ggplot-sf-categorical, results='show'}
philly_crimes_sf %>% 
  mutate(homic_rate_cat = cut(homic_rate, breaks_qt$brks)) %>% 
  ggplot() + 
    geom_sf(aes(fill=homic_rate_cat)) +
    scale_fill_brewer(palette = "OrRd") 
```

## Raster and ggplot

To visualize raster data using `ggplot2`, we will use the raster with the values for the digital terrain model (DTM). 

Before using ggplot we need to convert it to a dataframe. The `terra` package has an built-in function for conversion to a plotable dataframe.

```{r rast-df, results='show'}
# If you need to read it in again:
# HARV_DTM <- rast("data/HARV_dtmCrop.tif")

HARV_DTM_df <- as.data.frame(HARV_DTM, xy = TRUE)
str(HARV_DTM_df)
```

We can now use `ggplot()` to plot this data frame. We will set the color scale to `scale_fill_viridis_c` which is a color-blindness friendly color scale. [Here is more about the viridis color maps.](https://cran.r-project.org/web/packages/viridis/vignettes/intro-to-viridis.html) We will also use the `coord_fixed()` function with the default, ratio = 1, which ensures that one unit on the x-axis is the same length as one unit on the y-axis.

```{r ggplot-rast, results='show'}
ggplot() +
    geom_raster(data = HARV_DTM_df , aes(x = x, y = y, fill = HARV_dtmCrop)) +
    scale_fill_viridis_c() +
    coord_fixed()
```


## Choropleth with `tmap`

`tmap` is specifically designed to make creation of thematic maps more convenient. It borrows from the ggplot syntax and takes care of a lot of the styling and aesthetics. This reduces our amount of code significantly. We only need: 

- `tm_shape()` where we provide 
    - the `sf` object 
- `tm_polygons()` where we set
    - the attribute variable to map, 
    - the break style, and 
    - a title.


```{r tmap-plot, results='show'}
library(tmap)
tm_shape(philly_crimes_sf) +
  tm_polygons("homic_rate", 
              style="quantile", 
              title="Philadelphia \nhomicide density \nper sqKm")
```


`tmap` has a very nice feature that allows us to give basic interactivity to the map. We can switch from "plot" mode into "view" mode and call the last plot, like so:

```{r tmap-plot-viewmode, results='show', message=FALSE}
tmap_mode("view")
tmap_last()
```

Cool huh?

The `tmap` library also includes functions for simple spatial operations, geocoding and reverse geocoding using OSM. For more check `vignette("tmap-getstarted")`. 


## Raster with `tmap`

`tmap` can also plot raster files natively, for example:

```{r}
tmap_mode("plot")
tm_shape(HARV_DTM)+
    tm_raster(style = "cont", palette = "viridis")+
    tm_layout(legend.outside = TRUE)
```

See [Elegant and informative maps with `tmap`](https://r-tmap.github.io/tmap-book/) for more options.

## Web mapping with `leaflet`

`leaflet` provides bindings to the ['Leaflet' JavaScript library](http://leafletjs.com), "the leading open-source JavaScript library for mobile-friendly interactive maps". We have already seen a simple use of leaflet in the `tmap` example. 

The good news is that the `leaflet` library gives us loads of options to customize the web look and feel of the map. 

The bad news is that the `leaflet` library gives us loads of options to customize the web look and feel of the map.

Let's build up the map step by step.

First we load the `leaflet` library. Use the `leaflet()` function with an `sp` or `Spatial*` object and pipe it to `addPolygons()` function. It is not required, but improves readability if you use [the pipe operator `%>%`](https://github.com/tidyverse/magrittr) to chain the elements together when building up a map with `leaflet`. 

And while `tmap` was tolerant about our AEA projection of `philly_crimes_sf`, `leaflet` does require us to explicitly reproject the `sf` object.

```{r leaflet-polys, results='show'}
library(leaflet) 

# reproject
philly_WGS84 <- st_transform(philly_crimes_sf, 4326)

leaflet(philly_WGS84) %>%
  addPolygons()
```

To map the homicide density we use `addPolygons()` and:

- remove stroke (polygon borders)  
- set a fillColor for each polygon based on `homic_rate` and make it look nice by adjusting fillOpacity and smoothFactor (how much to simplify the polyline on each zoom level). The fill color is generated using `leaflet`'s `colorQuantile()` function, which takes the color scheme and the desired number of classes. To constuct the color scheme `colorQuantile()` returns a function that we supply to `addPolygons()` together with the name of the attribute variable to map.    
- add a popup with the `homic_rate` values. We will create as a vector of strings, that we then supply to `addPolygons()`.


```{r leaflet-popups, results='show'}
pal_fun <- colorQuantile("YlOrRd", NULL, n = 5)

p_popup <- paste0("<strong>Homicide Rate: </strong>", philly_WGS84$homic_rate)

leaflet(philly_WGS84) %>%
  addPolygons(
    stroke = FALSE, # remove polygon borders
    fillColor = ~pal_fun(homic_rate), # set fill color with function from above and value
    fillOpacity = 0.8, smoothFactor = 0.5, # make it nicer
    popup = p_popup)  # add popup
``` 


Here we add a basemap, which defaults to OSM, with `addTiles()`


```{r leaflet-basemap, results='show'}
leaflet(philly_WGS84) %>%
  addPolygons(
    stroke = FALSE, 
    fillColor = ~pal_fun(homic_rate),
    fillOpacity = 0.8, smoothFactor = 0.5,
    popup = p_popup) %>%
  addTiles()
```


Lastly, we add a legend. We will provide the `addLegend()` function with:

- the location of the legend on the map  
- the function that creates the color palette  
- the value we want the palette function to use  
- a title

```{r leaflet-legend, results='show'}
leaflet(philly_WGS84) %>%
  addPolygons(
    stroke = FALSE, 
    fillColor = ~pal_fun(homic_rate),
    fillOpacity = 0.8, smoothFactor = 0.5,
    popup = p_popup) %>%
  addTiles() %>%
  addLegend("bottomright",  # location
            pal=pal_fun,    # palette function
            values=~homic_rate,  # value to be passed to palette function
            title = 'Philadelphia homicide density per sqkm') # legend title
```

The labels of the legend show percentages instead of the actual value breaks[^20]. 

[^20]: The formatting is set with `labFormat()` and in the [documentation](https://cran.r-project.org/web/packages/leaflet/leaflet.pdf) we discover that: "By default, `labFormat` is basically `format(scientific = FALSE,big.mark = ',')` for the numeric palette, `as.character()` for the factor palette, and a function to return labels of the form `x[i] - x[i + 1]` for bin and quantile palettes (__in the case of quantile palettes, x is the probabilities instead of the values of breaks__)."

To set the labels for our breaks manually we replace the `pal` and `values` with the `colors` and `labels` arguments and set those directly using `brewer.pal()` and `breaks_qt` from an earlier section above.
    
```{r leaflet-labels, results='show'}
leaflet(philly_WGS84) %>%
  addPolygons(
    stroke = FALSE, 
    fillColor = ~pal_fun(homic_rate),
    fillOpacity = 0.8, smoothFactor = 0.5,
    popup = p_popup) %>%
  addTiles() %>%
  addLegend("bottomright", 
            colors = brewer.pal(7, "YlOrRd"), 
            labels = paste0("up to ", format(breaks_qt$brks[-1], digits = 2)),
            title =  'Philadelphia homicide density per sqkm')
```


That's more like it. Finally, I have added for you a control to switch to another basemap and turn the philly polygon off and on. Take a look at the changes in the code below.


```{r leaflet-control, results='show'}
leaflet(philly_WGS84) %>%
  addPolygons(
    stroke = FALSE, 
    fillColor = ~pal_fun(homic_rate),
    fillOpacity = 0.8, smoothFactor = 0.5,
    popup = p_popup,
    group = "philly") %>%
  addTiles(group = "OSM") %>%
  addProviderTiles("CartoDB.DarkMatter", group = "Carto") %>%
  addLegend("bottomright", 
            colors = brewer.pal(7, "YlOrRd"), 
            labels = paste0("up to ", format(breaks_qt$brks[-1], digits = 2)),
            title = 'Philadelphia homicide density per sqkm') %>%
  addLayersControl(baseGroups = c("OSM", "Carto"), 
                   overlayGroups = c("philly"))  
```


If you'd like to take this further here are a few pointers.

- [Leaflet for R](http://rstudio.github.io/leaflet/)
- [rayshader: Create Maps and Visualize Data in 2D and 3D](https://CRAN.R-project.org/package=rayshader )


[Here is an example](https://cengel.shinyapps.io/RioSlaveMarket/) using `ggplot`, `leaflet`, `shiny`, and [RStudio's flexdashboard](http://rmarkdown.rstudio.com/flexdashboard/) template to bring it all together.