
---
title: "In-class Exercises: Mapping"
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(janitor) # cleaning
library(leaflet)
library(macleish)

# Libraries for Mapping
library(sf) # shape files
library(maps) # we get the maps of VT and MA from here. 
library(tidygeocoder) # geocoding. 
library(tidycensus)

```

![](https://xkcd.com/1138)

# 3. Choropleth maps using census data

```{r}
# Get MA county sf object, that includes population estimate column in data
mass_pop_orig <- 
  tidycensus::get_acs(
    geography = "county", 
    variables = "B01003_001", 
    state = "MA",
    geometry = TRUE
  ) |> 
  # Add centroids to each region using purrr package
  mutate(
    lng = purrr::map_dbl(geometry, ~st_centroid(.x)[[1]]),
    lat = purrr::map_dbl(geometry, ~st_centroid(.x)[[2]])
  )  
mass_pop_orig

# Unfortunately the original NAME variable above is very long and redundant.
# Create a new data frame with the just the county name using the separate()
# function from the tidyr package. 
mass_pop <- mass_pop_orig |> 
  tidyr::separate(NAME, c("County", "State"), sep = ", ") |> 
  tidyr::separate(County, c("County", "Fluff"), sep = " ") 

# Plot gradient choropleth map based on numerical variable estimate
ggplot() +
  geom_sf(data = mass_pop, aes(fill = estimate)) +
  geom_text(data = mass_pop, aes(x = lng, y = lat, label = County)) +
  labs(fill = "Population\nEstimate")+
  scale_fill_continuous(type= "viridis", direction= -1) 
```

Note: You just made a population map based on a continuous scale. We'll use the graph you made above for reference. 

Below we bin our population variable, this will make three equally sized groups. 

```{r}
# What if instead we want a bin based choropleth map. Use the cut_number()
# function to convert the numerical variable estimate 
# to a categorical variable with 3 levels!
mass_pop <- mass_pop |>
  mutate(estimate_discrete = cut_number(estimate, 3))
mass_pop

# Plot!
ggplot() +
  geom_sf(data = mass_pop, aes(fill = estimate_discrete)) +
  geom_text(data = mass_pop, aes(x = lng, y = lat, label = County)) +
  labs(fill = "Population\nEstimate")+
  scale_fill_brewer(palette = "Set2")
```

Check out this site for more on [tidycensus](https://walker-data.com/tidycensus/articles/basic-usage.html)

Specifically the `load_variables` function.
```{r}
v20 <- load_variables(2020, "acs5", cache = TRUE)

View(v20)
```


# What counties have the highes proportion of under 5s?

We need to weight our data in some way. 
```{r}
# Where do all the under 5 year olds live?

mass_pop_under_5 <- 
  tidycensus::get_acs(
    geography = "county", 
    variables = c(population = "B01003_001", under_5m="B01001_003", under_5f= "B01001_027"),
    state = "MA",
    # We'll need to pivot wider, but that doesn't work with simple features. 
    geometry = FALSE
  )|>
  # Add centroids to each region using purrr package
 tidyr::separate(NAME, c("County", "State"), sep = ", ") |> 
 tidyr::separate(County, c("County", "Fluff"), sep = " ") 

# Just count them we get a similar map to before. 

props <- mass_pop_under_5 |> 
  select(County,variable,estimate)|>
  pivot_wider(names_from = variable, values_from = estimate)|>
  reframe(County, number_under_5 = under_5m+under_5f, prop_under_5 = (under_5m+under_5f)/population)

# Do a left join, a right join will lose the multipolygon-ness
mass_prop <- mass_pop |>
  left_join(props, by = "County")

ggplot() +
  geom_sf(data = mass_prop, aes(fill = prop_under_5)) +
  scale_fill_distiller(direction = 1)
```

What is wrong with this map? 

```{r}

# Here is the map
ggplot() +
  geom_sf(data = mass_prop, aes(fill = number_under_5)) +
  geom_text(data = mass_pop, aes(x = lng, y = lat, label = County)) +
  labs(fill = "Number of children under 5 years old")+
  scale_fill_continuous(type= "viridis", direction= -1) 
  
```

## Exercise

Use the `tidycensus` package to conduct a spatial analysis of the American Community Survey data it contains for Massachusetts. With a map answer the question: Which counties in Massachusetts have the highest number of retirees? What about the highest proportion of retirees? Which of these will look just like the population map? 

Hint: you want to load the variables, the look at the data frame and search "Not in labor force" Pick the ages for those who are likely to be retirees. When I did it for 2020 data I found six variables, pay close attention to the actual variable name. To find the proportion, include a seventh variable, population, which we grabbed above.  

[Instructions for using the tidycensus package](https://walker-data.com/tidycensus/articles/basic-usage.html)

```{r}
# Write your code below


```


# 4 Interactive maps using leaflet. If time. 

## Color in leaflet

The MacLeish package. https://github.com/beanumber/macleish

```{r}
library(macleish)
names(macleish_layers)
```

```{r}
leaflet() |>
  addTiles() |>
  addPolygons(
    data = macleish_layers[["boundary"]], 
    weight = 1, fillOpacity = 0.1
  )%>%
  addPolygons(
    data = macleish_layers[["buildings"]], 
    weight = 1, popup = ~ name
  ) %>%
  addPolylines(
    data = macleish_layers[["trails"]], 
    weight = 1, color = "brown",
    popup = ~ name
  ) %>%
  addPolylines(
    data = macleish_layers[["streams"]], 
    weight = 2
  ) %>%
  addMarkers(
    data = filter(macleish_layers[["landmarks"]], grepl("Met", Label)), 
    popup = ~Label
  )

```

```{r}
wetlands <- macleish_layers |>
  pluck("wetlands")

leaflet() |>
  addTiles()|>
  addPolygons(data = wetlands, 
              color = ~SHAPE_LEN, 
              popup = ~IT_VALDESC)

wetland_pal <- colorNumeric(
  palette = "viridis", 
  domain = wetlands |>
    pull(SHAPE_LEN)
)

leaflet() |>
  addTiles() |>
  addPolygons(data = wetlands, 
              color = ~wetland_pal(SHAPE_LEN), 
              fillOpacity = 0.8, weight = 0.1, 
              popup = ~IT_VALDESC) |>
  addLegend(position = "bottomleft", 
            pal = wetland_pal, 
            values = wetlands$SHAPE_LEN,
            title = "Wetland Area")
```



