```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Overview
“In February 2021, the state of Texas suffered a major power crisis, which came about as a result of three severe winter storms sweeping across the United States on February 10–11, 13–17, and 15–20.”[^wiki] For more background, check out these [engineering](https://www.youtube.com/watch?v=08mwXICY4JM&ab_channel=PracticalEngineering) and [political](https://www.youtube.com/watch?v=Zcrsgdl_hP0&ab_channel=Vox) perspectives. 

For this assignment, you are tasked with:  
-   estimating the number of homes in Houston that lost power as a result of the first two storms  
-   investigating if socioeconomic factors are predictors of communities recovery from a power outage  

Your analysis will be based on remotely-sensed night lights data, acquired from the [Visible Infrared Imaging Radiometer Suite (VIIRS)](https://en.wikipedia.org/wiki/Visible_Infrared_Imaging_Radiometer_Suite) onboard the Suomi satellite. In particular, you will use the VNP46A1 to detect differences in night lights before and after the storm to identify areas that lost electric power.

To determine the number of homes that lost power, you link (spatially join) these areas with [OpenStreetMap](https://www.openstreetmap.org/#map=4/38.01/-95.84) data on buildings and roads. 

To investigate potential socioeconomic factors that influenced recovery, you will link your analysis with data from the US Census Bureau. 

##### Learning objectives:
-   load vector/raster data  
-   simple raster operations  
-   simple vector operations  
-   spatial joins  

### Data
#### Night lights
Use NASA's Worldview to explore the data around the day of the storm. There are several days with too much cloud cover to be useful, but 2021-02-07 and 2021-02-16 provide two clear, contrasting images to visualize the extent of the power outage in Texas.

VIIRS data is distributed through NASA's [Level-1 and Atmospheric Archive & Distribution System Distributed Active Archive Center (LAADS DAAC)](https://ladsweb.modaps.eosdis.nasa.gov/). Many NASA Earth data products are distributed in 10x10 degree tiles in sinusoidal equal-area projection. Tiles are identified by their horizontal and vertical position in the grid. Houston lies on the border of tiles h08v05 and h08v06. We therefore need to download two tiles per date.

As you're learning in EDS 220, accessing, downloading, and preparing remote sensing data is a skill in it's own right! To prevent this assignment from being a large data wrangling challenge, we have downloaded and prepped the following files for you to work with, stored in the `VNP46A1` folder.\

-   `VNP46A1.A2021038.h08v05.001.2021039064328.h5.tif`: tile h08v05, collected on 2021-02-07  
-   `VNP46A1.A2021038.h08v06.001.2021039064329.h5.tif`: tile h08v06, collected on 2021-02-07  
-   `VNP46A1.A2021047.h08v05.001.2021048091106.h5.tif`: tile h08v05, collected on 2021-02-16  
-   `VNP46A1.A2021047.h08v06.001.2021048091105.h5.tif`: tile h08v06, collected on 2021-02-16  

#### Roads
Typically highways account for a large portion of the night lights observable from space (see Google's [Earth at Night](https://earth.google.com/web/@27.44405464,-84.7693044,206.63660162a,8916361.52264659d,35y,0h,0t,0r/data=CiQSIhIgMGY3ZTJkYzdlOGExMTFlNjk5MGQ2ZjgxOGQ2OWE2ZTc)). To minimize falsely identifying areas with reduced traffic as areas without power, we will ignore areas near highways.

[OpenStreetMap (OSM)](https://planet.openstreetmap.org/) is a collaborative project which creates publicly available geographic data of the world. Ingesting this data into a database where it can be subsetted and processed is a large undertaking. Fortunately, third party companies redistribute OSM data. We used [Geofabrik's download sites](https://download.geofabrik.de/) to retrieve a shapefile of all highways in Texas and prepared a Geopackage (`.gpkg` file) containing just the subset of roads that intersect the Houston metropolitan area.\  

- `gis_osm_roads_free_1.gpkg`

[^wiki]:Wikipedia. 2021. “2021 Texas power crisis.” Last modified October 2, 2021. https://en.wikipedia.org/wiki/2021_Texas_power_crisis.

#### Houses
We can also obtain building data from OpenStreetMap. We again downloaded from Geofabrick and prepared a GeoPackage containing only houses in the Houston metropolitan area.\

- `gis_osm_buildings_a_free_1.gpkg`

#### Socioeconomic
We cannot readily get socioeconomic information for every home, so instead we obtained data from the [U.S. Census Bureau's American Community Survey](https://www.census.gov/programs-surveys/acs) for census tracts in 2019. The *folder* `ACS_2019_5YR_TRACT_48.gdb` is an ArcGIS ["file geodatabase"](https://desktop.arcgis.com/en/arcmap/latest/manage-data/administer-file-gdbs/file-geodatabases.htm), a multi-file proprietary format that's roughly analogous to a GeoPackage file.\

You can use `st_layers()` to explore the contents of the geodatabase. Each layer contains a subset of the fields documents in the [ACS metadata](https://www2.census.gov/geo/docs/maps-data/data/tiger/prejoined/ACSMetadata2011.txt).\

The geodatabase contains a layer holding the geometry information, separate from the layers holding the ACS attributes. You have to combine the geometry with the attributes to get a feature layer that `sf` can use. 

## Assignment

Below is an outline of the steps you should consider taking to achieve the assignment tasks. 

#### Find locations of blackouts
For improved computational efficiency and easier interoperability with `sf`, I recommend using the `stars` package for raster handling.\

##### combine the data (5 points)
-   read in night lights tiles\
-    combine tiles into a single `stars` object for each date (2021-02-07 and 2021-02-16)\
    -   hint: use `st_mosaic`\
    
```{r include=TRUE,eval=FALSE}
library(sf)
library(tidyverse)
library(terra)
library(here)
library(tmap)
library(stars)
rm(list = ls())
setwd(here())

