diff --git a/terrainlib_lua/rivermapper.lua b/terrainlib_lua/rivermapper.lua
index bdeca47..a2d8663 100644
--- a/terrainlib_lua/rivermapper.lua
+++ b/terrainlib_lua/rivermapper.lua
@@ -98,30 +98,20 @@ local function flow_routing(dem, dirs, lakes) -- 'dirs' and 'lakes' are optional
 		-- i1, i2: coordinates of two nodes
 		-- b1: basin that contains i1
 		-- isY: whether the link is in Y direction
-		local b2
 		-- 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.
-		if i2 == 0 then -- If outside the map
-			b2 = 0
-		else
-			b2 = basin_id[i2]
-			if b2 == 0 then -- If basin of i2 is not already computed, skip
-				return
-			end
+		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 = {}
+			basin_links[b2] = l2
 		end
-		if b2 ~= b1 then -- If these two nodes don't belong to the same basin, we have found a link between two adjacent basins
-			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 = {}
-				basin_links[b2] = l2
-			end
-			if not l2.elev or 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
+		if not l2.elev or 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
 
@@ -149,7 +139,9 @@ local function flow_routing(dem, dirs, lakes) -- 'dirs' and 'lakes' are optional
 				singular[cur] = i+1
 			-- 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
-				add_link(i, i+1, ib, false)
+				if basin_id[i+1] ~= ib and basin_id[i+1] ~= 0 then
+					add_link(i, i+1, ib, false)
+				end
 			else -- If the eastern neighbour is outside the map
 				add_link(i, 0, ib, false)
 			end
@@ -159,7 +151,9 @@ local function flow_routing(dem, dirs, lakes) -- 'dirs' and 'lakes' are optional
 				cur = cur + 1
 				singular[cur] = i+X
 			elseif i <= X*(Y-1) then
-				add_link(i, i+X, ib, true)
+				if basin_id[i+X] ~= ib and basin_id[i+X] ~= 0 then
+					add_link(i, i+X, ib, true)
+				end
 			else
 				add_link(i, 0, ib, true)
 			end
@@ -169,7 +163,9 @@ local function flow_routing(dem, dirs, lakes) -- 'dirs' and 'lakes' are optional
 				cur = cur + 1
 				singular[cur] = i-1
 			elseif i%X ~= 1 then
-				add_link(i, i-1, ib, false)
+				if basin_id[i-1] ~= ib and basin_id[i-1] ~= 0 then
+					add_link(i, i-1, ib, false)
+				end
 			else
 				add_link(i, 0, ib, false)
 			end
@@ -178,7 +174,9 @@ local function flow_routing(dem, dirs, lakes) -- 'dirs' and 'lakes' are optional
 				cur = cur + 1
 				singular[cur] = i-X
 			elseif i > X then
-				add_link(i, i-X, ib, true)
+				if basin_id[i-X] ~= ib and basin_id[i-X] ~= 0 then
+					add_link(i, i-X, ib, true)
+				end
 			else
 				add_link(i, 0, ib, true)
 			end