2021 Texas Blackout

Analyzing affected areas & socioeconomic factors that influenced recovery
Author

Flora Hamilton

Published

October 26, 2023

Texas Blackouts

How did socioeconomic factors influence community recovery from 2021 power outages in Houston?

📦 repo link | https://github.com/floraham/texas-blackout

Context

The 2021 Texas blackout was a result of an arctic blast that led to the freezing of a significant portion of the state’s power generation capacity. The blackout was not solely due to one type of energy source, as every electricity source—natural gas, coal, nuclear, wind, and solar—fell short during the extreme weather conditions. Many components of energy plants froze with the vast majority of these plants being powered by natural gas, leading to widespread outage. The blackout highlighted the vulnerability of the U.S. electric grid to extreme weather events and the need for prioritizing energy infrastructure to prevent similar crises in the future [1][2][3].

The blackout also disproportionately affected communities of color, with approximately 70% of Texans losing electrical power at some point between February 14th and 20th. The recovery from the loss of health and property damage became less likely for those most severely affected, making them more vulnerable to future crises. This highlighted the need for the prioritization of these communities through direct investment and focused attention from local and state governments in coordination with utility providers [4].

Project objectives and methods

This project aims to estimate the number of households in Houston that experienced power loss due to the initial two storms. To determine the count of homes affected by power outages, we will spatially join these areas with data on buildings and roads.

This project also investigates whether socioeconomic factors serve as predictors of community recovery from a power outage. For the exploration of potential socioeconomic factors influencing recovery, the analysis will be linked with data from the US Census Bureau.

This analysis will rely on remotely-sensed night lights data acquired from the Visible Infrared Imaging Radiometer Suite (VIIRS) aboard the Suomi satellite. Specifically, we will use the VNP46A1 to identify variations in night lights before and after the storms as a method to identify areas that lost electrical power.

Data

Night lights ✨

We 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). 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.

Roads 🛣️ gis_osm_roads_free_1.gpkg

Typically highways account for a large portion of the night lights observable from space (see Google’s Earth at Night). To minimize falsely identifying areas with reduced traffic as areas without power, we will ignore areas near highways.

OpenStreetMap (OSM) is a collaborative project which creates publicly available geographic data of the world. We used Geofabrik’s download sites 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. 

Houses 🏠 gis_osm_buildings_a_free_1.gpkg

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

Socioeconomic 👨‍👩‍👧‍👦 ACS_2019_5YR_TRACT_48.gdb

We cannot readily get socioeconomic information for every home, so instead we obtained data from the U.S. Census Bureau’s American Community Survey for census tracts in 2019. The folder ACS_2019_5YR_TRACT_48.gdb is an ArcGIS “file geodatabase”, a multi-file proprietary format that’s roughly analogous to a GeoPackage file. We can use st_layers() to explore the contents of the geodatabase. Each layer contains a subset of the fields documents in the ACS metadata. 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.

Project Roadmap

Below is an outline of the steps taken to achieve the objectives.

1) Find locations of blackouts

2) Find homes impacted by blackouts

3) Investigate Socioeconomic Factors

4) Leverage the results for a summary and discussion

Importing necessary libraries and functions

Code 
library(sf)
library(stars)
library(tmap)
library(raster)
library(terra)
library(dplyr)
library(ggplot2)

Find locations of blackouts

For improved computational efficiency and easier inter-operability with sf, we use the stars package for raster handling.

Combining the data

We first read in night lights tiles, then combine tiles into a single stars object for each date (2021-02-07 and 2021-02-16) using st_mosaic.

The following files had been subsetted & stored in the VNP46A1 folder.

  • VNP46A1.A2021038.h08v05.001.2021039064328.tif: tile h08v05, collected on 2021-02-07
  • VNP46A1.A2021038.h08v06.001.2021039064329.tif: tile h08v06, collected on 2021-02-07
  • VNP46A1.A2021047.h08v05.001.2021048091106.tif: tile h08v05, collected on 2021-02-16
  • VNP46A1.A2021047.h08v06.001.2021048091105.tif: tile h08v06, collected on 2021-02-16
