visualization.qmd

# Visualization {#sec-visualization}
We'll explore a few of the primary libraries in R for mapping (both static and dynamic).

## Goals and Outcomes
* Gain familiarity with plotting and visualizing spatial data in R
* Work with four specific visualization and plotting libraries:
  + `ggplot2`
  + `leaflet`
  + `mapview`
  + `tmap`

Throughout this section we'll use the following packages:
```{r}
#| warning: false
library(tigris)
library(tidycensus)
library(ggplot2)
library(cowplot)
library(tmap)
library(tmaptools)
library(sf)
library(ggspatial)
library(leaflet)
library(Rspatialworkshop)
library(mapview)
mapviewOptions(fgb=FALSE)
library(readr)
library(dplyr)
library(readxl)
library(terra)
```
R is fantastic for making publication quality static maps, and for generating repetitive graphics through scripts;  we've already seen examples of how to make simple maps using base plotting,`ggplot`, and `tmap`. There are also several packages in R that link R code to plotting libraries developed in Javascript (or other languages) for interactive plotting and web integration.

It can be hard to decide which mapping packages to learn and use - some nice advice from Martin Tennekes who created `tmap`:

- If you know some `ggplot`, don't care about interactive maps, and don't want to spend a lot of time learning new packages, use `ggplot`
- If you want interactive maps as flexible as possible, use `leaflet`
- If you want to simply explore spatial objects ineractively as easily as possible, use `mapview`
- Otherwise, use `tmap`!
  
