From 23dca3d61cc4e9cd53b6652f297f9a5f521abe01 Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Mon, 1 Jul 2024 21:13:33 +0200 Subject: [PATCH 01/19] add waterflow algorithm that flows water over topography --- Project.toml | 3 +- src/GeophysicalModelGenerator.jl | 1 + src/WaterFlow.jl | 80 ++++++++++++++++++++++++++++++++ test/runtests.jl | 4 ++ test/test_WaterFlow.jl | 15 ++++++ 5 files changed, 102 insertions(+), 1 deletion(-) create mode 100644 src/WaterFlow.jl create mode 100644 test/test_WaterFlow.jl diff --git a/Project.toml b/Project.toml index 24b5a30c..dfee1050 100644 --- a/Project.toml +++ b/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"] diff --git a/src/GeophysicalModelGenerator.jl b/src/GeophysicalModelGenerator.jl index a5fe4891..9baf1652 100644 --- a/src/GeophysicalModelGenerator.jl +++ b/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) diff --git a/src/WaterFlow.jl b/src/WaterFlow.jl new file mode 100644 index 00000000..7eaefcd2 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 03b52fbb..96be2999 100644 --- a/test/runtests.jl +++ b/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())) diff --git a/test/test_WaterFlow.jl b/test/test_WaterFlow.jl new file mode 100644 index 00000000..cbeb632b --- /dev/null +++ b/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 + + From 8e36c35efe8d2721f82d8dc19e3a3bafc16800e5 Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Mon, 1 Jul 2024 21:19:14 +0200 Subject: [PATCH 02/19] update spellcheck --- .typos.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/.typos.toml b/.typos.toml index 54be0588..353f06a0 100644 --- a/.typos.toml +++ b/.typos.toml @@ -2,6 +2,7 @@ MOR = "MOR" dum = "dum" Shepard = "Shepard" +nin = "nin" [files] extend-exclude = ["tutorials/*.pvsm"] From 798af1f50d699a6c114a9517b1526a88edaa74f9 Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Mon, 1 Jul 2024 21:23:42 +0200 Subject: [PATCH 03/19] bump version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index dfee1050..9140a900 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "GeophysicalModelGenerator" uuid = "3700c31b-fa53-48a6-808a-ef22d5a84742" authors = ["Boris Kaus", "Marcel Thielmann"] -version = "0.7.4" +version = "0.7.5" [deps] Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" From 62834051cb13d939a9a50791c79182c109432544 Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Mon, 1 Jul 2024 21:58:51 +0200 Subject: [PATCH 04/19] restruict GMT version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 9140a900..5f47d7cb 100644 --- a/Project.toml +++ b/Project.toml @@ -45,7 +45,7 @@ Downloads = "1" FFMPEG = "0.4" FileIO = "1" GLMakie = "0.8, 0.9, 0.10" -GMT = "1" +GMT = "1.0 - 1.14" GeoParams = "0.2 - 0.6" Geodesy = "1" GeometryBasics = "0.1 - 0.4" From 6232aba7edc05ea616742f7293fd6f88c231624b Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Mon, 1 Jul 2024 22:20:08 +0200 Subject: [PATCH 05/19] allow specifying rainfall on the grid; also returns the area of each grid cell in m2 --- src/WaterFlow.jl | 18 ++++++++++++++---- test/test_WaterFlow.jl | 5 +++++ 2 files changed, 19 insertions(+), 4 deletions(-) diff --git a/src/WaterFlow.jl b/src/WaterFlow.jl index 7eaefcd2..bdbb63be 100644 --- a/src/WaterFlow.jl +++ b/src/WaterFlow.jl @@ -37,9 +37,13 @@ end """ Topo_water, sinks, pits, bnds = waterflows(Topo::GeoData; - flowdir_fn=WhereTheWaterFlows.d8dir_feature, feedback_fn=nothing, drain_pits=true, bnd_as_sink=true) + flowdir_fn=WhereTheWaterFlows.d8dir_feature, feedback_fn=nothing, drain_pits=true, bnd_as_sink=true, + rainfall = nothing) -Takes a GMG GeoData object of a topographic map and routes water through the grid. We add a number of +Takes a GMG GeoData object of a topographic map and routes water through the grid. Optionally, +you can specify `rainfall` in which case we accumulate the rain as specified in this 2D array instead of the cellarea. +This allows you to, for example, sum, up water if you have variable rainfall in the area. +The other options are as in the `waterflows` function of the package `WhereTheWaterFlows`. Example === @@ -60,9 +64,15 @@ GeoData ``` """ -function waterflows(Topo::GeoData, flowdir_fn= WhereTheWaterFlows.d8dir_feature; feedback_fn=nothing, drain_pits=true, bnd_as_sink=true) +function waterflows(Topo::GeoData, flowdir_fn= WhereTheWaterFlows.d8dir_feature; feedback_fn=nothing, drain_pits=true, bnd_as_sink=true, rainfall=nothing) cellarea = cell_area(Topo) + cellarea_m2 = cellarea + if !isnothing(rainfall) + @assert typeof(rainfall) == Array{Float64,2} + cellarea = rainfall + end + dem = Topo.depth.val[:,:,1] area = zeros(Float64,size(Topo.depth.val)) @@ -75,6 +85,6 @@ function waterflows(Topo::GeoData, flowdir_fn= WhereTheWaterFlows.d8dir_feature; 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 )) + Topo_water = addfield(Topo,(;area, slen, dir, nout, nin, c, cellarea_m2 )) return Topo_water, sinks, pits, bnds end \ No newline at end of file diff --git a/test/test_WaterFlow.jl b/test/test_WaterFlow.jl index cbeb632b..eed5703d 100644 --- a/test/test_WaterFlow.jl +++ b/test/test_WaterFlow.jl @@ -12,4 +12,9 @@ Topo_water, sinks, pits, bnds = waterflows(Topo) @test sum(Topo_water.fields.nin) == 459361 @test sum(Topo_water.fields.dir) == 2412566 +# With rain in m3/s per cell +rainfall = ones(size(Topo.lon.val[:,:,1]))*1e-3 # 2D array with rainfall per cell area +Topo_water1, sinks, pits, bnds = waterflows(Topo, rainfall=rainfall) + +@test maximum(Topo_water1.fields.area) ≈ 169.79800000000208 From aaacb02602e12b48cceab25e7fae1d70ddcdf8c5 Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Tue, 2 Jul 2024 12:36:43 +0200 Subject: [PATCH 06/19] fix GDAL_jll version --- Project.toml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Project.toml b/Project.toml index 5f47d7cb..a7425d40 100644 --- a/Project.toml +++ b/Project.toml @@ -9,6 +9,7 @@ DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" FFMPEG = "c87230d0-a227-11e9-1b43-d7ebe4e7570a" FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" +GDAL_jll = "a7073274-a066-55f0-b90d-d619367d196c" GeoParams = "e018b62d-d9de-4a26-8697-af89c310ae38" Geodesy = "0ef565a4-170c-5f04-8de2-149903a85f3d" GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326" @@ -44,6 +45,7 @@ Colors = "0.9 - 0.12" Downloads = "1" FFMPEG = "0.4" FileIO = "1" +GDAL_jll = "300.900.0 - 301.900.0" GLMakie = "0.8, 0.9, 0.10" GMT = "1.0 - 1.14" GeoParams = "0.2 - 0.6" From af993c4d76a7707ccad02f53f67f5447be3abacf Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Tue, 2 Jul 2024 12:37:02 +0200 Subject: [PATCH 07/19] add large catchment areas --- src/WaterFlow.jl | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/WaterFlow.jl b/src/WaterFlow.jl index bdbb63be..a888be9e 100644 --- a/src/WaterFlow.jl +++ b/src/WaterFlow.jl @@ -38,7 +38,8 @@ end """ Topo_water, sinks, pits, bnds = waterflows(Topo::GeoData; flowdir_fn=WhereTheWaterFlows.d8dir_feature, feedback_fn=nothing, drain_pits=true, bnd_as_sink=true, - rainfall = nothing) + rainfall = nothing, + minsize=300) Takes a GMG GeoData object of a topographic map and routes water through the grid. Optionally, you can specify `rainfall` in which case we accumulate the rain as specified in this 2D array instead of the cellarea. @@ -64,7 +65,7 @@ GeoData ``` """ -function waterflows(Topo::GeoData, flowdir_fn= WhereTheWaterFlows.d8dir_feature; feedback_fn=nothing, drain_pits=true, bnd_as_sink=true, rainfall=nothing) +function waterflows(Topo::GeoData, flowdir_fn= WhereTheWaterFlows.d8dir_feature; feedback_fn=nothing, drain_pits=true, bnd_as_sink=true, rainfall=nothing, minsize=300) cellarea = cell_area(Topo) cellarea_m2 = cellarea @@ -85,6 +86,8 @@ function waterflows(Topo::GeoData, flowdir_fn= WhereTheWaterFlows.d8dir_feature; 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, cellarea_m2 )) + catchment_large = prune_catchments(c, minsize; val=0) + + Topo_water = addfield(Topo,(;area, slen, dir, nout, nin, c, cellarea_m2, catchment_large)) return Topo_water, sinks, pits, bnds end \ No newline at end of file From c17239c1d023ef3a2d983e900675c89e44817f83 Mon Sep 17 00:00:00 2001 From: Boris Kaus <61824822+boriskaus@users.noreply.github.com> Date: Tue, 2 Jul 2024 14:14:54 +0200 Subject: [PATCH 08/19] Update src/WaterFlow.jl Co-authored-by: Albert de Montserrat <58044444+albert-de-montserrat@users.noreply.github.com> --- src/WaterFlow.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/WaterFlow.jl b/src/WaterFlow.jl index a888be9e..bdb13e21 100644 --- a/src/WaterFlow.jl +++ b/src/WaterFlow.jl @@ -11,8 +11,8 @@ export waterflows Computes the spacing with central differences """ function spacing(lon,lat) - dlon = zero(lon.val[:,:,1]) - dlat = zero(lat.val[:,:,1]) + dlon = zero(size(lon.val)[1:2]) + dlat = zero(size(lat.val)[1:2]) 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 From 07deceb7de864451b1134edea2a2688776a9ae0f Mon Sep 17 00:00:00 2001 From: Boris Kaus <61824822+boriskaus@users.noreply.github.com> Date: Tue, 2 Jul 2024 14:14:59 +0200 Subject: [PATCH 09/19] Update src/WaterFlow.jl Co-authored-by: Albert de Montserrat <58044444+albert-de-montserrat@users.noreply.github.com> --- src/WaterFlow.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/WaterFlow.jl b/src/WaterFlow.jl index bdb13e21..d35bd599 100644 --- a/src/WaterFlow.jl +++ b/src/WaterFlow.jl @@ -13,10 +13,10 @@ Computes the spacing with central differences function spacing(lon,lat) dlon = zero(size(lon.val)[1:2]) dlat = zero(size(lat.val)[1:2]) - 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] + @views dlon[2:end-1,:] = (lon.val[3:end,:,1] - lon.val[1:end-2,:,1])/2 + @views dlon[1,:], dlon[end,:] = dlon[2,:], dlon[end-1,:] + @views dlat[:,2:end-1] = (lat.val[:,3:end,1] - lat.val[:,1:end-2,1])/2 + @views dlat[:,1], dlat[:,end] = dlat[:,2], dlat[:,end-1] return dlon, dlat end From 9324c16af3f85293f576e65168d47c7303ab30ad Mon Sep 17 00:00:00 2001 From: Boris Kaus <61824822+boriskaus@users.noreply.github.com> Date: Tue, 2 Jul 2024 14:15:05 +0200 Subject: [PATCH 10/19] Update src/WaterFlow.jl Co-authored-by: Albert de Montserrat <58044444+albert-de-montserrat@users.noreply.github.com> --- src/WaterFlow.jl | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/WaterFlow.jl b/src/WaterFlow.jl index d35bd599..2cd0a379 100644 --- a/src/WaterFlow.jl +++ b/src/WaterFlow.jl @@ -76,12 +76,13 @@ function waterflows(Topo::GeoData, flowdir_fn= WhereTheWaterFlows.d8dir_feature; 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)) + ni = size(Topo.depth.val) + area = zeros(ni) + slen = zeros(Int64, ni) + dir = zeros(Int8, ni) + nout = zeros(Int8, ni) + nin = zeros(Int8, ni) + c = zeros(Int64, ni) 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) From 8069bd960cbddaf34ab6598054a43833fe84b70f Mon Sep 17 00:00:00 2001 From: Boris Kaus <61824822+boriskaus@users.noreply.github.com> Date: Tue, 2 Jul 2024 14:15:12 +0200 Subject: [PATCH 11/19] Update src/WaterFlow.jl Co-authored-by: Albert de Montserrat <58044444+albert-de-montserrat@users.noreply.github.com> --- src/WaterFlow.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/WaterFlow.jl b/src/WaterFlow.jl index 2cd0a379..ae95be2a 100644 --- a/src/WaterFlow.jl +++ b/src/WaterFlow.jl @@ -84,7 +84,7 @@ function waterflows(Topo::GeoData, flowdir_fn= WhereTheWaterFlows.d8dir_feature; nin = zeros(Int8, ni) c = zeros(Int64, ni) - area[:,:,1], slen[:,:,1], dir[:,:,1], nout[:,:,1], nin[:,:,1], sinks, pits, c[:,:,1], bnds = waterflows(dem, cellarea, flowdir_fn; + @views 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) catchment_large = prune_catchments(c, minsize; val=0) From c6c28b214e5a35b70bf810b97274db105e331e81 Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Wed, 3 Jul 2024 08:54:18 +0200 Subject: [PATCH 12/19] fix CI --- src/WaterFlow.jl | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/src/WaterFlow.jl b/src/WaterFlow.jl index ae95be2a..17ec243e 100644 --- a/src/WaterFlow.jl +++ b/src/WaterFlow.jl @@ -14,9 +14,12 @@ function spacing(lon,lat) dlon = zero(size(lon.val)[1:2]) dlat = zero(size(lat.val)[1:2]) @views dlon[2:end-1,:] = (lon.val[3:end,:,1] - lon.val[1:end-2,:,1])/2 - @views dlon[1,:], dlon[end,:] = dlon[2,:], dlon[end-1,:] - @views dlat[:,2:end-1] = (lat.val[:,3:end,1] - lat.val[:,1:end-2,1])/2 - @views dlat[:,1], dlat[:,end] = dlat[:,2], dlat[:,end-1] + dlon[1,:] = dlon[2,:] + dlon[end,:] = dlon[end-1,:] + dlat[:,2:end-1] = (lat.val[:,3:end,1] - lat.val[:,1:end-2,1])/2 + dlat[:,1] = dlat[:,2] + dlat[:,end] = dlat[:,end-1] + return dlon, dlat end @@ -84,7 +87,7 @@ function waterflows(Topo::GeoData, flowdir_fn= WhereTheWaterFlows.d8dir_feature; nin = zeros(Int8, ni) c = zeros(Int64, ni) - @views area[:,:,1], slen[:,:,1], dir[:,:,1], nout[:,:,1], nin[:,:,1], sinks, pits, c[:,:,1], bnds = waterflows(dem, cellarea, flowdir_fn; + 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) catchment_large = prune_catchments(c, minsize; val=0) From ed2b3c55d4b0be6de799bc65ba420b53befb55bb Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Wed, 3 Jul 2024 09:01:59 +0200 Subject: [PATCH 13/19] another fix --- src/WaterFlow.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/WaterFlow.jl b/src/WaterFlow.jl index 17ec243e..e05aaf6a 100644 --- a/src/WaterFlow.jl +++ b/src/WaterFlow.jl @@ -11,8 +11,8 @@ export waterflows Computes the spacing with central differences """ function spacing(lon,lat) - dlon = zero(size(lon.val)[1:2]) - dlat = zero(size(lat.val)[1:2]) + dlon = zeros(size(lon.val)[1:2]) + dlat = zeros(size(lat.val)[1:2]) @views dlon[2:end-1,:] = (lon.val[3:end,:,1] - lon.val[1:end-2,:,1])/2 dlon[1,:] = dlon[2,:] dlon[end,:] = dlon[end-1,:] From 004c293c5ec1669b08ef8c0f117be9fd59bdd628 Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Wed, 3 Jul 2024 09:09:20 +0200 Subject: [PATCH 14/19] add GLMakie heatmap and heatmap! multiple dispatch for surfaces --- ext/GLMakie_Visualisation.jl | 48 ++++++++++++++++++++++++++++++++++-- 1 file changed, 46 insertions(+), 2 deletions(-) diff --git a/ext/GLMakie_Visualisation.jl b/ext/GLMakie_Visualisation.jl index 517dff76..6d7f0110 100644 --- a/ext/GLMakie_Visualisation.jl +++ b/ext/GLMakie_Visualisation.jl @@ -3,7 +3,7 @@ module GLMakie_Visualisation using Statistics using GeophysicalModelGenerator: lonlatdepth_grid, GeoData, CartData, km, AbstractGeneralGrid -import GeophysicalModelGenerator: visualise +import GeophysicalModelGenerator: visualise, ustrip # We do not check `isdefined(Base, :get_extension)` as recommended since # Julia v1.9.0 does not load package extensions when their dependency is @@ -14,7 +14,9 @@ else using ..GLMakie end -export visualise +import GLMakie: heatmap!, heatmap + +export visualise, heatmap, heatmap! println("Loading GLMakie extensions for GMG") @@ -271,5 +273,47 @@ function visualise(Data::AbstractGeneralGrid; Topography=nothing, Topo_range=not return nothing end +""" + heatmap(x::GeoData, args...; field=:Topography, kwargs...) +heatmap for a 2D GeoData object (surface) +""" +function heatmap(x::GeoData, args...; field=:Topography, kwargs...) + @assert size(x.depth.val,3)==1 + + heatmap(x.lon.val[:,1], x.lat.val[1,:], ustrip.(x.fields[field][:,:,1]), args...; kwargs...) + +end + +""" + heatmap(x::CartData, args...; field=:Topography, kwargs...) +heatmap for a 2D CartData object (surface) +""" +function heatmap(x::CartData, args...; field=:Topography, kwargs...) + @assert size(x.z.val,3)==1 + + heatmap(x.x.val[:,1], x.y.val[1,:], ustrip.(x.fields[field][:,:,1]), args...; kwargs...) + +end + +""" + heatmap!(x::GeoData, args...; field=:Topography, kwargs...) +in-place heatmap for a 2D GeoData object (surface), +""" +function heatmap!(x::GeoData, args...; field=:Topography, kwargs...) + @assert size(x.z.val,3)==1 + + heatmap!(x.lon.val[:,1], x.lat.val[1,:], ustrip.(x.fields[field][:,:,1]), args...; kwargs...) + +end + +""" + heatmap!(x::CartData, args...; field=:Topography, kwargs...) +in-place heatmap for a 2D CartData object (surface) +""" +function heatmap!(x::CartData, args...; field=:Topography, kwargs...) + @assert size(x.z.val,3)==1 + heatmap!(x.x.val[:,1], x.y.val[1,:], ustrip.(x.fields[field][:,:,1]), args...; kwargs...) +end + end \ No newline at end of file From af532011ad45d536e765ac813fc53ab004a1523c Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Wed, 3 Jul 2024 10:32:43 +0200 Subject: [PATCH 15/19] add a few more fields to simplify plotting the largest catchment --- src/WaterFlow.jl | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/src/WaterFlow.jl b/src/WaterFlow.jl index e05aaf6a..e68de700 100644 --- a/src/WaterFlow.jl +++ b/src/WaterFlow.jl @@ -91,7 +91,13 @@ function waterflows(Topo::GeoData, flowdir_fn= WhereTheWaterFlows.d8dir_feature; feedback_fn=feedback_fn, drain_pits=drain_pits, bnd_as_sink=bnd_as_sink) catchment_large = prune_catchments(c, minsize; val=0) - - Topo_water = addfield(Topo,(;area, slen, dir, nout, nin, c, cellarea_m2, catchment_large)) + largest_catchment = catchment_large .== maximum(catchment_large) + largest_area = copy(area) + largest_area[.!largest_catchment] .= NaN + + log10_area = log10.(area) + log10_largest_area = log10.(largest_area) + + Topo_water = addfield(Topo,(;area, slen, dir, nout, nin, c, cellarea_m2, catchment_large, log10_area, largest_catchment, largest_area, log10_largest_area)) return Topo_water, sinks, pits, bnds end \ No newline at end of file From ac669710c30907d13ddc7ddaf5124b94e021af58 Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Wed, 3 Jul 2024 11:59:44 +0200 Subject: [PATCH 16/19] allow plotting an array computed on top of a surface --- ext/GLMakie_Visualisation.jl | 67 +++++++++++++++++++++++++++++++++++- 1 file changed, 66 insertions(+), 1 deletion(-) diff --git a/ext/GLMakie_Visualisation.jl b/ext/GLMakie_Visualisation.jl index 6d7f0110..5e07e998 100644 --- a/ext/GLMakie_Visualisation.jl +++ b/ext/GLMakie_Visualisation.jl @@ -279,30 +279,79 @@ heatmap for a 2D GeoData object (surface) """ function heatmap(x::GeoData, args...; field=:Topography, kwargs...) @assert size(x.depth.val,3)==1 + @assert hasfield(typeof(x.fields), field) heatmap(x.lon.val[:,1], x.lat.val[1,:], ustrip.(x.fields[field][:,:,1]), args...; kwargs...) end +""" + heatmap(x::GeoData, a::Array{_T,N}, args...; kwargs...) +in-place heatmap for a 2D GeoData object (surface) with an array `a`. +""" +function heatmap(x::GeoData, a::Array{_T,N}, args...; kwargs...) where{_T,N} + @assert size(x.depth.val,3)==1 + + if N==3 + heatmap(x.lon.val[:,1], x.lat.val[1,:], ustrip.(a[:,:,1]), args...; kwargs...) + elseif N==2 + heatmap(x.lon.val[:,1], x.lat.val[1,:], ustrip.(a), args...; kwargs...) + end + +end + """ heatmap(x::CartData, args...; field=:Topography, kwargs...) heatmap for a 2D CartData object (surface) """ function heatmap(x::CartData, args...; field=:Topography, kwargs...) @assert size(x.z.val,3)==1 + @assert hasfield(typeof(x.fields), field) heatmap(x.x.val[:,1], x.y.val[1,:], ustrip.(x.fields[field][:,:,1]), args...; kwargs...) end +""" + heatmap(x::CartData, a::Array{_T,N}, args...; kwargs...) +in-place heatmap for a 2D CartData object (surface) with an array `a`. +""" +function heatmap(x::CartData, a::Array{_T,N}, args...; kwargs...) where{_T,N} + @assert size(x.z.val,3)==1 + + if N==3 + heatmap(x.x.val[:,1], x.y.val[1,:], ustrip.(a[:,:,1]), args...; kwargs...) + elseif N==2 + heatmap(x.x.val[:,1], x.y.val[1,:], ustrip.(a), args...; kwargs...) + end + +end + + """ heatmap!(x::GeoData, args...; field=:Topography, kwargs...) in-place heatmap for a 2D GeoData object (surface), """ function heatmap!(x::GeoData, args...; field=:Topography, kwargs...) @assert size(x.z.val,3)==1 + @assert hasfield(typeof(x.fields), field) - heatmap!(x.lon.val[:,1], x.lat.val[1,:], ustrip.(x.fields[field][:,:,1]), args...; kwargs...) + heatmap(x.lon.val[:,1], x.lat.val[1,:], ustrip.(x.fields[field][:,:,1]), args...; kwargs...) + +end + +""" + heatmap!(x::GeoData, a::Array{_T,N}, args...; kwargs...) +in-place heatmap for a 2D GeoData object (surface) with an array `a`. +""" +function heatmap!(x::GeoData, a::Array{_T,N}, args...; kwargs...) where{_T,N} + @assert size(x.z.val,3)==1 + + if N==3 + heatmap!(x.lon.val[:,1], x.lat.val[1,:], ustrip.(a[:,:,1]), args...; kwargs...) + elseif N==2 + heatmap!(x.lon.val[:,1], x.lat.val[1,:], ustrip.(a), args...; kwargs...) + end end @@ -312,8 +361,24 @@ in-place heatmap for a 2D CartData object (surface) """ function heatmap!(x::CartData, args...; field=:Topography, kwargs...) @assert size(x.z.val,3)==1 + @assert hasfield(typeof(x.fields), field) + heatmap!(x.x.val[:,1], x.y.val[1,:], ustrip.(x.fields[field][:,:,1]), args...; kwargs...) end +""" + heatmap!(x::CartData, a::Array{_T,N}, args...; kwargs...) +in-place heatmap for a 2D CartData object (surface) with an array `a`. +""" +function heatmap!(x::CartData, a::Array{_T,N}, args...; kwargs...) where{_T,N} + @assert size(x.z.val,3)==1 + + if N==3 + heatmap!(x.x.val[:,1], x.y.val[1,:], ustrip.(a[:,:,1]), args...; kwargs...) + elseif N==2 + heatmap!(x.x.val[:,1], x.y.val[1,:], ustrip.(a[:,:]), args...; kwargs...) + end + +end end \ No newline at end of file From 64f4b9a7d474608c2dd3b8525797e1d2778360db Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Wed, 3 Jul 2024 12:00:08 +0200 Subject: [PATCH 17/19] remove log10 of largest area --- src/WaterFlow.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/WaterFlow.jl b/src/WaterFlow.jl index e68de700..37830828 100644 --- a/src/WaterFlow.jl +++ b/src/WaterFlow.jl @@ -96,8 +96,7 @@ function waterflows(Topo::GeoData, flowdir_fn= WhereTheWaterFlows.d8dir_feature; largest_area[.!largest_catchment] .= NaN log10_area = log10.(area) - log10_largest_area = log10.(largest_area) - Topo_water = addfield(Topo,(;area, slen, dir, nout, nin, c, cellarea_m2, catchment_large, log10_area, largest_catchment, largest_area, log10_largest_area)) + Topo_water = addfield(Topo,(;area, slen, dir, nout, nin, c, cellarea_m2, catchment_large, log10_area, largest_catchment, largest_area)) return Topo_water, sinks, pits, bnds end \ No newline at end of file From 1821f8737e0baf13e876237bf878e75c7136dc31 Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Wed, 3 Jul 2024 17:31:12 +0200 Subject: [PATCH 18/19] typo --- src/WaterFlow.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/WaterFlow.jl b/src/WaterFlow.jl index 37830828..4ed6e427 100644 --- a/src/WaterFlow.jl +++ b/src/WaterFlow.jl @@ -33,7 +33,7 @@ function cell_area(Topo::GeoData) Topo_cart = convert2CartData(Topo, proj) dx, dy = spacing(Topo_cart.x, Topo_cart.y) - area_m2 = dx.*dy*1e6 + area_m2 = dx*dy*1e6 return area_m2 end From 087a27352205c14811fceb332c10083d7afe502a Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Wed, 3 Jul 2024 18:05:15 +0200 Subject: [PATCH 19/19] undo last change. dot is required --- src/WaterFlow.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/WaterFlow.jl b/src/WaterFlow.jl index 4ed6e427..37830828 100644 --- a/src/WaterFlow.jl +++ b/src/WaterFlow.jl @@ -33,7 +33,7 @@ function cell_area(Topo::GeoData) Topo_cart = convert2CartData(Topo, proj) dx, dy = spacing(Topo_cart.x, Topo_cart.y) - area_m2 = dx*dy*1e6 + area_m2 = dx.*dy*1e6 return area_m2 end