Identifying the Impacts of Extreme Weather

EDS 223 Homework #3

Identifying the impacts of extreme weather in Houston, Texas, 2021

Homework Assignment 3

Background

Climate change is increasing the frequency and intensity of extreme weather events, with devastating impacts. “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.”

Trhis project will 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 not these impacts were disproportionately felt.

This analysis will be based on remotely-sensed night lights data, acquired from the Visible Infrared Imaging Radiometer Suite (VIIRS) onboard the Suomi satellite. In particular, it 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, I will link (spatially join) these areas with OpenStreetMap data on buildings and roads.

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

For the first step of this analysis, we will examine the raster maps that show the lights available at night in Houston, Texas, both before and after the major storm event. The result is a set of maps comparing night light intensities before and after the first two storms.

Note: In order to refine this analysis to the Houston area, I will crop (spatially subset) the blackout mask to the Houston area as defined by the following coordinates: (-96.5, 29), (-96.5, 30.5), (-94.5, 30.5), (-94.5, 29).

# Define the bounding box for the Houston area
houston_bbox <- st_bbox(c(xmin = -96.5, ymin = 29, xmax = -94.5, ymax = 30.5), crs = st_crs(nightlight_feb7))

# Convert the bounding box into an sf object
houston_sf <- st_as_sfc(houston_bbox) %>% 
  st_transform(st_crs(nightlight_feb7))

# Crop nightlight rasters to include just Houston
nightlight_feb7 <- st_crop(nightlight_feb7, houston_bbox)
nightlight_feb16 <- st_crop(nightlight_feb16, houston_bbox)

map_07 <- tm_shape(nightlight_feb7) + 
  tm_raster(palette = c("#010224", "#fffbf4"),
            breaks =c(51, 100, 200, 300, 500, 1000, 2000, 4000, 6000, 8000, Inf)) + 
  tm_layout(bg.color = "#010224",
            legend.show = FALSE,
            main.title = "Nightlights Feb. 7th",
            main.title.position = "center")

map_16 <- tm_shape(nightlight_feb16) + 
  tm_raster(palette = c("#010224", "#fffbf4"),
            breaks =c(51, 100, 200, 300, 500, 1000, 2000, 4000, 6000, 8000, Inf)) + 
  tm_layout(bg.color = "#010224",
            legend.show = FALSE,
            main.title = "Nightlights Feb. 16th",
            main.title.position = "center")

tmap_arrange(map_07, map_16, ncol = 2, asp = 1)

As we can see from these maps, there seems to be a drop in luminescence of the city at night after the storm event.

Create blackout mask

To identify places that experienced a blackout, I will create a “mask” that indicates for each cell of the raster whether or not it experienced a blackout.

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

The result of this analysis will be a map that show a vectorized raster of blackouts in the Houston area.

# Take the Feb 16th data subtract it from Feb 7th nightlight data. The 7th will be brighter...map algebra
nightlight_difference <- nightlight_feb7 - nightlight_feb16

# Assign NAs to all locations with a drop of less than 200
nightlight_difference[nightlight_difference <= 200] <- NA

# Eliminate invalid geometries and convert to efficient polygons
houston_blackouts <- as.polygons(rast(nightlight_difference)) %>%
  st_as_sf() %>% 
  st_make_valid()

#re-project the cropped blackout dataset to EPSG:3083 (NAD83 / Texas Centric Albers Equal Area)
houston_blackouts <- st_transform(houston_blackouts, crs = 3083)

# Make a map of the findings
tm_shape(houston_blackouts) +
  tm_polygons(col = "black",
              border.col = "black") +
  tm_layout(bg.color = "grey25",
  main.title = "Vectorized Raster of Blackouts",
            main.title.position = "center")

Exclude highways from the cropped blackout mask

Highways may have experienced changes in their night light intensities that are unrelated to the storm. Therefore, I will seek to exclude any locations within 200 meters of all highways in the Houston area.

In order to do so, I must first identify areas within 200m of all highways.

NULL

[1] "m"

# First, ensure that CRS of highways matches that of the Houston blackouts data
if (st_crs(roads) != st_crs(houston_blackouts)) {
  warning("Uh oh! Coordinate reference systems need to be transformed to match.")
  roads <- st_transform(roads, crs = st_crs(houston_blackouts))  
  message("Transformation successful! CRSs are now aligned.")
} else {
  message("Huzzah! CRSs already match.")
}

