Irregular Outages

Introduction

In February 2021, a historic winter storm swept across Texas, bringing record-low temperatures and overwhelming the state’s power grid. At the peak of the crisis, more than 4.5 million homes and businesses lost power, leaving residents without electricity for heating during freezing temperatures. The event resulted in at least 246 deaths and an estimated $195 billion in economic damages, making it one of the costliest natural disasters in U.S. history.

Research Question

This analysis addresses a critical question: Which communities in the Houston metropolitan area were most affected by the blackouts, and were there disparities in impacts across different income levels?

Understanding the spatial distribution of power outages and their relationship to socioeconomic factors is essential for:

• Improving power grid resilience and emergency response planning
• Identifying vulnerable communities that may need additional support during future crises
• Informing equitable infrastructure investment decisions

Background

The 2021 Texas Winter Storm Crisis

The February 2021 winter storm, known as Winter Storm Uri, exposed critical vulnerabilities in Texas’s isolated power grid. Unlike other states, Texas operates its own grid (ERCOT), which is not connected to neighboring states, limiting backup power options during emergencies. The unprecedented cold caused natural gas pipelines to freeze, wind turbines to ice over, and power plants to fail, creating a cascading crisis.

While the blackouts affected millions of Texans, their impacts were not uniformly distributed. Previous research on climate disasters has shown that low-income communities and communities of color often experience disproportionate impacts from extreme weather events. These communities may face longer recovery times, limited access to emergency resources, and inadequate infrastructure.

Remote Sensing for Disaster Assessment

Satellite-based night lights data provides a powerful tool for assessing power outages at scale. NASA’s Visible Infrared Imaging Radiometer Suite (VIIRS) captures daily nighttime light emissions, allowing researchers to detect changes in electrical grid function. By comparing light intensity before and after disasters, we can identify affected areas without relying solely on utility company reports, which may be incomplete during major crises.

This analysis combines VIIRS night lights data with OpenStreetMap building footprints and U.S. Census socioeconomic data to map blackout impacts and explore potential environmental justice dimensions of the 2021 Texas blackouts.

Data and Methods

Data Sources

This analysis integrates four primary datasets:

VIIRS Night Lights Data - NASA’s VIIRS instrument provides daily measurements of nighttime light intensity. I used two dates of imagery: February 7, 2021 (before the storm) and February 16, 2021 (during the height of the blackouts). The data is measured in nanowatts per square centimeter per steradian (nW cm⁻²sr⁻¹) (NASA LAADS DAAC 2021).

OpenStreetMap Building Footprints - Residential building polygons for the Houston metropolitan area, including houses, apartments, and other residential structures (OpenStreetMap Contributors 2021).

OpenStreetMap Road Network - Highway and major road data used to filter out areas where reduced light may be due to decreased traffic rather than power outages (OpenStreetMap Contributors 2021).

U.S. Census American Community Survey - 2015-2019 five-year estimates at the census tract level, including median household income data (U.S. Census Bureau 2019).

Analytical Approach

# Load required packages
library(tidyverse)      # Data manipulation and visualization
library(sf)             # Spatial data handling
library(here)           # File path management
library(terra)          # Raster data processing
library(stars)          # Spatiotemporal arrays
library(tmap)           # Thematic mapping
library(kableExtra)     # Table formatting

The analysis follows these steps:

1. Load and prepare spatial data - Read VIIRS night lights tiles, building footprints, road networks, and census data
2. Calculate light intensity changes - Compare pre- and post-storm imagery to identify areas with significant light reduction
3. Define blackout threshold - Classify areas with light intensity drops >200 nW cm⁻²sr⁻¹ as experiencing blackouts
4. Filter highway corridors - Exclude areas within 200m of major highways to avoid misidentifying reduced traffic as power outages
5. Identify impacted buildings - Intersect blackout areas with residential building footprints
6. Analyze socioeconomic patterns - Compare median household income distributions between census tracts with and without blackouts

Loading and Preparing the Data

