Next Article in Journal
Comparison of Post-fire Patterns in Brazilian Savanna and Tropical Forest from Remote Sensing Time Series
Previous Article in Journal
A Continuous, Semi-Automated Workflow: From 3D City Models with Geometric Optimization and CFD Simulations to Visualization of Wind in an Urban Environment
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Orinoco: Retrieving a River Delta Network with the Fast Marching Method and Python

Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA 91109, USA
*
Author to whom correspondence should be addressed.
ISPRS Int. J. Geo-Inf. 2020, 9(11), 658; https://0-doi-org.brum.beds.ac.uk/10.3390/ijgi9110658
Submission received: 1 September 2020 / Revised: 16 October 2020 / Accepted: 27 October 2020 / Published: 31 October 2020

Abstract

:
We present Orinoco, an open-source Python toolkit that applies the fast-marching method to derive a river delta channel network from a water mask and ocean delineation. We are able to estimate flow direction, along-channel distance, channel width, and network-related metrics for deltaic analyses including the steady-state fluxes. To demonstrate the capabilities of the toolkit, we apply our software to the Wax Lake and Atchafalaya River Deltas using water masks derived from Open Street Map (OSM) and Google Maps. We validate our width estimates using the Global River Width from Landsat (GRWL) database over the Mackenzie Delta as well as in situ width measurements from the National Water Information System (NWIS) in the southeastern United States. We also compare the stream flow direction estimates using products from RivGraph, a related Python package with similar functionality. With the exciting opportunities afforded with forthcoming surface water and topography (SWOT) data, we envision Orinoco as a tool to support the characterization of the complex structure of river deltas worldwide and to make such analyses easily accessible within a Python remote sensing workflow. To support that end, all the data, analyses, and figures in this paper can be found within Jupyter notebooks at Orinoco’s GitHub repository.

1. Introduction

Deltas are complex coastal landforms that constitute a vital link between the terrestrial and marine environment [1]. Their ecological value and potential for economic development has attracted millions of people over the course of history, rendering them among the most densely populated regions globally [2]. However, the implications of human activities along with accelerated subsidence and climate-driven sea level rise are threatening the future existence of these low-lying coastal environments and the economic, environmental, and ecological services they provide [3,4,5,6,7,8]. Additionally, the construction of dams has altered the sediment load of rivers needed to build and maintain deltaic landforms [9,10], whereas groundwater withdrawal, mining activities, and the extraction of oil and gas result in seawater intrusion and increased soil compaction [11].
The geophysical characterization of deltas can improve the representation of hydrological processes in Earth systems models [12] and of their role in biogeochemical cycles between continents, oceans, and atmosphere [13,14,15,16,17]. Parameters such as river width are used as input for modeling hydraulic geometry relationships [18,19], estimating discharge [20,21,22], and mapping river and stream surface area [17,23,24,25,26]. Within the deltaic habitats, the fluvial connectivity influences vegetation productivity and nutrients delivery [27]. Furthermore, understanding the evolution of delta systems is highly relevant for their management and protection [28], including flood hazard prediction [29], and requires accurate data on the geometry of their channel networks [30].
NASA’s upcoming surface water and topography (SWOT) satellite mission, planned for launch in 2022 [31], will provide unprecedented coverage over deltas and will provide insight into many of these deltaic processes. SWOT will generate water surface elevation raster products with a spatial resolution of 100 m and also estimate water surface slopes with an accuracy of 1.7 cm/km along a 10 km stretch for rivers wider than 100 m. However, SWOT measurements along a water channel may be discontinuous due to mixed land–water pixels over thin river reaches, “dark water” returns caused by low radar backscatter, and layover from radar shadows cast by high banks and tall vegetation [32]. Thus, an open source tool to facilitate the computation of along-channel distance between sparse measurements is required to estimate water surface slope within complex and often narrow deltaic channels, which may be left out of a SWOT water mask. Such slope estimates can in turn be used to quantify discharge [33]. Although the Global River Width from Landsat (GRWL) database [26] was designed specifically to support such SWOT hydrology applications [34], it excludes numerous river deltas because they lie near mean sea level [26].
In this paper, we introduce a Python toolkit Orinoco that extracts a channel network [35] as a graph determined by the channel geometry to support the analysis and processing of SWOT data over the world’s deltas. We extract this network by efficiently leveraging the fast-marching method to segment the delta’s channels. Orinoco is named after the Venezuelan delta which translates from Guarauno to “navigable place”. A graph representation of a channel network can support SWOT characterization of discharge distribution [21,22] across an entire delta. The graph G = ( V , E ) is defined with nodes (or vertices) V and links (or edges) E, modeling the movement of water through the deltaic channel network [35]. Ideally, links should lie at the center of each channel and run nearly parallel to the stream flow direction [25]. The links of a channel network can be used to estimate channel centerlines in a delta [34], but not conversely, at least not without significant ad hoc processing as centerline geometries may have gaps within the channels or their geometries may not easily be translated into a network, e.g., at a river junction, an endpoint of a centerline may intersect another in its interior. A channel network can shed insight into a delta’s geomorphology and connectivity [36], which are highly relevant properties to manage deltas [28]. With Orinoco, we can easily compute along-channel distances which are invaluable for approximating water surface slope using remotely sensed height data such as will be provided by the SWOT mission. We also illustrate the computation of steady-state fluxes in the network, which can be used to understand the distribution of discharge at the river mouth assuming the network has been properly constructed.
There are also numerous other tools that extract the hydraulic geometry of the channels within deltas, each with its merits and drawbacks [37,38,39,40,41]. RivGraph, to our knowledge, is the only other Python library under development that extracts a channel network according to the graphical definition above [42,43]. Orinoco introduces the use of the fast-marching method [44] to numerically navigate and segment deltaic channels to generate a network. The segmentation, in particular, provides a useful product for investigating measurements within deltaic channels.
While the channel networks obtained by Orinoco and RivGraph are similar, the methodologies for doing so are very different: where we use the fast-marching method and a channel segmentation, RivGraph uses the skeletonization procedure from [45].

2. The Sites and Data

2.1. Test and Validation Sites

We demonstrate Orinoco and extract the relevant data products over the Wax Lake and Atchafalaya River Deltas in the Southern United States using Open Street Map (OSM) derived water masks. OSM tiles have coverage over deltas globally at approximately 30 m resolution and require very modest processing to obtain a water mask for analysis. The Wax Lake and Atchafalaya Deltas are located about 175 km south-southeast of Louisiana and each covers an area of about 100 km 2 . These deltas are the few deltas within the Mississippi deltaic floodplain that are prograding [27,46], even as most of coastal Louisiana is experiencing land loss at the alarmingly rapid rate of approximately 57 km 2 /year between 1932 and 2016 [47]. In the case of the Wax Lake Delta, sediments accumulate at the mouth of the Wax Lake Outlet due to a man-made channel dredged by the U.S. Army Corps of Engineers in the early 1940s to reduce flooding in Morgan City [48].
Validation of the width estimates, as well as a comparison of the centerlines and stream flow directions, is performed over the Mackenzie Delta in Northern Canada using the Global River Width from Landsat (GRWL) database [26] and products from RivGraph [42]. With an area of about 13,000 km 2 , the Mackenzie River delta is North America’s largest arctic delta and the world’s second largest after the Lena Delta in Russia [49]. Permafrost underlies approximately 75% of the delta’s drainage basin, and due to a disproportionately high rate of warming in the Arctics [50], this permafrost is thawing, changing the hydrology [51] and sediment dynamics [52] in the region.
Orinoco requires only a binary water mask and ocean delineation as input, and thus a user has the flexibility to create their own water mask for the application. This is particularly important for deltaic analyses where flooding, seasonal vegetation, thawing permafrost, and sediment deposition can quickly alter the geometry and topology of the channels [48,52]. Otherwise, with the same publicly available OSM tiles as used here, our channel networks can be extracted over numerous deltas globally.

2.2. Water Masks and Ocean Delineation