Creating a blackout mask

We then find the change in night lights intensity (presumably) caused by the storm, and reclassify the difference raster, assuming that any location that experienced a drop of more than 200 nW cm-2sr-1 experienced a blackout. Then we assign NA to all locations that experienced a drop of less than 200 nW cm-2sr.-1

Code 
#find the change in night lights intensity (presumably) caused by the storm by calculating the light difference between the two dates 
difference <- combined_02_07 - combined_02_16


#check out difference in night lights intensity 
plot(difference)
downsample set to 6

Assigning all the houses that experienced a drop of 200 and more to a blackout mask.

Code 
mask_blackout <- difference > 200 
mask_blackout[mask_blackout == FALSE]  <- NA
Vectorize the mask

We use st_as_sf() to vectorize the blackout mask and fix any invalid geometries using st_make_valid.

Code 
blackout_mask_vector <- st_as_sf(mask_blackout) %>% st_make_valid() ##can plot() to check out what the blackout mask looks like
Cropping the vectorized map to our region of interest
  • We first define the Houston metropolitan area with the following coordinates : (-96.5, 29), (-96.5, 30.5), (-94.5, 30.5), (-94.5, 29) and then turn these coordinates into a polygon using st_polygon.
  • Then we convert the polygon into a simple feature collection using st_sfc() and assign a CRS, making sure that the polygon is in the same CRS and crop (spatially subset) the blackout mask to our region of interest. Lastly we re-project the cropped blackout dataset to EPSG:3083 (NAD83 / Texas Centric Albers Equal Area).
Code 
#crop the vectorized map to our region of interest using polygon of Houston's coordinates 
houston_polygon <- st_polygon(list(rbind(c(-96.5,29), c(-96.5,30.5), c(-94.5, 30.5), c(-94.5,29), c(-96.5,29))))


#convert the polygon into a simple feature collection using st_sfc() and assign the CRS 4326, which is the same as the night lights data
houston_border_sf <- st_sfc(houston_polygon, crs = st_crs(mask_blackout))


#inspect the polygon crs to make sure the crs is 4326, nightlights dataset
#st_crs(houston_border_sf) == st_crs(mask_blackout)

Spatially subsetting blackout mask to Houston region, then reprojecting to the correct CRS.

Excluding highways from blackout mask

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. First, we load just highway data from geopackage using st_read andreproject data to EPSG:3083, identifying areas within 200m of all highways using st_buffer. st_buffer produces undissolved buffers, so we use st_union to dissolve them. Then we find areas that experienced blackouts that are further than 200m from a highway.

Read in highway data from roads using SQL query :

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

Code 
query <- "SELECT * FROM gis_osm_roads_free_1 WHERE fclass='motorway'"
highways <- st_read("~/dev/eds223/assignment-3-floraham/data/gis_osm_roads_free_1.gpkg", query = query)

Then reproject and subset using st_difference()

Code 
# reproject data to EPSG:3083 
highways_3038 <- st_transform(highways, crs="EPSG:3083")

#verify it's the right crs
#st_crs(highways_3038)

# identify areas within 200m of all highways using st_buffer. Use st_union to dissolve overlapping polygons. 
a_near_hwys <- st_buffer(highways_3038, dist = 200) %>% st_union()

# find areas that experienced blackouts that are further than 200m from a highway by spacially subsetting houston blackout areas to ones near highways, and taking st_difference of that. 
blackouts_far_hwy_houston <- final_houston_blackout[a_near_hwys, , op = st_difference] 


#plot to see what it looks like 
#plot(blackouts_far_hwy_houston)

Finding homes impacted by blackouts

Load buildings data

We load buildings dataset using st_read and the following SQL query to select only residential buildings. We first reproject data to EPSG:3083.