# Load VIIRS night lights raster tiles for Houston area
# February 7, 2021 (pre-storm)

VIIRS_07a <- rast(here("posts/eds223-final/data/VNP46A1/VNP46A1.A2021038.h08v05.001.2021039064328.tif"))
VIIRS_07b <- rast(here("posts/eds223-final/data/VNP46A1/VNP46A1.A2021038.h08v06.001.2021039064329.tif"))

# February 16, 2021 (during blackouts)
VIIRS_16a <- rast(here("posts/eds223-final/data/VNP46A1/VNP46A1.A2021047.h08v05.001.2021048091106.tif"))
VIIRS_16b <- rast(here("posts/eds223-final/data/VNP46A1/VNP46A1.A2021047.h08v06.001.2021048091105.tif"))

# Load OpenStreetMap highways (motorways only)
highways <- st_read(here("posts/eds223-final/data/gis_osm_roads_free_1.gpkg/gis_osm_roads_free_1.gpkg"),
                    query = "SELECT * FROM gis_osm_roads_free_1 WHERE fclass='motorway'",
                    quiet = TRUE)

# Load OpenStreetMap residential buildings
houses <- st_read(here("posts/eds223-final/data/gis_osm_buildings_a_free_1.gpkg/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')",
                  quiet = TRUE)

# Load Census ACS data - requires joining geometry with income attributes
socio_geom <- st_read(here("posts/eds223-final/data/ACS_2019_5YR_TRACT_48_TEXAS.gdb/"),
                      layer = "ACS_2019_5YR_TRACT_48_TEXAS", quiet = TRUE)

socio_income <- st_read(here("posts/eds223-final/data/ACS_2019_5YR_TRACT_48_TEXAS.gdb/"),
                        layer = "X19_INCOME", quiet = TRUE)

# Join income data to geometries and reproject to match building data
socio <- socio_geom %>%
  rename("GEOID_Data" = "GEOID_Data") %>%
  left_join(socio_income %>% rename("GEOID_Data" = "GEOID"), by = "GEOID_Data") %>%
  st_transform(st_crs(houses))

All datasets use different coordinate reference systems initially, so I standardized them to EPSG:3083 (NAD83 / Texas Centric Albers Equal Area), which preserves area measurements needed for accurate spatial analysis.

Results

Visualizing the Blackouts: Before and After

The stark difference in nighttime light emissions reveals the massive scale of the power outages across Houston.

# Merge VIIRS tiles for each date
VIIRS_07 <- merge(VIIRS_07a, VIIRS_07b)  # Pre-storm
VIIRS_16 <- merge(VIIRS_16a, VIIRS_16b)  # During blackouts

# Define Houston bounding box
HOUSTON_BBOX <- st_bbox(c(xmin = -96.5, ymin = 29,
                          xmax = -94.5, ymax = 30.5),
                        crs = st_crs(VIIRS_07))

# Crop to Houston area
crop_07 <- crop(VIIRS_07, HOUSTON_BBOX)
crop_16 <- crop(VIIRS_16, HOUSTON_BBOX)

# Mask unrealistic values (>1000 or ≤0)
masked_07 <- crop_07
masked_07[(masked_07 > 1000) | (masked_07 <= 0)] <- NA

masked_16 <- crop_16
masked_16[(masked_16 > 1000) | (masked_16 <= 0)] <- NA

# Combine for side-by-side comparison
combined_crop <- c(masked_07, masked_16)
names(combined_crop) <- c("February 7, 2021 (Before)", "February 16, 2021 (After)")

tm_shape(combined_crop) +
  tm_raster(col.scale = tm_scale(values = "inferno"),
            col.legend = tm_legend(title = "Light Intensity\n(nW cm⁻²sr⁻¹)",
                                   title.size = 0.9,
                                   text.size = 0.6)) +
  tm_facets(ncol = 2, free.scales = FALSE) +
  tm_title("Houston Night Lights: Before and After the Storm", size = 1.2) +
  tm_layout(legend.outside = TRUE,
            legend.outside.position = "right",
            legend.outside.size = 0.12,
            panel.labels = c("February 7 (Before)", "February 16 (After)"),
            fontfamily = "sans") +
  tm_graticules(lwd = 0.2, col = "white", alpha = 0.3, labels.size = 0.5)

Night lights imagery showing Houston before (left) and during (right) the February 2021 blackouts. Darker areas in the right panel indicate regions that lost power. The widespread darkening across the metropolitan area illustrates the extensive scale of the grid failure.

The comparison reveals significant light reduction across much of the Houston metropolitan area. The bright urban core visible on February 7 dims substantially by February 16, indicating widespread power failures.

Identifying Blackout Areas

To systematically identify blackout zones, I calculated the difference in light intensity and classified areas with drops greater than 200 nW cm⁻²sr⁻¹ as experiencing blackouts.

# Calculate light intensity difference
VIIRS_houston_diff <- VIIRS_07 - VIIRS_16

# Classify blackout areas (drops >200 as blackouts)
blackout_mask <- classify(VIIRS_houston_diff,
                         matrix(c(-Inf, 200, NA,  # Not a blackout
                                  200, Inf, 1),   # Blackout
                                ncol = 3, byrow = TRUE))

# Convert raster to vector polygons
blackout_vector <- st_as_stars(blackout_mask) %>%
  st_as_sf() %>%
  st_make_valid()

# Crop to Houston and reproject
blackout_houston <- st_crop(blackout_vector, HOUSTON_BBOX) %>%
  st_transform(crs = 3083)

road_color <- "#FDB863"
fill_color <- "#d62728"
border_color <- "#8b0000"

blackout_houston_union <- st_union(blackout_houston)

tm_basemap("CartoDB.VoyagerNoLabels") +
  tm_shape(blackout_houston_union) +
  tm_polygons(fill = fill_color,
              col = border_color,
              lwd = 0.5) +
  tm_title("Detected Blackout Areas in Houston", size = 1.2) +
  tm_add_legend(type = "polygons",
                fill = fill_color,
                col = border_color,
                border.lwd = 0.5,
                labels = "Blackout Areas",
                title = "") +
  tm_layout(legend.position = c("right", "bottom"),
            legend.bg.color = "white",
            legend.bg.alpha = 0.8,
            text.fontfamily = "sans") +
  tm_scalebar(position = c("left", "bottom"), text.size = 0.6) +
  tm_compass(position = c("right", "top"), size = 1.5)

Initial identification of blackout areas based on night lights analysis. Red areas experienced light intensity drops exceeding 200 nW cm⁻²sr⁻¹, indicating likely power outages.

Refining the Analysis: Excluding Highways

Major highways may show reduced light intensity due to decreased traffic rather than residential power outages. To improve accuracy, I excluded areas within 200 meters of motorways.

# Reproject highways and create 200m buffer
highways <- st_transform(highways, crs = 3083)
highway_buffer <- highways %>%
  st_union() %>%
  st_buffer(dist = 200)

# Remove highway buffers from blackout areas
blackout_final <- st_difference(blackout_houston, highway_buffer)

# Visualize the highway exclusion process
tm_basemap("CartoDB.VoyagerNoLabels") +
  tm_shape(highways) +
  tm_lines(col = road_color,
           lwd = 2) +
  tm_shape(st_union(blackout_final)) +
  tm_polygons(fill = fill_color,
              col = "#4d0000",
              lwd = 0.5) +
  tm_add_legend(type = "polygons",
                labels = "Refined blackout areas (highways excluded)",
                fill = fill_color) +
  tm_add_legend(type = "lines",
                labels = "Major highways",
                col = road_color,
                lwd = 2) +
  tm_title("Highway Corridor Exclusion", size = 1.2) +
  tm_layout(legend.position = c("left", "top"),
            legend.bg.color = "white",
            legend.bg.alpha = 0.85,
            legend.text.size = 0.65,
            text.fontfamily = "sans") +
  tm_scalebar(position = c("right", "bottom"), text.size = 0.6) +
  tm_compass(position = c("right", "top"), size = 1.5)

Refinement of blackout areas by excluding highway corridors. The yellow lines show major highways with 200m buffers applied. This step removes areas where reduced light may be due to decreased traffic rather than residential power outages.

This refinement helps ensure we’re identifying genuine residential power outages rather than changes in traffic patterns.

Quantifying Residential Impacts

By intersecting the blackout areas with residential building footprints, we can estimate how many homes lost power.

# Reproject buildings and identify those in blackout zones
houses <- st_transform(houses, crs = 3083)
blackout_homes <- st_intersection(houses, blackout_final)

# Calculate statistics
n_blackout_homes <- length(unique(blackout_homes$osm_id))
n_total_homes <- nrow(houses)
pct_impacted <- round((n_blackout_homes / n_total_homes) * 100, 1)

Finding: Approximately 157,410 residential buildings (33.1% of structures in the dataset) experienced blackouts during the February 2021 storms.

# Create point representations for mapping
all_homes_points <- st_centroid(houses)
not_impacted_homes <- all_homes_points[!all_homes_points$osm_id %in% blackout_homes$osm_id, ]
impacted_homes_points <- st_centroid(blackout_homes)

tm_basemap("CartoDB.VoyagerNoLabels") +
  tm_shape(blackout_final) +
  tm_polygons(fill = "grey85", col = "grey70", lwd = 0.3) +
  tm_shape(not_impacted_homes) +
  tm_dots(fill = "#FFB3B3", size = 0.003, fill_alpha = 0.1) +
  tm_shape(impacted_homes_points) +
  tm_dots(fill = "#8B0000", size = 0.005, fill_alpha = 0.7) +
  tm_title("Residential Buildings Impacted by Blackouts", size = 1.3) +
  tm_layout(legend.outside = FALSE,
            legend.position = c("left", "top"),
            legend.bg.color = "white",
            legend.bg.alpha = 0.85,
            legend.text.size = 0.7,
            legend.title.size = 0.9,
            fontfamily = "sans") +
  tm_add_legend(type = "fill",
                labels = "Blackout zone",
                fill = "grey85",
                col = "grey70") +
  tm_add_legend(type = "symbol",
                labels = c("Power lost", "Power maintained"),
                fill = c("#8B0000", "#FFB3B3"),
                shape = 19,
                size = 0.5) +
  tm_scalebar(position = c("right", "bottom"), text.size = 0.6) +
  tm_compass(position = c("right", "top"), size = 1.5)

Spatial distribution of impacted and non-impacted residential buildings. Dark red points represent homes that experienced blackouts, while light pink points indicate buildings that maintained power throughout the crisis.

The map reveals that blackouts were widespread but not uniform. Some neighborhoods maintained power while adjacent areas went dark, suggesting the complexity of the grid failures.

Socioeconomic Analysis: Income and Blackout Exposure

Understanding where blackouts occurred is only part of the story. To assess whether impacts were distributed equitably, I aggregated affected buildings by census tract and examined their relationship to median household income. This tract-level analysis reveals whether certain communities bore a disproportionate burden during the crisis.

# Transform census data to analysis CRS
socio <- st_transform(socio, crs = 3083)

# Crop to Houston area
houston_tracts <- st_crop(socio, st_bbox(blackout_houston))

# Count impacted homes per tract
houston_tracts$n_homes <- as.numeric(lengths(st_intersects(houston_tracts, blackout_homes)))

# Create impact categories for mapping
houston_tracts$impact_category <- cut(houston_tracts$n_homes,
                                       breaks = c(0, 1, 10, 50, 150, 500, 1500, Inf),
                                       labels = c("None",
                                                 "Very Low (1-10)",
                                                 "Low (11-50)",
                                                 "Medium (51-150)",
                                                 "High (151-500)",
                                                 "Very High (501-1,500)",
                                                 "Extreme (>1,500)"),
                                       include.lowest = TRUE)

# Create binary blackout indicator
houston_tracts$blackout <- factor(houston_tracts$n_homes > 0,
                                  levels = c(TRUE, FALSE),
                                  labels = c("Blackout", "No Blackout"))

# Transform tracts to match crs of basemap (this is done automatically by tmap, code included for transparency)
houston_tracts_proj <- st_transform(houston_tracts, crs = 3857)

# Plot houses in census tracts affected by blackouts
tm_basemap("CartoDB.VoyagerNoLabels") +
  tm_shape(houston_tracts_proj) +
  tm_fill(fill = "impact_category",
          fill.scale = tm_scale(values = c("white", "#feedde", "#fdd0a2", "#fdae6b",
                                          "#fd8d3c", "#e6550d", "#a63603"),
                                value.na = "transparent"),
          fill.legend = tm_legend(title = "Impacted Homes\nper Census Tract",
                                  title.size = 0.9,
                                  text.size = 0.65)) +
  tm_borders(col = "grey50", lwd = 0.3) +
  tm_title("Blackout Severity by Census Tract", size = 1.3) +
  tm_scalebar(position = c("left", "bottom"), text.size = 0.6) +
  tm_compass(position = c("left", "top"), size = 1.5) +
  tm_layout(legend.outside = FALSE,
            legend.position = c("right", "top"),
            legend.bg.color = "white",
            legend.bg.alpha = 0.85,
            text.fontfamily = "sans")

Number of residential buildings that lost power by census tract. Darker colors indicate higher numbers of impacted homes. The spatial pattern shows that some tracts experienced extreme impacts while neighboring areas had minimal disruption.

The census tract analysis reveals substantial variation in blackout severity. Some tracts lost power for more than 1,500 residential buildings, while others experienced minimal disruption. This raises a critical question: were certain communities disproportionately affected?

Previous research on climate disasters shows that low-income communities and communities of color often experience greater exposure and more severe impacts. To investigate whether similar patterns emerged during the Texas blackouts, I analyzed the relationship between median household income and blackout exposure.

houston_tracts_no_na <- houston_tracts %>%
  filter(!is.na(B19013e1))

ggplot(houston_tracts_no_na, aes(x = blackout, y = B19013e1, fill = blackout)) +
  geom_boxplot(alpha = 0.75, outlier.alpha = 0.5) +
  scale_fill_manual(name = "Blackout Status",
                    values = c("Blackout" = "#d62728",
                               "No Blackout" = "#7f7f7f")) +
  scale_y_continuous(labels = scales::dollar_format(),
                    breaks = seq(0, 250000, 50000)) +
  labs(title = "Median Household Income by Blackout Status",
       subtitle = "Houston Metropolitan Area Census Tracts (2015-2019 ACS)",
       x = "",
       y = "Median Household Income",
       caption = "Source: U.S. Census Bureau American Community Survey") +
  theme_minimal(base_size = 12, base_family = "sans") +
  theme(legend.position = "none",
        plot.title = element_text(size = 14, face = "bold", hjust = 0),
        plot.subtitle = element_text(size = 11, hjust = 0),
        panel.grid.minor = element_blank())

Median household income comparison between census tracts that experienced blackouts versus those that maintained power. The similar distributions suggest that income was not a strong predictor of blackout exposure during this crisis.

income_stats <- houston_tracts %>%
  st_drop_geometry() %>%
  filter(!is.na(B19013e1)) %>%
  group_by(blackout) %>%
  summarise(
    `Census Tracts` = n(),
    `Median Income` = scales::dollar(median(B19013e1)),
    `Mean Income` = scales::dollar(mean(B19013e1)),
    `Standard Deviation` = scales::dollar(sd(B19013e1))
  )

kable(income_stats,
      caption = "Table 1: Income Statistics by Blackout Status",
      align = c("l", "r", "r", "r", "r")) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE,
                position = "left")

