class: center, middle ### Mapping with `leaflet` and `sf` <img src="img/hero_wall_pink.png" width="800px"/> ### Kelly McConville .large[Math 241 | Week 10 | Spring 2021] --- ## Announcements/Reminders * Lab 6 is due on Thursday. * Still time to sign up for the [Ugly Chart Competition](https://docs.google.com/forms/d/e/1FAIpQLSeQf3d87pfUKwXU7x-oFzegKWe1ptbibVQFPsXLdQmmlFbYWA/viewform) happening on Wednesday 4:45 - 6:00pm <img src="img/uglychart.png" width="640" style="display: block; margin: auto;" /> --- ## Goals for Today Visualizing Spatial Data * Choropleth maps * Raster maps * Interactive maps -- **Disclaimer**: * I know we are in COVID-19 overload and so throughout the semester I have intentionally minimized the number of COVID-19 related examples. * But, we are going to look at some COVID-19 related maps today. --- ## Choropleth Maps <img src="img/vaccine_rates.png" width="60%" style="display: block; margin: auto;" /> * Source: [The Guardian](https://www.theguardian.com/us-news/2021/feb/18/us-vaccine-distribution-tracker-by-state) --- ## Choropleth Maps <img src="img/choroplethCDC.png" width="100%" style="display: block; margin: auto;" /> * Source: [CDC](https://www.cdc.gov/coronavirus/2019-ncov/cases-updates/cases-in-us.html) --- ## Choropleth Maps <img src="img/choropleth1.png" width="80%" style="display: block; margin: auto;" /> * Source: [NYtimes.com](https://www.nytimes.com/interactive/2020/04/02/us/coronavirus-social-distancing.html) --- ## Choropleth Maps <img src="img/choropleth2.png" width="80%" style="display: block; margin: auto;" /> * Source: [NYtimes.com](https://www.nytimes.com/interactive/2020/04/02/us/coronavirus-social-distancing.html) --- ## Static, Raster Maps <img src="img/static_map.png" width="80%" style="display: block; margin: auto;" /> * Source: [JHU](https://coronavirus.jhu.edu/map.html) --- ## Interactive Maps <img src="img/static_map.png" width="80%" style="display: block; margin: auto;" /> * Source: [JHU](https://coronavirus.jhu.edu/map.html) --- ## Spatial Data Structures * Often stored in special data structures + Ex: Shapefile(s) + Consist of geometric objects like points, lines, and polygons -- * Visualizing Spatial Data + Need to draw the map + Add metadata on top of the map + Color in regions + Dots --- ## Polygon Maps * Want to draw various boundaries + EXs: Countries, states, counties, etc * polygon = closed area including the boundaries making up the area <img src="slidesWk10Tu_files/figure-html/unnamed-chunk-8-1.png" width="360" style="display: block; margin: auto;" /> --- ## Polygon Maps * Let's get a polygon map from the `maps` package ```r or_counties <- map_data("county", "oregon") glimpse(or_counties) ``` ``` ## Rows: 1,633 ## Columns: 6 ## $ long <dbl> -117, -117, -117, -118, -118, -118, -118, -118, -118, -118, … ## $ lat <dbl> 44, 44, 44, 44, 44, 44, 44, 44, 44, 44, 44, 44, 44, 44, 44, … ## $ group <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, … ## $ order <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 1… ## $ region <chr> "oregon", "oregon", "oregon", "oregon", "oregon", "oregon", … ## $ subregion <chr> "baker", "baker", "baker", "baker", "baker", "baker", "baker… ``` --- ## Polygon Maps * `lat`: latitude of a corner of the polygon * `long`: longitude of a corner of the polygon * `group`: county ```r glimpse(or_counties) ``` ``` ## Rows: 1,633 ## Columns: 6 ## $ long <dbl> -117, -117, -117, -118, -118, -118, -118, -118, -118, -118, … ## $ lat <dbl> 44, 44, 44, 44, 44, 44, 44, 44, 44, 44, 44, 44, 44, 44, 44, … ## $ group <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, … ## $ order <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 1… ## $ region <chr> "oregon", "oregon", "oregon", "oregon", "oregon", "oregon", … ## $ subregion <chr> "baker", "baker", "baker", "baker", "baker", "baker", "baker… ``` --- ## Polygon Maps ```r or_counties <- map_data("county", "oregon") ggplot(or_counties, aes(long, lat, color = factor(group))) + geom_point(show.legend = FALSE) + coord_quickmap() ``` <img src="slidesWk10Tu_files/figure-html/unnamed-chunk-11-1.png" width="360" style="display: block; margin: auto;" /> --- ## Polygon Maps ```r ggplot(or_counties, aes(long, lat, group = group)) + geom_polygon(fill = "white", colour = "olivedrab") + coord_quickmap() + theme_void() ``` <img src="slidesWk10Tu_files/figure-html/unnamed-chunk-12-1.png" width="360" style="display: block; margin: auto;" /> --- ## Simple Features Maps * Typically map data doesn't come in the lat-long format given in the maps package. * Often use the "simple features" standard produced by the Open Geospatial Consortium + R package: `sf` handles this type of data + Within `ggplot2`: `geom_sf()` and `coord_sf()` work with `sf` --- ## [`tidycensus`](https://walkerke.github.io/tidycensus/articles/spatial-data.html) * Data from the US Census or American Community Survey + Comes in simple features format (using `tigris`) + Need to obtain an [API key](https://api.census.gov/data/key_signup.html) ```r api_key <- "insert key" ``` ```r library(tidycensus) options(tigris_use_cache = TRUE) or <- get_acs(state = "OR", geography = "county", variables = "B19013_001", geometry = TRUE, key = api_key) ``` --- ## [`tidycensus`](https://walkerke.github.io/tidycensus/articles/spatial-data.html) ```r vars <- load_variables(2018, "acs5", cache = FALSE) vars ``` ``` ## # A tibble: 26,997 x 3 ## name label concept ## <chr> <chr> <chr> ## 1 B00001_001 Estimate!!Total UNWEIGHTED SAMPLE COUNT OF THE PO… ## 2 B00002_001 Estimate!!Total UNWEIGHTED SAMPLE HOUSING UNITS ## 3 B01001_001 Estimate!!Total SEX BY AGE ## 4 B01001_002 Estimate!!Total!!Male SEX BY AGE ## 5 B01001_003 Estimate!!Total!!Male!!Under 5… SEX BY AGE ## 6 B01001_004 Estimate!!Total!!Male!!5 to 9 … SEX BY AGE ## 7 B01001_005 Estimate!!Total!!Male!!10 to 1… SEX BY AGE ## 8 B01001_006 Estimate!!Total!!Male!!15 to 1… SEX BY AGE ## 9 B01001_007 Estimate!!Total!!Male!!18 and … SEX BY AGE ## 10 B01001_008 Estimate!!Total!!Male!!20 years SEX BY AGE ## # … with 26,987 more rows ``` --- ## Simple Features Object * Row for each polygon * `geometry` contains the polygon object (borders of the county) ```r or ``` ``` ## Simple feature collection with 36 features and 5 fields ## Geometry type: MULTIPOLYGON ## Dimension: XY ## Bounding box: xmin: -120 ymin: 42 xmax: -120 ymax: 46 ## CRS: 4269 ## First 10 features: ## GEOID NAME variable estimate moe ## 1 41017 Deschutes County, Oregon B19013_001 67043 1740 ## 2 41003 Benton County, Oregon B19013_001 62077 2483 ## 3 41015 Curry County, Oregon B19013_001 48440 4060 ## 4 41061 Union County, Oregon B19013_001 52171 3923 ## 5 41055 Sherman County, Oregon B19013_001 51071 4755 ## 6 41051 Multnomah County, Oregon B19013_001 69176 929 ## 7 41007 Clatsop County, Oregon B19013_001 54886 2779 ## 8 41033 Josephine County, Oregon B19013_001 45616 2104 ## 9 41031 Jefferson County, Oregon B19013_001 53277 3582 ## 10 41039 Lane County, Oregon B19013_001 52426 898 ## geometry ## 1 MULTIPOLYGON (((-122 44, -1... ## 2 MULTIPOLYGON (((-124 44, -1... ## 3 MULTIPOLYGON (((-124 42, -1... ## 4 MULTIPOLYGON (((-119 45, -1... ## 5 MULTIPOLYGON (((-121 45, -1... ## 6 MULTIPOLYGON (((-123 46, -1... ## 7 MULTIPOLYGON (((-124 46, -1... ## 8 MULTIPOLYGON (((-124 42, -1... ## 9 MULTIPOLYGON (((-122 44, -1... ## 10 MULTIPOLYGON (((-124 44, -1... ``` --- ## Polygon Map * `geom_sf()`: creates map from the `geometry` column * `coord_sf()`: controls the projection ```r ggplot(data = or, mapping = aes(geometry = geometry)) + geom_sf() + coord_sf() + theme_void() ``` <img src="slidesWk10Tu_files/figure-html/unnamed-chunk-18-1.png" width="360" style="display: block; margin: auto;" /> --- ## Choropleth Map ```r ggplot(data = or, mapping = aes(geometry = geometry, fill = estimate)) + geom_sf() + coord_sf() + scale_fill_viridis_c() + theme_void() ``` <img src="slidesWk10Tu_files/figure-html/unnamed-chunk-19-1.png" width="360" style="display: block; margin: auto;" /> --- ## Choropleth Map ```r ggplot(data = or, mapping = aes(geometry = geometry)) + geom_sf(aes(fill = estimate)) + coord_sf() + geom_sf_label(aes(label = NAME)) + scale_fill_viridis_c() + theme_void() ``` <img src="slidesWk10Tu_files/figure-html/unnamed-chunk-20-1.png" width="360" style="display: block; margin: auto;" /> --- ```r #devtools::install_github("yutannihilation/ggsflabel") library(ggsflabel) or <- or %>% mutate(NAME = str_replace_all(NAME, "County, Oregon", "")) ggplot(data = or, mapping = aes(geometry = geometry)) + geom_sf(aes(fill = estimate)) + coord_sf() + geom_sf_label_repel(aes(label = NAME), size = 3) + scale_fill_viridis_c() + theme_void() ``` <img src="slidesWk10Tu_files/figure-html/unnamed-chunk-21-1.png" width="576" style="display: block; margin: auto;" /> --- Can add data from another dataset ```r or_cities <- data_frame(city = c("Portland, OR", "Salem, OR", "Eugene, OR", "Gresham, OR"), lat = c(45.523, 44.943, 44.052, 45.498), long = c(-122.676, -123.035, -123.087, -122.431)) ggplot() + geom_sf(data = or, mapping = aes(geometry = geometry, fill = estimate)) + coord_sf() + geom_point(data = or_cities, mapping = aes(x = long, y = lat), color = "white") + scale_fill_viridis_c() + theme_void() ``` <img src="slidesWk10Tu_files/figure-html/unnamed-chunk-22-1.png" width="360" style="display: block; margin: auto;" /> --- ## Coordinate Reference System * Before covering maps, we plotted latitude and longitude on the Cartesian plane. * But the Earth is round (though not a perfect sphere). * **geodetic datum**: Set of assumptions about the shape of the Earth + North American Datum (NAD83) + World Geodetic System (WGS84) * Map projections: Convert from 3-D to 2-D + Area preserving + Shape preserving * Together give us a **coordinate reference system** --- ## Coordinate Reference System ```r or ``` ``` ## Simple feature collection with 36 features and 5 fields ## Geometry type: MULTIPOLYGON ## Dimension: XY ## Bounding box: xmin: -120 ymin: 42 xmax: -120 ymax: 46 ## CRS: 4269 ## First 10 features: ## GEOID NAME variable estimate moe geometry ## 1 41017 Deschutes B19013_001 67043 1740 MULTIPOLYGON (((-122 44, -1... ## 2 41003 Benton B19013_001 62077 2483 MULTIPOLYGON (((-124 44, -1... ## 3 41015 Curry B19013_001 48440 4060 MULTIPOLYGON (((-124 42, -1... ## 4 41061 Union B19013_001 52171 3923 MULTIPOLYGON (((-119 45, -1... ## 5 41055 Sherman B19013_001 51071 4755 MULTIPOLYGON (((-121 45, -1... ## 6 41051 Multnomah B19013_001 69176 929 MULTIPOLYGON (((-123 46, -1... ## 7 41007 Clatsop B19013_001 54886 2779 MULTIPOLYGON (((-124 46, -1... ## 8 41033 Josephine B19013_001 45616 2104 MULTIPOLYGON (((-124 42, -1... ## 9 41031 Jefferson B19013_001 53277 3582 MULTIPOLYGON (((-122 44, -1... ## 10 41039 Lane B19013_001 52426 898 MULTIPOLYGON (((-124 44, -1... ``` --- ## New Data ```r Eruptions <- read_csv("/home/courses/math241s21/Data/GVP_Eruption_Results.csv") select(Eruptions, VolcanoName, StartYear, Latitude, Longitude) %>% glimpse() ``` ``` ## Rows: 11,019 ## Columns: 4 ## $ VolcanoName <chr> "Bulusan", "Alaid", "Turrialba", "San Miguel", "Chill\xe0n… ## $ StartYear <dbl> 2016, 2016, 2016, 2016, 2016, 2016, 2015, 2015, 2015, 2015… ## $ Latitude <dbl> 12.77, 50.86, 10.03, 13.43, -36.86, 1.11, 11.98, 37.73, 12… ## $ Longitude <dbl> 124, 156, -84, -88, -71, 125, -86, 15, -87, 123, 100, 113,… ``` ```r world <- map_data("world") ``` --- ## Static Maps ```r ggplot(data = world, mapping = aes(x = long, y = lat, group = group)) + geom_polygon(fill = "grey", color = "blue") + coord_fixed(1.3) + geom_point(data = Eruptions, aes(Longitude, Latitude), inherit.aes = FALSE, color = "red") + theme_void() ``` <img src="slidesWk10Tu_files/figure-html/unnamed-chunk-25-1.png" width="360" style="display: block; margin: auto;" /> --- ## Raster Maps * Want a specific base layer map: `ggmap` package ```r library(ggmap) aleutian_box <- c(bottom = 50.977, left = -172.164, top = 59.5617, right = -157.1507) aleutian <- get_stamenmap(aleutian_box, maptype = "terrain-background", zoom = 5) aleutian %>% ggmap() ``` <img src="slidesWk10Tu_files/figure-html/unnamed-chunk-26-1.png" width="360" style="display: block; margin: auto;" /> --- ## Raster Maps * Can take a long time to download so not a bad idea to save: ```r save(aleutian, file = "aleutian.RData") ``` * And then load it when you need to: ```r load("aleutian.RData") ``` --- ## Raster Maps ```r aleutian %>% ggmap() + geom_point(data = Eruptions, aes(Longitude, Latitude), inherit.aes = FALSE, color = "red") + theme_void() ``` <img src="slidesWk10Tu_files/figure-html/unnamed-chunk-29-1.png" width="360" style="display: block; margin: auto;" /> --- ## Raster Maps ```r eruption_count <- count(Eruptions, VolcanoName, Latitude, Longitude) %>% filter(Longitude > -172.164, Longitude < -157.1507, Latitude > 50.977, Latitude < 59.5617) eruption_count ``` ``` ## # A tibble: 24 x 4 ## VolcanoName Latitude Longitude n ## <chr> <dbl> <dbl> <int> ## 1 Akutan 54.1 -166. 47 ## 2 Amak 55.4 -163. 3 ## 3 Amukta 52.5 -171. 8 ## 4 Aniakchak 56.9 -158. 15 ## 5 Black Peak 56.6 -159. 1 ## 6 Bogoslof 53.9 -168. 11 ## 7 Carlisle 52.9 -170. 4 ## 8 Cleveland 52.8 -170. 29 ## 9 Dana 55.6 -161. 1 ## 10 Fisher 54.6 -164. 6 ## # … with 14 more rows ``` --- ## Raster Maps ```r aleutian %>% ggmap() + geom_point(data = eruption_count, aes(Longitude, Latitude, size = n), inherit.aes = FALSE, color = "red", alpha = 0.5) + theme_void() ``` <img src="slidesWk10Tu_files/figure-html/unnamed-chunk-31-1.png" width="360" style="display: block; margin: auto;" /> --- ## Raster Maps ```r aleutian %>% ggmap() + geom_point(data = eruption_count, aes(Longitude, Latitude, color = n), inherit.aes = FALSE) + theme_void() + scale_color_viridis_c(direction = -1, option = "A") ``` <img src="slidesWk10Tu_files/figure-html/unnamed-chunk-32-1.png" width="360" style="display: block; margin: auto;" /> --- * But what about dynamic zooming? * What about clickable labels? ```r aleutian %>% ggmap() + geom_point(data = eruption_count, aes(Longitude, Latitude, color = n), inherit.aes = FALSE) + theme_void() + scale_color_viridis_c(direction = -1, option = "A") ``` <img src="slidesWk10Tu_files/figure-html/unnamed-chunk-33-1.png" width="360" style="display: block; margin: auto;" /> --- ## Interactive Maps with `leaflet` * `leaflet`: Built on a JavaScript library * Syntax similar to `ggplot2` and `dplyr` * Interactive zooming * Embed interactive map into RMarkdown documents (HTML output)
--- ## Interactive Maps with `leaflet` * First Layer: `leaflet()` + Creates the map widget * Second layer: `addTiles()` + Creates a base map using [OpenStreetMap](https://www.openstreetmap.org) * Additional layers: `addMarkers()`, `addPolygons()` ```r library(leaflet) leaflet() %>% addTiles() %>% addCircleMarkers(lng = ~Longitude, lat = ~Latitude, data = eruption_count, radius = ~sqrt(n), stroke = FALSE, fillOpacity = 0.9) ``` --- ## Interactive Maps with `leaflet` ```r leaflet() %>% addTiles() %>% addCircleMarkers(lng = ~Longitude, lat = ~Latitude, data = eruption_count, radius = ~sqrt(n), stroke = FALSE, fillOpacity = 0.9) ```
--- ## Clusters ```r leaflet() %>% setView(lng = -160, lat = 55, zoom = 4) %>% addTiles() %>% addCircleMarkers(lng = ~Longitude, lat = ~Latitude, data = eruption_count, clusterOptions = markerClusterOptions()) ```
--- ## Zooming Constraints ```r leaflet(options = leafletOptions(minZoom = 3, maxZoom = 7)) %>% addTiles() %>% addCircleMarkers(lng = ~Longitude, lat = ~Latitude, data = eruption_count, radius = ~sqrt(n)) ```
--- ## Other Markers ```r leaflet() %>% setView(lng = -160, lat = 55, zoom = 4) %>% addTiles() %>% addMarkers(lng = ~Longitude, lat = ~Latitude, data = eruption_count, label = ~as.character(n)) ```
--- ## Custom icons ```r volcano <- makeIcon( iconUrl = "https://openclipart.org/image/800px/svg_to_png/168263/volcano-mountain.png", iconWidth = 10, iconHeight = 20) volcano ``` ``` ## $iconUrl ## [1] "https://openclipart.org/image/800px/svg_to_png/168263/volcano-mountain.png" ## ## $iconWidth ## [1] 10 ## ## $iconHeight ## [1] 20 ## ## attr(,"class") ## [1] "leaflet_icon" ``` --- ```r leaflet() %>% setView(lng = -160, lat = 55, zoom = 4) %>% addTiles() %>% addMarkers(lng = ~Longitude, lat = ~Latitude, data = eruption_count, label = ~as.character(n), icon = volcano) ```
--- ## Pop-Ups ```r content <- paste("<b>", eruption_count$VolcanoName, "</b></br>", "Number of eruptions:", eruption_count$n) ``` --- ```r leaflet() %>% setView(lng = -160, lat = 55, zoom = 4) %>% addTiles() %>% addMarkers(lng = ~Longitude, lat = ~Latitude, data = eruption_count, popup = content) ```
--- ## Other Tiles ```r library(leaflet.extras) leaflet() %>% setView(lng = -122.630752, lat = 45.480542, zoom = 15) %>% addProviderTiles(providers$Stamen.Watercolor) ```
--- ## Other Tiles ```r leaflet() %>% setView(lng = -122.630752, lat = 45.480542, zoom = 15) %>% addProviderTiles(providers$CartoDB.Positron) %>% addMiniMap() ```
--- ## Other Tiles ```r leaflet() %>% setView(lng = -98.58, lat = 39.82, zoom = 3) %>% addProviderTiles("NASAGIBS.ViirsEarthAtNight2012") ```
--- ## Interactive Choropleth Map ```r pal <- colorNumeric(palette = "viridis", domain = or$estimate) or %>% leaflet() %>% addTiles() %>% addPolygons(popup = ~NAME, color = ~pal(estimate), stroke = FALSE, fillOpacity = 0.9) ```
--- ## Interactive Choropleth Map ```r or %>% sf::st_transform(crs = "+init=epsg:4326") %>% leaflet() %>% addTiles() %>% addPolygons(popup = ~NAME, fillColor = ~pal(estimate), stroke = FALSE, fillOpacity = 0.9) %>% addLegend("bottomright", pal = pal, values = ~estimate, title = "Median Income", opacity = 1) ```
--- ## Want More User Interactivity? -- We are getting there! Have to learn `shiny`. --- ## Final Project Idea * Dig into manipulating spatial data with the `sf` package: https://github.com/rstudio/cheatsheets/blob/master/sf.pdf * All things spatial data in R: https://geocompr.robinlovelace.net/index.html