The SQL query we will use is:

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')

Code 
#load buildings dataset using st_read and the following SQL query to select only residential buildings

query_buildings <- "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("~/dev/eds223/assignment-3-floraham/data/gis_osm_buildings_a_free_1.gpkg", query = query_buildings)
finding homes in blackout areas

We filter to homes within blackout areas, then we count number of impacted homes. First, checking for CRS compatibility. Since they aren’t, we need to transform and then subset.

Code 
##check CRS's are the same. No, they're not!
#st_crs(buildings) == st_crs(blackouts_far_hwy_houston)

##transform buildings crs 
buildings <- st_transform(buildings, crs = st_crs(blackouts_far_hwy_houston))

# filtering the houses data with the blackout mask
buildings_blackout <- buildings[blackouts_far_hwy_houston, drop = FALSE]

## View the data frame for inspection 
#View(buildings_blackout)

#count number of impacted homes
print(paste0(nrow(buildings_blackout), " Houston homes experienced power outage due to the Texas storms in Feburary, 2021"))
[1] "164867 Houston homes experienced power outage due to the Texas storms in Feburary, 2021"

Investigate socioeconomic factors

load ACS data

To set up for loading, we use st_read() to load the geodatabase layers and geometries are stored in the ACS_2019_5YR_TRACT_48_TEXAS layer. The income data is stored in the X19_INCOME layer, so we select the median income field B19013e1 andreproject data to EPSG:3083.

Determine which census tracts experienced blackouts

We join the income data to the census tract geometries by joining by geometry ID and spatially join census tract data with buildings determined to be impacted by blackouts, then find which census tracts had blackouts.

Code 
#determine which census tracts experienced blackouts 
#Join the income data to the census tract geometries. Join by geometry ID
median_income_df <- median_income %>% rename(GEOID_Data = GEOID,
         median_income = B19013e1) 


##check types to varify that they are the same class 
#class(median_income_df$GEOID_Data) == class(census$GEOID_Data)

joined_income_census <- left_join(census, 
                         median_income_df, 
                         by = "GEOID_Data")


## filter census data by buildings determined to be impacted by blackouts. 

# transforming both objects to the correct crs
joined_income_census <- st_transform(joined_income_census, crs = st_crs(census))
buildings_blackout <- st_transform(buildings_blackout, crs = st_crs(census))

#check that they are the same crs 
#st_crs(joined_income_census) == st_crs(buildings_blackout)

## filtering and mutating to add a blackout column for affected tracts 
census_blackout <- joined_income_census[buildings_blackout,]%>% mutate(blackout = 'Y')

print(paste0(nrow(census_blackout), " census tracts had blackouts"))
[1] "778 census tracts had blackouts"

Since we have blackout status of Y indicated, we will want to indicate “N” for no-blackout areas as well, for full dataset.

Comparing incomes of impacted tracts to unimpacted tracts

To do this, we create a map using tmap of median income by census tract, designating which tracts had blackouts. We then plot the distribution of income in impacted and unimpacted tracts.

Code 
#join income data to the census tract geometries and crop to region of interest
census_income_geom <- left_join(census, median_income_df, 
                         by = "GEOID_Data") %>% st_transform(3083) #this is texas#
houston_border_sf <- houston_border_sf %>% st_transform(3083)


#spatial crop joined census data to houston border 
census_income_geom_cropped <- census_income_geom[houston_border_sf, op = st_intersects] 


## use tmap to plot income distribution and affected areas. Make graph aesthetics. 
tmap_mode("plot")
tm_shape(census_income_geom_cropped) + 
  tm_graticules() + 
  tm_polygons("median_income", palette= "-viridis", title = "Median income") + 
  tm_shape(census_blackout) + tm_dots(size = 0.4, alpha = 0.5) + 
  tm_title( "Distribution of Med. Income & blackout status in Houston") +
  tm_scalebar(position = c("RIGHT", "BOTTOM")) + 
  tm_compass(position = c("LEFT", "BOTTOM")) + 
  tm_xlab("Longitude", size = 0.5) + 
  tm_ylab("Latitude", rotation = 90, size = 0.5)