Table 1: Income Statistics by Blackout Status

blackout	Census Tracts	Median Income	Mean Income	Standard Deviation
Blackout	749	$60,435	$71,262.24	$39,649.34
No Blackout	355	$57,385	$67,352.78	$35,900.53

Key Finding: Census tracts that experienced blackouts had a median household income of 60,435, compared to 57,385 for tracts without blackouts. The difference of approximately $3,000 represents only about 5% variation, suggesting that income level was not a strong predictor of blackout exposure in this event.

Discussion

Main Findings

This analysis identified approximately 157,410 residential buildings (33.1% of the study area) that lost power during the February 2021 storms. Unlike typical climate disasters, blackout exposure showed minimal income-based disparities-tracts with and without blackouts had median incomes differing by only ~$3,000 (5%). This likely reflects the systemic nature of the grid failure: when the entire system collapsed, infrastructure vulnerabilities crossed neighborhood boundaries.

However, exposure differs from impact. Low-income communities faced greater difficulty evacuating, accessing backup heating, and recovering-even when blackout rates were similar (Gronlund 2014).

Limitations

Temporal resolution: Two snapshots cannot capture outage duration, which may vary by neighborhood.

Threshold sensitivity: The 200 nW cm⁻²sr⁻¹ cutoff is somewhat arbitrary and may miss some outages.

