COMPASS-Flooding-Stats/Goodwin watershed delineation.R at main · Smithsonian/COMPASS-Flooding-Stats · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
# from the following lesson: https://vt-hydroinformatics.github.io/rgeowatersheds.html
library(tidyverse)
library(raster)
library(sf)
library(whitebox)
library(tmap)
library(stars)
library(rayshader)
library(rgl)
library(leaflet)
# library(terra)

whitebox::wbt_init()
theme_set(theme_classic())


dem <- raster("data/GIS/GoodwinIsland/Goodwin_elevation_3m.tif")

pal <- colorNumeric(palette = "PuOr", values(dem),
                    na.color = "transparent")

# leaflet() %>% addTiles() %>%
#   addRasterImage(dem, colors = pal, opacity = 0.8) %>%
#   addLegend(pal = pal, values = values(dem),
#             title = "Elv (m)")

# wbt_hillshade(dem = "data/GIS/GCREW_elevation.tif",
#               output = "data/GIS/GCREW_hillshade.tif",
#               azimuth = 115)
#
# hillshade <- raster("data/GIS/GCREW_hillshade.tif")
#
# leaflet() %>% addTiles() %>%
#   addRasterImage(hillshade, opacity = 0.8)
#
# Fill depressions
wbt_breach_depressions_least_cost(
  dem = "data/GIS/GoodwinIsland/Goodwin_elevation_3m.tif",
  output = "data/GIS/GoodwinIsland/Goodwin_breached.tif",
  dist = 5,
  fill = TRUE)

wbt_fill_depressions_wang_and_liu(
  dem = "data/GIS/GoodwinIsland/Goodwin_breached.tif",
  output = "data/GIS/GoodwinIsland/Goodwin_breached_filled.tif"
)

wbt_d8_flow_accumulation(input = "data/GIS/GoodwinIsland/Goodwin_breached_filled.tif",
                         output = "data/GIS/GoodwinIsland/Goodwin_D8FA.tif")

wbt_d8_pointer(dem = "data/GIS/GoodwinIsland/Goodwin_breached_filled.tif",
               output = "data/GIS/GoodwinIsland/Goodwin_pointer.tif")

ppoints <- tribble(
  ~Lon, ~Lat,
  -76.404642, 37.215855
)

ppointsSP <- SpatialPoints(ppoints, proj4string = CRS("+proj=longlat +datum=WGS84"))

shapefile(ppointsSP, filename = "data/GIS/GoodwinIsland/Goodwin_flume.shp", overwrite = TRUE)

wbt_extract_streams(flow_accum = "data/GIS/GoodwinIsland/Goodwin_D8FA.tif",
                    output = "data/GIS/GoodwinIsland/Goodwin_streams.tif",
                    threshold = 100)

wbt_jenson_snap_pour_points(pour_pts = "data/GIS/GoodwinIsland/Goodwin_flume.shp",
                            streams = "data/GIS/GoodwinIsland/Goodwin_streams.tif",
                            output = "data/GIS/GoodwinIsland/Goodwin_flume_snapped.shp",
                            snap_dist = 0.0005) #careful with this! Know the units of your data

pp <- shapefile("data/GIS/GoodwinIsland/Goodwin_flume_snapped.shp")
streams <- raster("data/GIS/GoodwinIsland/Goodwin_streams.tif")

# leaflet() %>% addProviderTiles("Esri.WorldImagery", group = "ESRI") %>%
#
#   addRasterImage(dem, colors = pal, opacity = 0.8) %>%
#   addRasterImage(streams, colors = "blue", opacity = 0.8) %>%
#   addCircleMarkers(color = "black", lat = pp@coords[,"y"], lng = pp@coords[,"x"])

wbt_watershed(d8_pntr = "data/GIS/GoodwinIsland/Goodwin_pointer.tif",
              pour_pts = "data/GIS/GoodwinIsland/Goodwin_flume_snapped.shp",
              output = "data/GIS/GoodwinIsland/Goodwin_flume_watershed.tif")

lil_shed <- raster("data/GIS/GoodwinIsland/Goodwin_flume_watershed.tif")

leaflet() %>% addProviderTiles("Esri.WorldImagery", group = "ESRI") %>%
  addRasterImage(lil_shed, colors = "white", opacity = 0.8) %>%
  addRasterImage(streams, colors = "blue", opacity = 0.8) %>%
  addCircleMarkers(color = "black", lat = pp@coords[,"y"], lng = pp@coords[,"x"])
# addRasterImage(streams, colors = "blue", opacity = 0.8)