5 Reach and catchment environmental metrics
I have calculated many environmental metrics for all reaches of the stream network. These metrics have many potential uses for stream management and research. The remaining chapters of this manual describe the derivation of the metrics and the structure of the tables holding them.
Two primary environmental metric tables hold the most commonly used metrics:
Table cat_env contains environmental variables for every reach that describe catchment attributes.
Table subc_env contains environmental variables for every reach that describe attributes of the reach or its immediate subcatchment.
The two tables contain site and reach fields matching those fields in the subcs table. The tables are broadly separated into variables that describe the reach or its local subcatchment (subc_env) and those that describe the catchment upstream of the reach (cat_env). Metrics fall into the following classes (some of which have fields in both subc_env and cat_env.). See Figure 1.1 for the full list of fields in these two tables.
There are two almost identical versions of subcatchment (and catchment) area (in m2) in the database. area_r_m2 in the subc_env table is the sum of gridcell areas for each subcatchment in the r_site raster table. scarea in the subcs table is the area of the polygons in the table. They differ because a small number of subcatchment polygons were manually altered to correct network relationships, resulting in an imperfect match with the DEM gridcells on which the subcatchments were originally calculated. scarea is a more accurate estimate of subcatchment area, but area_r_m2 was used for calculating subcatchment statistics (e.g. impervious cover and tree cover fractions) from raster data. cat_area_m2_r in the cat_env table sums area_r_m2 for all upstream subcatchments, and carea_km2 in the subcs table similarly sums and converts to km2 scarea.
5.1 Impervious cover
Impervious area as a fraction of the upstream catchment (total imperviousness, TI) is a primary measure of urban density. Effective imperviousness (EI) is an important measure of urban stormwater drainage impacts (Walsh et al. 2022). I have calculated these metrics for 2005, 2018, and 2022 (with 2012 and 204 to follow). The fields ti_XXXX and ei_XXXX in the cat_env table record these two metrics for year XXXX. The fields tia_XXXX_m2 and eia_XXXX_m2 in the subc_env table record the total and effect area of impervious surfaces in m2 within each subcatchment. See Chapter 6 for full details.
5.2 Tree canopy cover
Tree canopy cover near streams is a strong predictor of in-stream biotic assemblages (Peterson et al. 2010; Walsh and Webb 2014), and catchment tree canopy cover is an important predictor of streamflow volumes (Zhang, Dawes, and Walker 2001). I thus calculated both tree canopy area as a fraction of stream catchment (TF), and tree canopy area weighted by flow distance to and along streams as a fraction of similarly weighted catchment area (Attenuated forest cover, AF). I have calculated these metrics for 1750 (estimated from ecological vegetation classes), 2006, 2018 and 2022 (with 2012 and 2014 to follow). The fields “tf_XXXX” and “af_XXXX” in the cat_env table, and “tfa_XXXX_m2” (total area of tree canopy cover per subcatchment) in the *subcat_env” table follow the same convention as for ti and ei. See Chapter 7 for full details.
5.3 Flow
Two estimates of mean annual runoff were calculated: post-1997, and for the period 1911-2009 for equivalence with the measure of mean discharge used in earlier studies of the region. See Chapter 8 for full details.
5.4 Flow stress: degree of regulation
A “degree of regulation” index (DoR) was calculated as a measure of flow stress. See Chapter 9 for full details.
5.5 Barriers to fish passage
Barriers to fish (and potentially other fauna) migration such as weirs, dams, and piped streams can limit the distributions of species. The number of barriers to instream migration downstream of each reach was calculated from a revised dataset of barriers included in the database. Two variants were calculated n_all_barriers_ds summing partial and full barriers, and n_full_barriers_ds summing just full barriers were added to the subc_env table. See Chapter 10 for full details.
5.6 Septic tank density
Septic tanks can be major contributors to nitrogen export from catchments (Bernhardt et al. 2008), are their density is a strong predictor of nitrogen concentrations in Melbourne’s streams (Walsh and Kunapo 2009). Sewer backlog programs around peri-urban Melbourne are aiming to replace septic systems with reticulated sewerage (Victorian Auditor-General’s Office 2006), but septic tanks are common throughout the region, even in some areas with reticulated sewerage.
I calculated catchment septic tank density (sep_per_km2) for a portion of the region and added to the cat_env table. Data were only available for the Yarra Ranges Council area (data source: Yarra Ranges Council from around the year 2000), with a small amount of further data for the Dobsons Creek catchment (data source: Yarra Valley Water 2017), and augmented estimates for parts of the Menzies Creek catchment that fall outside the Yarra Ranges Council (Figure 5.1).
5.7 Catchment geology
Many measures of catchment geology or soil type are possible, using Government spatial data such as the Geomorphological map units of Victoria (Rees et al. (2017)) and Victorian Soil type mapping (Anon. (2018)). Rather than generate a large number of geological or soil metrics, I calculated two example metrics from the geormophological map: proportion of catchment area with basalt geology (c_basalt), and granitic geology (c_granitic), both in the cat_env table, and provide the code I used to calculate them as a template for other similar calculations (including any polygon land-cover data for the region). The resulting distributions of catchment geology are show in Figure 5.2.
# GMU250.shp downloaded from data.vic.gov.au and stored in geol_dir
system.time({
geol <- terra::vect(paste0(geol_dir,"GMU250.shp"))
geol <- terra::project(geol, "epsg:7855")
geol <- terra::makeValid(geol)
# unique(geol$LITHOLOGY) #15 classes including NAs...
# sum(is.na(geol$LITHOLOGY)) #12 - all in MW region are waterbodies, so..
geol$LITHOLOGY[is.na(geol$LITHOLOGY)]<-'Water Body'
# combine...
geol$LITHOLOGY[grep('Granite|Granitic',geol$LITHOLOGY)]<-'Granitic'
geol$LITHOLOGY[grep('Basalt|Volcanics',geol$LITHOLOGY)]<-'Basalt'
lith_groups <- unique(geol$LITHOLOGY)
geol$lith_id <- match(geol$LITHOLOGY,lith_groups) # Basalt = 1; Granitic = 5.
r_site <- terra::rast("directory_of_choice/mwstr_v13/r_site.tif")
geol_r <- terra::rasterize(geol, r_site, field = "lith_id") # ~ 10 s
}) # ~30 s
#### Note the next zonal calculations require ~ 70 Gb RAM
system.time({
basalt_site <- terra::zonal(geol_r, r_site,
fun = function(x){sum(x == 1, na.rm = TRUE)})
}) # ~ 4.5 min
gc()
system.time({
granite_site <- terra::zonal(geol_r, r_site,
fun = function(x){sum(x == 5, na.rm = TRUE)}) # ~ 4.5 min
})
gc()
# Assume an active connection to mwstr called db_m
subcs_geol <- DBI::dbGetQuery(db_m,
paste0("SELECT subcs.site, subcs.nextds, subcs.nus, ",
"subcs.agg_order,subcs.scarea, cat_env.cat_area_m2_r ",
"FROM subcs JOIN cat_env ON subcs.site = cat_env.site;"))
DBI::dbDisconnect(db_m)
system.time({
subcs_geol <- subcs_geol[order(subcs_geol$nus, subcs_geol$agg_order),]
subcs_geol$sc_basalt <-
basalt_site$lith_id[match(subcs_geol$site, basalt_site$site)]
subcs_geol$sc_granite <-
granite_site$lith_id[match(subcs_geol$site, basalt_site$site)]
subcs_geol$c_granite <- subcs_geol$c_basalt <- NA
# All subcatchments of agg_order 1 have no subcatchments upstream. Their
# catchment area (in km2) = their subcatchment area (in m2 * 1e-6)
subcs_geol$c_basalt[subcs_geol$agg_order == 1] <-
subcs_geol$sc_basalt[subcs_geol$agg_order == 1]
subcs_geol$c_granite[subcs_geol$agg_order == 1] <-
subcs_geol$sc_granite[subcs_geol$agg_order == 1]
for(i in min(which(subcs_geol$agg_order == 2)):nrow(subcs_geol)){
# Loop through all remaining subcatchments, adding its subcatchment area to the
# (already calculated) catchment areas of subcatchments immediately upstream
subcs_geol$c_basalt[i] <-
sum(subcs_geol$c_basalt[subcs_geol$nextds == subcs_geol$site[i]],
subcs_geol$sc_basalt[i])
subcs_geol$c_granite[i] <-
sum(subcs_geol$c_granite[subcs_geol$nextds == subcs_geol$site[i]],
subcs_geol$sc_granite[i])
if(i %% 10000 == 0) cat(i,"\n")
}
}) # 1.5 min
# multiply by 25 because each gridcell is 25 m2
subcs_geol$c_granite <- subcs_geol$c_granite*25/subcs_geol$cat_area_m2_r
subcs_geol$c_basalt <- subcs_geol$c_basalt*25/subcs_geol$cat_area_m2_r
5.8 Elevation, slope and aspect
I calculated elevation (above sea-level in m, subc_env.elev_asl_m, Figure 5.3 A) and slope (percentage, subc_env.slope_perc, Figure 5.3 B) from the filled 5-m digital elevation model (table r_dem). I calculated aspect (as northness, Figure 5.3 C) from the stream lines. These three reach attributes are stored in the subc_env table.
The stream lines are accurately delineated throughout the region, having been derived from a 1-m-resolution DEM, and a more accurate estimate of slope could be achieved using that DEM. I took a hybrid approach to estimating slope of each reach by calculating the difference in elevatoin between the 5-m raster pourpoints for each subcatchment and its upstream subcatchment(s), and used the length of the stream line through the subcatchment as the horizontal distance.
Slope was calculated as a percentage as:
\[ 100 \cdot \frac{\Delta E}{\sqrt{{L}^2 + {\Delta E}^2}} \tag{5.1}\]
where \(\Delta E\) = the difference in elevation between upstream and downstream points and \(L\) is the total length of the stream lines (of all types) in the reach. Slopes were not calculated for headwater reaches, most of which were very short line segments. Negative slopes, of which there were 1672, almost all very small, were set to zero.
Aspect can influence how shaded a reach is with a reach flowing north-south likely to be less shaded by its banks and hillslopes than a reach flowing east-west. A more complete measure of such an effect could integrate shading calculated from a DEM and daily patterns of solar radiation, but northness is offered as an initial approximation. Northness was calculated from the angle between the most downstream and most upstream points of the line(s) of each reach. It was calculated as
\[ abs(cos(\theta)) \tag{5.2}\]
where \(\theta\) is the angle in radians from North between the two points. Thus for a reach flowing north to south (or vice-versa) northness = 1, and for a reach flowing east to west (or vice versa) northness = 0.
db_m <- RPostgres::dbConnect(RPostgres::Postgres(), dbname = "mwstr")
str_slope_elev <- sf::st_read(db_m, query =
"SELECT streams.geom, subc_env.elev_asl_m, subc_env.slope_perc,
subc_env.northness, streams.reach, streams.sampleable, subcs.agg_order
FROM streams JOIN subc_env ON streams.site = subc_env.site
JOIN subcs on subc_env.site = subcs.site;")
DBI::dbDisconnect(db_m)
# Order by agg_order to place largest streams on top of plot
str_slope_elev <- str_slope_elev[order(str_slope_elev$agg_order),]
# Set colors for maps
str_slope_elev$elev_class <-
as.numeric(cut(str_slope_elev$elev_asl_m,
breaks = c(-10,10,20,40,80,160,320,640,1340)))
str_slope_elev$elev_col <-
RColorBrewer::brewer.pal(8,"Spectral")[match(str_slope_elev$elev_class,1:8)]
str_slope_elev$slope_class <-
as.numeric(cut(str_slope_elev$slope_perc, breaks = c(0,1,2,4,8,16,32,75)))
str_slope_elev$slope_col <-
RColorBrewer::brewer.pal(7,"Spectral")[match(str_slope_elev$slope_class,1:7)]
str_slope_elev$north_class <-
as.numeric(cut(str_slope_elev$northness,
breaks = c(-0.1,0.2,0.4,0.6,0.8,1.1)))
str_slope_elev$north_col <-
RColorBrewer::brewer.pal(5,"Spectral")[match(str_slope_elev$north_class,1:5)]
# Plot maps
par(mfrow = c(2,2), mar = c(0,0,0,0))
plot(str_slope_elev$geom, col = str_slope_elev$elev_col)
title(main = "A.", adj = 0, line = -1)
legend("bottomleft", lty = 1, col = RColorBrewer::brewer.pal(8, "Spectral"),
legend = c("<10","10-20","20-40","40-80","80-160","160-320","320-640",">640"),
cex = 0.7, title = "Elevation (m)")
plot(str_slope_elev$geom, col = str_slope_elev$slope_col)
title(main = "B.", adj = 0, line = -1)
legend("bottomleft", lty = 1, col = RColorBrewer::brewer.pal(7, "Spectral"),
legend = c("<1","1-2","2-4","4-8","8-16","16-32",">32"),
cex = 0.7, title = "Slope (%)")
plot(str_slope_elev$geom, col = str_slope_elev$north_col)
title(main = "C.", adj = 0, line = -1)
legend("bottomleft", lty = 1, col = RColorBrewer::brewer.pal(5, "Spectral"),
legend = c("0-0.2","0.2-0.4","0.4-0.6","0.6-0.8","0.8-1,0"),
cex = 0.7, title = "Northness")
5.9 Strahler stream order
I calculated Strahler order (variable strahler in subc_env, a numerical measure of branching complexity) for the reaches of the network using the following code, which applies the algorithm described at https://en.wikipedia.org/wiki/Strahler_number.
subcs_s <- subcs[c("site","nextds","nus","agg_order")]
subcs_s <- subcs_s[order(subcs_s$nus, subcs$agg_order),]
subcs_s$strahler <- NA
subcs_s$strahler[subcs_s$nus %in% c(0)] <- 1
# If there are no upstream lines strahler = 1
for(i in min(which(is.na(subcs_s$strahler))):dim(subcs_s)[1]){
if(sum(subcs_s$nextds == subcs_s$site[i]) == 1){
# If the site is a split in a stream line not at a confluence make strahler the
# value of the site upstream
subcs_s$strahler[i] <- subcs_s$strahler[subcs_s$nextds == subcs_s$site[i]]
}else{
usi <- subcs_s[subcs_s$nextds == subcs_s$site[i],]
subcs_s$strahler[i] <- ifelse(sum(usi$strahler == max(usi$strahler)) == 1,
# If the node has one child with Strahler number i, and all other children have
# Strahler numbers less than i, then the Strahler number of the node is i again.
max(usi$strahler),
# If the node has two or more children with Strahler number i, and no children
# with greater number, then the Strahler number of the node is i + 1.
max(usi$strahler) + 1)
}
} # took ~3 minutesTo plot a map showing spatial patterns of Strahler order, run the following code (or a variant thereof as explained in (Appendix C)[app3-db_assembly.qmd])
db_m <- RSQLite::dbConnect(RSQLite::SQLite(),
paste0(mwstr_dir,"mwstr_v13.gpkg"))
str_strahler_subcat <- sf::st_read(db_m, query =
paste0("SELECT streams.geom, subc_env.strahler ",
"FROM streams JOIN subc_env ON streams.site = subc_env.site ",
"WHERE streams.sampleable = 1;"))
DBI::dbDisconnect(db_m)
plot(str_strahler_subcat["strahler"])5.10 Melbourne Water’s ‘sub-catchments’
An administrative measure rather than an environmental one, Melbourne Water’s ‘subcatchments’ are spatial management units used for stream management reporting, including environmental condition. I thus included them as a field in subc_env, linked to the table mw_subcatchments
The Melbourne Water Healthy Waterways Strategy (HWS) 2018 divides the region into 69 management units, which MW call “subcatchments” (although many are not hydrological subcatchments). To permit reporting for the HWS to these spatial units, each site was allocated to a “subcatchment” (the mw_subcat field) and a “catchment” (the mw_cat field). To allocate each site, I identified the most downstream sites for MW subcatchments in QGIS (saved in terminal_segs_mw_scs).
The numeric identifiers mw_subcat and mw_cat correspond to 69 ‘subcatchment’ names, and 8 ‘catchment’ names. These ‘subcatchments’ and ‘catchments’ are defined by polygons that roughly correspond to the stream network (from 1). The table mw_subcatchments contains those polygons, with fields mw_subcat, mw_cat, mw_subcat_nameand mw_cat_name, to permit matching of numeric mw_subcat and mw_cat codes to names.
To plot a map showing spatial patterns of mw_subcat and mw_subcat_name, run the following code (or a variant thereof as explained in (Appendix C)[app3-db_assembly.qmd])
db_m <- RSQLite::dbConnect(RSQLite::SQLite(),
paste0(mwstr_dir,"mwstr_v13.gpkg"))
str_mwcat_subcat <- sf::st_read(db_m, query =
paste0("SELECT streams.geom, subc_env.mw_subcat, subc_env.mw_cat ",
"FROM streams JOIN subc_env ON streams.site = subc_env.site ",
"WHERE streams.sampleable = 1;"))
DBI::dbDisconnect(db_m)
plot(str_mwcat_subcat["mw_subcat"])
plot(str_mwcat_subcat["mw_cat"])5.11 Australian postcodes
A second administrative measure that may be used to narrow searches of the database, a postcode was assigned to (the bottom of) each reach by joining the polygons of postcode areas from Australia Post with table subc_pps. Twenty reaches (mainly estuaries, pipe outlets and a few coastal tributaries on French Island) were not assigned a postcode.
To plot a map showing spatial patterns of postcodes, run the following code (or a variant thereof as explained in (Appendix C)[app3-db_assembly.qmd])
db_m <- RSQLite::dbConnect(RSQLite::SQLite(),
paste0(mwstr_dir,"mwstr_v13.gpkg"))
str_postcode_subcat <- sf::st_read(db_m, query =
paste0("SELECT streams.geom, subc_env.postcode ",
"FROM streams JOIN subc_env ON streams.site = subc_env.site ",
"WHERE streams.sampleable = 1;"))
DBI::dbDisconnect(db_m)
plot(str_postcode_subcat["postcode"])see Grant, A. (2018) Subcatchments - Healthy Waterways Strategy 2018 -Data guide and metadata. A copy of this document and the data it describes resides on the unimelb drive “wergSpatial/Projects/MW_StreamNetwork/env_compilation”↩︎