Creating a water mask over a region of interest is very difficult [37], and creating a robust methodology to do so is beyond the scope of this work. Orinoco requires only a binary water mask and ocean delineation, which can be constructed from remotely sensed data or freely available maps. The binary channel mask for all subsequent discussion will be a map such that land pixels have value 0 and water pixels have value 1.
The water masks used to demonstrate the capabilities of Orinoco in Section 3 are derived from Stamen Terrain tiles, which are a part of the Open Street Map (OSM) project [53], and downloaded using the open-source software stitch [54]. We also demonstrate our product using a high resolution water mask derived from Google Map tiles by selecting colors within the map that were water. Such water masks can be obtained over any delta using these publicly available tiles. We did very minimal pre-processing of the water masks prior to extracting the network: we applied a binary dilation followed by an erosion using [55]. This ensures that bridges and small islands that are less than 2 grid cells in width do not impact connectivity. Potentially spurious connections could be introduced with this pre-processing and high curvature channels may be erroneously simplified though we did not observe any such issues in our study areas as a majority of channels were separated by more than 2 grid cells. This dilation-erosion scheme is one of the simplest pre-processing procedures that can be applied to improve the quality of a water mask and we refer to [37] for numerous additional strategies. The ocean mask was created by hand to ensure the mouth of the river channels nicely agreed with the water mask, though one could automate this ocean delineation using [56], which we do in Section 4.
We also note that within a water mask, there may be missing channels and gaps (for example due to large bridges or narrow channels) that can impact the resulting network and products. We observed such issues in the OSM and GRWL-derived water masks used here, and frequently such challenges will likely have to be overcome on a case-by-case basis.
The ocean delineation is also required for Orinoco. For our demonstrations in Section 3, we manually draw such a delineation. For our validation with GRWL, USGS, and RivGraph, we use the data from [56]. If the ocean delineation is too far from the outlet of the channels, the resulting products may need to be adjusted manually. Generally, we found that if the delta’s ocean and channel interface did not match with [56], it was easier to manually draw such a polygon rather than update the final products of Orinoco.

3. Methodology and Software Capabilities

We generate a channel network using a segmentation of the channels. The segments are derived according to the distance from the ocean using the fast-marching method. This approach of segmentation is efficient because our network structure is determined by spatially contiguous areas within the channels rather than individual pixels. The centroid of each segment is identified as a node in V and connected to its neighboring segments with links in E. The width of each segment, which spans the width of the channel and is measured perpendicular to the direction of the estimated stream flow direction, is stored as a node attribute. Orinoco exports the channel network to NetworkX [57], a popular, well-maintained Python library for studying networks and their properties. This allows the channel network to be merged into a Python remote sensing workflow.

3.1. The Fast-Marching Method to Extract a Network

The distance to the ocean and the resulting segmentation is obtained using the fast-marching method, which efficiently solves a special case of the Eikenol equation [44]. Given a domain Ω R 2 and an interface Γ Ω described by a level set f ( x ) = 0 with f : Ω R , let a signed distance function φ : Ω R from the interface Γ be described via the partial differential equation [44]:
| φ ( x ) | = 1 , x Ω φ ( x ) = 0 , x Γ .
Above, the solution φ has a sign ambiguity, but is easily resolved assuming φ > 0 away from Γ . We assume that all the water masks are in a Universal Transverse Mercator (UTM) coordinate system (or a coordinate reference system whose resolution cells are uniform measurements of distance and not degrees) so that φ represents distance. The fast-marching method is analogous to Dijkstra’s algorithm on a 4-connected grid induced by the image pixels [58]. We use the open source Python package scikit-fmm [59]. In the future, such distance may be computed on the 8-connected image grid, but currently, this is not implemented in scikit-fmm [60].
We apply the fast-marching method to the following scenario. Let Ω be the inland channel water mask. Let Γ specify the interface between the ocean and these channels. We obtain the approximate distance to this interface using the fast-marching method.
Using φ , we then segment the river channels according to some fixed distance D. Specifically, the segment label ( x ) at a point x Ω is defined as ( x ) : = φ ( x ) / D where · is the numeric floor. In our toolkit, we specify the threshold in terms of the number of pixels T pixel to determine the segmentation setting D = T pixel · r , where r is the resolution of the underlying mask. We ensure that all segments are contiguous and relabel accordingly using scipy [55].
These segments, independent of the subsequent network construction, are valuable for analysis within the deltaic channels because we can aggregate statistics and measurements within localized spatial areas of the deltaic channels. In Section 4, we take advantage of this to compare Orinoco measurements with other available products and measurements. This observation, too, would easily allow us to combine the approaches from other related software packages such as [39,42] into this workflow as each segment corresponds uniquely to a node.
A schematic of our setup is shown in Figure 1 over the Wax Lake and Atchafalaya Deltas using our OSM-derived water mask. The ocean, land, and channel areas are shown in Figure 1a. The associated distance function φ is shown in Figure 1b and subsequent segments in Figure 1c using T pixel = 5 pixels over a water mask with resolution 30 m (i.e., D = 150 m). At this point, we observe that the gradient of φ will be parallel to the fronts of the segments within the channel. We will use φ to estimate the stream flow direction later in Section 3.3.
From the segments, we obtain the so-called region adjacency graph (RAG) [61] and determine connectivity between the derived segments. The nodes are the centroids of each segment. The links are determined according to the RAG adjacency, i.e., which segments share a common boundary. These links can be used to approximate the channel centerlines.

3.2. Estimating Stream Flow Direction along Edges

We direct the network’s links and determine the approximate stream flow direction using the channel network structure. First, we partition the links into contiguous groups between nodes that are either (a) a network junction, i.e., having degree larger than 2, (b) a node with only one connection, i.e., having degree 1, or (c) a node whose associated segment is adjacent to the interface Γ . Using this partition of links, we obtain paths in the network whose endpoints are either inlets, outlets, or a junction. We assume that water flows in one direction along any such path within the network. We direct each contiguous path based on which endpoint is closer to the interface Γ with respect to φ . In the rare event that two endpoints are equidistant to Γ , we direct the link in both directions and therefore create a directed multi-graph to indicate the uncertainty of stream flow direction.
We caution that this simple graph-based approach to estimate stream flow direction along each edge can assign unrealistic flow directions for intricate channel networks over large areas. There are more detailed approaches to this such as the excellent approach taken by RivGraph, which incorporates numerous auxiliary data sets and empirical rules [42]. In [62], the direction of stream flow is estimated using high resolution Digital Elevation Models (DEMs) over basins. In [63], multiscale strategies are used to ensure estimated stream flow direction is consistent across a braided delta network. Although such estimation procedures are currently beyond the capabilities of Orinoco, such approaches and products should be integrated for more robust analyses. As mentioned earlier, stream flow direction estimates from other products and tools could be integrated into the network’s attributes using the channel segmentation. We emphasize estimating stream flow direction, even using a more sophisticated method, is extremely difficult particularly in deltas where most of the deltaic area lies at mean sea level, where a DEM may not provide much information about the direction of stream flow.
As in [43], we prune “dangling” paths that are likely not informative about the channel connectivity. Such paths have a degree 1 endpoint, are not connected to Γ and are composed of less than or equal to T prune links, where T prune is a user-defined threshold. This pruning can be done iteratively. In all the examples presented in this paper, we remove “dangling” paths that contain less than or equal to 3 links and apply 1 iteration of this pruning procedure.
After we have directed our network and pruned it, we recompute each node’s distance to the ocean using the network structure. Each edge in our network has a length determined by its straight-line distance in the associated UTM coordinate system. The graph distance between two nodes in our network is then the sum of edge lengths along the shortest path between them. We define a node’s distance to the ocean as the minimum distance for a node to reach another node whose segment is adjacent to Γ in our segmentation. The graph distance to the ocean can thus be seen to approximate the distance along the channel’s centerlines to the ocean. The spatial resolution of this computed graph distance is approximately the product of the resolution of the water mask r and T pixel . Indeed, we will be underestimating the distance of high curvature channels that require finer spatial resolution with respect to this approach.
In Figure 2, we display such products over the Atchafalaya River Delta. Figure 2a illustrates the grouping of links between junctions, inlets, and outlets to direct the network. The resulting directed network is shown in Figure 2b with the estimated stream flow direction along edges going towards the ocean. We also show the along-channel distance obtained via the network structure in Figure 2c. We note that the distances to the ocean using the channel network are larger than those determined by φ because a distance determined via the channel network structure is constrained to travel in the approximate center of the channels, which is not the case for φ . In addition to our OSM-derived water mask, we demonstrate our methodology on the much more computationally demanding 2 m Google Maps derived water mask over the Wax Lake and Atchafalaya deltas. We show a subset of this network over the Wax Lake and the computed distance to the ocean in Figure 3. The estimated stream flow direction along edges using our graphical strategy is less reliable because there are numerous fine channels with numerous paths to the ocean and many junctions between them. However, the intricately mapped channels of such a high-resolution commercial product make this example potentially useful when considering frequent floating vegetation makes such channel networks harder to determine from remotely sensed data alone. We again manually made an ocean delineation for this network as we did for the OSM water mask.

