-
Notifications
You must be signed in to change notification settings - Fork 0
/
03_gis.qmd
248 lines (195 loc) · 8.78 KB
/
03_gis.qmd
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
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
# GIS in Hydrology
Spatial Data Science\
Geocomputation with R
## WhiteboxTools
In this session we use the **WhiteboxTools** (WBT) a modern and advanced geospatial
package, tools collection, which contains \~450 functions. The WBT has an interface
to r and python.
The tool can be downloaded from [http://whiteboxgeo.com/](http://whiteboxgeo.com/).
## Raster and vector data
**Raster** data are represented by a matrix of pixels (cells) with values. Raster
is used for data which display continuous information across an area which cannot
be easily divided into vector features. For the purpose of watershed delineation
the raster input of Digital Elevation Model is used.
## Watershed delineation
The process of delineation is the first step in basin description. One simply has
to delineate the domain
The step-by-step process involves:\
- Acquiring digital elevation model of area\
- Pit and sink removal\
- Flow accumulation calculation \
- Flow direction calculation \
- Outlet identification \
- Delineation towards specified outlet\
Let's start with `whitebox` package that contains an API to the WhiteboxTools
executable binary. \
We need to reference path to the executable
Except for the `whitebox` package, some other packages for general work with spatial
data are necessary. The packages `terra` and `sf` are needed for working with the
raster and vector data. They also provide access to **PROJ**, **GEOS** and **GDAL**
which are open source libraries that handle projections, format standards and provide
geoscientific calculations. And the package `tmap` makes plotting both raster and
vector layers very easy.
## Watershed delineation workflow
```{r, collapse=TRUE}
library(whitebox) # <1>
library(terra) # <2>
library(sf) # <3>
library(tmap) # <4>
whitebox::wbt_init(exe_path = "~/Downloads/WBT/whitebox_tools") # <5>
whitebox::check_whitebox_binary() # <6>
```
1. Load the Whitebox API package
2. Load the `terra` package for raster use
3. Load the `sf` package for vector use
4. Load the `tmap` plotting functions for layers
5. Whitebox needs the information where the data executable is stored
6. Binary check of functionality
#### Sample data
```{r}
dem <- rast(whitebox::sample_dem_data()) # <1>
tm_layout() + # <2>
tm_shape(dem) + # <2>
tm_raster(palette = "-Greys", n = 50, legend.show = FALSE) # <2>
```
1. Use the `rast()` function to load the data
2. Plot via the `tmap` workflow
#### DEM workflow
The digital elevation model has to be adjusted for the watershed delineation
algorithm to be able to run successfully.
We need to specifiy a path to working directory. It has to be somewhere where you
as a user have access to write files. Be sure the folder exists.
```{r, collapse=TRUE, eval=TRUE}
wd_path <- "~/Desktop/GIS"
```
```{r, collapse=TRUE, eval=FALSE}
wbt_fill_depressions_wang_and_liu(dem = paste(wd_path, # <1>
"dem.tif", # <1>
sep = "/"), # <1>
output = paste(wd_path, # <1>
"filled_depresions.tif", # <1>
sep = "/")) # <1>
wbt_d8_pointer(dem = paste(wd_path, # <2>
"filled_depresions.tif", # <2>
sep = "/"), # <2>
output = paste(wd_path, # <2>
"d8pointer.tif", # <2>
sep = "/")) # <2>
wbt_d8_flow_accumulation(input = paste(wd_path, # <3>
"filled_depresions.tif", # <3>
sep = "/"), # <3>
output = paste(wd_path, # <3>
"flow_accu.tif", # <3>
sep = "/"), # <3>
out_type = "cells") # <3>
wbt_extract_streams(flow_accum = paste(wd_path, # <4>
"flow_accu.tif", # <4>
sep = "/"), # <4>
output = paste(wd_path, # <4>
"streams.tif", # <4>
sep = "/"), # <4>
threshold = 200) # <4>
```
1. This algorithm involves removing flat areas and filling depressions, thus producing hydrologically corrected DEM
2. Pointer is 
## Gauge
We have the raster prepared for the delineation, now we need to provide a point
layer with the gauge, to which the watershed should be delineated. The point
has to be placed at the stream network. We will create the layer from scratch.
```{r, eval=TRUE}
gauge <- st_sfc(st_point(x = c(671035, 4885783),
dim = "XY"),
crs = st_crs(26918))
st_write(obj = gauge,
dsn = paste(wd_path, "gauge", sep = "/"),
driver = "ESRI Shapefile",
delete_dsn = TRUE)
```
```{r, eval=TRUE}
tm_layout() +
tm_shape(dem) +
tm_raster(palette = "-Greys", n = 50, legend.show = FALSE) +
tm_shape(gauge) +
tm_dots(col = "red", size = 1)
```
If the point is not located directly in the stream, it could cause troubles.
The **Jenson snap pour** makes sure to move the point to the nearest stream pixel.
```{r, eval=FALSE}
wbt_jenson_snap_pour_points(pour_pts = paste(wd_path,
"gauge/gauge.shp",
sep = "/"),
streams = paste(wd_path,
"streams.tif",
sep = "/"),
output = paste(wd_path,
"gauge_snapped.shp",
sep = "/"),
snap_dist = 1000)
```
Now everything is set to delineate the watershed with using the gauge and D8 pointer.
```{r, eval=FALSE}
wbt_watershed(d8_pntr = paste(wd_path, "d8pointer.tif", sep = "/"),
pour_pts = paste(wd_path, "gauge_snapped.shp", sep = "/"),
output = paste(wd_path, "watershed.tif", sep = "/"))
```
The watershed is usually used in vector format, but now it is in raster. Let's
finish with the conversion.
```{r}
wtshd <- rast(paste(wd_path, "watershed.tif", sep = "/"))
watershed <- terra::as.polygons(wtshd)
```
## River morphology
Under the term **river morphology** we understand the description of the shape of
river channels. Hydrologists use indices such as **stream length** or **Strahler order**.
```{r}
wbt_strahler_stream_order(d8_pntr = paste(wd_path, "d8pointer.tif", sep = "/"),
streams = paste(wd_path, "streams.tif", sep = "/"),
output = paste(wd_path, "strahler.tif", sep = "/"))
```
```{r, echo=FALSE, eval=FALSE}
library(ggplot2)
library(tidyterra)
results <- list.files(wd_path, pattern = "tif$", full.names = TRUE)
r <- lapply(results, rast)
r1 <- rast(results[2])
r2 <- rast(results[1])
r3 <- rast(results[4])
r4 <- rast(results[6])
r5 <- watershed
r6 <- rast(results[5])
p1 <- ggplot() +
ggtitle(label = "a) Corrected DEM") +
geom_spatraster(data = r1, show.legend = FALSE) +
scale_fill_viridis_c() +
theme_minimal()
p2 <- ggplot() +
ggtitle(label = "b) D8 Pointer") +
geom_spatraster(data = r2, show.legend = FALSE) +
scale_fill_viridis_c() +
theme_minimal()
p3 <- ggplot() +
ggtitle(label = "c) Flow Accumulation") +
geom_spatraster(data = r3, show.legend = FALSE) +
scale_fill_viridis_b(na.value = "white") +
theme_minimal()
p4 <- ggplot() +
ggtitle(label = "d) Identified Streams") +
geom_spatraster(data = r4, show.legend = FALSE) +
scale_fill_viridis_c(na.value = "white") +
theme_minimal()
p5 <- ggplot() +
ggtitle(label = "e) Watershed") +
geom_spatraster(data = r1, show.legend = FALSE) +
scale_fill_viridis_c() +
geom_spatvector(data = r5, show.legend = FALSE) +
theme_minimal()
p6 <- ggplot() +
ggtitle(label = "f) Strahler Order") +
geom_spatraster(data = r6, show.legend = FALSE) +
scale_fill_viridis_b(na.value = "white") +
theme_minimal()
cowplot::plot_grid(p1, p2, p3, p4, p5, p6, ncol = 2, greedy = TRUE)
ggsave("~/Documents/github_projects/Hydrology_ZVX108E/data/output.png")
```
## Results
