---
title: "City Compare"
output:
  html_document:
    toc: true
    toc_depth: 5
    toc_float:
      collapsed: false
      smooth_scroll: true
editor_options:
  chunk_output_type: console
---

# Input Data & Configuration

## Libraries

```{r libs, eval = TRUE, echo = TRUE, results = "show", warning = FALSE, error = TRUE, message = FALSE}
date()
rm(list=ls())
library(tidyverse)
library(tidycensus)
library(sf)
library(openmeteo)
library(maps)
library(scales)

```

## API keys
```{r api_keys, eval = TRUE, echo = TRUE, results = "show", warning = FALSE, error = TRUE, message = FALSE}
# load census api key
census_api_key(key = substr(read_file(file = "api_keys/census_api_key"), 1, 40))

```

## Date ranges
```{r date_range, eval = TRUE, echo = TRUE, results = "show", warning = FALSE, error = TRUE, message = FALSE}
date_start <- "2010-01-01"
date_end <- "2024-12-31"

```
## Cities to compare

```{r cities, eval = TRUE, echo = TRUE, results = "show", warning = FALSE, error = TRUE, message = FALSE}
cities <- read_csv(file = "cities.csv")
cities <- cities %>%
  mutate(city_name = paste0(City, " ", Type))

```

# Get data

## Census data
```{r census, eval = TRUE, echo = TRUE, results = "show", warning = FALSE, error = TRUE, message = FALSE}
populations <- list(NULL)
for(city in cities$city_name){
  state <- cities %>% filter(city_name == city) %>% pull(State)
  populations[[city]] <- get_acs(
    geography = "place",
    variables = "B01003_001",
    state = state,
    year = 2023,
    geometry = TRUE
  ) %>%
  filter(str_detect(NAME, city))
}

populations <- bind_rows(populations)

city_center <- populations %>%
  st_centroid() %>%
  st_transform(4326) %>%  # Convert to WGS84 (standard lat/lon)
  mutate(
    lon = st_coordinates(.)[,1],
    lat = st_coordinates(.)[,2]
  ) %>%
  st_drop_geometry() %>%
  select(lat, lon)

cities <- bind_cols(cities, populations, city_center)
cities <- cities %>% mutate(
  density = estimate/as.double(st_area(geometry))*1000000
)

ggplot(cities) +
  geom_col(aes(x = City,
               y = estimate)) +
  labs(title = "City Population",
       x = "City",
       y = NULL) +
  scale_y_continuous(label = unit_format(unit = "K", scale = 1e-3, sep = ""))

ggplot(cities) +
  geom_col(aes(x = City,
               y = density)) +
  labs(title = "City Density",
       x = "City",
       y = "people/km^2") +
  scale_y_continuous()
```

## Map cities
```{r cities_map, eval = TRUE, echo = TRUE, results = "show", warning = FALSE, error = TRUE, message = FALSE}
ggplot(data = cities) +
  geom_polygon(data = map_data(map = "state"),
               aes(long, lat, group = group),
               fill = "white", colour = "grey50") +
  geom_point(aes(x = lon,
                 y = lat),
             shape = 21,
             fill = "lightgreen",
             color = "black",
             size = 4)


```

## weather
```{r weather, eval = TRUE, echo = TRUE, results = "show", warning = FALSE, error = TRUE, message = FALSE}
weather <- list(NULL)
for(city in cities$City){
  city_info <- cities %>% filter(City == city)
  city_run <- weather_history(
    location = c(city_info %>% pull(lat), city_info %>% pull(lon)),
    start = date_start,
    end = date_end,
    daily = c("apparent_temperature_max", "apparent_temperature_min"),
    response_units = list(temperature_unit = "fahrenheit")
  )
  city_run$city <- city
  weather[[city]] <- city_run
}

weather <- bind_rows(weather)

weather_summary <- weather %>%
  mutate(year = year(ymd(date)),
         month = month(ymd(date))) %>%
  group_by(year, city) %>%
  summarise(days_above_80 = sum(daily_apparent_temperature_max > 80)) %>%
  group_by(city) %>%
  summarise(median_days_above_80 = median(days_above_80))

ggplot(data = weather_summary) +
  geom_col(aes(x = city,
               y = median_days_above_80)) +
  labs(title = "Days above 80°F",
       y = "Median days per year",
       x = "City",
       fill = NULL)

ggplot(data = weather %>% pivot_longer(cols = starts_with("daily"), names_to = "max_min", values_to = "temperature") %>% filter(max_min %in% c("daily_apparent_temperature_min", "daily_apparent_temperature_max"))) +
  geom_violin(aes(x = city,
                  y = temperature,
                  fill = max_min)) +
  scale_fill_manual(labels = c("daily max", "daily min"), 
                    values = c("firebrick", "dodgerblue")) +
  labs(title = "Apparent Temperature",
       y = "°F",
       x = "City",
       fill = NULL)
```