3.3. Estimating Channel Width

We finally estimate the width at each node in our network. To do so, we first estimate the stream flow direction within each pixel of the channel mask and use this to approximate the stream flow direction at a particular node. Note that in the previous section we estimated stream flow along an edge, but we need a different approach to determine the stream flow direction at node, which may have numerous edges emanating from it. We model the local stream flow direction within the pixels of the channel according to φ , namely the direction of steepest descent relative to our distance function φ . We then associate the vector φ of the pixel at which the node is contained with that node’s stream flow direction. We define the candidate width geometry at a given node to be the line segment perpendicular to φ and whose midpoint is this node. We set the length of this candidate width geometry prior to intersection to be the perimeter of the node’s corresponding segment. We intersect this candidate line geometry with its corresponding segment and neighbors using [64] to obtain the width of the channel at the given node (a larger k-hop neighborhood can also be used thanks to the functionality of NetworkX). We provide an example of the estimated widths and their orientations in Figure 2d. Intersecting a candidate width with a neighborhood of segments rather than the full channel mask ensures that our widths are reasonably bounded within a neighborhood of the node which is particularly important near junctions. To this last point, if the width geometry is oriented in a direction almost perpendicular to the upstream stream flow direction, this bound prevents a width being significantly overestimated as seen in Figure 4c (we will discuss this more shortly). Moreover, having a smaller and simpler geometry for geometric intersection makes this process more efficient for the shapely algorithm; intersecting each candidate width geometry with the full complex geometry associated with the water mask is much more computationally expensive. We also note that some nodes may lie outside of the channel mask, e.g., when our segment is concave, and so φ is not defined and we cannot define a candidate width geometry. To ensure that all nodes have an associated width, we use the network structure to fill in missing widths. We simply take the average direction of a node’s neighboring links (once we reorient each link so they are within π / 2 of one another). In this average, we weight each link according to the inverse of its distance to that node; we assume closer neighboring nodes within the channel network are more important to determining the stream flow direction of the node under consideration. We note that this procedure is done entirely independent of φ because the distances are computed using the straight-line distance between nodes and therefore every node in our network will have a stream flow direction estimate using this approach. However, we found that this estimated stream flow direction was less reasonable visually in part because the network’s links direction could change abruptly due to channel islands or at junctions. Therefore, this approach is used to fill in those nodes without widths.
Near channel junctions, a node’s stream flow direction and hence the measured width is more challenging to determine. Firstly, our algorithm will tend to orient stream flow direction according to the incoming channel whose distance to the ocean is smaller as seen in Figure 2d. We note that the width geometries shown in Figure 2d are perpendicular to φ evaluated at the location of the node. Additional challenges occur if the difference in width of the two intersecting channels is large because the segments within this junction may be irregularly shaped making the placement of the node’s centroid not necessarily in the expected center of the channel. Lastly and importantly, intersecting channels introduce high curvature into the channel system and therefore φ is biased towards the downstream channels because that is the direction towards the ocean. We highlight some of these challenges in a subset over the Atchafalaya Delta in Figure 4. We see that at the junction at near 29.47 latitude in the center of the figure, the rightmost channel reaches this junction faster than the center-most channel. Therefore, the estimated stream flow direction is measured east-to-west rather than north-to-south, as we would expect; this can be seen inspecting the width geometries in Figure 4c. Furthermore, the difference in widths of the intersecting channels causes the segments to be more irregularly shaped, deviating from more rectangular segments seen away from junctions. At the junction at latitude 29.51 in the top-center of Figure 4c, we see that the width geometries and stream flow direction are biased towards the downstream channel and “lag” behind the expected stream flow direction as φ is pointing in the direction along the shortest path to the ocean in those areas. We also observe at this junction that even though the width geometries are roughly perpendicular to the larger channel it intersects, the geometric width is bounded. This occurs because we intersect the candidate widths geometry at these nodes near the junction only with neighboring segments (and not the full water mask) and this prevents a huge overestimate. As we move away from junctions, however, the width geometries become oriented as we might expect, particularly along those channels with low curvature as seen in Figure 2d.
Orinoco outputs the width geometries (in addition to the width estimates) as a final product. This allows users to not only inspect the final width measurement, but to inspect along what line the width is measured. This geometry output is valuable for doing detailed analysis over a delta’s channel and understanding why certain widths may deviate from other measurements and estimates of the channel’s widths. We summarize all of our steps for obtaining the deltaic channel network and related products in Algorithm 1.
Algorithm 1: A Summary of our Methodology.
Result: Channel network G = ( V , E ) and channel width w i for each i V .
Input: Channel Mask M, Interface Γ, pixel threshold T pixel , pruning threshold T prune , resolution r of
     water mask.
     Distance function φ ( M , Γ ) with fast-marching method [44].
     Segment Labels { i } i V ( φ , T pixel ) with i ( x ) : = x / ( T pixel · r ) .
     Undirected channel network G undirected = ( V , E ) { i } using RAG [61]
     Directed channel network G directed = ( V , E ) ( G undirected , φ ) from the link partition.
     Prune G directed as in [43] removing paths with T prune or less links.
     Channel widths { w i } i V ( G directed , φ , { i } ) .
     return G directed , { w i } i V

3.4. Network-Related Analysis and Computing the Steady-State Flux

The network structure affords many computations that are not possible with centerlines alone. Most straightforwardly, we can compute the along-channel distance using this network structure. Ideally, we would have the network’s edges aligned in the same direction as the stream flow direction so that a path between two nodes exists if and only if water from one node travels to another. However, even with an undirected network such a distance computation can still be performed and may be useful for estimating this along-channel distance. This distance can, in turn, be used to compute approximate slope over large areas when measurements of surface water heights are available within the delta’s channel. The SWOT mission will provide such measurements globally.
Additionally, the network structure and the channel’s width estimates allows us to compute the steady-state flux within the network. This can be used to compute the relative flux (or discharge) within the delta [65], and estimate how discharge is distributed at the delta’s outlets. This of course requires very careful accounting of the stream flow direction as well as all the sources and sinks within the network. Providing tools to assist with more efficient analysis may be valuable for detailed regional studies. Additionally, we compute the normalized entropy rate described in [35,36,65] that can be used to determine how evenly distributed this flux is throughout the delta.
We demonstrate these capabilities in a highly simplified setting using the OSM water mask. As noted in [65], for the computation to have physical meaning, we must have: all the sources and sinks accounted for (or assume unaccounted sources and sinks are negligible to the flux dynamics); have an accurate estimate of stream flow direction modeled by our network; and either assume the bathymetry is correlated throughout the network to width or incorporate the channel profile data into the edge attributes [65,66]. Such a careful set-up is beyond the scope of this work, though we mention this to illustrate this tool must be used in careful tandem with additional data sources and expert input. Our focus is to illustrate how the network structure can make such computations related to discharge much more readily accessible.
We concisely review the definition of network-related definitions following [35,65]. We first consider a subnetwork determined by a single source according to our estimated stream flow direction; if there are multiple channels where water is coming into the delta, then we would artificially connect them to a single additional source in the network. We select this source just south of the Intracoastal Waterway and the subnetwork is simply all those nodes reachable according to our estimated stream flow direction. We assume that this is the only source flowing downstream in our system. We merge our directed graph along links between junctions, outlets, and the source into a directed multigraph G m = ( V m , E m ) , meaning that two nodes may have multiple links between them. This merging is done relative to our subnetwork as junctions within the larger network may no longer be junctions in our subnetwork. This single-source multigraph over the Wax Lake Delta is shown in the second panel of Figure 5.
We then assign the width within each link of our multigraph according to the average width of the links contained in the original subnetwork. Using the widths to determine the flux at junctions, we set up the transition matrix P = ( p i j ) in which p i j is the probability that water at the junction i travels to j after one time step. This probability is determined according to the relative widths of the downstream links [35,65]. We then connect each outlet to the source so that the total water within the system is conserved after each time step. Computing the left eigenvector of P , we determine the steady state flux π, i.e., π T P = π T , with respect to the subnetwork, where π i is the percentage of water found at junction i as t . The steady state flux along the outlets is shown in the topmost panel of Figure 5; nodes are colored according to the relative flux along all outlets in this subnetwork. Assuming our network was constructed with the proper physical parameters, this would indicate how discharge would be distributed throughout each outlet in the delta.
The normalized entropy rate (nER) can then be computed from π . The nER is defined as
n E R = i V m π i j : ( i , j ) E m p i j log ( p i j ) log 1 d i
where d i is out-degree at node i within our subnetwork. The subnetwork in Figure 5 has nER = 0.532 , lower than the nER of ≈0.8 computed in [36] which was computed from a Landsat water mask. We note that our OSM derived water mask is missing some connections at the mouth of the Wax Lake Delta, though further investigation of the differences is beyond the scope of this work. The related notebook can be found at the Orinoco repository [67].

