add waterflow algorithm that flows water over topography

Author: boriskaus

Project.toml

@@ -26,6 +26,7 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+WhereTheWaterFlows = "ea906314-1493-4d22-a0af-f886a20c9fba"
 WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
 
 [weakdeps]
@@ -69,4 +70,4 @@ StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "GMT", "StableRNGs","GridapGmsh"]
+test = ["Test", "GMT", "StableRNGs", "GridapGmsh"]

src/GeophysicalModelGenerator.jl

@@ -47,6 +47,7 @@ include("IO.jl")
 include("event_counts.jl")
 include("surface_functions.jl")
 include("movies_from_pics.jl")
+include("WaterFlow.jl")
 
 # Add optional routines (only activated when the packages are loaded)
 

src/WaterFlow.jl

@@ -0,0 +1,80 @@
+# These are routes to perform waterflow calculations (upstream area etc.) on a DEM
+
+using WhereTheWaterFlows
+import WhereTheWaterFlows: waterflows
+
+export waterflows
+
+"""
+    dlon, dlat = spacing(lon,lat)
+
+Computes the spacing with central differences
+"""
+function spacing(lon,lat)
+    dlon = zero(lon.val[:,:,1])
+    dlat = zero(lat.val[:,:,1])
+    dlon[2:end-1,:] = (lon.val[3:end,:,1] - lon.val[1:end-2,:,1])/2
+    dlon[1,:], dlon[end,:] = dlon[2,:], dlon[end-1,:]
+    dlat[:,2:end-1] = (lat.val[:,3:end,1] - lat.val[:,1:end-2,1])/2
+    dlat[:,1], dlat[:,end] = dlat[:,2], dlat[:,end-1]
+    return dlon, dlat
+end
+
+"""
+    area_m2 = cell_area(Topo::GeoData)
+Returns the cell area for a Topographic dataset in m² (required for upstream area calculation)
+"""
+function cell_area(Topo::GeoData)
+
+    proj = ProjectionPoint(Lon=mean(Topo.lon.val[:]), Lat=mean(Topo.lat.val[:]))
+    Topo_cart = convert2CartData(Topo, proj)
+    dx, dy = spacing(Topo_cart.x, Topo_cart.y)
+
+    area_m2 = dx.*dy*1e6
+    return area_m2
+end
+
+
+"""
+    Topo_water, sinks, pits, bnds  = waterflows(Topo::GeoData; 
+        flowdir_fn=WhereTheWaterFlows.d8dir_feature, feedback_fn=nothing, drain_pits=true, bnd_as_sink=true)
+
+Takes a GMG GeoData object of a topographic map and routes water through the grid. We add a number of 
+
+Example
+===
+```julia
+# Download some topographic data
+julia> Topo = import_topo([6.5,7.3,50.2,50.6], file="@earth_relief_03s");
+
+# Flow the water through the area:
+julia> Topo_water, sinks, pits, bnds  = waterflows(Topo)
+julia> Topo_water
+GeoData 
+  size      : (961, 481, 1)
+  lon       ϵ [ 6.5 : 7.3]
+  lat       ϵ [ 50.2 : 50.59999999999999]
+  depth     ϵ [ 0.045 : 0.724]
+  fields    : (:Topography, :area, :slen, :dir, :nout, :nin, :c)
+
+```
+
+"""
+function waterflows(Topo::GeoData, flowdir_fn= WhereTheWaterFlows.d8dir_feature; feedback_fn=nothing, drain_pits=true, bnd_as_sink=true) 
+
+    cellarea = cell_area(Topo)
+    dem = Topo.depth.val[:,:,1]
+
+    area    = zeros(Float64,size(Topo.depth.val))
+    slen    = zeros(Int64,size(Topo.depth.val))
+    dir     = zeros(Int8,size(Topo.depth.val))
+    nout    = zeros(Int8,size(Topo.depth.val))
+    nin     = zeros(Int8,size(Topo.depth.val))
+    c       = zeros(Int64,size(Topo.depth.val))
+
+    area[:,:,1], slen[:,:,1], dir[:,:,1], nout[:,:,1], nin[:,:,1], sinks, pits, c[:,:,1], bnds = waterflows(dem, cellarea, flowdir_fn;
+                        feedback_fn=feedback_fn, drain_pits=drain_pits, bnd_as_sink=bnd_as_sink)
+
+    Topo_water =  addfield(Topo,(;area, slen, dir, nout, nin, c ))                        
+    return Topo_water, sinks, pits, bnds 
+end

test/runtests.jl

@@ -68,6 +68,10 @@ end
     include("test_create_movie.jl")
 end
 
+@testset "Waterflow" begin
+    include("test_WaterFlow.jl")
+end
+
 # Cleanup
 foreach(rm, filter(endswith(".vts"), readdir()))
 foreach(rm, filter(endswith(".vtu"), readdir()))

test/test_WaterFlow.jl

@@ -0,0 +1,15 @@
+using Test, GMT
+
+
+# Download some topographic data
+Topo = import_topo([6.5,7.3,50.2,50.6], file="@earth_relief_03s");
+
+# Flow the water through the area:
+Topo_water, sinks, pits, bnds  = waterflows(Topo)
+
+@test maximum(Topo_water.fields.area) ≈ 9.309204547276944e8
+@test sum(Topo_water.fields.c) == 834501044
+@test sum(Topo_water.fields.nin) == 459361
+@test sum(Topo_water.fields.dir) == 2412566
+
+