Evaluate the distribution of median incomes of homes by blackout status

Plot a histogram of median income for blackout or non-blackout regions.

Code 
#Plot from one dataset so that ggplot can plot both datasets on one axis & autogenerate legend
census_Y_N_houston <- census_Y_N %>% st_transform(3083) 
houston_border_sf <- houston_border_sf %>% st_transform(3083)
census_Y_N_houston <- census_Y_N_houston[houston_border_sf, op = st_intersects] 


ggplot(census_Y_N_houston, aes(median_income)) + geom_histogram(aes(color = blackout, fill = blackout), bins = 50) + ggtitle("Distribution of median income of homes in Houston by blackout status")
Warning: Removed 10 rows containing non-finite values (`stat_bin()`).

I decided that a box and whisker plot might be helpful to compare the distribution as well.

Code 
# plotting the comparison data via geom jitter plot
ggplot(census_Y_N_houston, aes(x = median_income, y = blackout)) +
  geom_boxplot(aes(color = blackout), alpha = 0.8) +
  labs(title = "Median Income of Houston Homes based on Blackout Status",
       y = "Blackout Status",
       x = "Median Income") +
  theme_minimal()
Warning: Removed 10 rows containing non-finite values (`stat_boxplot()`).

And included some summary statistics.

Code 
### Summarizing statistics 

census_Y_N_houston %>% 
  group_by(blackout) %>% 
  summarize(mean = mean(median_income, na.rm = TRUE),
            median = median(median_income, na.rm = TRUE),
            sd = sd(median_income, na.rm = TRUE))
Simple feature collection with 2 features and 4 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 1808144 ymin: 7178529 xmax: 2053896 ymax: 7423349
Projected CRS: NAD83 / Texas Centric Albers Equal Area
# A tibble: 2 Ă— 5
  blackout   mean median     sd                                            Shape
  <chr>     <dbl>  <dbl>  <dbl>                               <MULTIPOLYGON [m]>
1 N        67859. 57858. 36526. (((1833927 7192789, 1833903 7193076, 1833766 71…
2 Y        70939. 59916  39311. (((1895641 7236524, 1896074 7237007, 1896410 72…

Discussion

We found that 778 census tracts had been affected by Texas’s 2021 energy crisis, amounting to 164,867 total homes experiencing power outage between the studied dates (Feb 07, 2021 to Feb 16, 2021).

The average median income for homes that experienced a blackout was $70,939, slightly higher than the average median income for homes that didn’t experience a blackout, which was $67,859. It is important to note that 1) our factor of interest was median income within census tracts and that other socioeconomic factors could reveal unexplored spatial patterns, and that 2) weather conditions were not normalized to conduct this exercise, variable that could cause a different outcome.

References

[1] Ball, J. (Feb, 2021). The Texas Blackout is the Story of a Disaster Foretold. Texas Monthly. URL: https://www.texasmonthly.com/news-politics/texas-blackout-preventable/

[2] Henson, Bob. (Feb, 2021). Why the power is out in Texas… and why other states are vulnerable too. Yale Climate Connections. URL: https://yaleclimateconnections.org/2021/02/why-the-power-is-out-in-texas-and-why-other-states-are-vulnerable-too/

[3] Irfan, U. (Mar, 2021). Why every state is vulnerable to a Texas-style power crisis. Vox. URL: https://www.vox.com/22308149/texas-blackout-power-outage-winter-uri-grid-ercot

[4] National Center for Disaster Preparedness (NCDP). (Mar, 2023). Disaster Response and Equity: Reflecting on the Racial Disparities in Texas Power Outages. URL: https://ncdp.columbia.edu/ncdp-perspectives/disaster-response-and-equity-texas-power-outrages/