5.2 Exposure data
The next step in the Stanford Urban Risk Framework is to identify the intersection between the hazard and physical assets or communities which we care about (since a natural hazard happening in the middle of nowhere is of little concern to us). Many of the techniques you have developed for reading population data and working with it geospatially will directly apply here. For example, the wildfire information we collected in the last section can be combined with a variety of census data, using the spatial subsetting approach from the intro textbook, to produce approximations of direct community exposure to wildfires, and the extreme heat data we retrieved is already in the form of CBG data, so an equity analysis can immediately be performed. In this section, we’ll demonstrate one new technique for combining raster data, as we just experienced with flood maps, with polygons, such as building footprints, to estimate physical exposure at a granular scale.
We’ll use building data as we’ve used from Section 2.4, focused in East Palo Alto. Recall that we don’t have publicly accessible Assessor-Recorder data from San Mateo County to work with, which could be used to assign “value” to buildings based on property value, but alternatively we can think of economic value as the “replacement costs” for some square feet of damaged building, given that the costs to replace carpet or drywall are more related to today’s material and labor costs than they are to whatever the value of your property was when you bought it in the past. That being said, knowing building use type (which would be available from the Assessor-Recorder) would be useful, since construction costs can vary significantly depending on whether the building is residential, commercial, or industrial use. East Palo Alto does have an online zoning map, where zoning does not definitively tell you what the current use of a building is, but gives you a better sense than no information at all (consider scraping this as practice, using techniques from the last section). Of course, as we think about long-term flood risk, we also shouldn’t assume that a building stays the same forever, so both building footprint and land use can fundamentally change. For our demonstration, we’ll focus on just the direct question of how much physical exposure occurs for individual buildings, but keep in mind the many other layers of information that could be attached to buildings and parcels, as we explored in Chapter 2.
library(tidyverse)
library(sf)
library(leaflet)
library(mapboxapi)
library(raster)
library(stars)osm_bldg <- st_read("gis_osm_buildings_a_free_1.shp")
epa_boundary <-
places("CA", cb = T, progress_bar = F) %>%
filter(NAME == "East Palo Alto") %>%
st_transform(4326)
epa_bldg <- osm_bldg[epa_boundary, ]Since each flood map is additive, we can use the most intense flood to filter to just the buildings that are directly exposed to any of our floods.
flood_max <-
raster("flood/SLR050_RP100_epa_flood.tif")
flood_max_extent <-
flood_max %>%
st_as_stars() %>%
mutate(SLR050_RP100_epa_flood = ifelse(
!is.na(SLR050_RP100_epa_flood),
1,
NA
)) %>%
st_as_sf(merge = T) %>%
st_set_crs(26910) %>%
st_make_valid() %>%
st_transform(4326)
epa_bldg_flooded_max <-
epa_bldg %>%
st_transform(4326) %>%
.[flood_max_extent,]The stars package is particularly nice because it allows you to use many dplyr functions, and also provides the fastest method for converting a raster to its polygon extent, specifically st_as_sf(merge = T). Note that sometimes the resultant polygons from a raster conversion can be particularly prone to invalid properties; if st_make_valid() doesn’t resolve the issue, sometimes st_buffer(0) does the trick.
flood_pal <- colorNumeric(
palette = "Blues",
domain = values(flood_max),
na.color = "transparent"
)
leaflet() %>%
addMapboxTiles(
style_id = "satellite-streets-v11",
username = "mapbox",
options = tileOptions(opacity = 0.5)
) %>%
addRasterImage(
flood_max,
colors = flood_pal,
opacity = 0.75
) %>%
addPolygons(
data = epa_bldg_flooded_max,
fill = F,
color = "red",
weight = 0.5
) %>%
addLegend(
pal = flood_pal,
values = values(flood_max),
title = "Flood depth, cm"
)Using this maximum set of buildings, we can now loop through our 9 flood maps and perform an intersection between raster and each individual building polygon called “zonal statistics”, where the information contained in the raster is summarized in much the same way as group_by() and summarize() works (meaning you can pick to do mean or max, etc.), but the group is the pixels within the polygon boundaries (partial pixels are handled proportionally, much like spatial subsetting from the intro textbook). The function to do this is called extract() from raster.
epa_bldg_exposure <- NULL
for(slr in c("000","025","050")){
for(rp in c("001","020","100")){
print(paste0("SLR",slr,"_RP",rp))
flood <- raster( paste0("flood/SLR",slr,"_RP",rp,"_epa_flood.tif"))
flood_extent <-
(flood > -Inf) %>%
st_as_stars() %>%
st_as_sf(merge = T) %>%
st_set_crs(26910) %>%
st_make_valid() %>%
st_transform(4326)
epa_bldg_flooded <-
epa_bldg_flooded_max[flood_extent,] %>%
st_transform(26910)
flood_crop <-
crop(flood, epa_bldg_flooded)
flood_crop[is.na(flood_crop)] <- 0
temp <-
extract(
flood_crop,
epa_bldg_flooded,
fun = mean
) %>%
as.data.frame() %>%
rename(avg_depth = V1) %>%
cbind(
epa_bldg_flooded %>%
st_drop_geometry() %>%
dplyr::select(osm_id)
) %>%
mutate(
SLR = slr,
RP = rp
)
epa_bldg_exposure <-
epa_bldg_exposure %>%
rbind(temp)
}
}
saveRDS(epa_bldg_exposure,"epa_bldg_exposure.rds")Notes:
(flood > -Inf)is a shortcut of a step in a previous chunk where we converted all non-NA values in the raster to a single value to streamline the conversion to a polygon extent. The result of this logical condition is effectively the same, where this condition gives you TRUE for all real numbers.- It is not necessary to create a subset of buildings flooded in just the given scenario, or to crop the flood raster to just the extent of those buildings, but these steps can speed up the
extract()process, especially when there’s a smaller number of buildings.
The end result is a tidy dataframe called epa_bldg_exposure, which holds average depths for each building in each of 15 scenarios. We can now take this exposure information and calculate actual damages, which depends on vulnerability.