---
title: "Lec 29: In-class Exercise: Static Maps"
output:
  html_document
---

```{r, message=FALSE, warning=FALSE, echo=FALSE}
# Set R chunk default options:
# -do not output all R messages in html report
# -do not output all R warnings in html report
# -set default figure widths and heights for all code blocks, unless you 
#  override them
knitr::opts_chunk$set(
  message = FALSE, warning = FALSE, fig.width = 16/2, fig.height = 9/2
)
```

```{r, echo=FALSE}
# Load all packages here. It's good practice to load all packages at the
# beginning of a qmd file so that collaborators can see in one shot all packages
# they'll need installed:
library(tidyverse)
library(sf)
library(leaflet)
library(maps)
library(tmaptools)

# USA specific data
library(USAboundaries)
library(tidycensus)
library(macleish)
```

# Converting data frames to sf objects

Draw a map of the US with the 2 universities I've worked at marked:

```{r}
# How did I get the longitude and latitude for Worcester Polytechnic Institute and Smith College? By tapping into
# the OpenStreetMap (OSM) database using the geocode_OSM() function from the
# tmaptools package. Be careful though, apparently you are limited to 2000
# queries a day.
ne_colleges <- tibble(
  name = c("Worcester Polytechnic Institute", "Smith College"), 
  lng = c(-71.808, -72.640),
  lat = c(42.274, 42.316)
)
```

```{r}
# Convert ne_colleges data frame to sf object using st_as_sf() from sf
# package. Note we set the coordinate reference system to WGS84
# https://epsg.io/4326
ne_colleges_sf <- ne_colleges |>
  st_as_sf(coords = c("lng", "lat"), crs = 4326)
```

```{r}
# Load maps package map object of US states and plot using default R plot
# function (not ggplot2). Plot is not pretty!
USA_map <-  maps::map("state", plot = FALSE, fill = TRUE)
plot(USA_map)

# Convert maps package map object to sf object using st_as_sf() again:
USA_sf <- USA_map |> 
  st_as_sf(coords = c("x", "y")) |>
  st_transform(crs = 4326)
```

```{r}
# Plot both layers
ggplot() +
  geom_sf(data = USA_sf) +
  geom_sf(data = ne_colleges_sf, size = 2)
```

```{r}
# Obtain a map of Massachusetts counties from the USAboundaries package
MA_counties <- 
  USAboundaries::us_counties(states = c("massachusetts")) |> 
  select(-state_name)
```

```{r}
# Plot both layers
ggplot() +
  geom_sf(data = MA_counties) +
  geom_sf(data = ne_colleges_sf, size = 2)
```

## Exercise 1

Draw a map of only Massachusetts with some notable cities marked. This will require you to manually create a data frame like `ne_colleges` above. 

```{r}
# Write your code below

```

# 2. Loading shapefiles into R

Let's load a shapefile of all bike trails in Massachusetts as provided by
[MassGIS (Bureau of Geographic Information)](https://www.mass.gov/info-details/massgis-data-bicycle-trails)

```{r}
# Load the "shapefile" contained in the biketrails_arc/ folder and load
# it into R using the sf package's read_sf() function.
bike_trails <- read_sf("biketrails_arc/") |>
  st_transform(crs = 4326)

# Plot this data!
ggplot() +
  geom_sf(data = bike_trails, color = "red")

# The above map lacks context however. Where are these bikepaths?
# Let's obtain a map of Massachusetts counties from the USAboundaries package
MA_counties <- 
  USAboundaries::us_counties(states = c("massachusetts")) |> 
  select(-state_name)
  
# Plot this data!
ggplot() +
  geom_sf(data = MA_counties)

# Plot both maps by adding both layers!
ggplot() +
  geom_sf(data = bike_trails, color = "red") + 
  geom_sf(data = MA_counties) 

# Where did our bike paths go? The order you plot in matters! Switch the two
# geom_sf()
ggplot() +
  geom_sf(data = MA_counties) +
  geom_sf(data = bike_trails, color = "red")

# Better! But what are these counties? Let's mark each counties' centroid
MA_counties <- MA_counties |>
  # Add centroids to each region using purrr package (covered in SDS 270 I believe)
  mutate(
    lng = purrr::map_dbl(geometry, ~st_centroid(.x)[[1]]),
    lat = purrr::map_dbl(geometry, ~st_centroid(.x)[[2]])
    )
MA_counties

# Plot! Notice how we define data = in the geom_point() instead of the ggplot().
# This sets the data set only for the geom_point() and not the other geoms.
# Any data = defined in the global ggplot() trickles down by default to the other
# geoms.
ggplot() +
  geom_sf(data = MA_counties) +
  geom_sf(data = bike_trails, color = "red") +
  geom_point(data = MA_counties, aes(x = lng, y = lat), size = 1) 

# Can we do better? Yes! Use geom_text() instead with the label aesthetic!
ggplot() +
  geom_sf(data = MA_counties) +
  geom_sf(data = bike_trails, color = "red") +
  geom_text(data = MA_counties, aes(x = lng, y = lat, label = name)) 
```

## Exercise 2

Create a map of Washington state showing all 39 counties using the `WA_counties` shapefile data with Seattle, Olympia, and Spokane marked with a dot. Use a **relative** filepath and not an **absolute** filepath to load the `WA_counties` shapefile.

> Hint: You need to set `crs` when loading the `WA_counties` shapefile data

```{r}
# Write your code below

```
