In mkkallio/hydrostreamer: A Package for Interpolating Distributed Runoff Products on to Explicit River Segments

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

This tutorial aims to explain the workflow in hydrostreamer and showcase some of it's capabilities. More specific tutorials will be added in time, which concentrate on different parts of hydrostreamer

Using hydrostreamer

Some example data have been included to the package for tutorial purposes. It includes

• 1 by 1 degree DEM located in Southeast Asia. The DEM is originally ALOS World 3D at 30 meter resolution (Tanado et al 2014), which has been resampled to 0.005 degree resolution.

• 1 by 1 degree runoff timeseries in the same area as the DEM. Runoff is sourced from the Linear Optimal Runoff Aggregate (LORA) at 0.5 degree resolution (see Hobeichi et al 2019). The unit of runoff is mm/s (kg/m2/s), and are provided with a monthly timestep.

• A river network derived from the provided DEM with 216 river segments.

• River segment specific catchment areas, delineated from the provided DEM. Catchments are provided for a subset (n = 41) of the river segments.

Let's first load the data and inspect it:

library(hydrostreamer)
library(raster)
library(lubridate)
library(dplyr)
library(sf)

data(example_rivers)
data(example_basins)
runoff <- brick(system.file("extdata", "runoff.tif", package = "hydrostreamer"))
dem <- brick(system.file("extdata", "dem.tif", package = "hydrostreamer")) 

plot(runoff[[1]]) 
plot(st_union(basins), add=TRUE)
plot(river, add=TRUE)

plot(dem)
plot(st_union(basins), add=TRUE)
plot(river, add=TRUE)

hydrostreamer main workflow occurs in three steps:

1. Areal interpolation of runoff to explicit river segments
2. Routing down the river network to estimate discharge
3. Model averaging, if streamflow observations are available.

1. Converting raster runoff to a polygon network

The raster layers are converted to polygons in order to do all the computations using only vector processing, and thus. Each cell of the raster is polygonized, and if an area of interest is provided, the polygons are clipped to it. This removes any unneeded grid cells. The resulting HS object is a standard 'sf' object with information about each raster cell. The runoff timeseries can be found in a named list column runoff_ts. The elements are named by the ID, and can be accessed with the '$' notation.

We use raster_to_HS to convert a raster timeseries to a HS object.

source_runoff <- raster_to_HS(runoff, 
                       unit = "mm/s",
                       date = ymd("1980-01-01"), 
                       timestep = "month", 
                       aoi = st_union(basins),
                       names = "LORA")
source_runoff
plot(source_runoff)

# access runoff timeseries of the element with ID `2` 
source_runoff$runoff_ts$`2`

HS objects can also be created from an 'sf' (polygon) object and a runoff timeseries using function create_HS.

2. Areal interpolation

The four areal interpolation methods shown here are explained in more detail in our recent (2019) conference paper here. The following figure is from the publication, visualizing the different methods:

Area Weighted Interpolation

Hydrostreamer implements several different areal interpolation methods, all of which can be accessed using the interpolate_runoff function. The simplest form implemented in hydrostreamer is Area Weighted Interpolation (AWI), which divides runoff from source zones to target river reaches .

AWI <- interpolate_runoff(source_runoff, basins, riverID = "SEGMENT_ID")
AWI

The output contains the columns from basins, added with riverID and runoff timeseries converted to volume (cubic meters per second).

AWI$runoff_ts$`56`

AWI$mean_runoff <- sapply(AWI$runoff_ts, function(x) mean(x$LORA))
plot(AWI[,"mean_runoff"])

Dasymetric Mapping

We can also refine the estimates using a dasymetric variable. Here I assume that more runoff is generated at higher elevations than lower, which may or may not be a good assumptions. It is, however, easy to extract from the DEM.

In dasymetric mapping (DM), the output of area weighted interpolation is further refined by scaling it using an ancillary variable - in this case elevation.

elevation_values <- raster::extract(dem, basins)
basins$elevation <- sapply(elevation_values, mean)

DM <- interpolate_runoff(source_runoff, basins, dasymetric = "elevation",
                         riverID = "SEGMENT_ID")

DM$mean_runoff <- sapply(DM$runoff_ts, function(x) mean(x$LORA))
plot(DM[,"mean_runoff"])

Pycnophylactic Interpolation

A third possibility is to use Pycnophylactic Interpolation (PP).

source_runoff$mean_runoff <- sapply(source_runoff$runoff_ts, function(x) {
  mean(x$LORA)
})

PP <- interpolate_runoff(source_runoff, basins, pycnophylactic = "mean_runoff",
                         riverID = "SEGMENT_ID")

PP$mean_runoff <- sapply(PP$runoff_ts, function(x) mean(x$LORA))
plot(PP[,"mean_runoff"])

Combined Pycnophylactic - Dasymetric Interpolation