Data gaps: OpenStreetMap building coverage varies geographically, potentially introducing spatial bias.

Highway exclusion: The 200m buffer may inadvertently exclude legitimate residential blackouts.

Future work should incorporate daily time series data, analyze additional socioeconomic indicators, integrate health outcomes, and ground-truth results with utility records.

Policy Implications

The widespread blackouts indicate systemic vulnerabilities requiring comprehensive infrastructure winterization. Even when disasters affect communities uniformly, emergency response must account for differential vulnerabilities. Texas’s grid isolation limited backup options during the crisis, suggesting value in regional interconnection. Remote sensing approaches can provide rapid, independent disaster assessment to complement utility reporting.

Conclusion

The February 2021 blackouts revealed catastrophic infrastructure vulnerabilities affecting millions. This analysis demonstrates how satellite remote sensing combined with socioeconomic data enables rapid disaster assessment. While blackouts affected Houston relatively uniformly across income levels-likely due to systemic grid failure-this doesn’t mean impacts were equitable. Understanding both exposure patterns and underlying vulnerabilities remains essential for building resilient, just infrastructure as climate change intensifies extreme weather events.

Data Sources and Acknowledgments

VIIRS Night Lights Data (VNP46A1) - NASA’s Visible Infrared Imaging Radiometer Suite, Suomi NPP satellite, tiles h08v05 and h08v06. Accessed via NASA Earthdata: https://worldview.earthdata.nasa.gov/

OpenStreetMap - Road network and building footprint data. Accessed via Geofabrik: https://download.geofabrik.de/

U.S. Census Bureau - American Community Survey 5-Year Estimates (2015-2019), Texas census tracts. Accessed via: https://www.census.gov/programs-surveys/acs

References

Gronlund, Carina J. 2014. “Racial and Socioeconomic Disparities in Heat-Related Health Effects and Their Mechanisms: A Review.” Current Epidemiology Reports 1 (3): 165–73. https://doi.org/10.1007/s40471-014-0014-4.

NASA LAADS DAAC. 2021. “VIIRS/NPP Daily Gridded Day Night Band 15 Arc-Second Linear Lat Lon Grid Night.” NASA Level-1; Atmosphere Archive & Distribution System Distributed Active Archive Center. https://ladsweb.modaps.eosdis.nasa.gov/missions-and-measurements/products/VNP46A1/.

OpenStreetMap Contributors. 2021. “OpenStreetMap Data Extracts.” Geofabrik GmbH. https://download.geofabrik.de/.

U.S. Census Bureau. 2019. “American Community Survey 5-Year Data (2015-2019).” U.S. Census Bureau. https://www.census.gov/programs-surveys/acs/technical-documentation/table-and-geography-changes/2019/5-year.html.