Completely unroll nested function in flow routing algorithm

This commit is contained in:
Gaël C 2024-01-22 00:55:51 +01:00
parent b54f2c4546
commit 4bce5fab77

View File

@ -93,27 +93,6 @@ local function flow_routing(dem, dirs, lakes) -- 'dirs' and 'lakes' are optional
local links = {}
local basin_links
-- Function to analyse a link between two nodes
local function add_link(i1, i2, b1, isY)
-- i1, i2: coordinates of two nodes
-- b1: basin that contains i1
-- isY: whether the link is in Y direction
-- Note that basin number #0 represents the outside of the map; or if the coordinate is inside the map, means that the basin number is uninitialized.
local b2 = i2 == 0 and 0 or basin_id[i2]
local elev = i2 == 0 and dem[i1] or mmax(dem[i1], dem[i2]) -- Elevation of the highest of the two sides of the link (or only i1 if b2 is map outside)
local l2 = basin_links[b2]
if not l2 then
l2 = {b2, b1, elev=elev, i=mmax(i1,i2), is_y=isY}
basin_links[b2] = l2
elseif l2.elev > elev then -- If this link is lower than the lowest registered link between these two basins, register it as the new lowest pass
l2.elev = elev
l2.i = mmax(i1,i2)
l2.is_y = isY
l2[1] = b2
l2[2] = b1
end
end
for i=1, X*Y do
basin_id[i] = 0
end
@ -139,10 +118,32 @@ local function flow_routing(dem, dirs, lakes) -- 'dirs' and 'lakes' are optional
-- If no river is coming from the East, we might be at the limit of two basins, thus we need to test adjacency.
elseif i%X > 0 then
if basin_id[i+1] ~= ib and basin_id[i+1] ~= 0 then
add_link(i, i+1, ib, false)
local b2 = basin_id[i+1]
local elev = mmax(dem[i], dem[i+1]) -- Elevation of the highest of the two sides of the link (or only i1 if b2 is map outside)
local l2 = basin_links[b2]
if not l2 then
l2 = {b2, ib, elev=elev, i=i+1, is_y=false}
basin_links[b2] = l2
elseif l2.elev > elev then -- If this link is lower than the lowest registered link between these two basins, register it as the new lowest pass
l2.elev = elev
l2.i = i+1
l2.is_y = false
l2[1] = b2
l2[2] = ib
end
end
else -- If the eastern neighbour is outside the map
add_link(i, 0, ib, false)
local l2 = basin_links[0]
if not l2 then
l2 = {0, ib, elev=dem[i], i=i, is_y=false}
basin_links[0] = l2
elseif l2.elev > dem[i] then
l2.elev = dem[i]
l2.i = i
l2.is_y = false
l2[1] = 0
l2[2] = ib
end
end
if d >= 4 then -- River coming from the South
@ -151,10 +152,32 @@ local function flow_routing(dem, dirs, lakes) -- 'dirs' and 'lakes' are optional
singular[cur] = i+X
elseif i <= X*(Y-1) then
if basin_id[i+X] ~= ib and basin_id[i+X] ~= 0 then
add_link(i, i+X, ib, true)
local b2 = basin_id[i+X]
local elev = mmax(dem[i], dem[i+X])
local l2 = basin_links[b2]
if not l2 then
l2 = {b2, ib, elev=elev, i=i+X, is_y=true}
basin_links[b2] = l2
elseif l2.elev > elev then
l2.elev = elev
l2.i = i+X
l2.is_y = true
l2[1] = b2
l2[2] = ib
end
end
else
add_link(i, 0, ib, true)
local l2 = basin_links[0]
if not l2 then
l2 = {0, ib, elev=dem[i], i=i, is_y=true}
basin_links[0] = l2
elseif l2.elev > dem[i] then
l2.elev = dem[i]
l2.i = i
l2.is_y = true
l2[1] = 0
l2[2] = ib
end
end
if d >= 2 then -- River coming from the West
@ -163,10 +186,32 @@ local function flow_routing(dem, dirs, lakes) -- 'dirs' and 'lakes' are optional
singular[cur] = i-1
elseif i%X ~= 1 then
if basin_id[i-1] ~= ib and basin_id[i-1] ~= 0 then
add_link(i, i-1, ib, false)
local b2 = basin_id[i-1]
local elev = mmax(dem[i], dem[i-1])
local l2 = basin_links[b2]
if not l2 then
l2 = {b2, ib, elev=elev, i=i, is_y=false}
basin_links[b2] = l2
elseif l2.elev > elev then
l2.elev = elev
l2.i = i
l2.is_y = false
l2[1] = b2
l2[2] = ib
end
end
else
add_link(i, 0, ib, false)
local l2 = basin_links[0]
if not l2 then
l2 = {0, ib, elev=dem[i], i=i, is_y=false}
basin_links[0] = l2
elseif l2.elev > dem[i] then
l2.elev = dem[i]
l2.i = i
l2.is_y = false
l2[1] = 0
l2[2] = ib
end
end
if d >= 1 then -- River coming from the North
@ -174,10 +219,32 @@ local function flow_routing(dem, dirs, lakes) -- 'dirs' and 'lakes' are optional
singular[cur] = i-X
elseif i > X then
if basin_id[i-X] ~= ib and basin_id[i-X] ~= 0 then
add_link(i, i-X, ib, true)
local b2 = basin_id[i-X]
local elev = mmax(dem[i], dem[i-X])
local l2 = basin_links[b2]
if not l2 then
l2 = {b2, ib, elev=elev, i=i, is_y=true}
basin_links[b2] = l2
elseif l2.elev > elev then
l2.elev = elev
l2.i = i
l2.is_y = true
l2[1] = b2
l2[2] = ib
end
end
else
add_link(i, 0, ib, true)
local l2 = basin_links[0]
if not l2 then
l2 = {0, ib, elev=dem[i], i=i, is_y=true}
basin_links[0] = l2
elseif l2.elev > dem[i] then
l2.elev = dem[i]
l2.i = i
l2.is_y = true
l2[1] = 0
l2[2] = ib
end
end
end
dirs2 = nil