# Create a map-able highway system using st_buffer to understand areas around highways
houston_highways <- st_buffer(roads, dist = 200) %>%
  st_union() %>%
  st_transform(crs = 3083) %>% 
  st_make_valid() %>% 
  st_as_sf()

# Map this highway data
tm_shape(houston_highways) +
  tm_polygons(col = "#fffbf4",
              border.col = "#fffbf4") +
  tm_layout(main.title = "Houston Highways",
            main.title.position = "center",
            bg.color = "#010224")

Now that I have a firm understanding of where the highways are located in Houston, I can create a map-able set of data that will allow me exclude highway areas from my blackouts data, to account for the light provided by cars on highways. The result is a map of the Houston area that show lights at night, excluding highway zones.

no_highways <- st_difference(houston_blackouts, houston_highways)

tm_shape(no_highways) +
  tm_polygons(col = "black",
              border.col = "black") +
  tm_layout(main.title = "Blackout Areas Excluding Highways",
            main.title.position = "center",
            bg.color = "grey25")

Identify homes likely impacted by blackouts

In order to understand what buildings were impacted by these blackouts, I will map out the buildings that fall within these darkened zones of our nighlights data. The resulting map will show us buildings in the areas that were likely impacted by these blackouts. We will see a map of the homes in in Houston that lost power.

# First, make sure the CRS of the houses data frame matches, and get rid of invalid polygons
houses <- houses %>% 
  filter(!is.na(type)) %>% 
  st_transform(houses, crs = 3083) %>% 
  st_make_valid()

# Find homes that intersect with blackout areas far from highways
blackout_buildings <- st_intersection(houses, no_highways)

# Make a map of these affected buildings
tm_shape(no_highways) + 
  tm_polygons(col = "black",
              border.col = "black") +  
  tm_shape(blackout_buildings) + 
  tm_dots(col = "coral", size = 0.1, alpha = 0.3, title = "Impacted Homes") + 
  tm_layout(main.title = "Homes Impacted by Blackout",
            main.title.position = "center")

How many homes in Houston lost power, exactly? Using the data from the above map, we can get an estimate of the number of homes in Houston that lost power.

# Count the number of rows in the blackout_buildings dataset
estimate <- nrow(blackout_buildings)

# Create a data frame from this analysis
summary_impacted <- data.frame(description = "Total Affected Homes", total = estimate)

# Make a fun table that shows this estimate
table1 <- summary_impacted |>
  gt() |>
  gt_theme_nytimes() |>
  tab_header(title = "Estimated Homes Impacted by Blackout") |>
  cols_label(
    description = "Description",
    total = "Total",
  ) |>
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_body()
  )

table1

Estimated Homes Impacted by Blackout
Description	Total
Total Affected Homes	23980

Identify the census tracts likely impacted by blackouts

Now that I have identified the areas showing likely blackouts, excluding possible bias from car data, and mapped out the buildings in this area, I can overlay these maps with census data, which will give us a map of the census tracts in Houston that lost power, and will help us understand if people of lower incomes were more effected by these blackouts.

# Join income df to socio geo data frame on the GEOID(income) and GEOID_Data(socio) columns using a left join
income_geom <- socio %>%
  left_join(income, by = c("GEOID_Data" = "GEOID"))

# Make sure that the blackout(no highways) data and the income/socio data have the same CRS! 
if (st_crs(no_highways) != st_crs(income_geom)) {
  warning("Uh oh! Coordinate reference systems need to be transformed to match.")
  income_geom <- st_transform(income_geom, crs = st_crs(no_highways))
  print("Transformation successful!")
} else {
  print("Huzzah! Our transformation worked!")
}

[1] "Transformation successful!"

# Identify census tracts that that contained homes that experienced blackouts
# Directly filter income_geom by intersections with no_highways
tracts_affected <- st_filter(income_geom, no_highways)

# Create the map
census_map <- tm_shape(tracts_affected) + 
  tm_polygons() +
  tm_shape(no_highways) + 
  tm_polygons(col = "black", alpha = 0.3) +  
  tm_shape(blackout_buildings) + 
  tm_dots(col = "coral", size = 0.1, title = "Impacted Homes") + 
  tm_layout(main.title = "Census Blocks Impacted by Blackout")

census_map

Here is a plot comparing the distributions of median household income for census tracts that did and did not experience blackouts

Reflection

As these maps show, the winter storms of 2021 in Houston Texas resulted in blackouts across the city. These blackouts effected different census blocks, with lower income census blocks being impacted more heavily.