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-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.1V224H384V416H344V224H280V416H232V224H168V416H128V224zM256 64a32 32 0 1 1 0 64 32 32 0 1 1 0-64z"/></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="M304 128a80 80 0 1 0 -160 0 80 80 0 1 0 160 0zM96 128a128 128 0 1 1 256 0A128 128 0 1 1 96 128zM49.3 464H398.7c-8.9-63.3-63.3-112-129-112H178.3c-65.7 0-120.1 48.7-129 112zM0 482.3C0 383.8 79.8 304 178.3 304h91.4C368.2 304 448 383.8 448 482.3c0 16.4-13.3 29.7-29.7 29.7H29.7C13.3 512 0 498.7 0 482.3z"/></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 24c0-13.3-10.7-24-24-24s-24 10.7-24 24V64H64C28.7 64 0 92.7 0 128v16 48V448c0 35.3 28.7 64 64 64H384c35.3 0 64-28.7 64-64V192 144 128c0-35.3-28.7-64-64-64H344V24c0-13.3-10.7-24-24-24s-24 10.7-24 24V64H152V24zM48 192H400V448c0 8.8-7.2 16-16 16H64c-8.8 0-16-7.2-16-16V192z"/></svg> November 08, 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")
```

---

# Quiz 4

### Download the template from the #class channel

### Make sure you unzip it!

### When done, submit your `quiz4.qmd` on Blackboard

]

]

---

# 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 package: {sf}

Art by<br>[Allison Horst](https://allisonhorst.com/)

]

---

### 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: -70.06611 ymin: -18.01973 xmax: 74.89131 ymax: 60.40581
#> Geodetic CRS:  +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
#>          name                       geometry
#> 0       Aruba MULTIPOLYGON (((-69.89912 1...
#> 1 Afghanistan MULTIPOLYGON (((74.89131 37...
#> 2      Angola MULTIPOLYGON (((14.19082 -5...
#> 3    Anguilla MULTIPOLYGON (((-63.00122 1...
#> 4     Albania MULTIPOLYGON (((20.06396 42...
#> 5       Aland MULTIPOLYGON (((20.61133 60...
```

]]

---

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

]]

]

---

## Inset Hawaii and Alaska

Get the shape file from {tigris} package

```r
library(tigris)

us_sf <- tigris::states(class = "sf", cb = TRUE) %>%
    shift_geometry() %>%
    filter(GEOID < 60)
```

Make the plot with `geom_sf()`

```r
us_sf %>%
    ggplot() +
    geom_sf()
```

]]

]

---

## 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/2023-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/2023-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?]

---

]

# Not all maps are useful...

]

---

]

# ...but some maps are

]

---

# 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
#> Geodetic CRS:  WGS 84
#>           name milk_produced                       geometry
#> 1   Washington         6.526 MULTIPOLYGON (((-122.753 48...
#> 2        Idaho        14.627 MULTIPOLYGON (((-117.0382 4...
#> 3      Montana         0.288 MULTIPOLYGON (((-116.0482 4...
#> 4 North Dakota         0.345 MULTIPOLYGON (((-104.0476 4...
#> 5    Minnesota         9.864 MULTIPOLYGON (((-97.22609 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

---

# 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
#> Geodetic CRS:  WGS 84
#> # 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>

]

---

# Remember: Mini Project 3 Due 11/21

## (That's 2 weeks)

### https://eda.seas.gwu.edu/2023-Fall/mini/3-redesign.html

---

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