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.
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:
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.
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.
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 datndefinedandefine
# 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 formattingThe analysis follows these steps:
# 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.
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)")map <- 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)
map
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.
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)
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)
This refinement helps ensure we’re identifying genuine residential power outages rather than changes in traffic patterns.
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)
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.
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")
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())
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")| 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.
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).
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.
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.
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.
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
@online{miller2025,
author = {Miller, Emily},
title = {Mapping the {Texas} {Blackouts:} {A} {Spatial} {Analysis} of
the 2021 {Winter} {Storm} {Impacts}},
date = {2025-12-03},
url = {https://rellimylime.github.io/posts/eds223-final/},
langid = {en}
}