Lab 6

Terrain Analysis

#1 DATA INGEST

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

basin_elevation = elevatr::get_elev_raster(boundary, z = 13) %>%
  crop(boundary) %>% 
  mask(boundary) 

basin_in_feet = basin_elevation *3.281

writeRaster(basin_in_feet, filename = "/Users/lfinn443/github/geog-176A-labs/data/basin_elevation.tif", overwrite = T)

#Buildings and river-network data

building = osmdata::add_osm_feature(opq(boundary), "building") %>% 
  osmdata_sf()

stream = osmdata::add_osm_feature(opq(boundary), "waterway", "stream") %>% osmdata_sf()
building = st_intersection(building$osm_polygons, boundary)

river = st_intersection(stream$osm_lines, boundary)

#TERRAIN ANALYSIS

#Hillshade

wbt_hillshade("/Users/lfinn443/github/geog-176A-labs/data/basin_elevation.tif", "/Users/lfinn443/github/geog-176A-labs/data/hillshade.tif")

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

hillshade = raster("/Users/lfinn443/github/geog-176A-labs/data/hillshade.tif")

#Plots

plot(hillshade, col = gray.colors(256, alpha = .5), legend = F)
plot(river$geometry, col = "blue3", add = T)
plot(boundary$geometry, add = T)

#QUESTION 2: HEIGHT ABOVE NEAREST DRAINAGE

#Creating the River Raster

river_raster = st_transform(river, 5070) %>% 
  st_buffer(10) %>% 
  st_transform(crs(basin_in_feet)) %>% 
  fasterize::fasterize(basin_in_feet) %>% 
  writeRaster("/Users/lfinn443/github/geog-176A-labs/data/river_raster.tif", overwrite = T)

#Creating the hydrologically corrected surface

wbt_breach_depressions("/Users/lfinn443/github/geog-176A-labs/data/basin_elevation.tif",
                       "/Users/lfinn443/github/geog-176A-labs//data/breach_depression.tif")

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

#Creating the Hand Raster

wbt_elevation_above_stream("/Users/lfinn443/github/geog-176A-labs/data/breach_depression.tif", "/Users/lfinn443/github/geog-176A-labs/data/river_raster.tif", "/Users/lfinn443/github/geog-176A-labs/data/ft_above_str.tif")

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

river_raster = raster("/Users/lfinn443/github/geog-176A-labs/data/river_raster.tif")
hand_raster = raster("/Users/lfinn443/github/geog-176A-labs/data/ft_above_str.tif") +3.69

hand_raster[river_raster ==1] = 0

writeRaster(hand_raster, "/Users/lfinn443/github/geog-176A-labs/data/correct_to_lrd.tif", overwrite = T)