Maps

]

# Week 11: .fancy[Maps]

### <svg aria-hidden="true" role="img" viewBox="0 0 512 512" style="height:1em;width:1em;vertical-align:-0.125em;margin-left:auto;margin-right:auto;font-size:inherit;fill:white;overflow:visible;position:relative;"><path d="M243.4 2.6l-224 96c-14 6-21.8 21-18.7 35.8S16.8 160 32 160v8c0 13.3 10.7 24 24 24H456c13.3 0 24-10.7 24-24v-8c15.2 0 28.3-10.7 31.3-25.6s-4.8-29.9-18.7-35.8l-224-96c-8.1-3.4-17.2-3.4-25.2 0zM128 224H64V420.3c-.6 .3-1.2 .7-1.8 1.1l-48 32c-11.7 7.8-17 22.4-12.9 35.9S17.9 512 32 512H480c14.1 0 26.5-9.2 30.6-22.7s-1.1-28.1-12.9-35.9l-48-32c-.6-.4-1.2-.7-1.8-1.1V224H384V416H344V224H280V416H232V224H168V416H128V224zm128-96c-17.7 0-32-14.3-32-32s14.3-32 32-32s32 14.3 32 32s-14.3 32-32 32z"/></svg> EMSE 4572: Exploratory Data Analysis
### <svg aria-hidden="true" role="img" viewBox="0 0 448 512" style="height:1em;width:0.88em;vertical-align:-0.125em;margin-left:auto;margin-right:auto;font-size:inherit;fill:white;overflow:visible;position:relative;"><path d="M272 304h-96C78.8 304 0 382.8 0 480c0 17.67 14.33 32 32 32h384c17.67 0 32-14.33 32-32C448 382.8 369.2 304 272 304zM48.99 464C56.89 400.9 110.8 352 176 352h96c65.16 0 119.1 48.95 127 112H48.99zM224 256c70.69 0 128-57.31 128-128c0-70.69-57.31-128-128-128S96 57.31 96 128C96 198.7 153.3 256 224 256zM224 48c44.11 0 80 35.89 80 80c0 44.11-35.89 80-80 80S144 172.1 144 128C144 83.89 179.9 48 224 48z"/></svg> John Paul Helveston
### <svg aria-hidden="true" role="img" viewBox="0 0 448 512" style="height:1em;width:0.88em;vertical-align:-0.125em;margin-left:auto;margin-right:auto;font-size:inherit;fill:white;overflow:visible;position:relative;"><path d="M152 64H296V24C296 10.75 306.7 0 320 0C333.3 0 344 10.75 344 24V64H384C419.3 64 448 92.65 448 128V448C448 483.3 419.3 512 384 512H64C28.65 512 0 483.3 0 448V128C0 92.65 28.65 64 64 64H104V24C104 10.75 114.7 0 128 0C141.3 0 152 10.75 152 24V64zM48 448C48 456.8 55.16 464 64 464H384C392.8 464 400 456.8 400 448V192H48V448z"/></svg> November 09, 2022

]

---

## Today's data

```r
milk_production <- read_csv(here::here('data', 'milk_production.csv'))
us_coffee_shops <- read_csv(here::here('data', 'us_coffee_shops.csv'))
```

## New packages:

```r
install.packages('maps')
install.packages('mapproj')
install.packages('sf')
install.packages('rgeos')
install.packages('rnaturalearth')
devtools::install_github("ropensci/rnaturalearthhires")
devtools::install_github("ropensci/rnaturalearthdata")
```

---

# Week 11: .fancy[Maps]

## 1. Plotting maps
## 2. Adding data to maps

### BREAK

## 3. Projections

---

# Week 11: .fancy[Maps]

## 1. .orange[Plotting maps]
## 2. Adding data to maps

### BREAK

## 3. Projections

---

# How to make a map

#### Step 1: Load a shape file

a. Use a library

b. Read in a shape file

#### Step 2: Plot the shape file

a. Polygon data: `geom_polygon()`

b. Simple Features data: `geom_sf()`

---

## Polygon maps

Get the "World" shape file

```r
library(ggplot2)

*world <- map_data("world")
head(world)
```

```
#>        long      lat group order region subregion
#> 1 -69.89912 12.45200     1     1  Aruba      <NA>
#> 2 -69.89571 12.42300     1     2  Aruba      <NA>
#> 3 -69.94219 12.43853     1     3  Aruba      <NA>
#> 4 -70.00415 12.50049     1     4  Aruba      <NA>
#> 5 -70.06612 12.54697     1     5  Aruba      <NA>
#> 6 -70.05088 12.59707     1     6  Aruba      <NA>
```

]]

---

## Polygon maps

Get the "World" shape file

```r
library(ggplot2)

*world <- map_data("world")
```

Make the plot with `geom_polygon()`

```r
ggplot(world) +
*   geom_polygon(aes(x = long, y = lat, group = group),
                 fill = "grey90", color = "grey60")
```

]]

]

---

## Polygon maps

Get the "US States" shape file

```r
library(ggplot2)

*us_states <- map_data("state")
```

Make the plot with `geom_polygon()`

```r
ggplot(us_states) +
*   geom_polygon(aes(x = long, y = lat, group = group),
                 fill = "grey90", color = "grey60")
```

]]

]

---

### Simple Features (sf) maps

Library data from [Natural Earth](http://www.naturalearthdata.com/downloads/50m-cultural-vectors/)

```r
library(rnaturalearth)
library(rnaturalearthdata)

*world <- ne_countries(
* scale = "medium",
* returnclass = "sf")

world %>%
    select(name, geometry) %>%
    head()
```

]]

```
#> Simple feature collection with 6 features and 1 field
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -73.36621 ymin: -22.40205 xmax: 109.4449 ymax: 41.9062
#> CRS:           +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
#>        name                       geometry
#> 0  Zimbabwe MULTIPOLYGON (((31.28789 -2...
#> 1    Zambia MULTIPOLYGON (((30.39609 -1...
#> 2     Yemen MULTIPOLYGON (((53.08564 16...
#> 3   Vietnam MULTIPOLYGON (((104.064 10....
#> 4 Venezuela MULTIPOLYGON (((-60.82119 9...
#> 5   Vatican MULTIPOLYGON (((12.43916 41...
```

]]

---

## Simple Features (sf) maps

Get the "World" shape file

```r
library(rnaturalearth)
library(rnaturalearthdata)

world <- ne_countries(
  scale = "medium",
  returnclass = "sf")
```

Make the plot with `geom_sf()`

```r
library(sf)

ggplot(data = world) +
*   geom_sf(fill = "grey90", color = "grey60")
```

]]

]

---

## Simple Features (sf) maps

Get the "US States" shape file

```r
library(rnaturalearth)
library(rnaturalearthdata)

*us_states <- ne_states(
*   country = 'united states of america',
    returnclass = 'sf')
```

]]

---

## Simple Features (sf) maps

Get the "US States" shape file

```r
library(rnaturalearth)
library(rnaturalearthdata)

us_states <- ne_states(
    country = 'united states of america',
    returnclass = 'sf')
```

Make the plot with `geom_sf()`

```r
library(sf)

ggplot(data = us_states) +
*   geom_sf(fill = "grey90", color = "grey60")
```

]]

]

---

## Simple Features (sf) maps

Get the **Continental** "US States" shape file

```r
library(rnaturalearth)
library(rnaturalearthdata)

us_states_cont <- ne_states(
    country = 'united states of america',
    returnclass = 'sf') %>%
*   filter(! name %in% c('Alaska', 'Hawaii'))
```

Make the plot with `geom_sf()`

```r
library(sf)

ggplot(data = us_states_cont) +
*   geom_sf(fill = "grey90", color = "grey60")
```

]]

]

---

## The `maps` package

Includes data on:
- World: world, world.cities, lakes
- US: states, county, state, usa
- France: france
- Italy: italy
- New zealand: nz

Example:

```r
library(maps)

*us_counties <- st_as_sf(
*   map("county", plot = FALSE, fill = TRUE))

ggplot(data = us_counties) +
*   geom_sf(fill = 'grey90', color = 'grey60')
```

]]

]

---

## Simple Features (sf) maps: `st_read()`

Read in the "World" shape file from [Natural Earth](http://www.naturalearthdata.com/downloads/50m-cultural-vectors/)

```r
library(sf)

*world <- st_read(here::here(
    'data', 'natural_earth_countries',
    'ne_50m_admin_0_countries.shp')) %>%
    clean_names()
```

```
#> Reading layer `ne_50m_admin_0_countries' from data source `/Users/jhelvy/gh/teaching/EDA/2022-Fall/class/11-maps/data/natural_earth_countries/ne_50m_admin_0_countries.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 241 features and 94 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -180 ymin: -89.99893 xmax: 180 ymax: 83.59961
#> Geodetic CRS:  WGS 84
```

]]

---

## Simple Features (sf) maps: `st_read()`

Read in the "World" shape file

```r
library(sf)

*world <- st_read(here::here(
    'data', 'natural_earth_countries',
    'ne_50m_admin_0_countries.shp')) %>%
    clean_names()
```

```r
ggplot(data = world) +
*   geom_sf(fill = "grey90", color = "grey60")
```

]]

]

---

## Simple Features (sf) maps: `st_read()`

Read in the "Central Park" shape file [[source]](https://github.com/malcolmbarrett/designing_ggplots)

```r
library(sf)

central_park <- st_read(here::here(
    'data', 'central_park', 'CentralPark.shp'))
```

```
#> Reading layer `CentralPark' from data source `/Users/jhelvy/gh/teaching/EDA/2022-Fall/class/11-maps/data/central_park/CentralPark.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 2550 features and 6 fields
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: -73.99249 ymin: 40.7625 xmax: -73.93922 ymax: 40.8051
#> Geodetic CRS:  WGS 84
```

```r
ggplot(data = central_park) +
    geom_sf(color = 'grey75')
```

]]

]

---

## Your turn

Use the **rnaturalearth** library to extract and plot the shape files for China and Africa:

<br>

]

]

---

# Week 11: .fancy[Maps]

## 1. Plotting maps
## 2. .orange[Adding data to maps]

### BREAK

## 3. Projections

<!--
Squirrels and dogs plots:
https://designing-ggplots.netlify.com/#71

Hexmaps:
- Use https://www.r-graph-gallery.com/hexbin-map.html
- https://www.r-graph-gallery.com/328-hexbin-map-of-the-usa.html
- Practice with wind data?
-->

---

# First rule of adding data to maps:

# .red[Do you need to make a map?]

---

---

# 1. .orange[Choropleth maps]
# 2. Point maps

---

background-color: #f7fcfd

[Choropleth](https://en.wikipedia.org/wiki/Choropleth_map) - from Greek:

- χῶρος "choros" (area/region)
- πλῆθος "plethos" (multitude)

.font70[From reddit: [r/dataisbeautiful by u/GaudyAudi](https://www.reddit.com/r/dataisbeautiful/comments/fgq3d3/oc_us_state_abbreviation_conventions/)]

---

# Choropleth maps are easily misleading

---

## Number of events != Number of events **per capita**

---

## Manipulating fill scale produces wildly different maps

---

## Manipulating fill scale produces wildly different maps

<br>

<br><br>
.left[.font70[Source: [New York Times](https://www.nytimes.com/interactive/2020/04/02/us/coronavirus-social-distancing.html)]]

]

]

---

### Land doesn't vote - people vote

By [David Zumbach](https://twitter.com/DavidZumbach/status/1344547411985911808?s=19)

---

# Land doesn't vote - people vote

Election maps from: http://www-personal.umich.edu/~mejn/election/2016/

]

]

]

---

---

## Easy to lie with fake news

Bottom map is actually [this map](http://www-personal.umich.edu/~mejn/election/2012/countymaprb1024.png) of the 2012 election

---

## (here is what actual crime rates look like)

2016 Election map [[source]](http://www-personal.umich.edu/~mejn/election/2016/)

<center>
<img src="images/election_2016_county.png">
</center>
]

2014 Crime map [[source]](https://www.washingtonpost.com/graphics/national/crime-rates-by-county/)

]

---

# A choropleth alternative: hex maps

1994 Simpson Diversity Index in US Schools

<center>
<img src="images/Hexagons_SchoolDiversity_hex_1994.png">
</center>
.font80[https://github.com/malcolmbarrett/designing_ggplots]

]

2016 Electoral College

<center>
<img src="images/election_2016_chartograph.png" width=450>
</center>
https://fivethirtyeight.com/

]

---

## How to make a choropleth map

---

## How to make a choropleth map

Get the "fill" data

```r
milk_2017 <- milk_production %>%
    filter(year == 2017) %>%
    select(name = state, milk_produced) %>%
    mutate(milk_produced = milk_produced / 10^9)
```

Get the "map" data

```r
us_states <- ne_states(
    country = 'united states of america',
    returnclass = 'sf') %>%
    filter(! name %in% c('Alaska', 'Hawaii')) %>%
*   left_join(milk_2017, by = 'name')
```

]]

```r
us_states %>%
    select(name, milk_produced) %>%
    head()
```

```
#> Simple feature collection with 6 features and 2 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -124.7346 ymin: 41.69681 xmax: -82.4146 ymax: 49.36949
#> CRS:           +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
#>           name milk_produced                       geometry
#> 1    Minnesota         9.864 MULTIPOLYGON (((-95.16057 4...
#> 2   Washington         6.526 MULTIPOLYGON (((-122.6533 4...
#> 3        Idaho        14.627 MULTIPOLYGON (((-117.0382 4...
#> 4      Montana         0.288 MULTIPOLYGON (((-116.0482 4...
#> 5 North Dakota         0.345 MULTIPOLYGON (((-104.0476 4...
#> 6     Michigan        11.231 MULTIPOLYGON (((-84.4913 46...
```

]

---

## How to make a choropleth map

```r
ggplot(us_states) +
* geom_sf(aes(fill = milk_produced)) +
  scale_fill_viridis(
    option = "plasma",
    limits = c(0, 40)) +
  theme_void(base_size = 15) +
  theme(legend.position = 'bottom') +
  labs(fill = 'Milk produced\n(billions lbs)',
       title = 'Milk Production by State in 2017')
```

]]

]

---

## How to make a choropleth map

```r
ggplot(us_states) +
  geom_sf(aes(fill = milk_produced)) +
  scale_fill_viridis(
*   trans = 'sqrt',
    option = "plasma",
    limits = c(0, 40)) +
  theme_void(base_size = 15) +
  theme(legend.position = 'bottom') +
  labs(fill = 'Milk produced\n(billions lbs)',
       title = 'Milk Production by State in 2017')
```

]]

Non-linear scale:
<center>
<img src="images/plots/sf_us_milk_2017_quad.png">
<center>

]

---

# 1. Choropleth maps
# 2. .orange[Point maps]

---

# Points as locations

---

# Points encoding a variable

---

# For point size, use **area**, not radius

# `\(Area = \pi r^2\)`

## Radius

]

## Area

]

---

## Radius

]

## Area

]

---

## How to add points to a map

---

## How to add points to a map

Load the continental US shape file

```r
us_states_cont <- ne_states(
    country = 'united states of america',
    returnclass = 'sf') %>%
    filter(! name %in% c('Alaska', 'Hawaii'))
```

Read in the coffee shop data

```r
us_coffee_shops  <- read_csv(here::here(
    'data', 'us_coffee_shops.csv'))

# Only keep data in continental US
us_coffee_shops <- us_coffee_shops %>%
    filter(lat > 22,    lat < 50,
           long > -150, long < -66)
```

]]

```r
head(us_coffee_shops)
```

```
#> # A tibble: 6 × 8
#>   name             lat  long unique_id city               state_abb zip   state    
#>   <chr>          <dbl> <dbl>     <dbl> <chr>              <chr>     <chr> <chr>    
#> 1 Baskin Robbins  40.8 -73.4   3304448 Huntington Station NY        11746 New York 
#> 2 Baskin Robbins  42.1 -88.0  11342048 Rolling Meadows    IL        60008 Illinois 
#> 3 Baskin Robbins  34.0 -84.5   3304169 Marietta           GA        30066 Georgia  
#> 4 Baskin Robbins  29.8 -95.6   3304006 Houston            TX        77079 Texas    
#> 5 Baskin Robbins  36.4 -89.5   3303959 Tiptonville        TN        38079 Tennessee
#> 6 Baskin Robbins  40.7 -73.6   3304507 Merrick            NY        11566 New York
```

]

---

## How to add points to a map

Plot coffee shop locations over map

```r
ggplot() +
  geom_sf(data = us_states_cont) +
* geom_point(
*   data = us_coffee_shops,
*   aes(x = long, y = lat, color = name),
*   size = 0.3) +
  theme_void(base_size = 15) +
  theme(legend.position = 'bottom') +
  guides(color = guide_legend(
    # Move legend title to top
    title.position = "top",
    # Increase legend point size
    override.aes = list(size = 3))) +
  labs(color = 'Coffee shop',
       title = 'Coffee Shops in the US')
```

]]

]

---

## Your turn

Create this map of squirrels in NYC's Central Park using this data from the [Squirrel Census](https://www.thesquirrelcensus.com/):

- The `CentralPark.shp` file in the `data/central_park` folder.
- The `nyc_squirrels.csv` file in the `data` folder.

Hint: The color is mapped to the `primary_fur_color` variable

(More about the Squirrel Census [here](https://www.nytimes.com/interactive/2020/01/08/nyregion/central-park-squirrel-census.html))

]

---

### Intermission

---

# Quiz 4

### Link is in the #class channel

]

]

---

# Week 11: .fancy[Maps]

## 1. Plotting maps
## 2. Adding data to maps

### BREAK

## 3. .orange[Projections]

<!--
Content from Andrew's slides:
https://datavizf17.classes.andrewheiss.com/class/07-class/

Centroids with labels:
https://www.r-spatial.org/r/2018/10/25/ggplot2-sf-2.html
-->

---

## What's a map projection?

---

# What is the best projection?...it depends

<br>

# 1. [Compare projections](http://metrocosm.com/compare-map-projections.html)

<br>

# 2. [Compare country sizes](https://thetruesize.com/)

---

## Using projections

To modify the projection of a map, use `coord_sf(crs = st_crs(XXXX))`

```r
world <- ne_countries(scale = "medium", returnclass = "sf")
```

```r
ggplot(data = world) +
    geom_sf()
```

]]

```r
ggplot(data = world) +
    geom_sf() +
*   coord_sf(crs = "ESRI:54030")
```

]]

```r
ggplot(data = world) +
    geom_sf() +
*   coord_sf(crs = "ESRI:54009")
```

]]

---

## Common Projections

ggplot layer:

```r
coord_sf(crs = "ESRI:XXXX")
```

**World**

Code  | Projection
------|---------------------------------------------------
**`"ESRI:54030"`** | **Robinson**
`"ESRI:54002"` | Equidistant cylindrical
`"ESRI:54004"` | Mercator
`"ESRI:54008"` | Sinusoidal
`"ESRI:54009"` | Mollweide

]

**United States**

Code   | Projection
-------|---------------------------------------------------
**`"ESRI:102003"`** | **Albers**
`"ESRI:102004"` | Lambert Conformal Conic
`4269`          | NAD 83

]

---

## US projections

```r
us_states_cont <- ne_states(country = 'united states of america',
    returnclass = 'sf') %>%
    filter(! name %in% c('Alaska', 'Hawaii'))
```

NAD 83 projection

```r
ggplot(data = world) +
    geom_sf() +
*   coord_sf(crs = 4269)
```

]]

Mercator

```r
ggplot(data = world) +
    geom_sf() +
*   coord_sf(crs = "ESRI:54004")
```

]]

Albers

```r
ggplot(data = us_states_cont) +
    geom_sf() +
*   coord_sf(crs = "ESRI:102003")
```

]]

---

## Mapping data to projections - choropleth map

```r
milk_2017 <- milk_production %>%
  filter(year == 2017) %>%
  select(name = state, milk_produced) %>%
  mutate(milk_produced = milk_produced / 10^9)

us_states <- ne_states(
  country = 'united states of america',
  returnclass = 'sf') %>%
  filter(! name %in% c('Alaska', 'Hawaii')) %>%
* left_join(milk_2017, by = 'name')

ggplot(us_states) +
* geom_sf(aes(fill = milk_produced)) +
  scale_fill_viridis(
    option = "plasma",
    limits = c(0, 40)) +
  theme_void(base_size = 15) +
  theme(legend.position = 'bottom') +
  labs(fill = 'Milk produced\n(billions lbs)',
       title = 'Milk Production by State in 2017')
```

]]

]

---

## Mapping data to projections - choropleth map

```r
milk_2017 <- milk_production %>%
  filter(year == 2017) %>%
  select(name = state, milk_produced) %>%
  mutate(milk_produced = milk_produced / 10^9)

us_states <- ne_states(
  country = 'united states of america',
  returnclass = 'sf') %>%
  filter(! name %in% c('Alaska', 'Hawaii')) %>%
  left_join(milk_2017, by = 'name')

ggplot(us_states) +
  geom_sf(aes(fill = milk_produced)) +
  scale_fill_viridis(
    option = "plasma",
    limits = c(0, 40)) +
  theme_void(base_size = 15) +
  theme(legend.position = 'bottom') +
  labs(fill = 'Milk produced\n(billions lbs)',
       title = 'Milk Production by State in 2017') +
* coord_sf(crs = "ESRI:102003")
```

]]

<center>
Albers Projection

]

---

**Mapping data to projections - points**

```r
us_states_cont <- ne_states(
  country = 'united states of america',
  returnclass = 'sf') %>%
  filter(! name %in% c('Alaska', 'Hawaii'))

us_coffee_shops <- us_coffee_shops %>%
  filter(lat > 22,    lat < 50,
         long > -150, long < -66)

ggplot() +
  geom_sf(data = us_states_cont) +
* geom_point(
*   data = us_coffee_shops,
*   aes(x = long, y = lat, color = name),
*   size = 0.3) +
  theme_void(base_size = 15) +
  theme(legend.position = 'bottom') +
  guides(color = guide_legend(
    # Move legend title to top
    title.position = "top",
    # Increase legend point size
    override.aes = list(size = 3))) +
  labs(color = 'Coffee shop',
       title = 'Coffee Shops in the US')
```

]]

]

---

**Mapping data to projections - points**

```r
us_states_cont <- ne_states(
  country = 'united states of america',
  returnclass = 'sf') %>%
  filter(! name %in% c('Alaska', 'Hawaii'))

us_coffee_shops <- us_coffee_shops %>%
  filter(lat > 22,    lat < 50,
         long > -150, long < -66)

ggplot() +
  geom_sf(data = us_states_cont) +
  geom_point(
    data = us_coffee_shops,
    aes(x = long, y = lat, color = name),
    size = 0.3) +
  theme_void(base_size = 15) +
  theme(legend.position = 'bottom') +
  guides(color = guide_legend(
    # Move legend title to top
    title.position = "top",
    # Increase legend point size
    override.aes = list(size = 3))) +
  labs(color = 'Coffee shop',
       title = 'Coffee Shops in the US') +
* coord_sf(crs = "ESRI:102003")
```

]]

### .center[.red[Fail!]]

]

---

## Mapping data to projections - points

First match `us_coffee_shops` crs to `us_states_cont`

```r
us_states_cont <- ne_states(
    country = 'united states of america',
    returnclass = 'sf') %>%
    filter(! name %in% c('Alaska', 'Hawaii'))

us_coffee_shops <- us_coffee_shops %>%
    filter(lat > 22,    lat < 50,
           long > -150, long < -66)

*us_coffee_shops_sf <- st_as_sf(us_coffee_shops,
*   coords = c("long", "lat"),
*   crs = st_crs(us_states_cont))
```

]]

```r
head(us_coffee_shops_sf)
```

```
#> Simple feature collection with 6 features and 6 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: -95.60337 ymin: 29.77006 xmax: -73.41174 ymax: 42.05709
#> CRS:           +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
#> # A tibble: 6 × 7
#>   name           unique_id city               state_abb zip   state                 geometry
#>   <chr>              <dbl> <chr>              <chr>     <chr> <chr>              <POINT [°]>
#> 1 Baskin Robbins   3304448 Huntington Station NY        11746 New York  (-73.41174 40.81772)
#> 2 Baskin Robbins  11342048 Rolling Meadows    IL        60008 Illinois   (-88.0044 42.05709)
#> 3 Baskin Robbins   3304169 Marietta           GA        30066 Georgia   (-84.52889 34.02124)
#> 4 Baskin Robbins   3304006 Houston            TX        77079 Texas     (-95.60337 29.77006)
#> 5 Baskin Robbins   3303959 Tiptonville        TN        38079 Tennessee (-89.46582 36.37781)
#> 6 Baskin Robbins   3304507 Merrick            NY        11566 New York  (-73.55707 40.67515)
```

]

---

## Mapping data to projections - points

Plot coffee shop locations over map with `geom_sf()`

```r
ggplot() +
  geom_sf(data = us_states_cont) +
* geom_sf(
*   data = us_coffee_shops_sf,
*   aes(color = name),
*   size = 0.3) +
  theme_void(base_size = 15) +
  theme(legend.position = 'bottom') +
  guides(color = guide_legend(
    # Move legend title to top
    title.position = "top",
    # Increase legend point size
    override.aes = list(size = 3))) +
  labs(fill = 'Coffee shop',
       title = 'Coffee Shops in the US')
```

]]

]

---

## Mapping data to projections - points

Plot coffee shop locations over map with `geom_sf()`

```r
ggplot() +
  geom_sf(data = us_states_cont) +
  geom_sf(
    data = us_coffee_shops_sf,
    aes(color = name),
    size = 0.3) +
  theme_void(base_size = 15) +
  theme(legend.position = 'bottom') +
  guides(color = guide_legend(
    # Move legend title to top
    title.position = "top",
    # Increase legend point size
    override.aes = list(size = 3))) +
  labs(fill = 'Coffee shop',
       title = 'Coffee Shops in the US') +
* coord_sf(crs = "ESRI:102003")
```

]]

<center>
Albers Projection

]

---

## Mapping data to projections - points

Plot coffee shop locations over map with `geom_sf()`

]]

<center>
LCC Projection

]

---

## .center[Your turn]

Use the `internet_users_country.csv` data and the `world` data frame from the **rnaturalearth** library to create these two versions of internet access by country in 2015.

Hints:

- The `iso_a3` variable in the `worlds` data frame corresponds with the `code` variable in the `internet_users_country.csv` data frame (use this for joining).
- Use `scale_fill_gradient()` to fill the color:

```r
scale_fill_gradient(
    low = "#e7e1ef",
    high = "#dd1c77",
    na.value = "grey70",
    limits = c(0, 100))
```

]]

<center>
Robinson Projection
<img src="images/plots/sf_world_internet_robinson.png">

Mercator Projection
<img src="images/plots/sf_world_internet_mercator.png">
<center>

]

---

# Extra practice

---

## Your turn

Use the `us_states_cont` data frame and the `state_abbs` data frame to create a labeled map of the U.S.: