Handling flat areas

This is a simple demonstration of how large flat surfaces (like large water bodies) can be recognized, and straightened or excluded from analysis if desired ^[1].

Imporing packages and loading the surface.

using SurfaceWaterIntegratedModeling
import CairoMakie, Images # for visualization and loading of textures
import ColorSchemes
using Pkg.Artifacts

The package with SWIM testdata is provided as a Julia artifact, which can be accessed using the function datapath_testdata. We subsequently load and display the synthetic grid.

datapath = joinpath(datapath_testdata(), "data", "small")
grid = loadgrid(joinpath(datapath, "bay.txt"));

# for ease of use, we create our own label of key colors found in the
# colorsheme `:Paired_12` used below
cmap = Dict(:blue => 2, :green => 4, :red => 6, :orange => 8,
            :lilac => 10, :bright => 11);

# Also define a preset viewpoint
view1 = (CairoMakie.Vec(1059, 783, 622), CairoMakie.Vec(280, 311, 58), 0.68);

 Downloading artifact: swim_testdata
 Downloading artifact: swim_testdata

The terrain used for demonstration is a varied surface that includes hills, a mountainside, a road, and ocean areas.

tex = fill(cmap[:green], size(grid))
sf, fig, sc = plotgrid(grid, texture = tex,
                       colormap=ColorSchemes.:Paired_12,
                       colorrange=(1,12))
fig
set_camerapos(fig, sc, view1...)

Identifying flat areas and running spill analysis

While the ocean part is clearly visible on the surface, its cells in the terrain grid have not been identified. For this purpose, the identify_flat_areas function can be used:

rel_tol = 1e-3 # relative tolerance for an angle between two cells to be considered zero
max_cluster_size = 1 # threshold for excluding too small areas
isflat = identify_flat_areas(grid, rel_tol, max_cluster_size)

tex[isflat] .= cmap[:blue]
drape_surface(sf, tex)
set_camerapos(fig, sc, view1...)

We can see that the ocean part was correctly identified as flat, but so were many parts that ought not be included. We can increase the cluster size limit to filter these parts out.

max_cluster_size=100
isflat = identify_flat_areas(grid, rel_tol, max_cluster_size)

tex[isflat] .= cmap[:blue]
tex[.!isflat] .= cmap[:green]
drape_surface(sf, tex)
set_camerapos(fig, sc, view1...)

With this value for the threshold, only the ocean remains identified as flat.

For the surface used in this example, the ocean surface is already exactly flat. However, it is sometimes the case with terrain input data that conceptually flat regions still contains small irregularities that may cause noisy output from the spill analysis. One way to handle this, is to use the flatten_grid! function, which ensures that flat parts are exactly that. Even if not strictly necessary for our sample terrain, we demonstrate the use below:

flatten_grid!(grid, isflat, :min)

Another way to exclude flat areas from the analysis is to declare them as "sinks". For oceans, this makes conceptual sense, since they can be considered traps with "infinite" capacity. For other large flat areas, e.g. parking lots, this approach does not work so well, and flatten_grid! should be used. For our terrain, declaring the flat areas as sinks should be fine. We pass along isflat as the sink map in our call to spillanalysis below:

tstruct = spillanalysis(grid, sinks=isflat)