4. Validation of Orinoco Estimates

We validate Orinoco width and stream flow direction estimates comparing our results with products from GRWL, in situ width measurements available from the NWIS [68,69], and RivGraph products from [42]. We utilize the segments from the fast-marching method to compare nearby measurements, computing the mean of the various products within these segments.

4.1. Comparing Orinoco Widths with GRWL

Though some deltas are not a part of the GRWL database, many are included such as the Mackenzie Delta in Canada. We use the Mackenzie Delta (Tile NR08) to validate our width estimates and compare our channel centerlines with those from GRWL. While the GRWL widths are also estimates, they are frequently used for width estimates in deltaic channels and related analysis, supporting the need for comparison.
To obtain our channel network, we used the GRWL water masks. We determine the ocean boundary at the delta with [56]. Because some channels in the GRWL water masks are only connected via diagonal adjacency, we add a 1-pixel buffer to the GRWL water mask when computing φ using 4-connectivity. When computing our segments, we restrict φ to the original GRWL mask and employ 8-connectivity [70] to determine adjacency to ensure our widths agree with those obtained from GRWL. Ordinarily, this buffer may cause spurious connections in the channel network over more complex deltaic regions, when the space between two distinct channels is approximately equal to the spatial resolution of the water mask product. This was not the case in the Mackenzie Delta. In future work, we hope to expand the capabilities of Orinoco and be able to compute φ using 8-connectivity as discussed in [60].
For our GRWL comparison, we consider the GRWL widths within each segment. More explicitly, we assign each pixel intersecting a GRWL channel centerline with its associated GRWL width. We then average the widths within a segment excluding those pixels without data. For validation, we exclude segments that correspond to a network junction (i.e., have in- or out-degree greater than or equal to 2). We also exclude those segments corresponding to nodes that are adjacent to such network junctions. We remove such segments because our width estimate may be less accurate near such junctions as elaborated in Section 3.3.
For additional comparison, we also compute widths adapting the distance transform from scipy [55]. Explicitly, we consider the maximum distance d i within a segment i and assign a width to the corresponding node to be 2 d i 1 . This approach is an adaptation of the width computations found in [38,42]. This approach assumes that width geometries are approximated best to the closest land point from the centerline irrespective of stream flow direction. A more detailed comparison of such approaches to width measurements is discussed in [25].
We collect these comparisons in the scatter plot in Figure 6. We see that our width determination has comparable Root Mean Square Deviation (RMSD) and Mean Absolute Deviation (MAD) to the distance transform approach, but has significantly lower bias relative to GRWL. For context, our bias is subpixel (the GRWL water masks have 30 m resolution) and the bias using the distance transform is slightly over 1 pixel. This decrease in bias relative to GRWL is because our widths are measured perpendicular to the estimated stream flow which resembles more closely the approach taken with RivWidth [37] used to generate GRWL [26].
We inspect the difference of widths and connectivity over specific geographic subsets in Figure 7 and Figure 8. The yellow areas are those areas that have been excluded because the segment either does not contain a GRWL centerline, is a junction within the network, or is adjacent to one. We plot the Orinoco lines along which widths are measured, though we cannot make any additional visual width comparison with GRWL outside of this numeric difference. We see that Orinoco passes through smaller channels in the delta and moves through lakes to ensure connectivity throughout the deltaic channels. The Orinoco and GRWL widths have the highest deviation near junctions or channels with high curvature (lower left of Figure 8c). At the junctions, the computation of width is difficult due to the challenges associated with estimating stream flow direction there. In strongly curved channels, Orinoco’s width measurements may not align with the perpendicular of the channel’s center as the stream flow determined by φ may have some lag as discussed in Section 3.3. Furthermore, the analyses here are based on the mean of GRWL widths in the distance segments derived here and this comparison additionally makes the comparison to GRWL coarser.

4.2. Comparing Orinoco Widths with USGS Measurements and the GRWL Water Mask

We provide a comparison of in situ measurements of bank-to-bank reaches using NWIS data collected by the USGS. We use a GRWL water mask, specifically Tile NH15, for our width comparison. This data is obtained using the so-called Water Quality Data Portal [69,71] and can easily be validated using different GRWL tiles. This area is selected as it includes the Atchafalaya Delta, which we considered earlier. We demonstrate that our widths are within range of the in situ measurements using the GRWL water mask; we use the same setup for Orinoco as discussed in Section 4.1.
The width data is collected by the USGS and is the result of numerous in situ measurements at each station [71,72]. Because these measurements have quite a bit of variability, we consider their range of measurements and the qualitative statistics. In Figure 9a, we show the available stations within GRWL tile NH15. We have removed Matagorda because it’s width measurements varied by over 1 km. We further note these stations are those that lie within 60 m (a 2-pixel buffer) of the GRWL water mask; there are numerous other stations (68 stations in total) with width estimates in this GRWL tile though they are not located within or nearby one of the GRWL channels. In Figure 9b, we compare the ranges of USGS widths to the Orinoco width estimate. We see that in 6 of the 11 stations, Orinoco’s width lies within 1 standard deviation of the mean width. Furthermore, the Orinoco width estimate at 8 of the 11 stations is within the range of USGS measured widths; Ananhuac had only 1 measurement and the error (≈24 m) is smaller than the resolution of the water mask (30 m). In Figure 10, we show stations with more than one sample and whose measurement is not within 1 standard deviation of the mean USGS mean width. We see at Morgan City that the station lies near a junction and the width orientation is skewed by the fact that φ is biased by the downstream channel affecting the subsequent width measurement as discussed in Section 3.3. At Bayou Boeuf and Rosharon, the water mask appeared to be degraded. In particular, Bayou Boeuf had numerous freeways and bridges erroneously incorporated into the water width. We assume that the Ruliff error is likely due to the narrow channel and the sub-resolution variation is difficult to more carefully explain.

4.3. Comparing Stream Flow Direction Estimates with RivGraph

In this section, we demonstrate that our simplified approach to estimate stream flow direction shows acceptable agreement over a large majority of the MacKenzie River Delta. In particular, we measure the angle θ between centerlines produced by Orinoco and RivGraph using the data from [42] and then take the average angle between edges within each segment. Within the study area, at least 86.7% of segments have an average θ below 90 degrees and 84% of the segments have this average below 45 degrees indicating that the stream flow direction is comparable for those segments. We show the entire study area in Figure 11a. In Figure 11c, we show a subset of the delta where there is quite a bit of disagreement. We see that this area is where two different sub-networks with very distant outlets converge. This convergence and the difference in distance between these outlets likely contributes to this disagreement. In Figure 11d, we see more agreement overall likely because the junctions are only connecting different channels close to one primary outlet area. In the bottom panels of Figure 11, we show the estimated stream flows of Orinoco and RivGraph for further comparison. Again, RivGraph is much more advanced in its estimate of stream flow direction using numerous empirical rules. As noted earlier, we could integrate RivGraph estimates into the network for subsequent analyses using our segmentation.

