---
title: "NY HEAT Analysis"
output:
  pdf_document:
    fig_caption: true
    fig_crop: true
    fig_retina: true
    keep_tex: true
    number_sections: true
    toc: false
    pandoc_args: ["--quiet"]
execute:
  echo: false
  warning: false
---

```{r set_variables}
state_code = "NY"
state_fips = "36"
heat_cutoff = 0.06
```

## NY HEAT and the Hudson Valley

```{r}
library(tidyverse)
library(sf)
library(gt)
library(tigris)
library(sf)
library(tidycensus)
options(tigris_use_cache=TRUE, tigris_year=2019)
```

```{r define-heat_cutoffv}
hudson_valley <- tribble(
  ~county, ~region, ~fips,
  "Westchester", "Lower Hudson Valley", 36119,
  "Rockland", "Lower Hudson Valley", 36087,
  "Putnam", "Lower Hudson Valley", 36115,
  "Orange", "Mid-Hudson Valley", 36121,
  "Dutchess", "Mid-Hudson Valley", 36027,
  "Ulster", "Mid-Hudson Valley", 36111,
  "Sullivan", "Upper Hudson Valley", 36053,
  "Columbia", "Upper Hudson Valley", 36021,
  "Greene", "Upper Hudson Valley", 36017,
  "Albany", "Upper Hudson Valley", 36001,
  "Rensselaer", "Upper Hudson Valley", 36083
)
```

```{r load_mapping}
puma_county_mapping <- read_csv("/mnt/data/buildings2/geocorr/puma_county_mapping.csv", skip=1) |>
  janitor::clean_names()
```

```{r load_pums}
vars <- c("BLD", "PUMA", "HFL", "TEN", "YRBLT", "HINCP", "NP", "MV", "BDSP", "ELEP", "FULP", "GASP", "MRGP", "RNTP", "RMSP", "ELEFP", "FULFP", "GASFP")

# Get data from the Census API
filename <- paste0("/mnt/data/pums_1_22_", state_fips, ".Rds")

# Get data from the Census API
if (file.exists(filename)) {
  raw_data <- readRDS(filename)
} else {
  raw_data <- get_pums(variables = vars,
                 state = state_fips,
                 year = 2022,
                 survey = "acs1",
                 rep_weights = "housing",
                 recode = TRUE) |>
          distinct(SERIALNO, .keep_all = TRUE) # housing units only, not people
  saveRDS(raw_data, filename)
}
```

```{r lmi-data}
county_lmi <- readxl::read_excel("/mnt/data/buildings2/HUD_2021_income_limits_by_NY_county.xlsx") |>
  pivot_longer(-County, names_to = "NP") |>
  rename(lmi_cutoff = value)
```

```{r clean_pums}
data <- raw_data |>
  filter(HINCP != -60000) |> # remove missing for now
  mutate(HINCP = HINCP+1,
   bill_in_rent = (ELEFP == "1") | (GASFP == "1") | (FULP == "1"),
    GASP = case_when(
    GASP <= 4 ~ 0,
    .default = GASP
  ),
   ELEP = case_when(
    ELEP <= 3 ~ 0,
    .default = ELEP
  ),
   FULP = case_when(
    FULP <= 3 ~ 0,
    .default = FULP
   ),
  yearly_energy_bill = (ELEP + GASP) * 12 + FULP,
  energy_burden_pct = yearly_energy_bill / (HINCP+1),
  monthly_energy_burden = (yearly_energy_bill - (heat_cutoff*HINCP))/12,
  is_energy_burdened = energy_burden_pct > heat_cutoff,
  # `utility` bill does not include delivered fuel
  yearly_utility_bill = (ELEP + GASP) * 12,
  utility_burden_pct = yearly_utility_bill / (HINCP+1),
  monthly_utility_burden = (yearly_utility_bill - (heat_cutoff*HINCP))/12,
  is_utility_burdened = utility_burden_pct > heat_cutoff) |>
  inner_join(puma_county_mapping, by = c("PUMA" = "puma_2012"),
            relationship = 'many-to-many') |>
  inner_join(hudson_valley, by = c('county_code' = 'fips')) |>
  mutate(wt = WGTP * puma12_to_county_allocation_factor) |>
  mutate(NP = as.character(NP)) |>
  inner_join(county_lmi, by = c("county" = "County", "NP")) |>
  mutate(is_lmi = HINCP <= lmi_cutoff)
```

```{r}
total_rollup <- data |>
  filter(!bill_in_rent) |>
  summarise(
    households_included = sum(wt),
    median_income = DescTools::Median(HINCP, weights = wt),
    pct_energy_burdened = weighted.mean(is_energy_burdened, w = wt),
    pct_utility_burdened = weighted.mean(is_utility_burdened, w =wt),
    avg_monthly_bill_of_burdened = weighted.mean(case_when(
      is_energy_burdened ~ yearly_energy_bill/12), w = wt,
      na.rm = TRUE),
    utility_burden_of_burdened = mean(case_when(
      is_utility_burdened ~ monthly_utility_burden), w = wt,
      na.rm = TRUE)) |>
mutate(county = "All",
      region = "Hudson Valley")
```

```{r county-calculation}
county_rollup <- data |>
  filter(!bill_in_rent) |>
  summarise(
    households_included = sum(wt),
    median_income = DescTools::Median(HINCP, weights = wt),
    pct_energy_burdened = weighted.mean(is_energy_burdened, w = wt),
    pct_utility_burdened = weighted.mean(is_utility_burdened, w =wt),
    avg_monthly_bill_of_burdened = weighted.mean(case_when(
      is_energy_burdened ~ yearly_energy_bill/12), w = wt,
      na.rm = TRUE),
    utility_burden_of_burdened = mean(case_when(
      is_utility_burdened ~ monthly_utility_burden), w = wt,
      na.rm = TRUE),
    .by = county) |>
  left_join(hudson_valley, by = "county") |>
  select(-fips)
```