TrapStructure{Float64}([114.64 113.62 … 30.184 29.403; 114.93 114.0 … 30.002 29.197; … ; 86.445 86.805 … 32.437 32.877; 86.009 86.429 … 32.419 32.88], Int8[3 3 … 3 3; 3 7 … 5 3; … ; 1 6 … 2 2; 1 2 … 2 2], [-1 -1 … -56418 -56418; -1 -1 … -56515 -56612; … ; -52 -52 … -56417 -56417; -52 -52 … -56514 -56514], Spillpoint[Spillpoint(134, 83048, 83049, 4.3057), Spillpoint(58, 51885, 51884, 37.931), Spillpoint(40, 22595, 22594, 75.583), Spillpoint(134, 80346, 80798, 4.7259), Spillpoint(4, 25778, 26228, 70.685), Spillpoint(137, 82168, 82169, 3.2237), Spillpoint(154, 84888, 85340, 2.555), Spillpoint(161, 87159, 87609, 1.7134), Spillpoint(11, 19498, 19499, 74.053), Spillpoint(44, 19054, 19504, 73.804)  …  Spillpoint(-9145, 83628, 83629, 5.1506), Spillpoint(148, 83940, 83939, 4.4768), Spillpoint(4, 80798, 80346, 4.7259), Spillpoint(121, 81239, 80789, 4.8328), Spillpoint(127, 77697, 77245, 4.9052), Spillpoint(110, 83031, 83481, 5.3576), Spillpoint(145, 83929, 83928, 5.9195), Spillpoint(149, 85267, 85717, 7.0289), Spillpoint(124, 87064, 87063, 7.8472), Spillpoint(-9144, 87947, 87946, 8.4201)], [0.15759999999999952, 0.01899999999999835, 0.04399999999999693, 0.04630000000000045, 0.04200000000000159, 1.1318000000000001, 0.1871000000000005, 1.2660000000000002, 0.2009999999999934, 0.0  …  2560.1657999999984, 1289.7652000000007, 1621.6218, 1780.6920000000011, 1894.2151000000003, 2713.6629999999977, 3961.6346999999982, 7065.913900000009, 9881.127499999986, 12192.254699999985], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  2479.6928, 1277.9733000000008, 1289.7728000000006, 1621.6680999999999, 1780.748600000001, 1894.2469000000003, 2713.989299999998, 3964.7893999999983, 7066.0303000000085, 9893.510299999985], [[83046, 83047, 83048], [51885], [22595, 22596], [80345, 80346], [25778], [81711, 81712, 81713, 81714, 81715, 81716, 82163, 82164, 82165, 82166, 82167, 82168, 82618, 82619], [83985, 84436, 84437, 84888], [86258, 86259, 86260, 86261, 86708, 86709, 86710, 86711, 87159, 87160, 87161, 87162], [19498, 19948, 19949], [19054]  …  [70542, 70544, 70545, 70546, 70991, 70992, 70993, 70994, 70995, 70996  …  112963, 113414, 113865, 113866, 114318, 114770, 115222, 115673, 115674, 116125], [78143, 78593, 78594, 79044, 79045, 79495, 79496, 79945, 79946, 79947  …  97080, 97081, 97082, 97526, 97527, 97528, 97529, 97530, 97531, 97532], [76790, 77241, 77691, 77692, 78141, 78142, 78143, 78591, 78592, 78593  …  97529, 97530, 97531, 97532, 97533, 97979, 97980, 97981, 97982, 97983], [76790, 77241, 77242, 77691, 77692, 78141, 78142, 78143, 78144, 78147  …  97533, 97978, 97979, 97980, 97981, 97982, 97983, 97984, 98429, 98430], [76339, 76790, 77240, 77241, 77242, 77691, 77692, 77693, 77697, 78140  …  97978, 97979, 97980, 97981, 97982, 97983, 97984, 98429, 98430, 98880], [75887, 75888, 75891, 75892, 75893, 75894, 76337, 76338, 76339, 76340  …  98879, 98880, 98881, 98882, 98883, 98884, 98885, 98886, 99331, 99332], [74988, 74989, 74990, 74991, 74992, 74993, 74994, 75436, 75437, 75438  …  99782, 99783, 99784, 99785, 99786, 100232, 100233, 100234, 100235, 100684], [73189, 73190, 73191, 73192, 73639, 73640, 73641, 73642, 73643, 74083  …  101589, 102035, 102036, 102037, 102486, 102487, 102488, 102938, 102939, 103389], [71838, 72287, 72288, 72289, 72290, 72735, 72736, 72737, 72738, 72739  …  102939, 102940, 103388, 103389, 103390, 103391, 103839, 103840, 103841, 104291], [70936, 70937, 71384, 71385, 71386, 71387, 71388, 71389, 71834, 71835  …  103849, 104289, 104290, 104291, 104292, 104741, 104742, 105192, 105193, 105644]], [[1], [2], [3], [4], [5], [6], [7], [8], [9], [10]  …  [207, 187, 170, 180, 177, 176, 140, 143, 155, 157  …  119, 135, 174, 203, 198, 200, 194, 192, 208, 223], [1, 160, 137, 128, 167, 7, 146, 161, 8, 154, 166, 162, 159, 6, 142, 134, 11, 147, 151], [1, 160, 137, 128, 167, 7, 146, 161, 8, 154, 166, 162, 159, 6, 142, 134, 11, 147, 151, 148], [1, 160, 137, 128, 167, 7, 146, 161, 8, 154  …  162, 159, 6, 142, 134, 11, 147, 151, 148, 4], [1, 160, 137, 128, 167, 7, 146, 161, 8, 154  …  159, 6, 142, 134, 11, 147, 151, 148, 4, 121], [1, 160, 137, 128, 167, 7, 146, 161, 8, 154  …  6, 142, 134, 11, 147, 151, 148, 4, 121, 127], [110, 1, 160, 137, 128, 167, 7, 146, 161, 8  …  6, 142, 134, 11, 147, 151, 148, 4, 121, 127], [110, 1, 160, 137, 128, 167, 7, 146, 161, 8  …  134, 11, 147, 151, 148, 4, 121, 127, 141, 145], [110, 1, 160, 137, 128, 167, 7, 146, 161, 8  …  11, 147, 151, 148, 4, 121, 127, 141, 145, 149], [110, 1, 160, 137, 128, 167, 7, 146, 161, 8  …  148, 4, 121, 127, 141, 145, 149, 153, 158, 124]], [[1, 1108, 1110, 1111, 1112, 1113, 1114, 1115, 1116, 1117, 1118], [2], [3], [4, 1112, 1113, 1114, 1115, 1116, 1117, 1118], [5], [6, 1098, 1101, 1103, 1105, 1106, 1108, 1110, 1111, 1112, 1113, 1114, 1115, 1116, 1117, 1118], [7, 1057, 1071, 1076, 1085, 1087, 1094, 1097, 1098, 1101  …  1108, 1110, 1111, 1112, 1113, 1114, 1115, 1116, 1117, 1118], [8, 897, 1002, 1038, 1057, 1071, 1076, 1085, 1087, 1094  …  1108, 1110, 1111, 1112, 1113, 1114, 1115, 1116, 1117, 1118], [9], [10]  …  [886], [887], [888], [889], [890, 991], [891, 993], [892, 993], [893], [894], [895]], Graphs.SimpleGraphs.SimpleDiGraph{Int64}(446, [[1108], Int64[], Int64[], [1112], Int64[], [1098], [1057], [897], Int64[], Int64[]  …  Int64[], [1111], [1112], [1113], [1114], [1115], [1116], [1117], [1118], Int64[]], [Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[]  …  [223, 1107], [151, 1108], [148, 1110], [4, 1111], [121, 1112], [127, 1113], [110, 1114], [911, 1115], [149, 1116], [1029, 1117]]), nothing, [(448, 113), (451, 113), (442, 114), (443, 114), (444, 114), (447, 114), (448, 114), (449, 114), (450, 114), (451, 114)  …  (161, 601), (162, 601), (163, 601), (164, 601), (165, 601), (166, 601), (167, 601), (168, 601), (169, 601), (170, 601)])