day_night_band_day38_tile5 <- read_stars("data/VNP46A1/VNP46A1.A2021038.h08v05.001.2021039064328.tif")
day_night_band_day38_tile6 <- read_stars("data/VNP46A1/VNP46A1.A2021038.h08v06.001.2021039064329.tif")
day_night_band_day47_tile5 <- read_stars("data/VNP46A1/VNP46A1.A2021047.h08v05.001.2021048091106.tif")
day_night_band_day47_tile6 <- read_stars("data/VNP46A1/VNP46A1.A2021047.h08v06.001.2021048091105.tif")


feb_7_comb <- st_mosaic(day_night_band_day38_tile5, day_night_band_day38_tile6)
feb_16_comb <- st_mosaic(day_night_band_day47_tile5, day_night_band_day47_tile6)
```
    
##### create a blackout mask (10 points)
-   find the change in night lights intensity (presumably) caused by the storm\
-   reclassify the difference raster, assuming that any location that experienced a drop of more than 200 nW cm^-2^sr^-1^ experienced a blackout\
-   assign `NA` to all locations that experienced a drop of *less* than 200 nW cm^-2^sr^-1^\

```{r include=TRUE,eval=FALSE}
difference <- (feb_7_comb - feb_16_comb) > 200
difference[difference == FALSE] <- NA 
```


##### vectorize the mask (5 points)
-   use `st_as_sf()` to vectorize the blackout mask\
-   fix any invalid geometries using `st_make_valid`

```{r include=TRUE,eval=FALSE}
vec_diff <- st_as_sf(difference)
vec_diff <- st_make_valid(vec_diff)
```


##### crop the vectorized map to our region of interest (10 points)
-   define the Houston metropolitan area with the following coordinates\
    -   (-96.5, 29), (-96.5, 30.5), (-94.5, 30.5), (-94.5, 29)\
-   turn these coordinates into a polygon using `st_polygon`\
-   convert the polygon into a simple feature collection using `st_sfc()` and assign a CRS\
    -   hint: because we are using this polygon to crop the night lights data it needs the same CRS\
-   crop (spatially subset) the blackout mask to our region of interest\ 
-   re-project the cropped blackout dataset to EPSG:3083 (NAD83 / Texas Centric Albers Equal Area)\

```{r include=TRUE,eval=FALSE}
pt1 <- st_point(c(-96.5, 29))
pt2 <- st_point(c(-96.5, 30.5))
pt3 <- st_point(c(-94.5, 30.5))
pt4 <- st_point(c(-94.5, 29))
coords <- list(rbind(pt1, pt2, pt3, pt4, pt1))
polygon <- st_polygon(x = coords)
houston <- st_sfc(polygon, crs = "EPSG:4326")

