-
Notifications
You must be signed in to change notification settings - Fork 0
/
hdem.Rmd
350 lines (271 loc) · 19.9 KB
/
hdem.Rmd
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
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
---
title: "Watershed analysis of the ORMGP region"
knit: (function(input_file, encoding) {
out_dir <- 'docs';
rmarkdown::render(input_file,
encoding=encoding,
output_file=file.path(dirname(input_file), out_dir, 'index.html'))})
author: "Oak Ridges Moraine Groundwater Program"
date: "13/12/2021"
output:
html_document:
css: "./styles.css"
toc: true
toc_float: true
toc_depth: 2
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
<!-- <br></br><br></br> -->
<!-- __*<mark>PLEASE NOTE THAT THE FOLLOWING PAGE IS IN DRAFT</mark>*__ -->
<!-- <br></br><br></br> -->
# Map
Presented up front is the __*final*__ product for which the rest of this document describes. It is a map of the ORMGP jurisdiction discretized into ~10km² sub-watersheds. *Clicking* at any sub-watershed will return a number of properties. In the figure below, sub-watersheds are colour-coded according to their degree of impervious cover.
```{r, echo=FALSE, message=FALSE, warning=FALSE, out.height='600px', out.width='100%', fig.cap="ORMGP v.2020 10km² sub-watershed map. Click a subwatershed to view spatial properties."}
library(leaflet)
library(rgdal)
geojson <- readOGR("shp/owrc20-50a_SWS10-final.geojson",verbose = FALSE)
geojson$iperc <- as.numeric(factor(geojson$perm))
leaflet(geojson) %>%
addProviderTiles( providers$Stamen.TonerLite, options = providerTileOptions(noWrap = TRUE) ) %>%
addPolygons(color = ~colorQuantile("YlGnBu", geojson$perimp)(perimp),
popup = ~paste0('<b>sub-watershed: ',mmID,'</b>',
'<br> area: ',round(Area/1000000,1),'km²',
'<br> permeability: ',perm,
'<br> impervious cover: ',round(perimp*100,0),'%',
'<br> canopy cover: ',round(percov*100,0),'%',
'<br> open water cover: ',round(perow*100,0),'%',
'<br> wetland cover: ',round(perwl*100,0),'%'
),
group = "Imperviousness"
) %>%
addPolygons(color = ~colorQuantile("YlGn", geojson$percov)(percov),
popup = ~paste0('<b>sub-watershed: ',mmID,'</b>',
'<br> area: ',round(Area/1000000,1),'km²',
'<br> permeability: ',perm,
'<br> impervious cover: ',round(perimp*100,0),'%',
'<br> canopy cover: ',round(percov*100,0),'%',
'<br> open water cover: ',round(perow*100,0),'%',
'<br> wetland cover: ',round(perwl*100,0),'%'
),
group = "Canopy coverage"
) %>%
addPolygons(color = ~colorQuantile("Blues", geojson$perow)(perow),
popup = ~paste0('<b>sub-watershed: ',mmID,'</b>',
'<br> area: ',round(Area/1000000,1),'km²',
'<br> permeability: ',perm,
'<br> impervious cover: ',round(perimp*100,0),'%',
'<br> canopy cover: ',round(percov*100,0),'%',
'<br> open water cover: ',round(perow*100,0),'%',
'<br> wetland cover: ',round(perwl*100,0),'%'
),
group = "Open water coverage"
) %>%
addPolygons(color = ~colorQuantile("GnBu", geojson$perwl)(perwl),
popup = ~paste0('<b>sub-watershed: ',mmID,'</b>',
'<br> area: ',round(Area/1000000,1),'km²',
'<br> permeability: ',perm,
'<br> impervious cover: ',round(perimp*100,0),'%',
'<br> canopy cover: ',round(percov*100,0),'%',
'<br> open water cover: ',round(perow*100,0),'%',
'<br> wetland cover: ',round(perwl*100,0),'%'
),
group = "Wetland coverage"
) %>%
addPolygons(color = ~colorFactor("Dark2", geojson$iperc)(iperc),
popup = ~paste0('<b>sub-watershed: ',mmID,'</b>',
'<br> area: ',round(Area/1000000,1),'km²',
'<br> permeability: ',perm,
'<br> impervious cover: ',round(perimp*100,0),'%',
'<br> canopy cover: ',round(percov*100,0),'%',
'<br> open water cover: ',round(perow*100,0),'%',
'<br> wetland cover: ',round(perwl*100,0),'%'
),
group = "Permeability"
) %>%
addLayersControl(
baseGroups = c("Imperviousness", "Canopy coverage","Wetland coverage","Open water coverage","Permeability"),
options = layersControlOptions(collapsed = FALSE)
)
```
# Introduction
The 3 million hectare ORMGP jurisdiction is subdivided to a number of 10km² sub-watersheds as a basis for hydrometeorological data analysis. Every sub-watershed has a defined topological order in which headwater sub-watersheds can easily be mapped to subsequent downstream sub-watersheds, and so on until feeding the great lakes. The intent here is to deem these sub-watersheds a "logical unit" for climatological and water budget analyses. Below is a description of the derivation the v.2020 OWRC 10km² sub-watershed map and its derivatives including:
- a catchment area delineation tool
- an interpolated real-time daily meteorological dataset dating back to the year 1900
- an overlay analysis performed to characterize sub-watershed:
+ impervious area
+ canopy coverage
+ water body coverage
+ wetland coverage
+ relative permeability/infiltrability
+ mean slope and dominant aspect
+ mean depth to water table.
Below first describes the processing of a Provincial Digital Elevation Model (DEM) yielding a "hydrologically-correct" model of the ORMGP ground surface. From this information, the ORMGP region is portioned into 2,813 ~10km² sub-watersheds.
Next, the hydrologically-corrected digital elevation model (HDEM) is further process to derive a number of metrics aggregated at the sub-watershed scale.
# Digital Elevation Model
Ground surface elevations were collected from the 2006 version of the Provincial (Ministry of Natural Resources) Digital Elevation Model (DEM). This is a 10x10m² DEM derived from 1:10,000 scale OBM base mapping based on photogrammetic elevation points and 5m contours where the photogrammetic elevation points did not exist. An up-scaled 50x50m² Digital Elevation Model (DEM) is produced by merging the tiles shown below:
```{r fig.asp=0.50, fig.align="center", fig.width=12, message=FALSE, warning=FALSE, echo=FALSE, fig.cap="MNR 2006 Provincial DEM tiles shown in green."}
library(ggplot2)
library(geojsonio)
library(sf)
spGL <- geojson_read("shp/greatLakes.geojson", what="sp")
spMNR <- geojson_read("shp/MNR_2006-tiles.geojson", what="sp")
spORMGP <- geojson_read("shp/ORMGP-region.geojson", what="sp")
spYPDTH <- geojson_read("shp/YPDTH.geojson", what="sp")
spMNR_points <- as.data.frame(do.call(rbind, lapply(spMNR@polygons, function(x) x@labpt))) # Convert nested list to data frame by column
names(spMNR_points) <- c('long','lat')
spMNR_points$num = spMNR$NUMBER
ggplot() +
geom_polygon(data=spGL, aes(x=long, y=lat, group=group), fill="#51a8ff", color="#51a8ff") +
geom_polygon(data=spYPDTH, aes(x=long, y=lat, group=group), fill=NA, color="brown", size=1) +
geom_polygon(data=spORMGP, aes(x=long, y=lat, group=group), fill=NA, color="orange4", size=1) +
geom_polygon(data=spMNR, aes(x=long, y=lat, group=group), fill=NA, color="darkgreen",size=1) +
geom_text(data=spMNR_points, aes(x=long, y=lat, label=num), color="darkgreen") +
theme_bw() + theme(axis.title = element_blank(),
title = element_blank(),
legend.box.spacing = unit(0, "mm")) +
# plot.background=element_rect(fill="red")) +
#labs(x='longitude',y='latitude', title="MNR 2006 Provincial DEM tiles shown in green.") +
coord_sf(xlim=c(-81, -77), ylim=c(43, 45.75), expand=FALSE)
```
Elevation data were up-scaled by taking the average of known elevations occurring within every 50x50m² cell. The resulting grid is a 5000x5000 50m-uniform cell grid, with an upper-left-corner origin at (E: 545,000, N: 5,035,000) NAD83 UTM zone 17.
## Hydrological "correction" {#hyd.corr}
An automated topological topological analysis is performed on the DEM using the following methodologies:
1. Automated depression filling (Wang and Liu, 2006) was applied to the DEM. This filtering of elevation data ensures that every grid cell has at least 1 neighbouring cell with an assigned elevation at or below the current cell's elevation. This ensures that "drainage" is never impeded.
2. While the above code works for most of the area, it will leave flat (zero-gradient) regions especially around lakes and wetlands. A fix by Garbrecht and Martz (1997) was added to ensure a consistent flow direction with negligible change to the corrected DEM.
3. Flow paths are then computed based on the "D8" algorithm (O'Callaghan and Mark, 1984).
Cell slopes and aspects are computed using a 9-point planar regression from the cell's elevation plus it's 8 neighbouring elevations.
## Manual adjustments
While automated hydrological correction is quite powerful when applied to the Provincial DEM, there are in rare places where the algorithm fails to capture mapped flow paths (usually in flatter rural regions close to embanked roads). Fortunately, these errors can be easily corrected by imposing flow directions using hand-drawn flow paths. Flow paths are save as polylines, where its vertices are ordered according to flow direction.
With the current version (v.2020), 10 flow corrections have been imposed and are saved in a set of shapefiles. There is at the moment 1 new flow path correction in queue and will be imposed for the next release. *This is to say that this layer and its derivatives are continually being updated.*
## Validation {#fig.hdem}
<center>
The ORMGP regional 50x50m² Hydrologically corrected Digital Elevation Model (HDEM, v.2020)
</center>
\
The benefit of computing the flow direction topology associated with the HDEM is that any user a could place a "virtual particle" on the landscape and follow its drainage path as it traverses toward its terminus; in this case, the Great Lakes.
More importantly however, given an HDEM, one can also efficiently compute the contributing area to any selected point on the landscape. (A live demonstration of this [is presented below](#cat.delin).) This feature provides the means to validate the processing presented in this document.
Locations of 96 current and historical stream flow gauges are used to compare drainage areas reported to the drainage areas computed using the HDEM. Sources of the station locations were collected from the Water Survey of Canada and the Toronto Region Conservation Authority. The table below lists the stations used, followed by a figure comparing their match to the HDEM. As mentioned earlier, this layer will continually be updated with input from users and it is expected that this validation can only improve.
<br></br>
<center>List of streamflow gauging stations with reported contributing area.
```{r message=FALSE, warning=FALSE, echo=FALSE}
library(knitr)
library(dplyr)
library(kableExtra)
library(rgdal)
library(ggplot2)
options(knitr.kable.NA = '')
df <- read.csv('shp/t2-check.csv')
sputm <- SpatialPoints(cbind(df$E,df$N), proj4string=CRS("+proj=utm +zone=17 +datum=WGS84"))
spgeo <- spTransform(sputm, CRS("+proj=longlat +datum=WGS84"))
lnlt <- coordinates(spgeo) # https://stackoverflow.com/questions/30018098/how-to-convert-utm-coordinates-to-lat-and-long-in-r
df %>%
mutate(DAdb=DAdb/1000/1000,DAcalc=DAcalc/1000/1000,Long=lnlt[,1],Lat=lnlt[,2]) %>%
select(2,6:9) %>%
relocate(LOC_NAME,Long,Lat,DAdb,DAcalc) %>%
mutate(diff=(DAcalc-DAdb)/DAdb*100) %>%
arrange(LOC_NAME) %>%
kbl(
col.names = c("Station ID", "Long", "Lat", "reported (km²)", "computed (km²)", "*difference* (%)"),
align = c("l","r","r","r","r","r"),
digits = 1
) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed")) %>%
add_header_above(c(" " = 3, "Contributing Area" = 3), align = "c") %>%
kableExtra::scroll_box(width = "90%", height = "300px")
df %>% ggplot(aes(DAdb/1000/1000,DAcalc/1000/1000)) +
geom_point() +
geom_abline(linetype="dotted") +
labs(title="\n\nGraphical comparison",x="reported contributing area (km²)",y="computed contributing area (km²)") +
coord_fixed()
```
</center>
# Derivative products
An overlay analysis is the process of overlaying 2 or more spatial layers and capturing statistics associated with their relative coverage. In this case, the sub-watershed layer is overlain by Provincial land-use and surficial geology layers to obtain information like percent impervious, relative permeability, etc.
Provincial layers discussed in more detail below have in all cases been re-sampled to the 50x50m² grid associated with the [hydrologically corrected DEM](#hyd.corr). It is from these rasters where the aggregation of watershed characteristics is computed.
## Land use
The Ministry of Natural Resources and Forestry (2019) SOLRIS version 3.0 provincial land use layer is employed to aggregate imperviousness and canopy coverage at the sub-watershed scale. In areas to the north, where the SOLRIS coverage discontinues, interpretation was applied by:
1. Using Provincial mapping of roads, wetlands and water bodies, areas outside of the SOLRIS data bounds *(typically up on the Canadian Shield)* are filled in with the appropriate SOLRIS land use class index (201, 150, 170, respectively--MNRF, 2019); and,
2. All remaining area not covered by SOLRIS is assumed Forest (SOLRIS land use class index of 90), as observed with satellite imagery.
The dominant SOLRIS land use class (by area) is assigned the Land use class index for every 50x50m² grid cell.
<center>
Final 50x50m SOLRIS mapping. *(For illustrative purposes only [see here](https://github.com/OWRC/subwatershed/blob/main/jupyter/OWRC-SWS.ipynb) to reproduce shown raster.)*
</center>
\
### Land use coverage
For any ~10km² sub-watershed and give a 50x50m² grid , there should be a set of roughly 4,000 SOLRIS land use class indices. Using a look-up system, the set of cells contained within a sub-watershed are assigned a value of imperviousness, water body, wetland and canopy coverage (according to their SOLRIS index) and accumulated to a sub-watershed sum.
<br></br>
<center>Percent impervious and canopy coverage as per SOLRIS v3.0 (MNRF, 2019) land use classification.
```{r message=FALSE, warning=FALSE, echo=FALSE}
library(knitr)
library(dplyr)
library(kableExtra)
options(knitr.kable.NA = '')
read.csv('shp/lookup_200731.csv') %>%
select(1:2,7:8) %>%
mutate(PerImp=PerImp*100,PerCov=PerCov*100) %>%
kbl(
col.names = c("Index", "Name", "Imperviousness (%)", "Canopy cover (%)"),
align = c("r","l","c","c"),
digits = 0
) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed")) %>%
add_header_above(c("SOLRIS Land use classification" = 2, " " = 2), align = "l") %>%
kableExtra::scroll_box(width = "90%", height = "300px")
```
</center>
<br></br>
<center>
*(For illustrative purposes only [see here](https://github.com/OWRC/subwatershed/blob/main/jupyter/OWRC-SWS.ipynb) to reproduce shown raster.)*
</center>
\
<center>
Final 50x50m canopy mapping. *(For illustrative purposes only [see here](https://github.com/OWRC/subwatershed/blob/main/jupyter/OWRC-SWS.ipynb) to reproduce shown raster.)*
</center>
\
## Surficial geology
The Ontario Geological Survey's 2010 Surficial geology of southern Ontario layer also assigns a 50x50m² grid by the dominant class.
<center>
Final 50x50m permeability mapping. *(For illustrative purposes only [see here](https://github.com/OWRC/subwatershed/blob/main/jupyter/OWRC-SWS.ipynb) to reproduce shown raster.)*
</center>
\
### Permeability
The OGS classes have been grouped according to the attribute "permeability" using a similar look-up table cross-referencing scheme. OGS (2010) adds: *"Permeability classification is a very generalized one, based purely on characteristics of material types."*
After assigning an assumed "effective" hydraulic conductivity to every permeability group, sub-watershed "permeability" is then calculated as the geometric mean of 50x50m² grid cells contained within a sub-watershed. Effective hydraulic conductivity value assumed for every permeability group is shown here:
<br></br>
<center>Permeability classifications (after OGS, 2010) and assumed effective hydraulic conductivities.
```{r message=FALSE, warning=FALSE, echo=FALSE}
library(knitr)
library(dplyr)
par <- c("Low","Low-medium","Medium","Medium-high","high","unknown/variable","fluvial","organics")
val <- format(c(1e-9,1e-8,1e-7,1e-6,1e-5,1e-8,1e-5,1e-6),digits=3)
data.frame(par,val) %>%
kbl(
col.names = c(" ", "K (m/s)"),
align = c("l","r")
) %>%
kable_styling(full_width = F, bootstrap_options = c("striped", "hover", "condensed"))
```
</center>
The resulting effective hydraulic conductivity is then reverted back to the nearest Low--High OGS (2010) classification.
## Contributing area delineation {#cat.delin}
One (of many) APIs (application programming interfaces) hosted by the ORMGP leverages the drainage topology computed in the area. Now, users have the ability to have returned a delineated catchment area polygon for any given point that lies within the HDEM extent. Try below:
```{r, echo=FALSE, message=FALSE, warning=FALSE, out.height='600px', out.width='100%', fig.cap="ORMGP v.2020 HDEM. Click anywhere (within our jurisdiction) to return its contributing area. *Hint: best to click along a watercourse*"}
knitr::include_app("https://owrc.shinyapps.io/CAdemo/", height = "600px")
```
# Source code
Processing discussed above has been documented using a [jupyter notebook](https://github.com/OWRC/subwatershed/blob/main/jupyter/OWRC-SWS.ipynb). Source data can be found [here](https://www.dropbox.com/sh/rkdu5bwn1xhm7mh/AAA1bdaplpZAZIYl0DE1e39Oa?dl=0) and additional outputs can be found [here](https://github.com/OWRC/subwatershed/tree/main/jupyter/output).
## Final shapefile
[owrc20-50a_SWS10-final.geojson (v.2020)](https://raw.githubusercontent.com/OWRC/subwatershed/main/shp/owrc20-50a_SWS10-final.geojson)
# References
Garbrecht Martz 1997 The assignment of drainage direction over flat surfaces in raster digital elevation models
O'Callaghan, J.F., and D.M. Mark, 1984. The extraction of drainage net-works from digital elevation data, Comput. Vision Graphics Image Process., 28, pp. 328-344
Ontario Geological Survey 2010. Surficial geology of southern Ontario; Ontario Geological Survey, Miscellaneous Release— Data 128 – Revised.
Ministry of Natural Resources and Forestry, 2019. Southern Ontario Land Resource Information System (SOLRIS) Version 3.0: Data Specifications. Science and Research Branch, April 2019
Wang, L., H. Liu, 2006. An efficient method for identifying and filling surface depressions in digital elevation models for hydrologic analysis and modelling. International Journal of Geographical Information Science 20(2): 193-213.