The spill anaysis was able to identify 895 spill regions and 1118 traps:

print("Number of spill regions identified: ", numregions(tstruct), '\n')
print("Number of traps identified: ", numtraps(tstruct), '\n')

Number of spill regions identified: 895
Number of traps identified: 1118

We can visualize the traps, spill regions, rivers and oceans on the terrain surface by creating a texture with show_region_selection:

tex = show_region_selection(tstruct,
                            region_color=cmap[:green],
                            river_color=cmap[:blue]-1,
                            trap_color=cmap[:blue]);

Terrain that spills directly to the ocean, or out of the domain is assigned a light green color, to distinguish it from terrain that spills to identified traps, while the ocean is set to blue:

tex[tex.==0] .= cmap[:green]-1;
tex[isflat] .= cmap[:blue];

drape_surface(sf, tex)
set_camerapos(fig, sc, view1...)

Adding more sinks

The ridge in front of the large lake in the upper left part of the plot is actually a railroad. The regions immediately in front of it and behind it should be considered drained, so water does not massively accumulated like the above plot indicates. By adding a drain near its lowest point, the lake goes away:

sinks = [(119,193), (180, 193)]; # two drain locations
for s in sinks
    # add the sinks to the map that we used to represent sinks (and flat areas)
    isflat[s...] = true;