```{r}
#| label: pct-hudson-valley
total_rollup$pct_energy_burdened[[1]] |> scales::percent(1)
```

```{r}
#| label: tbl-hudson-valley-table
bind_rows(total_rollup, county_rollup) |>
  select(region, county, pct_energy_burdened, avg_monthly_bill_of_burdened, utility_burden_of_burdened) |>
  group_by(region) |>
  gt() |>
  row_group_order(c("Hudson Valley", "Upper Hudson Valley", "Mid-Hudson Valley", "Lower Hudson Valley")) |>
  fmt_currency(c(avg_monthly_bill_of_burdened,
    utility_burden_of_burdened), decimals = 0) |>
  fmt_percent(c(pct_energy_burdened), decimals = 0) |>
  cols_label(
    county = "County",
    pct_energy_burdened = "Energy Burdened Households (%)",
    avg_monthly_bill_of_burdened = "Avg. Monthly Energy Bills of Burdened Households",
    utility_burden_of_burdened = "Avg. Monthly Savings Under NY HEAT for Burdened Households",
  ) |>
  tab_header("NY HEAT in the Hudson Valley", subtitle = "Household energy burdens and savings by county") |>
  tab_footnote("Excluding households with utility bills in rent, or missing fuel or income information.\nSource: 2022 ACS. Win Climate, 2024")
```

```{r}
#| label: total-hudson-valley-savings
county_rollup |>
  summarise(sum(utility_burden_of_burdened * households_included))
```

```{r}
#| label: hudson-valley-burden-income
data |>
  filter(!is.na(is_lmi)) |>
  group_by(is_lmi) |>
  summarise(pct_energy_burdened = weighted.mean(is_energy_burdened, w = wt))
```

```{r}
#| label: fig-hudson-valley-map-pct
county_boundaries <- counties("NY")

county_boundaries |>
  right_join(county_rollup, by = c("NAME" = "county")) |>
  mutate(label = paste0(NAME," (", scales::percent(pct_energy_burdened, 1) ,")")) |>
  ggplot(aes(fill=pct_energy_burdened, label=label)) +
  geom_sf(color = "black") +
  ggrepel::geom_label_repel(
    aes(geometry = geometry),
    fill = "white",
    stat = "sf_coordinates",
    min.segment.length = 0
    ) +
  scale_fill_distiller(palette = "Spectral",
                       name = "Percentage",
                       limits = c(0, .5),
                       labels = scales::percent_format(scale=100)) +
  theme_void() +
  labs(title = "Energy Burden by County",
       subtitle = "Win Climate, 2024\n(Source: 2022 ACS)",
       pct = "Percent")
```

```{r}
#| label: fig-hudson-valley-map-amt
county_boundaries |>
  right_join(county_rollup, by = c("NAME" = "county")) |>
  mutate(label = paste0(NAME," (", scales::dollar(utility_burden_of_burdened, 1) ,")")) |>
  ggplot(aes(fill=utility_burden_of_burdened, label=label)) +
  geom_sf(color = "black") +
  ggrepel::geom_label_repel(
    aes(geometry = geometry),
    fill = "white",
    stat = "sf_coordinates",
    min.segment.length = 0
    ) +
  scale_fill_distiller(palette = "Greens",
                       direction = 1,
                       name = "NY HEAT Savings",
                       labels = scales::dollar_format(accuracy = 1),
                       limits = c(0, 250)) +
  theme_void() +
  labs(title = "Average Monthly NY HEAT Savings by County",
       subtitle = "Win Climate, 2024\n(Source: 2022 ACS)")
```

```{r}
#| label: hudson-valley-map-puma
#| fig-height: 8
puma_boundaries <- pumas("NY", 2019)

puma_rollup <- data |>
  filter(!bill_in_rent) |>
  summarise(
    households_included = sum(wt),
    median_income = DescTools::Median(HINCP, weights = wt),
    pct_energy_burdened = weighted.mean(is_energy_burdened, w = wt),
    avg_monthly_bill_of_burdened = weighted.mean(case_when(
      is_energy_burdened ~ yearly_energy_bill/12), w = wt,
      na.rm = TRUE),
    utility_burden_of_burdened = mean(case_when(
      is_energy_burdened ~ monthly_energy_burden), w = wt,
      na.rm = TRUE),
    .by = c(puma12_name, PUMA))  |>
  mutate(GEOID10 = paste0(state_fips, PUMA))

puma_boundaries |>
  right_join(puma_rollup, by = c("GEOID10")) |>
  #mutate(label = paste0(NAME," (", scales::percent(pct_energy_burdened, 1) ,")")) |>
  ggplot(aes(fill=pct_energy_burdened)) +#, label=label)) +
  geom_sf(color = "black", alpha=.5) +
  #ggrepel::geom_label_repel(
  #  aes(geometry = geometry),
  #  fill = "white",
  #  stat = "sf_coordinates",
  #  min.segment.length = 0
  #  ) +
  scale_fill_distiller(palette = "Spectral",
                       name = "Percentage",
                       limits = c(0, .5),
                       labels = scales::percent_format(scale=100)) +
  theme_void() +
  theme(legend.position = "top") +
  labs(title = "Energy Burden by County",
       subtitle = "Win Climate for Spring Street Climate Fund, 2024\n(Source: 2022 ACS)",
       pct = "Percent")
```

```{r}
#| label: hudson-valley-post-burden

```
