Learning objectives
- Load vector data
- Load raster data
- Simple raster operations
- Simple geoprocessing / vector operations
- Spatial joins
- Integrating SQL into an RMD
Introduction
During February 2021, Texas experienced three extreme winter storms on February 10-11, February 13-17, and Februrary 15-20 that caused power outages in Houston’s residential areas. We are interested in quantifying the severity of those power outages by calculating the number of Houston homes that lost power due to the first two storms.
To do so, we are used NASA’s VNP46A1 data product from the Visible Infrared Imaging Radiometer Suite (VIIRS) to examine night lights image data from February 7-16. We used four VIIRS data tiles to study the Houston area.
We also incorporated OpenStreetMap roads and building data to account for light emitted from traffic and buildings.
We defined an area as having experienced a blackout if we calculate a difference in night lights intensity greater than 200 nW cm-2 sr-1.
Load necessary packages
Load and explore data
Define a function to load the DNB dataset from VNP46A1 granules
show
read_dnb <- function(file_name) {
dataset_name <- "//HDFEOS/GRIDS/VNP_Grid_DNB/Data_Fields/DNB_At_Sensor_Radiance_500m"
h_string <- gdal_metadata(file_name)[199]
v_string <- gdal_metadata(file_name)[219]
tile_h <- as.integer(str_split(h_string, "=", simplify = TRUE)[[2]])
tile_v <- as.integer(str_split(v_string, "=", simplify = TRUE)[[2]])
west <- (10 * tile_h) - 180
north <- 90 - (10 * tile_v)
east <- west + 10
south <- north - 10
delta <- 10 / 2400
dnb <- read_stars(file_name, sub = dataset_name, quiet = TRUE)
st_crs(dnb) <- st_crs(4326)
st_dimensions(dnb)$x$delta <- delta
st_dimensions(dnb)$x$offset <- west
st_dimensions(dnb)$y$delta <- -delta
st_dimensions(dnb)$y$offset <- north
return(dnb)
}
Invoke the function to read in all 4 datasets
show
file_name <- "data/VNP46A1.A2021038.h08v05.001.2021039064328.h5"
dnb_38_v05 <- read_dnb(file_name = file_name)
file_name <- "data/VNP46A1.A2021038.h08v06.001.2021039064329.h5"
dnb_38_v06 <- read_dnb(file_name = file_name)
file_name <- "data/VNP46A1.A2021047.h08v05.001.2021048091106.h5"
dnb_47_v05 <- read_dnb(file_name = file_name)
file_name <- "data/VNP46A1.A2021047.h08v06.001.2021048091105.h5"
dnb_47_v06 <- read_dnb(file_name = file_name)
Combine the data into 2 datasets
dnb_38is data from before the stormdnb_47is data from after the storm
Create exploratory plots for night lights intensity before and after the storm
Note: This is a crude data analysis of the “DNB_At_Sensor_Radiance_500m” dataset, which contains values as seen from space. There is no correction for stray light, moonshine, cloud presence, or any other confounding variables As a result, it is possible we have surprising outliers and regions that systemically actually appear to be brighter from space after the storm.
show
downsample set to c(9,9)
Figure 1: Night lights intensity in Houston, TX before the storm.
show
downsample set to c(9,9)
Figure 2: Night lights intensity in Houston, TX after the storm.
Create a blackout mask
- To find the difference, we subtract night lights intensity after the storm (which is less bright) from night lights intensity before the storm (which is brighter)
- Next, we mark all observations that did not experience a blackout as NA
show
difference <- (dnb_38 - dnb_47) > 200
difference[difference == FALSE] <- NA
Vectorize the mask and fix invalid geometries
show
diff_vector <- st_as_sf(difference)
diff_vector <- st_make_valid(diff_vector)
Crop the vectorized map to our ROI
- Create a bounding box for the Houston area
- Convert the bounding box to a
sf - Project the bounding box to WGS 84
- Crop the vectorized map to the Houston bounding box
- Project the cropped, vectorized map to NAD 83
show
houston <- 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))))
houston <- st_sfc(houston)
houston <- houston %>%
st_set_crs(4326)
diff_cropped <- diff_vector[houston, op = st_intersects]
diff_cropped <- diff_cropped %>%
st_transform(3083)
Read in roads Data
- Use SQL to read in the highways data
- Project
highwaysto NAD 83 - Create and dissolve 200m buffers around the highways
show
query <- "SELECT * FROM gis_osm_roads_free_1 WHERE fclass='motorway'"
highways <- st_read("data/gis_osm_roads_free_1.gpkg", query = query, quiet = TRUE)
highways <- highways %>%
st_transform(3083)
hwy_buff_200m <- st_buffer(highways, dist = 200)
hwy_buff_200m <- st_union(hwy_buff_200m)
Read in buildings data
- Use SQL to read in the buildings data
- Project
buildingsto NAD 83 - Note:
buildingscontains residential, apartments, house, static_caravan, and detached types
show
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')"
buildings <- st_read("data/gis_osm_buildings_a_free_1.gpkg", query = query, quiet = TRUE)
buildings <- buildings %>%
st_transform(3083)
Read in census tract data
- Read in both geometries and income data
- Extract median income field and GEOID field from
acs_income - Join the datasets
- Project the joined dataset to NAD 83
- Project Houston bounding box to NAD 83
- Crop the joined data to the Houston bounding box
show
acs_geoms <- st_read("data/ACS_2019_5YR_TRACT_48_TEXAS.gdb",
layer = "ACS_2019_5YR_TRACT_48_TEXAS",
quiet = TRUE)
acs_income <- st_read("data/ACS_2019_5YR_TRACT_48_TEXAS.gdb",
layer = "X19_INCOME",
quiet = TRUE)
median_income <- acs_income[c("GEOID", "B19013e1")]
geo_id <- acs_income[c("GEOID")]
geoms_income <- left_join(acs_geoms, acs_income, by = c("GEOID_Data" = "GEOID"))
geoms_income <- geoms_income %>%
st_transform(3083)
houston <- houston %>%
st_transform(3083)
geoms_income_cropped <- geoms_income[houston, op = st_intersects]
Merge the datasets
- Remove highway 200m buffer from the vectorized blackout mask of Houston so our area of interest ignores highways
show
no_hwy <- st_difference(diff_cropped, hwy_buff_200m)
Question 1: How many residential buildings were without power on 2021-02-16?
- Use spatial subsetting to find all residential buildings in the blackout areas
- Crop all residential buildings based on
no_hwy(which is our area of interest) - Extract the number of rows (aka the number of buildings that experienced a blackout )
show
buildings_blackout <- buildings[no_hwy, op = st_intersects]
q1_ans <- nrow(buildings_blackout)
There were 157411 residential buildings without power on 2021-02-16.
Visualize buildings where the blackout occurred
- Created a map of residential buildings that experienced a blackout (
res) - Added
reson top of a base map of the Houston area
show
res <- tm_shape(buildings_blackout) +
tm_borders()
pt1 <- st_point(c(-96.5, 29))
pt2 <- st_point(c(-96.5, 30.5))
pt3 <- st_point(c(-94.5, 30.5))
pt4 <- st_point(c(-94.5, 29))
coords <- list(rbind(pt1, pt2, pt3, pt4, pt1))
polygon <- st_polygon(x = coords)
houston <- st_sfc(polygon, crs = "EPSG:4326")
houston_bbox <- st_bbox(houston)
houston_map <- osm.raster(houston_bbox)
tm_shape(houston_map) +
tm_rgb() +
res
Figure 3: Approximately 157,411 residential buildings in Houston, TX were affected by the February 2021 storm.
Conclusion
In conclusion, we identified 157411 residential buildings in Houston without power due to the severe February 2021 storms. While night lights intensity may just be one way to measure the impact of these events, this method illustrates the vast potential of remotely sensed data from sources such as the Visible Infrared Imaging Radiometer Suite (VIIRS) aboard the Suomi satellite.
GitHub
Full code for this assignment can be found here.