Identifying the Impacts of Extreme Weather

Research Question

To what extent did the February 2021 Texas winter storm blackouts disproportionately affect low-income communities in the Houston metropolitan area?

Background

Climate change is increasing the frequency and intensity of extreme weather events, with devastating impacts on infrastructure and communities. In February 2021, Texas experienced a major power crisis resulting from three severe winter storms that swept across the United States between February 10-20, leaving millions without power for days in freezing temperatures¹.

In this blog post, I identify the impacts of these series of extreme winter storms by estimating the number of homes in the Houston metropolitan area that lost power and investigate whether these impacts were disproportionately felt by low-income populations.

Data

Sources

NIGHT LIGHTS DATA

Imagery from the Visible Infrared Imaging Radiometer Suite (VIIRS) distributed in 10x10 degree tiles in sinusoidal equal-area projection. Houston is included in two tiles: h08v05 and h08v06. We accessed tiles through NASA’s Level-1 and Atmospheric Archive & Distribution System Distributed Active Archive Center (LAADS DAAC)

ROADS DATA

OpenStreetMap highway data for Texas, accessed via Geofabrik’s download site, was used to retrieve a shapefile of all highways in Texas, and a Geopackage (.gpkg file) was prepared containing only the subset of roads that intersect the Houston metropolitan area.

BUILDING DATA

OpenStreetMap building data for Texas, accessed via Geofabrik’s download site, was used to retrieve a shapefile of all houses in Texas, and a Geopackage (.gpkg file) was prepared containing only the subset of houses that intersect the Houston metropolitan area.

SOCIOECONOMIC DATA

From the U.S. Census Bureau’s American Community Survey for census tracts in 2019, which contains social, economic, housing, and demographic information across the 50 states of the U.S., the District of Columbia, and Puerto Rico. The dataset can be found here: Census

Analysis

To investigate the impact of the February 2021 winter storms on Houston, I analyzed satellite imagery and socioeconomic data through three main steps: identifying blackout locations, estimating affected homes, and examining socioeconomic patterns.

Step 1: Identifying Blackout Locations

I began by comparing nighttime light intensity before and after the storm using NASA’s satellite data.

Load Packages and Import the datasets

# Load libraries 
library(tidyverse)
library(stars)
library(sf)
library(here)
library(dplyr)
library(raster) 
library(terra)
library(tmap)
library(paletteer)

# Read raster tiles for both dates
rt_0705 <- read_stars(here("data", "VNP46A1",
                           "VNP46A1.A2021038.h08v05.001.2021039064328.tif"))
rt_0706 <- read_stars(here("data", "VNP46A1", 
                           "VNP46A1.A2021038.h08v06.001.2021039064329.tif"))
rt_1605 <- read_stars(here("data", "VNP46A1", 
                           "VNP46A1.A2021047.h08v05.001.2021048091106.tif"))
rt_1606 <- read_stars(here("data", "VNP46A1", 
                           "VNP46A1.A2021047.h08v06.001.2021048091105.tif"))

# Read roads data
roads <- st_read(here("data",
       "gis_osm_roads_free_1.gpkg"),
  query = "SELECT * FROM gis_osm_roads_free_1 WHERE fclass='motorway'") %>% 
  st_transform(crs = "epsg:3083")

