Add river-constrained D8 flow direction

Project.toml

@@ -15,6 +15,9 @@ SpatialRasterLite = "25a4c9b9-c310-4d6e-b9a7-fa79274760a1"
 ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
+[extensions]
+RiverGraphsArchGDALExt = "ArchGDAL"
+
 [sources]
 SpatialRasterLite = {url = "https://github.com/jl-pkgs/SpatialRasterLite.jl"}
 

ext/RiverGraphsArchGDALExt.jl

@@ -0,0 +1,89 @@
+module RiverGraphsArchGDALExt
+
+import ArchGDAL as AG
+import RiverGraphs
+import RiverGraphs: read_river_lines, hydro_enforced_flowdir
+
+"""
+    read_river_lines(path::AbstractString; layer=0, strict=true)
+
+Read LineString/MultiLineString features from a GDAL-supported vector file and
+return a `Vector{Vector{Tuple{Float64,Float64}}}`. MultiLineString parts are
+returned as separate river lines.
+
+Coordinates are not reprojected. The vector data must use the same CRS as the
+raster axes passed later to `rasterize_flowpath` / `hydro_enforced_flowdir`.
+
+With `strict=true`, non-line geometries raise an error. With `strict=false`,
+unsupported or empty geometries are skipped. GeometryCollection is traversed
+recursively, so mixed collections can be handled with `strict=false`.
+"""
+function read_river_lines(path::AbstractString;
+  layer::Integer=0, strict::Bool=true)
+  lines = Vector{Vector{Tuple{Float64,Float64}}}()
+
+  AG.read(path) do dataset
+    nlayers = AG.nlayer(dataset)
+    0 <= layer < nlayers || throw(ArgumentError(
+      "layer index $layer is outside 0:$(nlayers - 1) for $path"))
+
+    source = AG.getlayer(dataset, layer)
+    for feature in source
+      geom = AG.getgeom(feature, 0)
+      _append_river_geometry!(lines, geom; strict)
+    end
+  end
+
+  lines
+end
+
+"""
+    hydro_enforced_flowdir(dem, path, lon, lat; layer=0, strict=true, kwargs...)
+
+Read a GIS river vector file with ArchGDAL and run RiverGraphs' complete
+hydro-enforcement pipeline on the supplied raster grid.
+"""
+function hydro_enforced_flowdir(dem::AbstractMatrix, path::AbstractString,
+  lon::AbstractVector, lat::AbstractVector;
+  layer::Integer=0, strict::Bool=true, kwargs...)
+  lines = read_river_lines(path; layer, strict)
+  RiverGraphs.hydro_enforced_flowdir(dem, lines, lon, lat; kwargs...)
+end
+
+function _append_river_geometry!(lines, geom; strict::Bool)
+  name = AG.geomname(geom)
+  if ismissing(name)
+    strict && throw(ArgumentError("river feature has an empty geometry"))
+    return lines
+  end
+
+  typename = uppercase(String(name))
+  if typename == "LINESTRING" || typename == "LINEARRING"
+    n = Int(AG.ngeom(geom))
+    if n < 2
+      strict && throw(ArgumentError("river LineString has fewer than two points"))
+      return lines
+    end
+
+    line = Vector{Tuple{Float64,Float64}}(undef, n)
+    @inbounds for i in 0:n-1
+      point = AG.getpoint(geom, i)
+      line[i+1] = (Float64(point[1]), Float64(point[2]))
+    end
+    push!(lines, line)
+
+  elseif typename == "MULTILINESTRING" || typename == "GEOMETRYCOLLECTION"
+    for i in 0:Int(AG.ngeom(geom))-1
+      part = AG.getgeom(geom, i)
+      _append_river_geometry!(lines, part; strict)
+    end
+
+  elseif strict
+    throw(ArgumentError(
+      "unsupported river geometry $typename; expected LineString or MultiLineString"))
+  end
+
+  lines
+end
+
+end # module RiverGraphsArchGDALExt

src/RiverGraphs.jl

@@ -8,6 +8,10 @@ import RTableTools: cbind, fwrite, fread
 
 include("IO.jl")
 include("RiverGraph.jl")
+include("validate_flowdir.jl")
+include("constrained_flowdir.jl")
+include("priority_flood.jl")
+include("flowdir_qc.jl")
 include("fillnodata_upbasin.jl")
 include("fillnodata_upriver.jl")
 
@@ -22,12 +26,23 @@ include("utils.jl")
 include("sf.jl")
 
 
+"""
+    read_river_lines(source; kwargs...)
+
+Read river centreline geometries as vectors of `(x, y)` coordinate tuples.
+Methods for GIS vector files are provided by the optional ArchGDAL extension;
+load ArchGDAL before calling this function on a file path.
+"""
+function read_river_lines end
+
+
 export RiverGraph
 export active_indices, reverse_index, pcr_dir,
   graph_flow,
   topological_sort_kahn,
   stream_order, stream_link, stream_network,
   fillnodata_upbasin, fillnodata_upriver
+export read_river_lines
 
 const path_flowdir_GuanShan = abspath("$(@__DIR__)/../data/GuanShan_flwdir.tif")