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

: 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 ﬂow direction, along-channel distance, channel width, and network-related metrics for deltaic analyses including the steady-state ﬂuxes. 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 ﬂow 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 workﬂow. To support that end, all the data, analyses, and ﬁgures in this paper can be found within Jupyter notebooks at Orinoco’s GitHub repository.


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 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.

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.

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.

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]: 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.
(a) Inland channels and ocean (b) Distance function (ϕ) in channels to interface Γ.
(c) Segmentation of channels using ϕ

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. (a) (b) 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.

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. 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 ∇ϕ.
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].
Prune G directed as in [43] removing paths with T prune or less links.

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. 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].
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 ij ) in which p ij 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 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].

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.

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 2d 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 Figures 7 and 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.
(a) (b) 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).

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.   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.

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.

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.  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.