# Read building data
buildings <- st_read(here("data",
                          "gis_osm_buildings_a_free_1.gpkg"),
  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')")

# Read the layers holding the Texas geometry information
texas_geom <- st_read(here("data/ACS_2019_5YR_TRACT_48_TEXAS.gdb"),
                layer = "ACS_2019_5YR_TRACT_48_TEXAS")

texas_income <- st_read(here("data/ACS_2019_5YR_TRACT_48_TEXAS.gdb"),
                         layer = "X19_INCOME",
                         query = "SELECT B19013e1, GEOID FROM X19_INCOME") %>%
  mutate(GEOID_Data = GEOID)

Visualization of Night Light Changes

Maps comparing night light intensities before and after the storms

# Merge tiles from February 7, 2021
mosaic_7 <- st_mosaic(rt_0705, rt_0706)
# Merge tiles from February 16, 2021
mosaic_16 <- st_mosaic(rt_1605, rt_1606)

# Create color palette for the night lights
pal_lights <- c("#000011", "#000066", "#2222AA", "#FFFF66", "#FFFFFF")

# Plot raster of night lights before the storms
map_before <- tm_shape(mosaic_7) +
  tm_raster(palette = pal_lights,
            style = "quantile", 
            title = "Light Intensity(nW cm⁻²sr⁻¹)") +
  tm_title(text = "Night Lights\nFebruary 07, 2021 (Before storm)", size = 0.7)+
  tm_graticules()+
  tm_layout(inner.margins = c(0.001, 0.001, 0.001, 0.001),
            component.autoscale = FALSE,
            legend.text.size = 0.5,
            legend.title.size = 0.5,)

# Plot raster of night lights after the storms
map_after <- tm_shape(mosaic_16) +
  tm_raster(palette = pal_lights,
            style = "quantile", 
            title = "Light Intensity(nW cm⁻²sr⁻¹)") +
  tm_title(text = "Night Lights\nFebruary 16, 2021 (After storm)", size = 0.7)+
  tm_graticules() +
  tm_layout(inner.margins = c(0.001, 0.001, 0.001, 0.001),
            component.autoscale = FALSE,
            legend.text.size = 0.5,
            legend.title.size = 0.5,)

tmap_arrange(map_before, map_after, nrow = 1)

To isolate areas that truly experienced blackouts, I calculated the change in light intensity between these dates. Areas showing a drop greater than 200 nW cm⁻²sr⁻¹ were classified as experiencing blackout.

Calculate light intensity difference

# Subtract light intensity (February 7 - February 16)
difference <- mosaic_7 - mosaic_16

#  keep only locations with light drops > 200 nW cm⁻²sr⁻¹
lights_mask_diff <- difference
lights_mask_diff[difference < 200] <- NA

# Vectorize mask and ensure valid geometry
mask_vec <- st_as_sf(lights_mask_diff, as_points = FALSE, merge = TRUE) %>% 
  st_make_valid() %>% 
  st_as_sf(crs = 'EPSG:3083')

I then focused the analysis on the Houston metropolitan area and excluded highways, which typically remain lit even during residential blackouts, to ensure the results reflected actual residential impacts.

Selecting Houston Area

# Define coordinates of Houston
lon <- c(-96.5, -96.5, -94.5, -94.5)
lat <- c(29, 30.5, 30.5, 29)

# Create bbox from min and max of lon and lat
houston_bbox <- st_bbox(c(
  xmin = min(lon),
  ymin = min(lat),
  xmax = max(lon),
  ymax = max(lat)), crs = st_crs(mask_vec))

# Crop the blackout mask to the Houston bounding box
cropped_blackout <-  st_crop(mask_vec, houston_bbox) %>% 
st_make_valid() %>% 
  st_transform(cropped_blackout, crs = "epsg:3083")

Removing highways

# create a 200m buffer around highways and dissolve overlapping areas
roads_buffer <- st_buffer(roads, dist = 200) %>% 
st_union

# Check if the datasets have the same coordinate reference system.
# If not, print a warning and transform the second dataset to match the first.
if(st_crs(cropped_blackout) != st_crs(roads_buffer)){
  warning("coordinate reference systems do not match")
  roads_buffer <- st_transform(roads_buffer, crs = st_crs(cropped_blackout))}

Extract blackout areas outside highways

# Ensure the blackout layer has valid geometries
cropped_blackout <- st_make_valid(cropped_blackout)

# Ensure the roads buffer layer has valid geometries
roads_buffer <- st_make_valid(roads_buffer)

# Subtract the roads buffer from the blackout
blackout_not_road <- st_difference(cropped_blackout, roads_buffer)

# Check if the datasets have the same coordinate reference system.
# If not, print a warning and transform the second dataset to match the first.
if(st_crs(blackout_not_road) != st_crs(buildings)){
  warning("coordinate reference systems do not match")
  buildings <- st_transform(buildings, crs = st_crs(blackout_not_road))}

# Select homes that overlap with blackout_not_road.
homes_blackout <- buildings[blackout_not_road, 
                     op = st_intersects, 
                     drop = FALSE] 

Map of the homes in Houston that lost power

# Plot a map of homes affected by the blackout
map_homes_blackout <- tm_shape(homes_blackout) +
  tm_polygons(col = "#480091") +
  tm_title(text = "Homes in Houston that lost power") +
  tm_graticules() +
  tm_basemap("OpenStreetMap", alpha = 0.7) +
  tm_layout(inner.margins = c(0.01, 0.01, 0.01, 0.1),
            component.autoscale = FALSE) +
  tm_add_legend(type = "fill",
                labels = "Homes Affected",
                col = "#480091",
                title = "Legend")

map_homes_blackout

The map above shows the spatial distribution of affected homes, with a notable concentration in central Houston. This visualization helps us understand not just how many homes were impacted, but where the impacts were most severe.

Step 2: Estimating number of Affected Homes

Using building data from OpenStreetMap, I identified residential structures within the blackout zones. The analysis revealed that approximately 157,970 homes lost power during the storm.

# Count the number of homes that lost power
length(homes_blackout$fclass)

[1] 157970

Step 3: Exploring Socioeconomic Patterns

I wanted to investigate whether the blackouts disproportionately affected lower-income communities. To explore this, I overlaid the blackout data with census tract information, examining median household income across affected and unaffected areas.

Prepare median income data

# Left join geometry and median income layers
inc_geom <- left_join(texas_geom, texas_income, by = "GEOID_Data") %>% 
  st_transform(crs = "epsg:3083")

# Transform CRS
houston_bbox <- houston_bbox %>% 
  st_transform(crs = "epsg:3083")

# Crop income layer to the Houston bounding box
houston_income <- st_crop(inc_geom, houston_bbox) %>% 
st_make_valid() %>% 
  st_transform(houston_income, crs = "epsg:3083")

Map of the census tracts in Houston that lost power

# Verify that both datasets use the same CRS and add a stop error if they don't
if(st_crs(houston_income) != st_crs(homes_blackout)){
    stop("coordinate reference systems must match to join datasets")}

# Intersect census tracts with blackout areas
homes_census_blackout <- houston_income[homes_blackout,
                               op = st_intersects]

# Create map showing median income for blackout-affected census tracts
tm_shape(homes_census_blackout) +
  tm_fill(
    col = "B19013e1",
    palette = "Reds",              
    title = "Median Income") + 
 tm_basemap(server = "OpenStreetMap", alpha = 0.7) +
 tm_graticules() +
 tm_title(text = "Median Income by Census Tract for Houses
          Affected by the Blackout in Houston" )

The map reveals the spatial distribution of median household income across census tracts affected by the blackout. The darker red areas indicate higher-income neighborhoods, while lighter shades represent lower-income areas.

To quantify whether income played a role in blackout vulnerability, I compared all Houston census tracts, categorizing them by whether they experienced blackouts or not.

Find census tracts that did not experience blackout

# Compare all Houston tracts against those identified with blackouts,
# Remove geometry and select only income and status variables
census_combined <- houston_income %>%
  mutate(blackout_status = if_else(GEOID_Data %in% homes_census_blackout$GEOID_Data, 
                                   "Blackout", "No Blackout")) %>%
  st_drop_geometry() %>%
  dplyr::select(B19013e1, blackout_status) %>%
  drop_na()

Visualize the distribution of median household income across census tracts to determine if the blackout disproportionately impacted certain socioeconomic groups.

# Plot income distribution by blackout status
ggplot() +
  geom_boxplot(data = census_combined,
           aes(x = blackout_status, 
               y = B19013e1,
               fill = blackout_status)) +
  scale_x_discrete(labels=c("Blackout", "No Blackout")) +
  scale_fill_manual(values = c("#2955AA", "#F9cF66"),
                    name = "Blackout Status") +
  labs(x = "Affected Status",
       y = "Median Income of Census Tract",
       title = "Distribution of Median Income in Areas Affected 
       and Unaffected by Blackout in Houston") +
  theme_minimal()

The boxplot reveals similar income distributions between areas that experienced blackouts and those that did not. The yellow box representing unaffected areas shows a slightly lower median income compared to the blue box representing blackout-affected census tracts.

Conclusion and Next Steps

The results show that median household incomes are similar between census tracts that experienced blackouts and those that did not, suggesting that income level was not a strong factor in determining which areas lost power during the storm. However, we need to consider some limitations of this result, as they relate to the spatial scale and join methods used in our analysis. Some possible bias in the results may be related to the concentration of buildings in the central part of Houston affected by blackouts, as shown in the Homes Map, and the reliance on tract-level overlaps that generalize patterns and hide smaller spatial or demographic variations.

Future research should employ finer spatial resolution (e.g., block group or individual building level) and could incorporate additional socioeconomic variables beyond income, including housing age, renter vs. owner occupancy, and racial demographics.

References

Home - LAADS DAAC. (2017, February 22). Nasa.gov. https://ladsweb.modaps.eosdis.nasa.gov/

Download OpenStreetMap for OpenStreetMap Data Extracts | Geofabrik Download Server. (2025). Geofabrik Download Server. https://download.geofabrik.de/

Bureau, U. C. (2025, September 29). American Community Survey (ACS). Census.gov. https://www.census.gov/programs-surveys/acs

DEITABASE. (2024, December 27). 2021 Texas Power Grid Failure – a preventable disaster. DEITABASE - a Community Resource for Merging DEI and CEE Education. https://limos.engin.umich.edu/deitabase/2024/12/27/2021-texas-power-grid-failure/

Footnotes

1. 2021 Texas Power Grid Failure – a preventable disaster. DEITABASE. Last modified 2024, December 27. https://limos.engin.umich.edu/deitabase/2024/12/27/2021-texas-power-grid-failure/.