The fourth option is to use a combination of PP and DM. In this case, PP is first performed instead of AWI, followed by the same scaling process with the dasymetric variable.

PPDM <- interpolate_runoff(source_runoff, basins, 
                         dasymetric = "elevation",
                         pycnophylactic = "mean_runoff",
                         riverID = "SEGMENT_ID")
PPDM$mean_runoff <- sapply(PPDM$runoff_ts, function(x) mean(x$LORA))
plot(PPDM[,"mean_runoff"])

Area-to-Line interpolation

All of the above examples are using catchment areas as the target units where runoff is estimated. However, existing vector-based river network data may not come specified with the DEM they were created with, or the processing steps are not clearly defined. In such cases, delineating the catchment areas for each segment may be difficult or outright impossible. For such cases, hydrostreamer provides the possibility to do interpolation from the source areas to linestrings. In such case, interpolation is based on the length of the linestring instead of the polygon area.

A2L <- interpolate_runoff(source_runoff, river,
                          riverID = "SEGMENT_ID")

A2L$mean_runoff <- sapply(A2L$runoff_ts, function(x) mean(x$LORA))
plot(A2L[,"mean_runoff"])

# DASYMETRIC MAPPING WITH LINES
river <- dplyr::filter(river, SEGMENT_ID %in% basins$SEGMENT_ID) %>% 
  dplyr::mutate(elevation = basins$elevation)
A2LDM <- interpolate_runoff(source_runoff, river,
                            dasymetric = "elevation",
                          riverID = "SEGMENT_ID")

A2LDM$mean_runoff <- sapply(A2LDM$runoff_ts, function(x) mean(x$LORA))
plot(A2LDM[,"mean_runoff"])

Using both lines and basins

The last example shown here is using both linestrings, and their respective catchment areas. The advantage here is that we can represent the rivers as they are in the network - as connected lines - and still use the catchment areas as a more accurate representation of the runoff produced in the catchment of the segment. In addition, using the lines allows us to use the constant river routing algorithm in the next step instead of only instantaneous routing.

RB <- interpolate_runoff(source_runoff,
                         river,
                         basins = basins,
                         riverID = "SEGMENT_ID")

RB$mean_runoff <- sapply(RB$runoff_ts, function(x) mean(x$LORA))
plot(RB[,"mean_runoff"])

3. Apply river routing

While the runoff generated at each river segment is already useful for many applications, knowing river discharge is often also desirable. hydrostreamer provides two simple river routing algorithms for this purpose: instantaneous routing, useful for e.g. estimating runoff in the entire upstream catchment of each river segment, and constant flow velocity routing.

Each routing method is accessible through the function accumulate_runoff(). Since the catchment provided here is small, and the timestep in runoff is one month, we'll just use instantaneous routing here. There would be negligible difference between constant velocity and instantaneous routing in this case.

However, we cannot use the catchments-only downscaled runoff here directly because there is no routing information. Since we can derive the routing information for the river lines, we can use that same routing info for the catchments also.

Note that running river_network() is not explicitly necessary, since the routing algorithm does it automatically if it has not been run in advance.

routed_river <- river_network(river, riverID = "SEGMENT_ID")
routed_river

PPDM$NEXT <- routed_river$NEXT
PPDM$PREVIOUS <- routed_river$PREVIOUS
PPDM$UP_SEGMENTS <- routed_river$UP_SEGMENTS

catchment_discharge <- accumulate_runoff(PPDM, routing_method = "instant")
river_discharge <- accumulate_runoff(RB, routing_method = "instant")

river_discharge

The algorithm adds a new list column discharge_ts, containing the routed discharge estimates.

We can also plot the estimated discharge from the two approaches at segment 200 by accessing the new list column discharge_ts.

plot(catchment_discharge$discharge_ts$`200`, type = 'l')
lines(river_discharge$discharge_ts$`200`, col='red')
title(main = "PPDM discharge in black, and Area-to-Line discharge in red.",
      sub = "Unit = m3/s, riverID = `200`")

As seen from the plots, in this catchment the difference between PPDM and Area-to-Line Interpolation is very small. In fact, we've found that as the basin size increases, the difference becomes increasingly small. This is particularly true for discharge at a monthly timestep. Larger difference can be expected with a daily timestep and with constant flow velocity routing.

References

T. Tadono, H. Ishida, F. Oda, S. Naito, K. Minakawa, H. Iwamoto : Precise Global DEM Generation By ALOS PRISM, ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Vol.II-4, pp.71-76, 2014.

Hobeichi, S., Abramowitz, G., Evans, J., and Beck, H. E.: Linear Optimal Runoff Aggregate (LORA): a global gridded synthesis runoff product, Hydrol. Earth Syst. Sci., 23, 851-870, https://doi.org/10.5194/hess-23-851-2019, 2019.

mkkallio/hydrostreamer documentation built on Oct. 14, 2023, 9:38 p.m.