cropped_diff <- vec_diff[houston, op = st_intersects]
cropped_diff <- st_transform(cropped_diff, crs = "EPSG:3083")
```


##### exclude highways from blackout mask (10 points)
The roads geopackage includes data on roads other than highways. However, we can avoid reading in data we don't need by taking advantage of `st_read`'s ability to subset using a SQL query.\

-   define SQL query\
-   load just highway data from geopackage using `st_read`\
  -   reproject data to EPSG:3083\
-   identify areas within 200m of all highways using `st_buffer`\
    -   hint: `st_buffer` produces undissolved buffers, use `st_union` to dissolve them\
-   find areas that experienced blackouts that are further than 200m from a highway

`query <- "SELECT * FROM gis_osm_roads_free_1 WHERE fclass='motorway'"`\
`highways <- st_read("data/gis_osm_roads_free_1.gpkg", query = query)`

```{r include=TRUE,eval=FALSE}
query <- "SELECT * FROM gis_osm_roads_free_1 WHERE fclass='motorway'"

highways <- st_read("data/gis_osm_roads_free_1.gpkg", 
                    query = query, 
                    quiet = TRUE) %>%
  st_transform(crs = "EPSG:3083")

highway_buffer <- st_buffer(highways, dist = 200)
highway_union <- st_union(highway_buffer)

blackout_nohighway <- st_difference(cropped_diff, highway_union)
```


#### Find homes impacted by blackouts

##### load buildings data (10 points)
-   load buildings dataset using `st_read` and the following SQL query to select only residential buildings\
  -   hint: reproject data to EPSG:3083\

`SELECT *`\ 
`FROM gis_osm_buildings_a_free_1`\
`WHERE (type IS NULL AND name IS NULL)`\
`OR type in ('residential', 'apartments', 'house', 'static_caravan', 'detached')`\

```{r include=TRUE,eval=FALSE,eval=FALSE}
query <- "SELECT * FROM gis_osm_buildings_a_free_1 WHERE (type IS NULL AND name IS NULL) OR type in ('residential', 'apartments', 'house', 'static_caravan', 'detached')"

buildings <- st_read("data/gis_osm_buildings_a_free_1.gpkg", 
                     query = query, 
                     quiet = TRUE) %>%
  st_transform(crs = "EPSG:3083")
```

##### find homes in blackout areas (20 points)
-   filter to homes within blackout areas\
-   count number of impacted homes\

```{r include=TRUE,eval=FALSE}
blackout_nohighway <- st_difference(cropped_diff, highway_union)

affected_buildings <- buildings[blackout_nohighway, op = st_intersects]
nrow(affected_buildings)
```


#### Investigate socioeconomic factors

##### load ACS data (10 points)
-   use `st_read()` to load the geodatabase layers\
-   geometries are stored in the `ACS_2019_5YR_TRACT_48_TEXAS` layer\
-   income data is stored in the `X19_INCOME` layer\
-   select the median income field `B19013e1`\
  -   hint: reproject data to EPSG:3083\

```{r include=TRUE,eval=FALSE}
acs_geoms <- st_read("data/ACS_2019_5YR_TRACT_48_TEXAS.gdb",
                     layer = "ACS_2019_5YR_TRACT_48_TEXAS", 
                     quiet = TRUE)
acs_income <- st_read("data/ACS_2019_5YR_TRACT_48_TEXAS.gdb",
                      layer = "X19_INCOME", 
                      quiet = TRUE)

median_income <- acs_income[c("GEOID", "B19013e1")]

```


##### determine which census tracts experienced blackouts (10 points)
-   join the income data to the census tract geometries\
  -   hint: make sure to join by geometry ID\
-   spatially join census tract data with buildings determined to be impacted by blackouts\
-   find which census tracts had blackouts\

```{r include=TRUE,eval=FALSE}

acs_data <- dplyr::left_join(acs_geoms,
                             median_income,
                             by = c("GEOID_Data" = "GEOID")) %>%
  st_transform(crs = "EPSG:3083")


tracts_with_blackout_info <- st_join(acs_data,
                                     affected_buildings,
                                     left = TRUE) %>%
  group_by(NAME) %>%
  summarise(blackout = any(!is.na(osm_id)), median_income = first(B19013e1))
```


##### compare incomes of impacted tracts to unimpacted tracts (10 points)
-   create a map of median income by census tract, designating which tracts had blackouts\
-   plot the distribution of income in impacted and unimpacted tracts
-   write approx. 100 words summarizing your results and discussing any limitations to this study