5. Conclusions

We introduced Orinoco, an open-source Python library for extracting a deltaic channel network from a water mask and ocean delineation using the fast-marching method. Applying the fast-marching method, we obtain a channel segmentation that allows us to efficiently extract a network structure and compare channel-related measurements. We demonstrated the application of our software to large areas including an entire GRWL tile (Section 4) as well as a high resolution water mask over the Wax Lake and Atchafalaya using Google Map tiles. To illustrate the capabilities and products of our software, we obtained products from binary water masks derived from OSM or Google map tiles which are available across deltas globally.
The extracted network and its related products aim to support the processing and analysis of forthcoming SWOT data over the world’s deltas, though these tools can just as easily be used in conjunction with in situ data or other remotely sensed data over deltas such as airborne lidar and AirSWOT so long as there is a binary water mask and delineated area which water flows to or from.
We validated width estimates using GRWL products and USGS in situ measurements. We demonstrated good agreement, specifically subpixel bias, of our width products with the GRWL database over the Mackenzie Delta in Northern Canada, noting challenges at junctions and higher curvature channels. We compared USGS width data from the Water Quality Data Portal across the GRWL tile NH15 and saw a majority of the Orinoco width measurements were in the range of USGS measurements. In both cases, we inspected Orinoco outputs width geometries to explain the possible sources of errors related to our estimated channel width.
We also validated the stream flow direction estimates comparing Orinoco’s with those from RivGraph. While RivGraph has numerous more rules for estimating the stream flow direction and is likely more accurate, we saw overall good agreement using our purely network-based approach. We discussed the challenges associated with estimating stream flow direction and width using this new distance-motivated approach. We also noted that width and direction products from related tools can easily be integrated into the Orinoco networks using the segmentation of the channels.
We also demonstrated the value of the directed network structure in computing along-channel distance as well as computing the steady-state flux along the outlets, which can help estimate how discharge is distributed through a deltaic network. Although this demonstration of steady state flux was more theoretical, the ability to compute fluxes once a suitable network is constructed will be valuable for future discharge-related analyses.
Orinoco can easily be integrated into Python workflows. Additionally, the channel network being exported as a NetworkX graph can easily be updated and modified. All the analyses and plots shown in this paper are available at Orinoco’s GitHub page [67] within Jupyter Notebooks. We hope that the tools presented here can be used in conjunction with the numerous tools for studying remotely sensed data over deltas to promote robust hydrological analysis of these important coastal ecosystems.

Author Contributions

Each author helped to conceptualize this work. Marc Simard provided supervision of the research and guided the team toward SWOT-related applications. Charlie Marshak and Michael Denbina wrote the software with regular input from Marc Simard and Tom Van der Stocken. Charlie Marshak worked on the early draft of the manuscript. Marc Simard and Tom Van der Stocken greatly expanded this early version relating it to SWOT and deltaic ecology. The team worked on the final draft. Marc Simard, Charlie Marshak, and Johan Nilsson did the validation and formal analysis related to GRWL and USGS data. All authors have read and agreed to the published version of the manuscript.

Funding

This project was funded under the NASA Physical Oceanography Program and M.S. is a member of the SWOT Science Team.

Acknowledgments

