commit 0bf351b2f67ea6d9fcd0e90b90d043209b61b7d1
Author: Gael-de-Sailly <gael.chretien@akeonet.com>
Date:   Thu Apr 9 21:13:38 2020 +0200

    Initial commit: working example using a basis of Simplex noise and implementing river flowing, lakes, and erosion

diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..e257aae
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,7 @@
+__pycache__/
+dem
+lakes
+links
+rivers
+size
+unused/
diff --git a/erosion.py b/erosion.py
new file mode 100644
index 0000000..ce113f5
--- /dev/null
+++ b/erosion.py
@@ -0,0 +1,77 @@
+import numpy as np
+import scipy.ndimage as im
+import rivermapper as rm
+
+def advection(dem, dirs, rivers, time, K=1, m=0.5, sea_level=0):
+    dirs = dirs.copy()
+    dirs[0,:] = 0
+    dirs[-1,:] = 0
+    dirs[:,0] = 0
+    dirs[:,-1] = 0
+
+    adv_time = 1 / (K*rivers**m)
+    dem = np.maximum(dem, sea_level)
+    dem_new = np.zeros(dem.shape)
+
+    for y in range(dirs.shape[0]):
+        for x in range(dirs.shape[1]):
+            x0, y0 = x, y
+            x1, y1 = x, y
+            remaining = time
+            while True:
+                flow_dir = dirs[y0,x0]
+                if flow_dir == 0:
+                    remaining = 0
+                    break
+                elif flow_dir == 1:
+                    y1 += 1
+                elif flow_dir == 2:
+                    x1 += 1
+                elif flow_dir == 3:
+                    y1 -= 1
+                elif flow_dir == 4:
+                    x1 -= 1
+
+                if remaining <= adv_time[y0,x0]:
+                    break
+                remaining -= adv_time[y0,x0]
+                x0, y0 = x1, y1
+
+            c = remaining / adv_time[y0,x0]
+            dem_new[y,x] = c*dem[y1,x1] + (1-c)*dem[y0,x0]
+
+    return np.minimum(dem, dem_new)
+
+def diffusion(dem, time, d=1):
+    radius = d * time**.5
+    return im.gaussian_filter(dem, radius, mode='reflect')
+
+class EvolutionModel:
+    def __init__(self, dem, K=1, m=0.5, d=1, sea_level=0, flow=False):
+        self.dem = dem
+        #self.bedrock = dem
+        self.K = K
+        self.m = m
+        self.d = d
+        self.sea_level = sea_level
+        #self.ref_isostasy = im.gaussian_filter(dem, radius, mode='reflect')
+        if flow:
+            self.calculate_flow()
+        else:
+            self.lakes = dem
+            self.dirs = np.zeros(dem.shape, dtype='u1')
+            self.rivers = np.zeros(dem.shape, dtype='u4')
+            self.flow_uptodate = False
+
+    def calculate_flow(self):
+        self.dirs, self.lakes, self.rivers = rm.flow(self.dem)
+        self.flow_uptodate = True
+
+    def advection(self, time):
+        dem = advection(self.lakes, self.dirs, self.rivers, time, K=self.K, m=self.m, sea_level=self.sea_level)
+        self.dem = np.minimum(dem, self.dem)
+        self.flow_uptodate = False
+
+    def diffusion(self, time):
+        self.dem = diffusion(self.dem, time, d=self.d)
+        self.flow_uptodate = False
diff --git a/init.lua b/init.lua
new file mode 100644
index 0000000..8f934b0
--- /dev/null
+++ b/init.lua
@@ -0,0 +1,178 @@
+local modpath = minetest.get_modpath(minetest.get_current_modname()) .. '/'
+local worldpath = minetest.get_worldpath() .. '/'
+local load_map = dofile(modpath .. 'load.lua')
+
+local function copy_if_needed(filename)
+	local wfilename = worldpath..filename
+	local wfile = io.open(wfilename, 'r')
+	if wfile then
+		wfile:close()
+		return
+	end
+	local mfilename = modpath..filename
+	local mfile = io.open(mfilename, 'r')
+	local wfile = io.open(wfilename, 'w')
+	wfile:write(mfile:read("*all"))
+	mfile:close()
+	wfile:close()
+end
+
+copy_if_needed('size')
+local sfile = io.open(worldpath..'size')
+local X = tonumber(sfile:read('*l'))
+local Z = tonumber(sfile:read('*l'))
+
+copy_if_needed('dem')
+local dem = load_map(worldpath..'dem', 2, true)
+copy_if_needed('lakes')
+local lakes = load_map(worldpath..'lakes', 2, true)
+copy_if_needed('links')
+local links = load_map(worldpath..'links', 1, false)
+copy_if_needed('rivers')
+local rivers = load_map(worldpath..'rivers', 4, false)
+
+local function index(x, z)
+	return z*X+x+1
+end
+
+local function interp(v00, v01, v10, v11, xf, zf)
+	v0 = v01*xf + v00*(1-xf)
+	v1 = v11*xf + v10*(1-xf)
+	return v1*zf + v0*(1-zf)
+end
+
+local data = {}
+
+local blocksize = 6
+local sea_level = 1
+local min_catchment = 25
+
+local storage = minetest.get_mod_storage()
+if storage:contains("blocksize") then
+	blocksize = storage:get_int("blocksize")
+else
+	storage:set_int("blocksize", blocksize)
+end
+if storage:contains("sea_level") then
+	sea_level = storage:get_int("sea_level")
+else
+	storage:set_int("sea_level", sea_level)
+end
+if storage:contains("min_catchment") then
+	min_catchment = storage:get_float("min_catchment")
+else
+	storage:set_float("min_catchment", min_catchment)
+end
+
+local function generate(minp, maxp, seed)
+	local c_stone = minetest.get_content_id("default:stone")
+	local c_dirt = minetest.get_content_id("default:dirt")
+	local c_lawn = minetest.get_content_id("default:dirt_with_grass")
+	local c_sand = minetest.get_content_id("default:sand")
+	local c_water = minetest.get_content_id("default:water_source")
+	local c_rwater = minetest.get_content_id("default:river_water_source")
+
+	local vm, emin, emax = minetest.get_mapgen_object("voxelmanip")
+	vm:get_data(data)
+
+	local a = VoxelArea:new({MinEdge = emin, MaxEdge = emax})
+	local ystride = a.ystride -- Tip : the ystride of a VoxelArea is the number to add to the array index to get the index of the position above. It's faster because it avoids to completely recalculate the index.
+
+	for x = minp.x, maxp.x do
+		for z = minp.z, maxp.z do
+			local xb = x/blocksize
+			local zb = z/blocksize
+
+			if xb >= 0 and xb < X-1 and zb >= 0 and zb < Z-1 then
+				local x0 = math.floor(xb)
+				local x1 = x0+1
+				local z0 = math.floor(zb)
+				local z1 = z0+1
+
+				local xf = xb - x0
+				local zf = zb - z0
+
+				local i00 = index(x0,z0)
+				local i01 = index(x1,z0)
+				local i10 = index(x0,z1)
+				local i11 = index(x1,z1)
+
+				local terrain_height = math.floor(interp(
+					dem[i00],
+					dem[i01],
+					dem[i10],
+					dem[i11],
+					xf, zf
+				))
+
+				local lake_height = math.floor(math.min(
+					lakes[i00],
+					lakes[i01],
+					lakes[i10],
+					lakes[i11]
+				))
+
+				local is_lake = lake_height > terrain_height
+
+				local is_river = false
+				if xf == 0 then
+					if links[i00] == 1 and rivers[i00] >= min_catchment then
+						is_river = true
+					elseif links[i10] == 3 and rivers[i10] >= min_catchment then
+						is_river = true
+					end
+				end
+
+				if zf == 0 then
+					if links[i00] == 2 and rivers[i00] >= min_catchment then
+						is_river = true
+					elseif links[i01] == 4 and rivers[i01] >= min_catchment then
+						is_river = true
+					end
+				end
+				
+				local ivm = a:index(x, minp.y-1, z)
+
+				if terrain_height >= minp.y then
+					for y=minp.y, math.min(maxp.y, terrain_height) do
+						if y == terrain_height then
+							if is_lake or y <= sea_level then
+								data[ivm] = c_sand
+							elseif is_river then
+								data[ivm] = c_rwater
+							else
+								data[ivm] = c_lawn
+							end
+						else
+							data[ivm] = c_stone
+						end
+						ivm = ivm + ystride
+					end
+				end
+
+				if lake_height > sea_level then
+					if is_lake and lake_height > minp.y then
+						for y=math.max(minp.y, terrain_height+1), math.min(maxp.y, lake_height) do
+							data[ivm] = c_rwater
+							ivm = ivm + ystride
+						end
+					end
+				else
+					for y=math.max(minp.y, terrain_height+1), math.min(maxp.y, sea_level) do
+						data[ivm] = c_water
+						ivm = ivm + ystride
+					end
+				end
+			end
+		end
+	end
+
+	vm:set_data(data)
+	minetest.generate_ores(vm, minp, maxp)
+	vm:set_lighting({day = 0, night = 0})
+	vm:calc_lighting()
+	vm:update_liquids()
+	vm:write_to_map()
+end
+
+minetest.register_on_generated(generate)
diff --git a/load.lua b/load.lua
new file mode 100644
index 0000000..597c356
--- /dev/null
+++ b/load.lua
@@ -0,0 +1,28 @@
+local function load_map(filename, bytes, signed)
+	local file = io.open(filename, 'r')
+	local data = file:read('*all')
+
+	local map = {}
+
+	local size = math.floor(#data/bytes)
+
+	for i=1, size do
+		local i0, i1 = (i-1)*bytes+1, i*bytes
+		local elements = {data:byte(i0, i1)}
+		local n = elements[1]
+		if signed and n >= 128 then
+			n = n - 256
+		end
+
+		for j=2, bytes do
+			n = n*256 + elements[j]
+		end
+
+		map[i] = n
+	end
+	file:close()
+
+	return map
+end
+
+return load_map
diff --git a/rivermapper.py b/rivermapper.py
new file mode 100644
index 0000000..52b3aa4
--- /dev/null
+++ b/rivermapper.py
@@ -0,0 +1,98 @@
+import numpy as np
+import heapq
+import sys
+
+# Conventions:
+# 1 = South (+Y)
+# 2 = East  (+X)
+# 3 = North (-Y)
+# 4 = West  (-X)
+
+sys.setrecursionlimit(65536)
+
+neighbours_dirs = np.array([
+    [0,1,0],
+    [2,0,4],
+    [0,3,0],
+], dtype='u1')
+
+neighbours_pattern = neighbours_dirs > 0
+
+def flow_dirs_lakes(dem):
+    (Y, X) = dem.shape
+
+    dem_margin = np.zeros((Y+2, X+2))
+    dem_margin[1:-1,1:-1] = dem
+
+    # Initialize: list map borders
+    borders = []
+
+    for x in range(1,X+1):
+        dem_north = dem_margin[1,x]
+        borders.append((dem_north, dem_north, 1, x))
+        dem_south = dem_margin[Y,x]
+        borders.append((dem_south, dem_south, Y, x))
+
+    for y in range(2,Y):
+        dem_west = dem_margin[y,1]
+        borders.append((dem_west, dem_west, y, 1))
+        dem_east = dem_margin[y,X]
+        borders.append((dem_east, dem_east, y, X))
+
+    heapq.heapify(borders)
+
+    dirs = np.zeros((Y+2, X+2), dtype='u1')
+    dirs[-2:,:] = 1
+    dirs[:,-2:] = 2
+    dirs[ :2,:] = 3
+    dirs[:, :2] = 4
+
+    lakes = np.zeros((Y, X))
+
+    def add_point(y, x, altmax):
+        alt = dem_margin[y, x]
+        heapq.heappush(borders, (alt, altmax, y, x))
+
+    while len(borders) > 0:
+        (alt, altmax, y, x) = heapq.heappop(borders)
+        neighbours = dirs[y-1:y+2, x-1:x+2]
+        empty_neighbours = (neighbours == 0) * neighbours_pattern
+        neighbours += empty_neighbours * neighbours_dirs
+
+        lake = max(alt, altmax)
+        lakes[y-1,x-1] = lake
+
+        coords = np.transpose(empty_neighbours.nonzero())
+        for (dy,dx) in coords-1:
+            add_point(y+dy, x+dx, lake)
+
+    return dirs[1:-1,1:-1], lakes
+
+def accumulate(dirs, dem=None):
+    (Y, X) = dirs.shape
+    dirs_margin = np.zeros((Y+2,X+2))-1
+    dirs_margin[1:-1,1:-1] = dirs
+    quantity = np.zeros((Y, X), dtype='u4')
+
+    def calculate_quantity(y, x):
+        if quantity[y,x] > 0:
+            return quantity[y,x]
+        q = 1
+        neighbours = dirs_margin[y:y+3, x:x+3]
+        donors = neighbours == neighbours_dirs
+
+        coords = np.transpose(donors.nonzero())
+        for (dy,dx) in coords-1:
+            q += calculate_quantity(y+dy, x+dx)
+        quantity[y, x] = q
+        return q
+
+    for x in range(X):
+        for y in range(Y):
+            calculate_quantity(y, x)
+
+    return quantity
+
+def flow(dem):
+    dirs, lakes = flow_dirs_lakes(dem)
+    return dirs, lakes, accumulate(dirs, dem)
diff --git a/save.py b/save.py
new file mode 100644
index 0000000..f1c356d
--- /dev/null
+++ b/save.py
@@ -0,0 +1,8 @@
+import numpy as np
+
+def save(data, fname, dtype=None):
+    if dtype is not None:
+        data = data.astype(dtype)
+
+    with open(fname, 'wb') as f:
+        f.write(data.tobytes())
diff --git a/terrain_rivers.py b/terrain_rivers.py
new file mode 100755
index 0000000..55ed5a1
--- /dev/null
+++ b/terrain_rivers.py
@@ -0,0 +1,101 @@
+#!/usr/bin/env python3
+
+import numpy as np
+import noise
+from save import save
+from erosion import EvolutionModel
+import os
+import sys
+
+# Always place in this script's parent directory
+os.chdir(os.path.dirname(sys.argv[0]))
+argc = len(sys.argv)
+
+if argc > 1:
+	mapsize = int(sys.argv[1])
+else:
+	mapsize = 400
+
+scale = mapsize / 2
+n = np.zeros((mapsize, mapsize))
+
+#micronoise_depth = 0.05
+
+params = {
+    "octaves" : 8,
+    "persistence" : 0.5,
+    "lacunarity" : 2.,
+}
+
+xbase = np.random.randint(65536)
+ybase = np.random.randint(65536)
+
+for x in range(mapsize):
+    for y in range(mapsize):
+        n[x,y] = noise.snoise2(x/scale + xbase, y/scale + ybase, **params)
+
+#micronoise = np.random.rand(mapsize, mapsize)
+#nn = np.exp(n*2) + micronoise*micronoise_depth
+nn = n*mapsize/5 + mapsize/20
+
+print('Initializing model')
+model = EvolutionModel(nn, K=1, m=0.35, d=1, sea_level=0)
+
+print('Flow calculation 1')
+model.calculate_flow()
+
+print('Advection 1')
+model.advection(2)
+
+print('Flow calculation 2')
+model.calculate_flow()
+
+print('Diffusion')
+model.diffusion(4)
+
+print('Advection 2')
+model.advection(2)
+
+print('Flow calculation 3')
+model.calculate_flow()
+
+print('Done')
+
+save(model.dem, 'dem', dtype='>i2')
+save(model.lakes, 'lakes', dtype='>i2')
+save(model.dirs, 'links', dtype='u1')
+save(model.rivers, 'rivers', dtype='>u4')
+
+with open('size', 'w') as sfile:
+    sfile.write('{:d}\n{:d}'.format(mapsize, mapsize))
+
+try:
+    import matplotlib.pyplot as plt
+
+    plt.subplot(2,2,1)
+    plt.pcolormesh(nn, cmap='viridis')
+    plt.gca().set_aspect('equal', 'box')
+    #plt.colorbar(orientation='horizontal')
+    plt.title('Raw elevation')
+
+    plt.subplot(2,2,2)
+    plt.pcolormesh(model.lakes, cmap='viridis')
+    plt.gca().set_aspect('equal', 'box')
+    #plt.colorbar(orientation='horizontal')
+    plt.title('Lake surface elevation')
+
+    plt.subplot(2,2,3)
+    plt.pcolormesh(model.dem, cmap='viridis')
+    plt.gca().set_aspect('equal', 'box')
+    #plt.colorbar(orientation='horizontal')
+    plt.title('Elevation after advection')
+
+    plt.subplot(2,2,4)
+    plt.pcolormesh(model.rivers, vmin=0, vmax=mapsize**2/25, cmap='Blues')
+    plt.gca().set_aspect('equal', 'box')
+    #plt.colorbar(orientation='horizontal')
+    plt.title('Rivers discharge')
+
+    plt.show()
+except:
+    pass