Skip to content

Commit 23dca3d

Browse files
committed
add waterflow algorithm that flows water over topography
1 parent f946bf4 commit 23dca3d

File tree

5 files changed

+102
-1
lines changed

5 files changed

+102
-1
lines changed

Project.toml

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -26,6 +26,7 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
2626
StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
2727
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
2828
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
29+
WhereTheWaterFlows = "ea906314-1493-4d22-a0af-f886a20c9fba"
2930
WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
3031

3132
[weakdeps]
@@ -69,4 +70,4 @@ StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
6970
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
7071

7172
[targets]
72-
test = ["Test", "GMT", "StableRNGs","GridapGmsh"]
73+
test = ["Test", "GMT", "StableRNGs", "GridapGmsh"]

src/GeophysicalModelGenerator.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -47,6 +47,7 @@ include("IO.jl")
4747
include("event_counts.jl")
4848
include("surface_functions.jl")
4949
include("movies_from_pics.jl")
50+
include("WaterFlow.jl")
5051

5152
# Add optional routines (only activated when the packages are loaded)
5253

src/WaterFlow.jl

Lines changed: 80 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,80 @@
1+
# These are routes to perform waterflow calculations (upstream area etc.) on a DEM
2+
3+
using WhereTheWaterFlows
4+
import WhereTheWaterFlows: waterflows
5+
6+
export waterflows
7+
8+
"""
9+
dlon, dlat = spacing(lon,lat)
10+
11+
Computes the spacing with central differences
12+
"""
13+
function spacing(lon,lat)
14+
dlon = zero(lon.val[:,:,1])
15+
dlat = zero(lat.val[:,:,1])
16+
dlon[2:end-1,:] = (lon.val[3:end,:,1] - lon.val[1:end-2,:,1])/2
17+
dlon[1,:], dlon[end,:] = dlon[2,:], dlon[end-1,:]
18+
dlat[:,2:end-1] = (lat.val[:,3:end,1] - lat.val[:,1:end-2,1])/2
19+
dlat[:,1], dlat[:,end] = dlat[:,2], dlat[:,end-1]
20+
return dlon, dlat
21+
end
22+
23+
"""
24+
area_m2 = cell_area(Topo::GeoData)
25+
Returns the cell area for a Topographic dataset in m² (required for upstream area calculation)
26+
"""
27+
function cell_area(Topo::GeoData)
28+
29+
proj = ProjectionPoint(Lon=mean(Topo.lon.val[:]), Lat=mean(Topo.lat.val[:]))
30+
Topo_cart = convert2CartData(Topo, proj)
31+
dx, dy = spacing(Topo_cart.x, Topo_cart.y)
32+
33+
area_m2 = dx.*dy*1e6
34+
return area_m2
35+
end
36+
37+
38+
"""
39+
Topo_water, sinks, pits, bnds = waterflows(Topo::GeoData;
40+
flowdir_fn=WhereTheWaterFlows.d8dir_feature, feedback_fn=nothing, drain_pits=true, bnd_as_sink=true)
41+
42+
Takes a GMG GeoData object of a topographic map and routes water through the grid. We add a number of
43+
44+
Example
45+
===
46+
```julia
47+
# Download some topographic data
48+
julia> Topo = import_topo([6.5,7.3,50.2,50.6], file="@earth_relief_03s");
49+
50+
# Flow the water through the area:
51+
julia> Topo_water, sinks, pits, bnds = waterflows(Topo)
52+
julia> Topo_water
53+
GeoData
54+
size : (961, 481, 1)
55+
lon ϵ [ 6.5 : 7.3]
56+
lat ϵ [ 50.2 : 50.59999999999999]
57+
depth ϵ [ 0.045 : 0.724]
58+
fields : (:Topography, :area, :slen, :dir, :nout, :nin, :c)
59+
60+
```
61+
62+
"""
63+
function waterflows(Topo::GeoData, flowdir_fn= WhereTheWaterFlows.d8dir_feature; feedback_fn=nothing, drain_pits=true, bnd_as_sink=true)
64+
65+
cellarea = cell_area(Topo)
66+
dem = Topo.depth.val[:,:,1]
67+
68+
area = zeros(Float64,size(Topo.depth.val))
69+
slen = zeros(Int64,size(Topo.depth.val))
70+
dir = zeros(Int8,size(Topo.depth.val))
71+
nout = zeros(Int8,size(Topo.depth.val))
72+
nin = zeros(Int8,size(Topo.depth.val))
73+
c = zeros(Int64,size(Topo.depth.val))
74+
75+
area[:,:,1], slen[:,:,1], dir[:,:,1], nout[:,:,1], nin[:,:,1], sinks, pits, c[:,:,1], bnds = waterflows(dem, cellarea, flowdir_fn;
76+
feedback_fn=feedback_fn, drain_pits=drain_pits, bnd_as_sink=bnd_as_sink)
77+
78+
Topo_water = addfield(Topo,(;area, slen, dir, nout, nin, c ))
79+
return Topo_water, sinks, pits, bnds
80+
end

test/runtests.jl

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -68,6 +68,10 @@ end
6868
include("test_create_movie.jl")
6969
end
7070

71+
@testset "Waterflow" begin
72+
include("test_WaterFlow.jl")
73+
end
74+
7175
# Cleanup
7276
foreach(rm, filter(endswith(".vts"), readdir()))
7377
foreach(rm, filter(endswith(".vtu"), readdir()))

test/test_WaterFlow.jl

Lines changed: 15 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,15 @@
1+
using Test, GMT
2+
3+
4+
# Download some topographic data
5+
Topo = import_topo([6.5,7.3,50.2,50.6], file="@earth_relief_03s");
6+
7+
# Flow the water through the area:
8+
Topo_water, sinks, pits, bnds = waterflows(Topo)
9+
10+
@test maximum(Topo_water.fields.area) 9.309204547276944e8
11+
@test sum(Topo_water.fields.c) == 834501044
12+
@test sum(Topo_water.fields.nin) == 459361
13+
@test sum(Topo_water.fields.dir) == 2412566
14+
15+

0 commit comments

Comments
 (0)