This work was performed at the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration (NASA). We are grateful for the conversations with Kyle Wright and Jessica Fayne during the early stages of this project. We also are grateful to the anonymous reviewers for the extensive feedback to improve our original submission.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Evans, G. Deltas: The Fertile Dustbins of the Continents. Proc. Geol. Assoc. 2012, 123, 397–418. [Google Scholar] [CrossRef]
  2. Ericson, J.P.; Vörösmarty, C.J.; Dingman, S.L.; Ward, L.G.; Meybeck, M. Effective Sea-level Rise and Deltas: Causes of Change and Human Dimension Implications. Glob. Planet. Chang. 2006, 50, 63–82. [Google Scholar] [CrossRef]
  3. Syvitski, J.P. Deltas at Risk. Sustain. Sci. 2008, 3, 23–32. [Google Scholar] [CrossRef]
  4. Tessler, Z.D.; Vörösmarty, C.J.; Grossberg, M.; Gladkova, I.; Aizenman, H. A Global Empirical Typology of Anthropogenic Drivers of Environmental Change in Deltas. Sustain. Sci. 2016, 11, 525–537. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  5. Allison, M.A.; Nittrouer, C.A.; Ogston, A.S.; Mullarney, J.C.; Nguyen, T.T. Sedimentation and Survival of the Mekong Delta: A Case Study of Decreased Sediment Supply and Accelerating Rates of Relative Sea Level Rise. Oceanography 2017, 30, 98–109. [Google Scholar] [CrossRef] [Green Version]
  6. Schmitt, R.; Rubin, Z.; Kondolf, G. Losing Ground-Scenarios of Land Loss as Consequence of Shifting Sediment Budgets in the Mekong Delta. Geomorphology 2017, 294, 58–69. [Google Scholar] [CrossRef]
  7. Voosen, P. Seas Are Rising Faster than Believed at Many River Deltas. Science 2019. [Google Scholar] [CrossRef]
  8. Grill, G.; Lehner, B.; Thieme, M.; Geenen, B.; Tickner, D.; Antonelli, F.; Babu, S.; Borrelli, P.; Cheng, L.; Crochetiere, H. Mapping the World’s Free-Flowing Rivers. Nature 2019, 569, 215. [Google Scholar] [CrossRef]
  9. Syvitski, J.P.M.; Vörösmarty, C.J.; Kettner, A.J.; Green, P. Impact of Humans on the Flux of Terrestrial Sediment to the Global Coastal Ocean. Science 2005, 308, 376–380. [Google Scholar] [CrossRef]
  10. Blum, M.D.; Roberts, H.H. Drowning of the Mississippi Delta due to Insufficient Sediment Supply and Global Sea-level Rise. Nat. Geosci. 2009, 2, 488–491. [Google Scholar] [CrossRef]
  11. Syvitski, J.P.M.; Kettner, A.J.; Overeem, I.; Hutton, E.W.H.; Hannon, M.T.; Brakenridge, G.R.; Day, J.; Vörösmarty, C.; Saito, Y.; Giosan, L.; et al. Sinking Deltas due to Human Activities. Nat. Geosci. 2009, 2, 681–686. [Google Scholar] [CrossRef]
  12. Clark, M.P.; Fan, Y.; Lawrence, D.M.; Adam, J.C.; Bolster, D.; Gochis, D.J.; Hooper, R.P.; Kumar, M.; Leung, L.R.; Mackay, D.S. Improving the Representation of Hydrologic Processes in Earth System Models. Water Resour. Res. 2015, 51, 5929–5956. [Google Scholar] [CrossRef]
  13. Aufdenkampe, A.K.; Mayorga, E.; Raymond, P.A.; Melack, J.M.; Doney, S.C.; Alin, S.R.; Aalto, R.E.; Yoo, K. Riverine Coupling of Biogeochemical Cycles between Land, Oceans, and Atmosphere. Front. Ecol. Environ. 2011, 9, 53–60. [Google Scholar] [CrossRef] [Green Version]
  14. Bastviken, D.; Tranvik, L.J.; Downing, J.A.; Crill, P.M.; Enrich-Prast, A. Freshwater Methane Emissions Offset the Continental Carbon Sink. Science 2011, 331, 50. [Google Scholar] [CrossRef] [Green Version]
  15. Butman, D.; Raymond, P.A. Significant Efflux of Carbon Dioxide from Streams and Rivers in the United States. Nat. Geosci. 2011, 4, 839–842. [Google Scholar] [CrossRef]
  16. Bauer, J.E.; Cai, W.J.; Raymond, P.A.; Bianchi, T.S.; Hopkinson, C.S.; Regnier, P.A.G. The Changing Carbon Cycle of the Coastal Ocean. Nature 2013, 504, 61–70. [Google Scholar] [CrossRef]
  17. Raymond, P.A.; Hartmann, J.; Lauerwald, R.; Sobek, S.; McDonald, C.; Hoover, M.; Butman, D.; Striegl, R.; Mayorga, E.; Humborg, C.; et al. Global Carbon Dioxide Emissions from Inland Waters. Nature 2013, 503, 355–359. [Google Scholar] [CrossRef] [Green Version]
  18. Harman, C.; Stewardson, M.; DeRose, R. Variability and Uncertainty in Reach Bankfull Hydraulic Geometry. J. Hydrol. 2008, 351, 13–25. [Google Scholar] [CrossRef]
  19. Gleason, C.J.; Smith, L.C. Toward Global Mapping Of River Discharge Using Satellite Images and At-Many-Stations Hydraulic Geometry. Proc. Natl. Acad. Sci. USA 2014, 111, 4788–4791. [Google Scholar] [CrossRef] [Green Version]
  20. Bolla Pittaluga, M.; Tambroni, N.; Canestrelli, A.; Slingerland, R.; Lanzoni, S.; Seminara, G. Where River and Tide Meet: The Morphodynamic Equilibrium of Alluvial Estuaries. J. Geophys. Res. Earth Surf. 2015, 120, 75–94. [Google Scholar] [CrossRef]
  21. Yoon, Y.; Beighley, E.; Lee, H.; Pavelsky, T.; Allen, G. Estimating Flood Discharges in Reservoir-Regulated River basins by Integrating Synthetic SWOT Satellite Observations and Hydrologic Modeling. J. Hydrol. Eng. 2016, 21, 05015030. [Google Scholar] [CrossRef]
  22. Tuozzolo, S.; Lind, G.; Overstreet, B.; Mangano, J.; Fonstad, M.; Hagemann, M.; Frasson, R.; Larnier, K.; Garambois, P.A.; Monnier, J. Estimating River Discharge with Swath Altimetry: A Proof of Concept using AirSWOT Observations. Geophys. Res. Lett. 2019, 46, 1459–1466. [Google Scholar] [CrossRef]
  23. Lehner, B.; Döll, P. Development and Validation of a Global Database of Lakes, Reservoirs and Wetlands. J. Hydrol. 2004, 296, 1–22. [Google Scholar] [CrossRef]
  24. Downing, J.; Cole, J.; Duarte, C.; Middelburg, J.; Melack, J.; Prairie, Y.; Kortelainen, P.; Striegl, R.; McDowell, W.; Tranvik, L. Global abundance and size distribution of streams and rivers. Inland Waters 2012, 2, 229–236. [Google Scholar] [CrossRef]
  25. Yamazaki, D.; O’Loughlin, F.; Trigg, M.A.; Miller, Z.F.; Pavelsky, T.M.; Bates, P.D. Development of the Global Width Database for Large Rivers. Water Resour. Res. 2014, 50, 3467–3480. [Google Scholar] [CrossRef]
  26. Allen, G.H.; Pavelsky, T.M. Global Extent of Rivers and Streams. Science 2018, 361, 585–588. [Google Scholar] [CrossRef] [Green Version]
  27. Twilley, R.; Day, J.; Bevington, A.; Castañeda-Moya, E.; Christensen, A.; Holm, G.; Heffner, L.; Lane, R.; McCall, A.; Aarons, A.; et al. Ecogeomorphology of Coastal Deltaic Floodplains and Estuaries in an Active Delta: Insights from the Atchafalaya Coastal Basin. Estuar. Coast. Shelf Sci. 2019, 227, 106341. [Google Scholar] [CrossRef]
  28. Giosan, L.; Syvitski, J.; Constantinescu, S.; Day, J. Climate Change: Protect the World’s Deltas. Nat. News 2014, 516, 31. [Google Scholar] [CrossRef] [Green Version]
  29. Francesch-Huidobro, M.; Dabrowski, M.; Tai, Y.; Chan, F.; Stead, D. Governance Challenges of Flood-Prone Delta Cities: Integrating Flood Risk Management and Climate Change in Spatial Planning. Prog. Plan. 2017, 114, 1–27. [Google Scholar] [CrossRef]
  30. Liang, M.; Van Dyk, C.; Passalacqua, P. Quantifying the Patterns and Dynamics of River Deltas under Conditions of Steady Forcing and Relative Sea Level Rise. J. Geophys. Res. Earth Surf. 2016, 121, 465–496. [Google Scholar] [CrossRef]
  31. Biancamaria, S.; Lettenmaier, D.P.; Pavelsky, T.M. The SWOT mission and Its Capabilities for Land Hydrology. Surv. Geophys. 2016, 37, 307–337. [Google Scholar] [CrossRef] [Green Version]
  32. Fu, L.L.; Alsdorf, D.; Morrow, R.; Rodriguez, E. SWOT Mission Science Document. Available online: https://swot.jpl.nasa.gov/system/documents/files/2179_2179_SWOT_MSD_1202012.pdf (accessed on 15 September 2020).
  33. Frasson, R.P.d.M.; Pavelsky, T.M.; Fonstad, M.A.; Durand, M.T.; Allen, G.H.; Schumann, G.; Lion, C.; Beighley, R.E.; Yang, X. Global Relationships Between River Width, Slope, Catchment Area, Meander Wavelength, Sinuosity, and Discharge. Geophys. Res. Lett. 2019, 46, 3252–3262. [Google Scholar] [CrossRef] [Green Version]
  34. Rodriguez, E.; Williams, B.; Fore, A. RiverObs. 2010. Available online: https://github.com/SWOTAlgorithms/RiverObs (accessed on 15 September 2020).
  35. Tejedor, A.; Longjas, A.; Zaliapin, I.; Foufoula-Georgiou, E. Delta Channel Networks: 1. A Graph-Theoretic Approach for Studying Connectivity and Steady State Transport on Deltaic Surfaces. Water Resour. Res. 2015, 51, 3998–4018. [Google Scholar] [CrossRef]
  36. Tejedor, A.; Longjas, A.; Edmonds, D.A.; Zaliapin, I.; Georgiou, T.T.; Rinaldo, A.; Foufoula-Georgiou, E. Entropy and Optimality in River Deltas. Proc. Natl. Acad. Sci. USA 2017, 114, 11651–11656. [Google Scholar] [CrossRef] [Green Version]
  37. Pavelsky, T.M.; Smith, L.C. RivWidth: A Software Tool for the Calculation of River Widths from Remotely Sensed Imagery. IEEE Geosci. Remote Sens. Lett. 2008, 5, 70–73. [Google Scholar] [CrossRef]
  38. Fisher, G.B.; Amos, C.B.; Bookhagen, B.; Burbank, D.W.; Godard, V. Channel Widths, Landslides, Faults, and Beyond: The New World Order of High-Spatial Resolution Google Earth imagery in the Study of Earth Surface Processes. Geol. Soc. Am. Spec. Pap. 2012, 492, 1–22. [Google Scholar]
  39. Isikdogan, F.; Bovik, A.; Passalacqua, P. RivaMap: An automated river analysis and mapping engine. Remote Sens. Environ. 2017, 202, 88–97. [Google Scholar] [CrossRef]
  40. Golly, A.; Turowski, J.M. Deriving Principal Channel Metrics from Bank and Long-profile Geometry with the R package CMGO. Earth Surf. Dyn. 2017, 5, 557–570. [Google Scholar] [CrossRef] [Green Version]
  41. Yang, X.; Pavelsky, T.M.; Allen, G.H.; Donchyts, G. RivWidthCloud: An Automated Google Earth Engine Algorithm for River Width Extraction From Remotely Sensed Imagery. IEEE Geosci. Remote Sens. Lett. 2020, 17, 217–221. [Google Scholar] [CrossRef]
  42. Schwenk, J.; Piliouras, A.; Rowland, J.C. Determining Flow Directions in River Channel Networks using Planform Morphology and Topology. Earth Surf. Dyn. 2020, 8, 87–102. [Google Scholar] [CrossRef] [Green Version]
  43. Schwenk, J. RivGraph. 2019. Available online: https://github.com/jonschwenk/RivGraph (accessed on 1 September 2020).
  44. Sethian, J.A. A Fast Marching Level Set Method for Monotonically Advancing Fronts. Proc. Natl. Acad. Sci. USA 1996, 93, 1591–1595. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  45. Zhang, T.; Suen, C.Y. A Fast Parallel Algorithm for Thinning Digital Patterns. Commun. ACM 1984, 27, 236–239. [Google Scholar] [CrossRef]
  46. Jensen, D.; Cavanaugh, K.C.; Simard, M.; Okin, G.S.; Castañeda-Moya, E.; McCall, A.; Twilley, R.R. Integrating Imaging Spectrometer and Synthetic Aperture Radar Data for Estimating Wetland Vegetation Aboveground Biomass in Coastal Louisiana. Remote Sens. 2019, 11, 2533. [Google Scholar] [CrossRef] [Green Version]
  47. Couvillion, B.R.; Beck, H.; Schoolmaster, D.; Fischer, M. Land area Change in Coastal Louisiana (1932 to 2016); Technical Report; US Geological Survey: Denver, CO, USA, 2017.
  48. Fisk, H.N. Geological Investigation of the Atchafalaya Basin and the Problem of Mississippi River Diversion; Waterways Experiment Station: Vicksburg, MS, USA, 1952.
  49. Burn, C.R.; Kokelj, S. The Environment and Permafrost of the Mackenzie Delta Area. Permafr. Periglac. Process. 2009, 20, 83–105. [Google Scholar] [CrossRef]
  50. Overland, J.E.; Wang, M.; Walsh, J.E.; Stroeve, J.C. Future Arctic Climate Changes: Adaptation and Mitigation Time Scales. Earth’s Future 2014, 2, 68–74. [Google Scholar] [CrossRef]
  51. Yang, D.; Shi, X.; Marsh, P. Variability and Extreme of Mackenzie River Daily Discharge during 1973–2011. Quat. Int. 2015, 380, 159–168. [Google Scholar] [CrossRef]
  52. Doxaran, D.; Devred, E.; Babin, M. A 50% increase in the Mass of Terrestrial Particles Delivered by the Mackenzie River into the Beaufort Sea (Canadian Arctic Ocean) over the Last 10 Years. Biogeosciences 2015, 12, 3551–3565. [Google Scholar] [CrossRef] [Green Version]
  53. Haklay, M.; Weber, P. Open Street Map: User-Generated Street Maps. IEEE Pervasive Comput. 2008, 7, 12–18. [Google Scholar] [CrossRef] [Green Version]
  54. Fisher, E. Stitch. 2014. Available online: https://github.com/ericfischer/tile-stitch (accessed on 15 September 2020).
  55. Virtanen, P.; Gommers, R.; Oliphant, T.E.; Haberland, M.; Reddy, T.; Cournapeau, D.; Burovski, E.; Peterson, P.; Weckesser, W.; Bright, J.; et al. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nat. Methods 2020, 17, 261–272. [Google Scholar] [CrossRef] [Green Version]
  56. Department of Geography, UCLA. World Water Bodies. 2015. Available online: https://apps.gis.ucla.edu/geodata/dataset/world_water_bodies/resource/a6b40af0-84cb-40ce-b1c5-b024527a6943 (accessed on 15 September 2020).
  57. Hagberg, A.; Swart, P.; S Chult, D. Exploring Network Structure, Dynamics, and Function Using NetworkX; Technical Report; Los Alamos National Lab. (LANL): Los Alamos, NM, USA, 2008.
  58. Dijkstra, E.W. A Note on Two Problems in Connexion with Graphs. Numer. Math. 1959, 1, 269–271. [Google Scholar] [CrossRef] [Green Version]
  59. Furtney, J. Scikit-Fmm. 2010. Available online: https://https://github.com/scikit-fmm/scikit-fmm (accessed on 15 September 2020).
  60. Marshak, C. Orinoco. 2019. Available online: https://github.com/scikit-fmm/scikit-fmm/issues/32 (accessed on 15 September 2020).
  61. Saarinen, K. Color Image Segmentation by a Watershed Algorithm and Region Adjacency Graph Processing. In Proceedings of the 1st International Conference on Image Processing, Austin, TX, USA, 13–16 November 1994; Volume 3, pp. 1021–1025. [Google Scholar]
  62. Yamazaki, D.; Ikeshima, D.; Sosa, J.; Bates, P.D.; Allen, G.H.; Pavelsky, T.M. MERIT Hydro: A High-Resolution Global Hydrography Map based on Latest Topography Dataset. Water Resour. Res. 2019, 55, 5053–5073. [Google Scholar] [CrossRef] [Green Version]
  63. O’Loughlin, F.; Trigg, M.A.; Schumann, G.J.P.; Bates, P.D. Hydraulic Characterization of the Middle Reach of the Congo River. Water Resour. Res. 2013, 49, 5059–5070. [Google Scholar] [CrossRef]
  64. Gillies, S.; Bierbaum, A.; Lautaportti, K.; Tonnhofer, O. Shapely: Manipulation and Analysis of Geometric Objects; 2007-Present. Available online: https://github.com/Toblerity/Shapely (accessed on 15 September 2020).
  65. Tejedor, A.; Longjas, A.; Zaliapin, I.; Foufoula-Georgiou, E. Delta channel networks: 2. Metrics of Topologic and Dynamic Complexity for Delta Comparison, Physical Inference, and Vulnerability Assessment. Water Resour. Res. 2015, 51, 4019–4045. [Google Scholar] [CrossRef]
  66. Andreadis, K.M.; Schumann, G.J.P.; Pavelsky, T. A simple global river bankfull width and depth database. Water Resour. Res. 2013, 49, 7164–7168. [Google Scholar] [CrossRef]
  67. Marshak, C.; Simard, M.; Van der Stocken, T.; Denbina, M.; Nilsson, J. Orinoco. 2020. Available online: https://github.com/simard-landscape-lab/orinoco (accessed on 15 September 2020).
  68. Goodall, J.L.; Horsburgh, J.S.; Whiteaker, T.L.; Maidment, D.R.; Zaslavsky, I. A First Approach to Web Services for the National Water Information System. Environ. Model. Softw. 2008, 23, 404–411. [Google Scholar] [CrossRef]
  69. Read, E.K.; Carr, L.; De Cicco, L.; Dugan, H.A.; Hanson, P.C.; Hart, J.A.; Kreft, J.; Read, J.S.; Winslow, L.A. Water Quality Data for National-scale Aquatic Research: The Water Quality Portal. Water Resour. Res. 2017, 53, 1735–1745. [Google Scholar] [CrossRef]
  70. Annadurai, S. Fundamentals of Digital Image Processing; Pearson Education: Bengaluru, India, 2007. [Google Scholar]
  71. Water Quality Portal. 2020. Available online: https://www.waterqualitydata.us/portal (accessed on 15 September 2020).
  72. Hirsch, R.M.; Costa, J.E. US Stream Flow Measurement and Data Dissemination Improve. Eos Trans. Am. Geophys. Union 2004, 85, 197–203. [Google Scholar] [CrossRef]
