In this lab:
- Finding and bringing data into R using Web APIs
- Aligning retrieved object and field data
- Preparing terrain data for analysis
- Building a partial Flood Inundation Map Library for the Mission Creek Basin
Libraries
data wrangling tools
- library(tidyverse)
raster tools
- library(raster)
- library(fasterize)
spatial and terrian tools
- library(whitebox)
- library(sf)
data used
- library(osmdata)
- library(elevatr)
- library(gifski)
Collect Data
Create Basin boundary
basindata <- sf::read_sf("https://labs.waterdata.usgs.gov/api/nldi/linked-data/nwissite/USGS-11119750/basin")
# convert values from meters to feet (1m = 3.281 feet)
convertmtoft <- function(meters){
return (meters * 3.281)
}
# crop/mask the elevation data to basin boundary
elevation <- get_elev_raster(basindata, z = 13) %>%
crop(basindata) %>%
mask(basindata) %>%
convertmtoft()
# write raster data folder as a tif file
writeRaster(elevation, "C:/Users/Catherine/Documents/github/geog-176A-labs/data/elevation_basin.tif", overwrite = TRUE)Get elevation, buildings and river-network data
elevationraster <- raster("C:/Users/Catherine/Documents/github/geog-176A-labs/data/elevation_basin.tif")
basinbox <- st_bbox(elevationraster) %>%
st_as_sfc() %>%
st_transform(4326)
buildings <- opq(basinbox) %>%
add_osm_feature(key = 'building') %>%
osmdata_sf()
pointbuildings <- buildings$osm_points %>%
st_intersection(basindata)
buildings <- buildings$osm_polygons %>%
st_intersection(basindata)
# extract the railway POINT from the buildings for latter plotting
railwaypoints <- opq(basinbox) %>%
add_osm_feature(key = 'railway', value = 'station' ) %>%
osmdata_sf()
railwaypoints <- railwaypoints$osm_points %>%
st_intersection(basindata)
streams <- opq(basinbox) %>%
add_osm_feature(key = 'waterway', value = 'stream') %>%
osmdata_sf()
streams <- streams$osm_lines %>%
st_intersection(basindata) Terrain Analysis
Create Hillshade
# a “Hillshade” raster for visualization
wbt_hillshade("C:/Users/Catherine/Documents/github/geog-176A-labs/data/elevation_basin.tif", 'C:/Users/Catherine/Documents/github/geog-176A-labs/data/basin_hillshade.tif')## [1] "hillshade - Elapsed Time (excluding I/O): 0.45s"basinhillshade <- raster('C:/Users/Catherine/Documents/github/geog-176A-labs/data/basin_hillshade.tif')
plot(basinhillshade, col = gray.colors(256, alpha = .5), box = F, main = '', legend = F)
plot(basindata, add = T)
plot(streams, add = T , col = 'blue')
***
Create Height Above Nearest Drainage (HAND) raster
# buffer streamlines using a distance of 10 meters
buffer_streamlines <- streams %>%
st_transform(5070) %>%
st_buffer(10) %>%
st_transform(4326)
# use fasterize to create river network raster using the elevation grid as the template
rasterrivers <- fasterize::fasterize(buffer_streamlines, elevationraster)
# write river raster to data folder
writeRaster(rasterrivers, "C:/Users/Catherine/Documents/github/geog-176A-labs/data/riverraster.tif", overwrite = T)
riverraster <- raster('C:/Users/Catherine/Documents/github/geog-176A-labs/data/riverraster.tif')Creating a hydrologically corrected surface
# use the whitebox tool wbt_breach_depressions to correct errors in our DEM/elevation grid.
wbt_breach_depressions("C:/Users/Catherine/Documents/github/geog-176A-labs/data/elevation_basin.tif", 'C:/Users/Catherine/Documents/github/geog-176A-labs/data/breachelevationraster.tif')## [1] "breach_depressions - Elapsed Time (excluding I/O): 0.357s"#
wbt_elevation_above_stream('C:/Users/Catherine/Documents/github/geog-176A-labs/data/breachelevationraster.tif', "data/riverraster.tif", "C:/Users/Catherine/Documents/github/geog-176A-labs/data/HANDelevationraster.tif" )## [1] "elevation_above_stream - Reading streams data..."# a “Height Above Nearest Drainage” (HAND) raster for rapid flood assessment
HANDraster <- raster('C:/Users/Catherine/Documents/github/geog-176A-labs/data/HANDelevationraster.tif')Add 3.69 foot offset to current HAND raster and then ensure river network still has a HAND value of 0
# Add the offset to the entire HAND raster
offsetHANDraster <- HANDraster + 3.69
# if the cell has a value of 1 in river raster then its HAND value is set to 0.
index <- which(values(riverraster) == 1)
values(offsetHANDraster)[index] <- 0
# save the new, offset raster to your data folder
writeRaster(offsetHANDraster, 'C:/Users/Catherine/Documents/github/geog-176A-labs/data/offsetHANDraster.tif', overwrite = T)2017 Impact Assessment
Map the flood
# Read in the corrected HAND raster
offsetHANDraster <- raster('C:/Users/Catherine/Documents/github/geog-176A-labs/data/offsetHANDraster.tif')
# Create a flood map where any cell with a HAND value greater then 10.02 is set to NA.
values(offsetHANDraster)[which(values(offsetHANDraster) > 10.02)] <- NA
plot(basinhillshade, box = F, main = '', legend = F, col = gray.colors(256, alpha = .5))
# Overlay the flood map using the rev(blues9) color palette and add = TRUE
plot(riverraster, add = T, box = F, legend = F, col = rev(blues9))
# Plot the railway station as a point of reference colored in green, with cex = 1, and pch = 16
plot(railwaypoints, add = T, box = F, cex = 1, pch = 16, col = 'green')
Estimate and visualize the impacts
# use the OSM building centroids to extract the flood depth at each structure from your flood raster
flooddepth <- raster::extract(offsetHANDraster, pointbuildings)
# count the number not equal to NA in order to determine the number of impacted structures.
pointbuildings$flooded <- as.factor(ifelse(!is.na(flooddepth), 1, 0))
plot(basinhillshade, col = gray.colors(256, alpha = .5), box = F, legend = F,
main = paste(sum(pointbuildings$flooded==1),'Buildings Impacted - 2017 Santa Barbara Basin Flood'))
# add flood map
plot(offsetHANDraster,add = T, box = F, legend = F, col = rev(blues9))
# add the buildings centroids to the map:
plot(pointbuildings, add = T, box = F, cex = .08, pch = 16, col = c('black','red')[pointbuildings$flooded])
# Add the railroad as a reference point
plot(railwaypoints, add = T, box = F, cex = 1, pch = 16, col = 'green')
legend(x = 'topright', legend = c('Flooded', 'Not Flooded'),
col = c('red', 'black'), cex = 1, pch = 16)