Also, as pointed out in [Spatial Data Science](https://r-spatial.org/book/08-Plotting.html) by Edzar Pebesma and Roger Bivand, 'Every plot is a projection' so it's essential to have an understanding of coordinate reference systems and projections when visualizing spatial data in R - as they point out 'any time we visualize, in any way, the world on a flat device, we project: we convert ellipsoidal coordinates into Cartesian coordinate'.


## ggplot2
[ggplot2](https://ggplot2.tidyverse.org/) now has support for `geom_sf` that was developed in conjunction with the development of `sf` and helps creating publication quality maps directly using `sf` objects. An introduction to this is found in [Moreno and Basille (2018)](https://www.r-spatial.org/r/2018/10/25/ggplot2-sf.html). 

Here we'll show some of the useful functionality of `ggplot2` with `geom_sf` objects pulling census and American Community Survey data using the [tidycensus](https://github.com/walkerke/tidycensus) package.  

:::{.callout-note}
Note that to use `tidycensus` you'll need to set your Census API key. A key can be obtained from [here](http://api.census.gov/data/key_signup.html).
:::

```{r message=FALSE, error=FALSE, warning=FALSE}
#| message: false
#| error: false
#| warning: false
#| echo: true
#| results: hide

mult_tracts <- tidycensus::get_acs(state='OR',county='Multnomah',geography='tract', variables=c('B19013_001','B16010_028','B01003_001'), geometry=TRUE)  

# tidy data
mult_wide <- mult_tracts |> 
  sf::st_transform(2153) |> #UTM 11N
  dplyr::filter(!is.na(estimate)) |> 
  tidyr::pivot_wider(
    names_from = c("variable"),
    values_from = c("estimate","moe")
  ) |> 
  dplyr::rename(Median_Income=estimate_B19013_001, College_Education=estimate_B16010_028,Total_Pop=estimate_B01003_001) |> 
  dplyr::mutate(Perc_College_Ed = (College_Education / Total_Pop) * 100) |> 
  dplyr::filter(!is.na(Median_Income) & !is.na(Perc_College_Ed))
  
    
# Median Income
mult_wide |> 
  ggplot() + geom_sf(aes(fill = Median_Income)) + 
  scale_y_continuous() +
  scale_fill_viridis_c(option = "mako") +
  theme_minimal_grid(12)
```

```{r message=FALSE, error=FALSE, warning=FALSE}
#| message: false
#| error: false
#| warning: false
mult_wide |> 
  ggplot() + geom_sf(aes(fill = Perc_College_Ed)) + 
  scale_y_continuous() +
  scale_fill_viridis_c(option = "mako") +
  theme_minimal_grid(12)
```

A great resource for mapping with `ggplot2` is [this blog post from 2018](https://r-spatial.org/r/2018/10/25/ggplot2-sf-2.html). 

::: {.callout-note appearance="simple" icon="false"}
### Exercise
Using the blog post link above, get counties for your state using `tigris` (or other!) package. Ignore the extra packages and some of the bells and whistles in the blog post, but try using a few of their ideas and adding a few of the ggplot parameters to make an interesting map using county data. You can try loading city data as well.
:::

::: {.callout-note appearance="simple" icon="false" collapse="true"}
### Solution
One idea
```{r}
#| message: false
#| error: false
#| warning: false
#| echo: true
#| results: hide
counties <- tigris::counties("Oregon", cb = TRUE)
counties <- counties |> 
  sf::st_transform(2991) |> 
  dplyr::mutate(area=as.numeric(st_area(counties)))

ggplot() +
    geom_sf(data = counties, fill = NA, color = gray(.5)) 
```

```{r message=FALSE, error=FALSE, warning=FALSE}
ggplot() +
    geom_sf(data = counties, aes(fill = area))  +
  scale_fill_viridis_c(trans = "sqrt", alpha = .4) +
    
    theme(panel.grid.major = element_line(color = gray(0.5), linetype = "dashed", 
        size = 0.5), panel.background = element_rect(fill = "aliceblue")) +
    annotation_scale(location = "bl", width_hint = 0.4) +
    annotation_north_arrow(location = "tr", which_north = "true",
        style = north_arrow_fancy_orienteering)
```
:::

## leaflet
[Leaflet](https://leafletjs.com/) is an extremely popular open-source javascript library for interactive web mapping, and the `leaflet` R package allows R users to create Leaflet maps from R. `Leaflet` can plot `sf` or `sp` objects, or x / y coordinates, and can plot points, lines or polygons. There are a number of [base layers you can choose from](http://leaflet-extras.github.io/leaflet-providers/preview/index.html).  It's worth spending some time exploring the excellent [Leaflet for R site](https://rstudio.github.io/leaflet/).

Here we make a simple leaflet map of our the location of the EPA Pacific Ecological Systems Division lab in Corvallis where I work with a custom popup we create:
```{r leaflet, message=FALSE, warning=FALSE, error=FALSE}
content <- paste(sep = "<br/>",
  "<b><a href='https://www.epa.gov/greeningepa/pacific-ecological-systems-division-pesd-laboratory'>EPA Lab Corvallis</a></b>",
  "200 S.W. 35th Street ",
  "Corvallis, OR 97333 "
)

leaflet()  |>  addTiles()  |> 
  addPopups(-123.290391, 44.565548,  content,
    options = popupOptions(closeButton = FALSE)
  )
```

## mapview
[Mapview](https://r-spatial.github.io/mapview/index.html) is a package designed for quick and easy interactive visualizations of spatial data - it makes use of `leaflet` but simplifies mapping functions compared to the `leaflet` package.  

It's easy to layer features with `mapview` - you can supply a list of objects to `mapview` or use + syntax as with `ggplot`.

Here we'll plot stream gages within Benton County:
```{r}
#| echo: true
#| results: hide
#| error: false
#| message: false
#| warning: false
counties <- counties("Oregon", cb = TRUE)
benton <- counties[counties$NAME=='Benton',]
fpath <- system.file("extdata", "Gages_flowdata.csv", package="Rspatialworkshop")

gages <- read_csv(fpath,show_col_types = FALSE)
gages_sf <- gages %>%
  st_as_sf(coords = c("LON_SITE", "LAT_SITE"), crs = 4269, remove = FALSE) 
st_crs(gages_sf)==st_crs(benton)
# remember spatial indexing from Geoprocessing section?
gages_benton <- gages_sf[benton,]
mapview(gages_benton) + benton
```

We can also use handy convenience packages like the [AOI](https://mikejohnson51.github.io/AOI/) package by Mike Johnson for flexible, term based geocoding to return areas of interest as `sf` objects and map them:
```{r}
AOI::aoi_get(list("Corvallis, OR", 10, 10)) |> mapview()
```

::: {.callout-note appearance="simple" icon="false"}
### Exercise
Glance through the [`mapview` basics](https://r-spatial.github.io/mapview/articles/articles/mapview_01-basics.html) and adjust legend and attributes.  Take a look at [`mapview` advanced controls](https://r-spatial.github.io/mapview/articles/articles/mapview_02-advanced.html) as well and try plotting stations and the county polygon together, this time modifying features such as thicker black outline and transparent fill for the county outline and colorizing gages by their average flow ('AVE').
:::

::: {.callout-note appearance="simple" icon="false" collapse="true"}
### Solution
```{r message=FALSE, error=FALSE, warning=FALSE}
mapview(gages_benton,zcol='AVE') + mapview(benton, alpha.regions=.07, color='black', lwd=2)
```
:::

### Adding Web Map services in `mapview`
We'll visualize data with `mapview` and load a web map service layers alongside using `mapview` and underlying `leaflet` functionality.

First we load load an excel file containing coordinate information in a known projection and promote to an `sf` spatial data frame.
```{r message=FALSE, error=FALSE, warning=FALSE}
fpath <- system.file("extdata", "Station_Locations.xlsx", package="Rspatialworkshop")
stations <- read_xlsx(fpath)
glimpse(stations)
summary(stations$x)
# common clean up steps for spatial data - we can't use data missing coordinates so drop those records
stations <- stations[complete.cases(stations),]
# often spatial data in projected coordinates will have missing negative values for x values - common thing to fix:
stations$x[stations$x > 0] <- 0 - stations$x[stations$x > 0]
stations <- stations  |>  
  st_as_sf(coords = c("x", "y"), remove = FALSE)

# in this case we know the particular Albers projection and have the information as a proj string
st_crs(stations) <- "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs" 
```

Basic interactive map of our spatial stations with `mapview`:
```{r message=FALSE, error=FALSE, warning=FALSE}
mapview(stations)
```

Here we'll load a web map servcice (WMS) for the  National Hydrography dataset. We're looking at stream stations so imagine we want to visualize how closely these sites match a known rivers and stream network:
```{r message=FALSE, error=FALSE, warning=FALSE}
# create a mapview object with our stations:
m <- mapview(stations, legend=FALSE)

# we configure the map attribute of our mapview object - try:
# 'attributes(m) 
# to see those attributes

#  The map attribute for mapview accepts leaflet methods - in this case we use addWMSTiles to add web map service tiles to the map
m@map <- m@map  |>  addWMSTiles(group = 'NHDPlus',
                              "https://watersgeo.epa.gov/arcgis/services/NHDPlus_NP21/NHDSnapshot_NP21/MapServer/WmsServer?",
                              layers  = 4,
                              options = WMSTileOptions(format = "image/png", transparent = TRUE),
                              attribution = "")  |>  addWMSTiles(group = 'NHDPlusHR',
                                                                "https://hydro.nationalmap.gov/arcgis/services/NHDPlus_HR/MapServer/WMSServer?",
                                                                layers  = 9,
                                                                options = WMSTileOptions(format = "image/png", transparent = TRUE),
                                                                attribution = "")   |>  mapview:::mapViewLayersControl(names = c("NHDPlus","NHDPlusHR"))
m
```

## tmap
[tamp](https://cran.r-project.org/web/packages/tmap/vignettes/tmap-getstarted.html) uses the same syntax as ggplot: the grammar of graphics - and it supports both static and interactive modes. Like `ggplot2`, each dataset added to a `tmap` plot can be mapped in a number of different ways such as location, color, size, etc. The [Making maps with R section of Geocomputation with R](https://r.geocompx.org/adv-map) has an excellent and in depth treatment of using `tmap`.

We can explore the basics using the counties of Oregon data we've been using in previous examples and exersices.
```{r}
# Add fill layer to Oregon counties
tm_shape(counties) +
  tm_fill() 
# Add border layer to Oregon counties
tm_shape(counties) +
  tm_borders() 
# Add fill and border layers to Oregon counties
tm_shape(counties) +
  tm_fill() +
  tm_borders() 
```

### Basemaps with `tmap`
```{r}
data(metro)
oregon_cnties <- st_transform(counties, 4326) 
metro |> 
  st_transform(4326) |> 
  st_crop(st_bbox(oregon_cnties))  -> or_metro
tmap_mode("view")
tm_basemap("OpenStreetMap") +
tm_shape(or_metro) + tm_bubbles(size = "pop2020", col = "red") +
tm_tiles(paste0("http://services.arcgisonline.com/arcgis/rest/services/Canvas/",
    "World_Light_Gray_Reference/MapServer/tile/{z}/{y}/{x}"), group = "Labels")
```

### Choropleths with `tmap`
```{r}
tm_shape(mult_wide) + tm_polygons("Median_Income")
```

```{r}
tm_shape(mult_wide) + tm_polygons("Perc_College_Ed")
```
That’s a pretty basic map - we can adjust a number of settings such as: 

- breaks: we can set different breaks for our map 
- bins: we can control the number of bins
- palette: we can change the color palette
- layout: put the legend outside of the map, increase legend width

These are just a few - let’s play with those to start with.
```{r}
breaks = c(0, 20, 40, 80)
t1 <- tm_shape(mult_wide) + tm_polygons("Perc_College_Ed", breaks=breaks,palette = "BuGn") + tm_layout(legend.outside=TRUE, legend.outside.position = "right", legend.outside.size=.5)
t2 <- tm_shape(mult_wide) + tm_polygons("Median_Income", n=3,palette = "BuGn") + tm_layout(legend.outside=TRUE, legend.outside.position = "right", legend.outside.size=.5)
tmap_arrange(t1, t2, nrow = 2)
```

### Faceting
```{r}
# Set mode to interactive
tmap_mode("view")

# Plot it out
tm_shape(mult_wide) + tm_polygons(c("Median_Income", "Perc_College_Ed")) + tm_facets(sync = TRUE, ncol = 2)
```

### Plotting rasters and vectors with tmap
Bring in boundary and elevation of Crater Lake NP (datasets in `Rspatialworkshop` package) and plot with tmap
```{r tmap_raster_plot, message=FALSE, warning=FALSE, error=FALSE}
# Set mode to plot
tmap_mode("plot")
data(CraterLake)
raster_filepath <- system.file("extdata", "elevation.tif", package = "Rspatialworkshop")
elevation <- rast(raster_filepath)

map_crlk <- tm_shape(CraterLake) + tm_polygons(lwd = 2)
map_crlkel = map_crlk +
  tm_shape(elevation) + tm_raster(alpha = 0.7,palette = terrain.colors(12)) + tm_layout(legend.position = c("left","bottom"),
          legend.width = 1)

map_crlkel
```

::: {.callout-note appearance="simple" icon="false"}
### Exercise
Take a few minutes to read through [Making maps with R in Geocomputation with R](https://r.geocompx.org/adv-map#adv-map) or the [tmap GitHub pages](https://r-tmap.github.io/tmap/) and try to find a few modifications to plotting the  elevation and the Crater Lake park boundary you can come up with to make a more interesting map with `tmap`.
:::

::: {.callout-note appearance="simple" icon="false" collapse="true"}
### Solution
```{r}
tm_shape(elevation) +
  tm_raster(title = "Elevation", 
            style = "cont",
            palette = "-Spectral") +
  tm_shape(CraterLake) +
  tm_borders(col = "black", 
             lwd = 1)+ 
  tm_scale_bar(breaks = c(0, 10, 20),
               text.size = .5,
               position = c("left", "bottom")) +
  tm_grid(n.x = 4, n.y = 3) +
  tm_compass(position = c("right", "top"), 
             type = "arrow", 
             size = 1.5)
```


:::