Geography 176A

library(sf)        # vector manipulation
library(raster)    # raster manipulation
library(fasterize) # "faster" raster
library(whitebox)  # terrain analysis

# Data libraries
library(osmdata)   # OSM API
library(elevatr)   # Elevation  Web Tiles

library(tidyverse)
library(fasterize) 
library(whitebox)  
library(gifski)
library(knitr)

basin = read_sf("https://labs.waterdata.usgs.gov/api/nldi/linked-data/nwissite/USGS-11119750/basin")

elev = elevatr::get_elev_raster(basin, z = 13) %>% 
  crop(basin) %>% 
  mask(basin)
elev = elev * 3.281
writeRaster(elev, "../data/basin-elevation.tif", overwrite = TRUE)

elev = raster("../data/basin-elevation.tif")

bb = st_bbox(basin) %>% 
  st_as_sfc() %>% 
  st_transform(4326)

building <- opq(bb) %>% 
  add_osm_feature(key = "building") %>%
  osmdata_sf()

building_points = building$osm_points %>% 
  st_intersection(basin) %>% 
  st_transform(crs(basin)) %>% 
  st_centroid()

building_poly = building$osm_polygons %>% 
  st_intersection(basin) %>% 
  st_transform(crs(basin))

railway = opq(bb) %>% 
  add_osm_feature(key = "railway", value = "station") %>% 
  osmdata_sf()

railway_points = railway$osm_points %>% 
  st_intersection(basin) %>% 
  st_transform(crs(basin))

stream = opq(bb) %>% 
  add_osm_feature(key = "waterway", value = "stream") %>% 
  osmdata_sf()

stream_lines = stream$osm_lines %>% 
  st_intersection(basin) %>% 
  st_transform(crs(basin))

wbt_hillshade("../data/basin-elevation.tif", '../data/basin-hillshade.tif')

## [1] "hillshade - Elapsed Time (excluding I/O): 0.36s"

hillshade = raster("../data/basin-hillshade.tif")

plot(hillshade, col = gray.colors(256, alpha = .5), box = FALSE, legend = FALSE, main = "Basin and Stream")
plot(stream_lines, col = "blue", add = TRUE)
plot(basin, add = TRUE)

stream_buffer = stream_lines %>% 
  st_transform(5070) %>% 
  st_buffer(10) %>% 
  st_transform(4326)

network_river = fasterize::fasterize(stream_buffer, elev)

writeRaster(network_river, "../data/network-river.tif", overwrite = TRUE)
wbt_breach_depressions("../data/basin-elevation.tif", "../data/breach-depressions.tif")

## [1] "breach_depressions - Elapsed Time (excluding I/O): 0.224s"

wbt_elevation_above_stream("../data/breach-depressions.tif", "../data/network-river.tif", "../data/basin-hand.tif")

## [1] "elevation_above_stream - Elapsed Time (excluding I/O): 0.79s"

basin_river = raster("../data/network-river.tif")
basin_hand = raster("../data/basin-hand.tif")

basin_offset = basin_hand + 3.69
basin_offset[basin_river == 1] = 0
writeRaster(basin_offset, "../data/basin-offset.tif", overwrite = TRUE)

basin_offset = raster("../data/basin-offset.tif")
basin_flood = basin_offset
basin_offset[basin_offset > 10.02] = NA

Geography 176A

Lab 06: Terrain Analysis

Yan Wang