end

# run the spill analysis again
tstruct = spillanalysis(grid, sinks=isflat)

TrapStructure{Float64}([114.64 113.62 … 30.184 29.403; 114.93 114.0 … 30.002 29.197; … ; 86.445 86.805 … 32.437 32.877; 86.009 86.429 … 32.419 32.88], Int8[3 3 … 3 3; 3 7 … 5 3; … ; 1 6 … 2 2; 1 2 … 2 2], [-1 -1 … -56420 -56420; -1 -1 … -56517 -56614; … ; -52 -52 … -56419 -56419; -52 -52 … -56516 -56516], Spillpoint[Spillpoint(134, 83048, 83049, 4.3057), Spillpoint(58, 51885, 51884, 37.931), Spillpoint(40, 22595, 22594, 75.583), Spillpoint(134, 80346, 80798, 4.7259), Spillpoint(4, 25778, 26228, 70.685), Spillpoint(137, 82168, 82169, 3.2237), Spillpoint(154, 84888, 85340, 2.555), Spillpoint(-9626, 86710, 86711, 1.569), Spillpoint(11, 19498, 19499, 74.053), Spillpoint(44, 19054, 19504, 73.804)  …  Spillpoint(264, 129688, 129236, 5.1194), Spillpoint(338, 138570, 139022, 19.647), Spillpoint(101, 58108, 58558, 6.7165), Spillpoint(876, 261659, 262109, 2.5455), Spillpoint(245, 137214, 136762, 19.654), Spillpoint(260, 131499, 131049, 5.1203), Spillpoint(792, 249963, 249513, 2.5764), Spillpoint(-183, 52280, 52731, 6.721), Spillpoint(879, 263907, 264357, 2.5931), Spillpoint(873, 264806, 265256, 2.6536)], [0.15759999999999952, 0.01899999999999835, 0.04399999999999693, 0.04630000000000045, 0.04200000000000159, 1.1318000000000001, 0.1871000000000005, 0.13719999999999977, 0.2009999999999934, 0.0  …  4.020999999999961, 23.639999999999706, 120.73149999999991, 11.514600000000007, 25.008999999999986, 4.166600000000044, 13.898600000000002, 122.85570000000013, 15.211200000000012, 20.7984], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  3.991800000000029, 22.356000000000098, 117.4655999999999, 8.557000000000006, 23.645999999999706, 4.03339999999996, 11.785000000000007, 120.9713999999999, 13.929200000000002, 15.247000000000012], [[83046, 83047, 83048], [51885], [22595, 22596], [80345, 80346], [25778], [81711, 81712, 81713, 81714, 81715, 81716, 82163, 82164, 82165, 82166, 82167, 82168, 82618, 82619], [83985, 84436, 84437, 84888], [86709, 86710, 87161, 87162], [19498, 19948, 19949], [19054]  …  [129688, 130140, 130141, 130593, 130594, 131044, 131045, 131046, 131047, 131495  …  137823, 137824, 137825, 137826, 138274, 138275, 138276, 138277, 138278, 138729], [137215, 137216, 137664, 137665, 137667, 137668, 137669, 138112, 138113, 138114  …  144429, 144430, 144873, 144874, 144875, 144876, 144877, 144878, 144879, 144880], [45072, 45073, 45074, 45075, 45076, 45077, 45078, 45081, 45498, 45499  …  57214, 57659, 57660, 57661, 57662, 57663, 58108, 58109, 58110, 58111], [250412, 250413, 250862, 250863, 251312, 251313, 251762, 251763, 251764, 252212  …  259410, 259860, 260309, 260310, 260759, 260760, 261209, 261210, 261659, 261660], [137214, 137215, 137216, 137664, 137665, 137666, 137667, 137668, 137669, 138112  …  144430, 144873, 144874, 144875, 144876, 144877, 144878, 144879, 144880, 145329], [129236, 129688, 130140, 130141, 130593, 130594, 131044, 131045, 131046, 131047  …  137823, 137824, 137825, 137826, 138274, 138275, 138276, 138277, 138278, 138729], [249963, 250412, 250413, 250862, 250863, 251312, 251313, 251314, 251762, 251763  …  261659, 261660, 262109, 262558, 262559, 263008, 263009, 263458, 263459, 263908], [45072, 45073, 45074, 45075, 45076, 45077, 45078, 45079, 45081, 45498  …  58109, 58110, 58111, 58555, 58556, 58557, 58558, 59004, 59005, 59006], [248163, 248613, 249063, 249513, 249963, 250412, 250413, 250862, 250863, 251312  …  261660, 262109, 262558, 262559, 263008, 263009, 263458, 263459, 263907, 263908], [246813, 247263, 247713, 248163, 248613, 249063, 249512, 249513, 249962, 249963  …  263009, 263458, 263459, 263907, 263908, 263909, 264356, 264357, 264358, 264806]], [[1], [2], [3], [4], [5], [6], [7], [8], [9], [10]  …  [270, 293, 272, 284, 304, 307, 312, 311, 317, 324, 305], [325, 329, 369, 345, 342, 328, 332, 357, 364, 349], [99, 96, 85, 65, 59, 28, 87, 90, 36, 86, 106, 105, 104, 102, 103, 112], [870, 855, 864, 837, 847, 834, 825, 867, 860], [325, 329, 369, 345, 342, 328, 332, 357, 364, 349, 338], [270, 293, 272, 284, 304, 307, 312, 311, 317, 324, 305, 264], [872, 876, 870, 855, 864, 837, 847, 834, 825, 867, 860], [99, 96, 85, 65, 59, 28, 87, 90, 36, 86, 106, 105, 104, 102, 103, 112, 100, 101], [872, 876, 870, 855, 864, 837, 847, 834, 825, 867, 860, 792], [872, 876, 870, 855, 864, 837, 847, 834, 825, 867, 860, 792, 879]], [[1], [2], [3], [4], [5], [6], [7], [8], [9], [10]  …  [886], [887], [888], [889], [890, 990], [891, 992], [892, 992], [893], [894], [895]], Graphs.SimpleGraphs.SimpleDiGraph{Int64}(366, [Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[]  …  [1074], [1073], [1076], [1075], Int64[], Int64[], [1077], Int64[], [1078], Int64[]], [Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[], Int64[]  …  [270, 1065], [933, 1067], [99, 1066], [860, 1068], [338, 1070], [264, 1069], [989, 1072], [905, 1071], [792, 1075], [879, 1077]]), nothing, [(448, 113), (451, 113), (442, 114), (443, 114), (444, 114), (447, 114), (448, 114), (449, 114), (450, 114), (451, 114)  …  (161, 601), (162, 601), (163, 601), (164, 601), (165, 601), (166, 601), (167, 601), (168, 601), (169, 601), (170, 601)])

The result of this analysis no longer indicate a large lake in this region:

tex = show_region_selection(tstruct,
                            region_color=4,
                            river_color=2,
                            trap_color=2);
tex[tex.==0] .= cmap[:green]-1
tex[isflat] .= cmap[:blue]

drape_surface(sf, tex)
set_camerapos(fig, sc, view1...)

This page was generated using Literate.jl.

1The data used in this example was originally obtained from Kartverket (the Norwegian Mapping Authority) under the Creative Commons Attribution 4.0 International (CC BY 4.0) license.