Figure 1. Illustration of our methodology to obtain segments for a channel network using Stamen Terrain-derived water masks over the Wax Lake and Atchaflaya Deltas. The resolution of the water mask is 30 m (zoom level 12). Panel (a) shows the initialization of areas including the water mask. Panel (b) shows the distance function φ. The segments in (c) are obtained according to φ / D with D   =   150 m with each segment randomly colored. We note in panel (a) the ocean is a partially transparent yellow, so land areas in the ocean mask appear as orange.
Figure 1. Illustration of our methodology to obtain segments for a channel network using Stamen Terrain-derived water masks over the Wax Lake and Atchaflaya Deltas. The resolution of the water mask is 30 m (zoom level 12). Panel (a) shows the initialization of areas including the water mask. Panel (b) shows the distance function φ. The segments in (c) are obtained according to φ / D with D   =   150 m with each segment randomly colored. We note in panel (a) the ocean is a partially transparent yellow, so land areas in the ocean mask appear as orange.
Ijgi 09 00658 g001
Figure 2. Example of Orinoco data products over the Atchafalaya River Delta. Panel (a) shows the partition of the links into contiguous paths as specified in the text. Using this partition and φ, panel (b) shows the estimated stream flow direction. Each node’s distance to the ocean within the channel network is shown in panel (c) and the width and their associated orientation in the channels, in panel (d).
Figure 2. Example of Orinoco data products over the Atchafalaya River Delta. Panel (a) shows the partition of the links into contiguous paths as specified in the text. Using this partition and φ, panel (b) shows the estimated stream flow direction. Each node’s distance to the ocean within the channel network is shown in panel (c) and the width and their associated orientation in the channels, in panel (d).
Ijgi 09 00658 g002
Figure 3. Using Google map tiles, we select water areas to create a 2 m water mask over the Wax Lake Delta. Panel (a) shows the channel centerlines and panel (b) are the nodes colored by their distance to the ocean according to the network structure.
Figure 3. Using Google map tiles, we select water areas to create a 2 m water mask over the Wax Lake Delta. Panel (a) shows the channel centerlines and panel (b) are the nodes colored by their distance to the ocean according to the network structure.
Ijgi 09 00658 g003
Figure 4. Panel (a) shows the distance function, (b) the segments (randomly colored), and (c) the width along and associated orientation over a subset of the Atchafalaya Delta using the OSM-derived water mask shown in Figure 2. The black dots in panel (c) are the network’s nodes. We note the width’s geometry should be approximately perpendicular to the front of the segment it belongs to as this will be parallel to φ .
Figure 4. Panel (a) shows the distance function, (b) the segments (randomly colored), and (c) the width along and associated orientation over a subset of the Atchafalaya Delta using the OSM-derived water mask shown in Figure 2. The black dots in panel (c) are the network’s nodes. We note the width’s geometry should be approximately perpendicular to the front of the segment it belongs to as this will be parallel to φ .
Ijgi 09 00658 g004
Figure 5. Using the channel network from the OSM water mask, we compute the steady state fluxes to estimate the percentage of discharge from each outlet in the Wax Lake considering a source chosen just south of the Intracoastal Waterway. The topmost panel shows the stationary fluxes normalized across the outlets. The bottom panel illustrates the directed multigraph to compute the flux and normalized entropy rate. In the bottom panel, the source (blue square) and the outlets (green circles) are connected in our computation of the stationary distribution π described in [36].
Figure 5. Using the channel network from the OSM water mask, we compute the steady state fluxes to estimate the percentage of discharge from each outlet in the Wax Lake considering a source chosen just south of the Intracoastal Waterway. The topmost panel shows the stationary fluxes normalized across the outlets. The bottom panel illustrates the directed multigraph to compute the flux and normalized entropy rate. In the bottom panel, the source (blue square) and the outlets (green circles) are connected in our computation of the stationary distribution π described in [36].
Ijgi 09 00658 g005
Figure 6. The comparison of GRWL widths over the Mackenzie (NR08) (a) with the Orinoco widths and (b) with widths derived using the distance transform. We have removed segments whose corresponding nodes are junctions or are adjacent to them. The deviation statistics are computed subtracting the values of the estimated width (x-axis) from the GRWL segment-level width (y-axis). The percent error statistics (%) are relative to the mean GRWL widths (y-axis).
Figure 6. The comparison of GRWL widths over the Mackenzie (NR08) (a) with the Orinoco widths and (b) with widths derived using the distance transform. We have removed segments whose corresponding nodes are junctions or are adjacent to them. The deviation statistics are computed subtracting the values of the estimated width (x-axis) from the GRWL segment-level width (y-axis). The percent error statistics (%) are relative to the mean GRWL widths (y-axis).
Ijgi 09 00658 g006
Figure 7. Comparison of GRWL and Orinoco over a subset of the Mackenzie (NR08). Top Row: comparison of channel centerlines and connectivity in (a) GRWL and (b) Orinoco. Bottom Row: (a) is the per segment difference in widths (when applicable) and (d) is the histogram of errors. Note that we also have removed segments whose corresponding node in the network is a junction or adjacent to one (see panel (c), “No Data”).
Figure 7. Comparison of GRWL and Orinoco over a subset of the Mackenzie (NR08). Top Row: comparison of channel centerlines and connectivity in (a) GRWL and (b) Orinoco. Bottom Row: (a) is the per segment difference in widths (when applicable) and (d) is the histogram of errors. Note that we also have removed segments whose corresponding node in the network is a junction or adjacent to one (see panel (c), “No Data”).
Ijgi 09 00658 g007aIjgi 09 00658 g007b
Figure 8. Comparison of GRWL and Orinoco over a subset of the Mackenzie (NR08). Top Row: comparison of channel centerlines and connectivity in (a) GRWL and (b) Orinoco. Bottom Row: (a) is the per segment difference in widths (when applicable) and (d) is the histogram of errors. Note that we also have removed segments whose corresponding node in the network is a junction or adjacent to one (see panel (c), “No Data”).
Figure 8. Comparison of GRWL and Orinoco over a subset of the Mackenzie (NR08). Top Row: comparison of channel centerlines and connectivity in (a) GRWL and (b) Orinoco. Bottom Row: (a) is the per segment difference in widths (when applicable) and (d) is the histogram of errors. Note that we also have removed segments whose corresponding node in the network is a junction or adjacent to one (see panel (c), “No Data”).
Ijgi 09 00658 g008aIjgi 09 00658 g008b
Figure 9. In the top panel (a), we show those USGS stations intersecting the GRWL water mask or those stations within 60 m. In the bottom panel (b), we compare the closest Orinoco width with the range of USGS widths as well as the the mean width µ and the range of widths within 1 standard deviation µ ± σ . The value n below the USGS station name indicates how many measurements were made.
Figure 9. In the top panel (a), we show those USGS stations intersecting the GRWL water mask or those stations within 60 m. In the bottom panel (b), we compare the closest Orinoco width with the range of USGS widths as well as the the mean width µ and the range of widths within 1 standard deviation µ ± σ . The value n below the USGS station name indicates how many measurements were made.
Ijgi 09 00658 g009
Figure 10. The USGS stations and closest Orinoco node along with the measured width. In the background is the GRWL water mask. The error is Orinoco width subtracted from the mean USGS measurement. The percent error is this error as a percent relative to the mean USGS measurement.
Figure 10. The USGS stations and closest Orinoco node along with the measured width. In the background is the GRWL water mask. The error is Orinoco width subtracted from the mean USGS measurement. The percent error is this error as a percent relative to the mean USGS measurement.
Ijgi 09 00658 g010
Figure 11. Estimated stream flow compared with RivGraph over the Mackenzie Delta.
Figure 11. Estimated stream flow compared with RivGraph over the Mackenzie Delta.
Ijgi 09 00658 g011
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Marshak, C.; Simard, M.; Denbina, M.; Nilsson, J.; Van der Stocken, T. Orinoco: Retrieving a River Delta Network with the Fast Marching Method and Python. ISPRS Int. J. Geo-Inf. 2020, 9, 658. https://0-doi-org.brum.beds.ac.uk/10.3390/ijgi9110658

AMA Style

Marshak C, Simard M, Denbina M, Nilsson J, Van der Stocken T. Orinoco: Retrieving a River Delta Network with the Fast Marching Method and Python. ISPRS International Journal of Geo-Information. 2020; 9(11):658. https://0-doi-org.brum.beds.ac.uk/10.3390/ijgi9110658

Chicago/Turabian Style

Marshak, Charlie, Marc Simard, Michael Denbina, Johan Nilsson, and Tom Van der Stocken. 2020. "Orinoco: Retrieving a River Delta Network with the Fast Marching Method and Python" ISPRS International Journal of Geo-Information 9, no. 11: 658. https://0-doi-org.brum.beds.ac.uk/10.3390/ijgi9110658

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop