Maps

.leftcol30[
<center>
<img src="https://raw.githubusercontent.com/emse-eda-gwu/2021-Spring/master/images/eda_hex_sticker.png" width=250>
</center>
]

### <svg style="height:0.8em;top:.04em;position:relative;fill:white;" viewBox="0 0 512 512"><path d="M496 128v16a8 8 0 0 1-8 8h-24v12c0 6.627-5.373 12-12 12H60c-6.627 0-12-5.373-12-12v-12H24a8 8 0 0 1-8-8v-16a8 8 0 0 1 4.941-7.392l232-88a7.996 7.996 0 0 1 6.118 0l232 88A8 8 0 0 1 496 128zm-24 304H40c-13.255 0-24 10.745-24 24v16a8 8 0 0 0 8 8h464a8 8 0 0 0 8-8v-16c0-13.255-10.745-24-24-24zM96 192v192H60c-6.627 0-12 5.373-12 12v20h416v-20c0-6.627-5.373-12-12-12h-36V192h-64v192h-64V192h-64v192h-64V192H96z"/></svg> EMSE 4575: Exploratory Data Analysis
### <svg style="height:0.8em;top:.04em;position:relative;fill:white;" viewBox="0 0 448 512"><path d="M224 256c70.7 0 128-57.3 128-128S294.7 0 224 0 96 57.3 96 128s57.3 128 128 128zm89.6 32h-16.7c-22.2 10.2-46.9 16-72.9 16s-50.6-5.8-72.9-16h-16.7C60.2 288 0 348.2 0 422.4V464c0 26.5 21.5 48 48 48h352c26.5 0 48-21.5 48-48v-41.6c0-74.2-60.2-134.4-134.4-134.4z"/></svg> John Paul Helveston
### <svg style="height:0.8em;top:.04em;position:relative;fill:white;" viewBox="0 0 448 512"><path d="M0 464c0 26.5 21.5 48 48 48h352c26.5 0 48-21.5 48-48V192H0v272zm320-196c0-6.6 5.4-12 12-12h40c6.6 0 12 5.4 12 12v40c0 6.6-5.4 12-12 12h-40c-6.6 0-12-5.4-12-12v-40zm0 128c0-6.6 5.4-12 12-12h40c6.6 0 12 5.4 12 12v40c0 6.6-5.4 12-12 12h-40c-6.6 0-12-5.4-12-12v-40zM192 268c0-6.6 5.4-12 12-12h40c6.6 0 12 5.4 12 12v40c0 6.6-5.4 12-12 12h-40c-6.6 0-12-5.4-12-12v-40zm0 128c0-6.6 5.4-12 12-12h40c6.6 0 12 5.4 12 12v40c0 6.6-5.4 12-12 12h-40c-6.6 0-12-5.4-12-12v-40zM64 268c0-6.6 5.4-12 12-12h40c6.6 0 12 5.4 12 12v40c0 6.6-5.4 12-12 12H76c-6.6 0-12-5.4-12-12v-40zm0 128c0-6.6 5.4-12 12-12h40c6.6 0 12 5.4 12 12v40c0 6.6-5.4 12-12 12H76c-6.6 0-12-5.4-12-12v-40zM400 64h-48V16c0-8.8-7.2-16-16-16h-32c-8.8 0-16 7.2-16 16v48H160V16c0-8.8-7.2-16-16-16h-32c-8.8 0-16 7.2-16 16v48H48C21.5 64 0 85.5 0 112v48h448v-48c0-26.5-21.5-48-48-48z"/></svg> March 31, 2021
]

---

# Quiz 4

---

## 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 12: .fancy[Maps]

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

### BREAK

## 3. Projections

---

# Week 12: .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

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

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

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

.leftcol40[.code70[
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

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

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

.leftcol60[.code70[
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()`

.leftcol60[.code70[
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/0gw/EDA/2021-Spring/class/12-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()`

```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()`

.leftcol60[.code70[
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/0gw/EDA/2021-Spring/class/12-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')
```
]]

---
class: inverse

## Your turn

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

<br>

---

# Week 12: .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)

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

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

<center>
<img src="images/crime_map.jpeg">
</center>
]

---

# A choropleth alternative: hex maps

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

<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

```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')
```
]]

---

# 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\)`

---

---

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

.font90[
Create the following map of squirrels in NYC's Central Park using the following 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

---

# Week 12: .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()
```

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

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

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

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

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

---

## Common Projections

ggplot layer:

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

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

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'))
```

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

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

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

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

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

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

---

## 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")
```
]]

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

---

**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>
<img src="images/plots/sf_us_coffee_albers_bad.png">
<center>
]

---

## 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 x 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")
```
]]

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

---

## Mapping data to projections - points

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

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

---

## .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))
```
]]

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