Added glaciers, and re-organized noise definitions

Uses noise parameters of builtin biomegen

.gitignore

@@ -1,7 +1,12 @@
 __pycache__/
 dem
 lakes
-links
-rivers
 size
+offset_x
+offset_y
+bounds_x
+bounds_y
+dirs
+rivers
 unused/
+data/

README.md

@@ -1,27 +1,36 @@
 mapgen_rivers
 =============
 
-Procedural map generator for Minetest 5.x. Still experimental and basic.
+Procedural map generator for Minetest 5.x. Focused on river networks, and features valley erosion and lakes.
 
 Contains two distinct programs: Python scripts for pre-processing, and Lua scripts to generate the map on Minetest.
 
+![Screenshot](https://user-images.githubusercontent.com/6905002/79541028-687b3000-8089-11ea-9209-c23c15d75383.png)
+
 # Installation
 This mod should be placed in the `/mods` directory like any other Minetest mod.
 
 The Python part relies on external libraries that you need to install:
-- `numpy`, a widely used library for numerical calculations
+- `numpy` and `scipy`, widely used libraries for numerical calculations
 - `noise`, doing Perlin/Simplex noises
 - optionally, `matplotlib` (for map preview)
 
 They are commonly found on `pip` or `conda` Python distributions.
 
 # Usage
+By default, the mod contains a demo 400x400 map (so you can start the game directly), but it is recommended that you run the pre-processing script to generate a new map before world creation, if you can.
+
 ## Pre-processing
 Run the script `terrain_rivers.py` via command line. You can optionally append the map size (by default 400). Example for a 1000x1000 map:
 ```
 ./terrain_rivers.py 1000
 ```
-For a default 400x400 map, it should take between 1 and 2 minutes. It will generate 5 files directly in the mod folder, containing the map data (1.4 MB for the default size).
+For a default 400x400 map, it should take between 1 and 2 minutes. It will generate 5 files directly in the mod folder, containing the map data.
+
+If you have `matplotlib` installed, the script `view_map.py` can be used to get a map preview. Example:
+```
+./view_map.py data/
+```
 
 ## Map generation
 Just create a Minetest world with `singlenode` mapgen, enable this mod and start the world. The data files are immediately copied in the world folder so you can re-generate them afterwards, it won't affect the old worlds.

bounds.py

@@ -0,0 +1,74 @@
+import numpy as np
+
+def make_bounds(dirs, rivers):
+    """
+    Give an array of all horizontal and vertical bounds
+    """
+
+    (Y, X) = dirs.shape
+    bounds_h = np.zeros((Y, X-1), dtype='i4')
+    bounds_v = np.zeros((Y-1, X), dtype='i4')
+
+    bounds_v += (rivers * (dirs==1))[:-1,:]
+    bounds_h += (rivers * (dirs==2))[:,:-1]
+    bounds_v -= (rivers * (dirs==3))[1:,:]
+    bounds_h -= (rivers * (dirs==4))[:,1:]
+
+    return bounds_h, bounds_v
+
+def get_fixed(dirs):
+    """
+    Give the list of points that should not be twisted
+    """
+
+    borders = np.zeros(dirs.shape, dtype='?')
+    borders[-1,:] |= dirs[-1,:]==1
+    borders[:,-1] |= dirs[:,-1]==2
+    borders[0,:] |= dirs[0,:]==3
+    borders[:,0] |= dirs[:,0]==4
+
+    donors = np.zeros(dirs.shape, dtype='?')
+    donors[1:,:] |= dirs[:-1,:]==1
+    donors[:,1:] |= dirs[:,:-1]==2
+    donors[:-1,:] |= dirs[1:,:]==3
+    donors[:,:-1] |= dirs[:,1:]==4
+    return borders | ~donors
+
+def twist(bounds_x, bounds_y, fixed, d=0.1, n=5):
+    """
+    Twist the grid (define an offset for every node). Model river bounds as if they were elastics.
+    Smoothes preferentially big rivers.
+    """
+
+    moveable = ~fixed
+
+    (Y, X) = fixed.shape
+    offset_x = np.zeros((Y, X))
+    offset_y = np.zeros((Y, X))
+
+    for i in range(n):
+        force_long = np.abs(bounds_x) * (1+np.diff(offset_x, axis=1))
+        force_trans = np.abs(bounds_y) * np.diff(offset_x, axis=0)
+
+        force_x = np.zeros((Y, X))
+        force_x[:,:-1] = force_long
+        force_x[:,1:] -= force_long
+        force_x[:-1,:]+= force_trans
+        force_x[1:,:] -= force_trans
+
+        force_long = np.abs(bounds_y) * (1+np.diff(offset_y, axis=0))
+        force_trans = np.abs(bounds_x) * np.diff(offset_y, axis=1)
+
+        force_y = np.zeros((Y, X))
+        force_y[:-1,:] = force_long
+        force_y[1:,:] -= force_long
+        force_y[:,:-1]+= force_trans
+        force_y[:,1:] -= force_trans
+
+        length = np.hypot(force_x, force_y)
+        length[length==0] = 1
+        coeff = d / length * moveable # Normalize, take into account the direction only
+        offset_x += force_x * coeff
+        offset_y += force_y * coeff
+
+    return offset_x, offset_y

demo_data/size

@@ -0,0 +1,2 @@
+401
+401

erosion.py

@@ -3,22 +3,33 @@ import scipy.ndimage as im
 import rivermapper as rm
 
 def advection(dem, dirs, rivers, time, K=1, m=0.5, sea_level=0):
+    """
+    Simulate erosion by rivers.
+    This models erosion as an upstream advection of elevations ("erosion waves").
+    Advection speed depends on water flux and parameters:
+
+    v = K * flux^m
+    """
+
     dirs = dirs.copy()
     dirs[0,:] = 0
     dirs[-1,:] = 0
     dirs[:,0] = 0
     dirs[:,-1] = 0
 
-    adv_time = 1 / (K*rivers**m)
+    adv_time = 1 / (K*rivers**m) # For every pixel, calculate the time an "erosion wave" will need to cross it.
     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]):
+            # Elevations propagate upstream, so for every pixel we seek the downstream pixel whose erosion wave just reached the current pixel.
+            # This means summing the advection times downstream until we reach the erosion time.
             x0, y0 = x, y
             x1, y1 = x, y
             remaining = time
             while True:
+                # Move one pixel downstream
                 flow_dir = dirs[y0,x0]
                 if flow_dir == 0:
                     remaining = 0
@@ -32,19 +43,19 @@ def advection(dem, dirs, rivers, time, K=1, m=0.5, sea_level=0):
                 elif flow_dir == 4:
                     x1 -= 1
 
-                if remaining <= adv_time[y0,x0]:
+                if remaining <= adv_time[y0,x0]: # Time is over, we found it.
                     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]
+            dem_new[y,x] = c*dem[y1,x1] + (1-c)*dem[y0,x0] # If between 2 pixels, perform linear interpolation.
 
     return np.minimum(dem, dem_new)
 
 def diffusion(dem, time, d=1):
     radius = d * time**.5
-    return im.gaussian_filter(dem, radius, mode='reflect')
+    return im.gaussian_filter(dem, radius, mode='reflect') # Diffusive erosion is a simple Gaussian blur
 
 class EvolutionModel:
     def __init__(self, dem, K=1, m=0.5, d=1, sea_level=0, flow=False, flex_radius=100):
@@ -78,9 +89,9 @@ class EvolutionModel:
         self.flow_uptodate = False
 
     def define_isostasy(self):
-        self.ref_isostasy = im.gaussian_filter(self.dem, self.flex_radius, mode='reflect')
+        self.ref_isostasy = im.gaussian_filter(self.dem, self.flex_radius, mode='reflect') # Define a blurred version of the DEM that will be considered as the reference isostatic elevation.
 
     def adjust_isostasy(self, rate=1):
-        isostasy = im.gaussian_filter(self.dem, self.flex_radius, mode='reflect')
+        isostasy = im.gaussian_filter(self.dem, self.flex_radius, mode='reflect') # Calculate blurred DEM
-        correction = (self.ref_isostasy - isostasy) * rate
+        correction = (self.ref_isostasy - isostasy) * rate # Compare it with the reference isostasy
-        self.dem = self.dem + correction
+        self.dem = self.dem + correction # Adjust

geometry.lua

@@ -0,0 +1,37 @@
+local function distance_to_segment(x1, y1, x2, y2, x, y)
+	-- get the distance between point (x,y) and segment (x1,y1)-(x2,y2)
+	local a = (x1-x2)^2 + (y1-y2)^2 -- square of distance
+	local b = (x1-x)^2 + (y1-y)^2
+	local c = (x2-x)^2 + (y2-y)^2
+	if a + b < c then
+		-- The closest point of the segment is the extremity 1
+		return math.sqrt(b)
+	elseif a + c < b then
+		-- The closest point of the segment is the extremity 2
+		return math.sqrt(c)
+	else
+		-- The closest point is on the segment
+		return math.abs(x1 * (y2-y) + x2 * (y-y1) + x * (y1-y2)) / math.sqrt(a)
+	end
+end
+
+local function transform_quadri(X, Y, x, y)
+	-- To index points in an irregular quadrilateral, giving x and y between 0 (one edge) and 1 (opposite edge)
+	-- X, Y 4-vectors giving the coordinates of the 4 vertices
+	-- x, y position to index.
+	local x1, x2, x3, x4 = unpack(X)
+	local y1, y2, y3, y4 = unpack(Y)
+
+	-- Compare distance to 2 opposite edges, they give the X coordinate
+	local d23 = distance_to_segment(x2,y2,x3,y3,x,y)
+	local d41 = distance_to_segment(x4,y4,x1,y1,x,y)
+	local xc = d41 / (d23+d41)
+
+	-- Same for the 2 other edges, they give the Y coordinate
+	local d12 = distance_to_segment(x1,y1,x2,y2,x,y)
+	local d34 = distance_to_segment(x3,y3,x4,y4,x,y)
+	local yc = d12 / (d12+d34)
+	return xc, yc
+end
+
+return transform_quadri

heightmap.lua

@@ -0,0 +1,119 @@
+local modpath = minetest.get_modpath(minetest.get_current_modname()) .. '/'
+
+local make_polygons = dofile(modpath .. 'polygons.lua')
+local transform_quadri = dofile(modpath .. 'geometry.lua')
+
+local blocksize = mapgen_rivers.blocksize
+local sea_level = mapgen_rivers.sea_level
+local riverbed_slope = mapgen_rivers.riverbed_slope
+
+-- Linear interpolation
+local function interp(v00, v01, v11, v10, xf, zf)
+	local v0 = v01*xf + v00*(1-xf)
+	local v1 = v11*xf + v10*(1-xf)
+	return v1*zf + v0*(1-zf)
+end
+
+local function heightmaps(minp, maxp)
+
+	local polygons = make_polygons(minp, maxp)
+	local incr = maxp.x-minp.x+1
+	local i0 = (minp.z-minp.z) * incr + (minp.x-minp.x) + 1
+
+	local terrain_height_map = {}
+	local lake_height_map = {}
+
+	local i = 1
+	for x=minp.x, maxp.x do
+		for z=minp.z, maxp.z do
+			local poly = polygons[i]
+
+			if poly then
+				local xf, zf = transform_quadri(poly.x, poly.z, x/blocksize, z/blocksize)
+				local i00, i01, i11, i10 = unpack(poly.i)
+
+				-- Load river width on 4 edges and corners
+				local r_west, r_north, r_east, r_south = unpack(poly.rivers)
+				local c_NW, c_NE, c_SE, c_SW = unpack(poly.river_corners)
+
+				-- Calculate the depth factor for each edge and corner.
+				-- Depth factor:
+				-- < 0: outside river
+				-- = 0: on riverbank
+				-- > 0: inside river
+				local depth_factors = {
+					r_west - xf,
+					r_north - zf,
+					xf - r_east,
+					zf - r_south,
+					c_NW-xf-zf,
+					xf-zf-c_NE,
+					xf+zf-c_SE,
+					zf-xf-c_SW,
+				}
+
+				-- Find the maximal depth factor and determine to which river it belongs
+				local depth_factor_max = 0
+				local imax = 0
+				for i=1, 8 do
+					if depth_factors[i] >= depth_factor_max then
+						depth_factor_max = depth_factors[i]
+						imax = i
+					end
+				end
+
+				-- Transform the coordinates to have xf and zf = 0 or 1 in rivers (to avoid rivers having lateral slope and to accomodate the surrounding smoothly)
+				if imax == 0 then
+					local x0 = math.max(r_west, c_NW-zf, zf-c_SW)
+					local x1 = math.min(r_east, c_NE+zf, c_SE-zf)
+					local z0 = math.max(r_north, c_NW-xf, xf-c_NE)
+					local z1 = math.min(r_south, c_SW+xf, c_SE-xf)
+					xf = (xf-x0) / (x1-x0)
+					zf = (zf-z0) / (z1-z0)
+				elseif imax == 1 then
+					xf = 0
+				elseif imax == 2 then
+					zf = 0
+				elseif imax == 3 then
+					xf = 1
+				elseif imax == 4 then
+					zf = 1
+				elseif imax == 5 then
+					xf, zf = 0, 0
+				elseif imax == 6 then
+					xf, zf = 1, 0
+				elseif imax == 7 then
+					xf, zf = 1, 1
+				elseif imax == 8 then
+					xf, zf = 0, 1
+				end
+
+				-- Determine elevation by interpolation
+				local vdem = poly.dem
+				local terrain_height = math.floor(0.5+interp(
+					vdem[1],
+					vdem[2],
+					vdem[3],
+					vdem[4],
+					xf, zf
+				))
+
+				local lake_height = math.max(math.floor(poly.lake), terrain_height)
+				if imax > 0 and depth_factor_max > 0 then
+					terrain_height = math.min(math.max(lake_height, sea_level) - math.floor(1+depth_factor_max*riverbed_slope), terrain_height)
+				end
+
+				terrain_height_map[i] = terrain_height
+				lake_height_map[i] = lake_height
+			else
+				terrain_height_map[i] = -31000
+				lake_height_map[i] = -31000
+			end
+			i = i + 1
+		end
+	end
+
+	return terrain_height_map, lake_height_map
+end
+
+return heightmaps

init.lua

@@ -1,76 +1,105 @@
+mapgen_rivers = {}
+
 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)
+dofile(modpath .. 'settings.lua')
-	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 blocksize = mapgen_rivers.blocksize
-local sfile = io.open(worldpath..'size')
+local sea_level = mapgen_rivers.sea_level
-local X = tonumber(sfile:read('*l'))
+local riverbed_slope = mapgen_rivers.riverbed_slope
-local Z = tonumber(sfile:read('*l'))
+local elevation_chill = mapgen_rivers.elevation_chill
 
-copy_if_needed('dem')
+dofile(modpath .. 'noises.lua')
-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)
+local make_polygons = dofile(modpath .. 'polygons.lua')
-	return z*X+x+1
-end
 
-local function interp(v00, v01, v10, v11, xf, zf)
+local transform_quadri = dofile(modpath .. 'geometry.lua')
-	v0 = v01*xf + v00*(1-xf)
+
-	v1 = v11*xf + v10*(1-xf)
+local heightmaps = dofile(modpath .. 'heightmap.lua')
+
+-- Linear interpolation
+local function interp(v00, v01, v11, v10, xf, zf)
+	local v0 = v01*xf + v00*(1-xf)
+	local v1 = v11*xf + v10*(1-xf)
 	return v1*zf + v0*(1-zf)
 end
 
 local data = {}
 
-local blocksize = 6
+local noise_x_obj, noise_z_obj, noise_distort_obj, noise_heat_obj, noise_heat_blend_obj
-local sea_level = 1
+local noise_x_map = {}
-local min_catchment = 25
+local noise_z_map = {}
-
+local noise_distort_map = {}
-local storage = minetest.get_mod_storage()
+local noise_heat_map = {}
-if storage:contains("blocksize") then
+local noise_heat_blend_map = {}
-	blocksize = storage:get_int("blocksize")
+local mapsize
-else
+local init = false
-	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 chulens = {
+		x = maxp.x-minp.x+1,
+		y = maxp.y-minp.y+1,
+		z = maxp.z-minp.z+1,
+	}
+
+	if not init then
+		mapsize = {
+			x = chulens.x,
+			y = chulens.y+1,
+			z = chulens.z,
+		}
+		noise_x_obj = minetest.get_perlin_map(mapgen_rivers.noise_params.distort_x, mapsize)
+		noise_z_obj = minetest.get_perlin_map(mapgen_rivers.noise_params.distort_z, mapsize)
+		noise_heat_obj = minetest.get_perlin_map(mapgen_rivers.noise_params.heat, chulens)
+		noise_heat_blend_obj = minetest.get_perlin_map(mapgen_rivers.noise_params.heat_blend, chulens)
+		noise_distort_obj = minetest.get_perlin_map(mapgen_rivers.noise_params.distort_amplitude, chulens)
+		init = true
+	end
+
+	local minp2d = {x=minp.x, y=minp.z}
+	noise_x_obj:get_3d_map_flat(minp, noise_x_map)
+	noise_z_obj:get_3d_map_flat(minp, noise_z_map)
+	noise_distort_obj:get_2d_map_flat(minp2d, noise_distort_map)
+	noise_heat_obj:get_2d_map_flat(minp2d, noise_heat_map)
+	noise_heat_blend_obj:get_2d_map_flat(minp2d, noise_heat_blend_map)
+
+	local xmin, xmax, zmin, zmax = minp.x, maxp.x, minp.z, maxp.z
+	local i = 0
+	local i2d = 0
+	for z=minp.z, maxp.z do
+		for y=minp.y, maxp.y+1 do
+			for x=minp.x, maxp.x do
+				i = i+1
+				i2d = i2d+1
+				local distort = noise_distort_map[i2d]
+				local xv = noise_x_map[i]*distort + x
+				if xv < xmin then xmin = xv end
+				if xv > xmax then xmax = xv end
+				noise_x_map[i] = xv
+				local zv = noise_z_map[i]*distort + z
+				if zv < zmin then zmin = zv end
+				if zv > zmax then zmax = zv end
+				noise_z_map[i] = zv
+			end
+			i2d = i2d-chulens.x
+		end
+	end
+
+	local pminp = {x=math.floor(xmin), z=math.floor(zmin)}
+	local pmaxp = {x=math.floor(xmax)+1, z=math.floor(zmax)+1}
+	local incr = pmaxp.z-pminp.z+1
+	local i_origin = 1 - pminp.x*incr - pminp.z
+	local terrain_map, lake_map = heightmaps(pminp, pmaxp)
+
 	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_dirtsnow = minetest.get_content_id("default:dirt_with_snow")
+	local c_snow = minetest.get_content_id("default:snowblock")
 	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 c_ice = minetest.get_content_id("default:ice")
 
 	local vm, emin, emax = minetest.get_mapgen_object("voxelmanip")
 	vm:get_data(data)
@@ -78,93 +107,71 @@ local function generate(minp, maxp, seed)
 	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
+	local nid = mapsize.x*(mapsize.y-1) + 1
+	local incrY = -mapsize.x
+	local incrX = 1 - mapsize.y*incrY
+	local incrZ = mapsize.x*mapsize.y - mapsize.x*incrX - mapsize.x*mapsize.y*incrY
+
+	local i2d = 1
+	
 	for z = minp.z, maxp.z do
-			local xb = x/blocksize
+		for x = minp.x, maxp.x do
-			local zb = z/blocksize
+			local ivm = a:index(x, minp.y, z)
+			local ground_above = false
+			local temperature = noise_heat_map[i2d]+noise_heat_blend_map[i2d]
+			for y = maxp.y+1, minp.y, -1 do
+				local xn = noise_x_map[nid]
+				local zn = noise_z_map[nid]
+				local x0 = math.floor(xn)
+				local z0 = math.floor(zn)
 
-			if xb >= 0 and xb < X-1 and zb >= 0 and zb < Z-1 then
+				local i0 = i_origin + x0*incr + z0
-				local x0 = math.floor(xb)
+				local i1 = i0+incr
-				local x1 = x0+1
+				local i2 = i1+1
-				local z0 = math.floor(zb)
+				local i3 = i0+1
-				local z1 = z0+1
 
-				local xf = xb - x0
+				local terrain = interp(terrain_map[i0], terrain_map[i1], terrain_map[i2], terrain_map[i3], xn-x0, zn-z0)
-				local zf = zb - z0
+				if y <= maxp.y then
+					local lake = math.min(lake_map[i0], lake_map[i1], lake_map[i2], lake_map[i3])
 
-				local i00 = index(x0,z0)
+					local is_lake = lake > terrain
-				local i01 = index(x1,z0)
+					local ivm = a:index(x, y, z)
-				local i10 = index(x0,z1)
+					if y <= terrain then
-				local i11 = index(x1,z1)
+						if y <= terrain-1 or ground_above then
-
-				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
+						elseif is_lake or y < sea_level then
+							data[ivm] = c_sand
+						else
+							local temperature_y = temperature - y*elevation_chill
+							if temperature_y >= 15 then
+								data[ivm] = c_lawn
+							elseif temperature_y >= 0 then
+								data[ivm] = c_dirtsnow
+							else
+								data[ivm] = c_snow
 							end
-						ivm = ivm + ystride
+						end
+					elseif y <= lake and lake > sea_level then
+						local temperature_y = temperature - y*elevation_chill
+						if temperature_y >= 0 then
+							data[ivm] = c_rwater
+						else
+							data[ivm] = c_ice
+						end
+					elseif y <= sea_level then
+						data[ivm] = c_water
 					end
 				end
 
-				if lake_height > sea_level then
+				ground_above = y <= terrain
-					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
+				nid = nid + incrY
 			end
+			nid = nid + incrX
+			i2d = i2d + 1
 		end
-				else
+		nid = nid + incrZ
-					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)

load.lua

@@ -1,11 +1,14 @@
-local function load_map(filename, bytes, signed)
+local worldpath = minetest.get_worldpath() .. "/river_data/"
-	local file = io.open(filename, 'r')
+
+local function load_map(filename, bytes, signed, size)
+	local file = io.open(worldpath .. filename, 'r')
 	local data = file:read('*all')
+	if #data < bytes*size then
+		data = minetest.decompress(data)
+	end
 
 	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)}

noises.lua

@@ -0,0 +1,37 @@
+mapgen_rivers.noise_params = {
+	distort_x = {
+		offset = 0,
+		scale = 1,
+		seed = -4574,
+		spread = {x=64, y=32, z=64},
+		octaves = 3,
+		persistence = 0.75,
+		lacunarity = 2,
+	},
+
+	distort_z = {
+		offset = 0,
+		scale = 1,
+		seed = -7940,
+		spread = {x=64, y=32, z=64},
+		octaves = 3,
+		persistence = 0.75,
+		lacunarity = 2,
+	},
+
+	distort_amplitude = {
+		offset = 0,
+		scale = 10,
+		seed = 676,
+		spread = {x=1024, y=1024, z=1024},
+		octaves = 5,
+		persistence = 0.5,
+		lacunarity = 2,
+		flags = "absvalue",
+	},
+
+	heat = minetest.get_mapgen_setting_noiseparams('mg_biome_np_heat'),
+	heat_blend = minetest.get_mapgen_setting_noiseparams('mg_biome_np_heat_blend'),
+}
+
+mapgen_rivers.noise_params.heat.offset = mapgen_rivers.noise_params.heat.offset + mapgen_rivers.sea_level*mapgen_rivers.elevation_chill

polygons.lua

@@ -0,0 +1,205 @@
+local modpath = minetest.get_modpath(minetest.get_current_modname()) .. '/'
+local mod_data_path = modpath .. 'data/'
+if not io.open(mod_data_path .. 'size', 'r') then
+	mod_data_path = modpath .. 'demo_data/'
+end
+
+local world_data_path = minetest.get_worldpath() .. '/river_data/'
+minetest.mkdir(world_data_path)
+
+local load_map = dofile(modpath .. 'load.lua')
+
+local function copy_if_needed(filename)
+	local wfilename = world_data_path..filename
+	local wfile = io.open(wfilename, 'r')
+	if wfile then
+		wfile:close()
+		return
+	end
+	local mfilename = mod_data_path..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(world_data_path..'size')
+local X = tonumber(sfile:read('*l'))
+local Z = tonumber(sfile:read('*l'))
+sfile:close()
+
+copy_if_needed('dem')
+local dem = load_map('dem', 2, true, X*Z)
+copy_if_needed('lakes')
+local lakes = load_map('lakes', 2, true, X*Z)
+copy_if_needed('dirs')
+local dirs = load_map('dirs', 1, false, X*Z)
+copy_if_needed('rivers')
+local rivers = load_map('rivers', 4, false, X*Z)
+
+copy_if_needed('offset_x')
+local offset_x = load_map('offset_x', 1, true, X*Z)
+for k, v in ipairs(offset_x) do
+	offset_x[k] = (v+0.5)/256
+end
+
+copy_if_needed('offset_y')
+local offset_z = load_map('offset_y', 1, true, X*Z)
+for k, v in ipairs(offset_z) do
+	offset_z[k] = (v+0.5)/256
+end
+
+-- To index a flat array representing a 2D map
+local function index(x, z)
+	return z*X+x+1
+end
+
+local blocksize = mapgen_rivers.blocksize
+local min_catchment = mapgen_rivers.min_catchment
+local max_catchment = mapgen_rivers.max_catchment
+
+-- Width coefficients: coefficients solving
+--   wfactor * min_catchment ^ wpower = 1/(2*blocksize)
+--   wfactor * max_catchment ^ wpower = 1
+local wpower = math.log(2*blocksize)/math.log(max_catchment/min_catchment)
+local wfactor = 1 / max_catchment ^ wpower
+local function river_width(flow)
+	flow = math.abs(flow)
+	if flow < min_catchment then
+		return 0
+	end
+
+	return math.min(wfactor * flow ^ wpower, 1)
+end
+
+local noise_heat -- Need a large-scale noise here so no heat blend
+local elevation_chill = mapgen_rivers.elevation_chill
+local function get_temperature(x, y, z)
+	local pos = {x=x, y=z}
+	return noise_heat:get2d(pos) - y*elevation_chill
+end
+
+local glaciers = mapgen_rivers.glaciers
+local glacier_factor = mapgen_rivers.glacier_factor
+
+local init = false
+
+-- On map generation, determine into which polygon every point (in 2D) will fall.
+-- Also store polygon-specific data
+local function make_polygons(minp, maxp)
+	if not init then
+		if glaciers then
+			noise_heat = minetest.get_perlin(mapgen_rivers.noise_params.heat)
+		end
+		init = true
+	end
+
+	local chulens = maxp.z - minp.z + 1
+
+	local polygons = {}
+	-- Determine the minimum and maximum coordinates of the polygons that could be on the chunk, knowing that they have an average size of 'blocksize' and a maximal offset of 0.5 blocksize.
+	local xpmin, xpmax = math.max(math.floor(minp.x/blocksize - 0.5), 0), math.min(math.ceil(maxp.x/blocksize), X-2)
+	local zpmin, zpmax = math.max(math.floor(minp.z/blocksize - 0.5), 0), math.min(math.ceil(maxp.z/blocksize), Z-2)
+
+	-- Iterate over the polygons
+	for xp = xpmin, xpmax do
+		for zp=zpmin, zpmax do
+			local iA = index(xp, zp)
+			local iB = index(xp+1, zp)
+			local iC = index(xp+1, zp+1)
+			local iD = index(xp, zp+1)
+			-- Extract the vertices of the polygon
+			local poly_x = {offset_x[iA]+xp, offset_x[iB]+xp+1, offset_x[iC]+xp+1, offset_x[iD]+xp}
+			local poly_z = {offset_z[iA]+zp, offset_z[iB]+zp, offset_z[iC]+zp+1, offset_z[iD]+zp+1}
+			local polygon = {x=poly_x, z=poly_z, i={iA, iB, iC, iD}}
+
+			local bounds = {} -- Will be a list of the intercepts of polygon edges for every X position (scanline algorithm)
+			-- Calculate the min and max X positions 
+			local xmin = math.max(math.floor(blocksize*math.min(unpack(poly_x)))+1, minp.x)
+			local xmax = math.min(math.floor(blocksize*math.max(unpack(poly_x))), maxp.x)
+			-- And initialize the arrays
+			for x=xmin, xmax do
+				bounds[x] = {}
+			end
+
+			local i1 = 4
+			for i2=1, 4 do -- Loop on 4 edges
+				local x1, x2 = poly_x[i1], poly_x[i2]
+				-- Calculate the integer X positions over which this edge spans
+				local lxmin = math.floor(blocksize*math.min(x1, x2))+1
+				local lxmax = math.floor(blocksize*math.max(x1, x2))
+				if lxmin <= lxmax then -- If there is at least one position in it
+					local z1, z2 = poly_z[i1], poly_z[i2]
+					-- Calculate coefficient of the equation defining the edge: Z=aX+b
+					local a = (z1-z2) / (x1-x2)
+					local b = blocksize*(z1 - a*x1)
+					for x=math.max(lxmin, minp.x), math.min(lxmax, maxp.x) do
+						-- For every X position involved, add the intercepted Z position in the table
+						table.insert(bounds[x], a*x+b)
+					end
+				end
+				i1 = i2
+			end
+			for x=xmin, xmax do
+				-- Now sort the bounds list
+				local xlist = bounds[x]
+				table.sort(xlist)
+				local c = math.floor(#xlist/2)
+				for l=1, c do
+					-- Take pairs of Z coordinates: all positions between them belong to the polygon.
+					local zmin = math.max(math.floor(xlist[l*2-1])+1, minp.z)
+					local zmax = math.min(math.floor(xlist[l*2]), maxp.z)
+					local i = (x-minp.x) * chulens + (zmin-minp.z) + 1
+					for z=zmin, zmax do
+						-- Fill the map at these places
+						polygons[i] = polygon
+						i = i + 1
+					end
+				end
+			end
+
+			local poly_dem = {dem[iA], dem[iB], dem[iC], dem[iD]}
+			polygon.dem = poly_dem
+			polygon.lake = math.min(lakes[iA], lakes[iB], lakes[iC], lakes[iD])
+
+			-- Now, rivers.
+			-- Load river flux values for the 4 corners
+			local riverA = river_width(rivers[iA])
+			local riverB = river_width(rivers[iB])
+			local riverC = river_width(rivers[iC])
+			local riverD = river_width(rivers[iD])
+			if glaciers then -- Widen the river
+				if get_temperature(poly_x[1]*blocksize, poly_dem[1], poly_z[1]*blocksize) < 0 then
+					riverA = math.min(riverA*glacier_factor, 1)
+				end
+				if get_temperature(poly_x[2]*blocksize, poly_dem[2], poly_z[2]*blocksize) < 0 then
+					riverB = math.min(riverB*glacier_factor, 1)
+				end
+				if get_temperature(poly_x[3]*blocksize, poly_dem[3], poly_z[3]*blocksize) < 0 then
+					riverC = math.min(riverC*glacier_factor, 1)
+				end
+				if get_temperature(poly_x[4]*blocksize, poly_dem[4], poly_z[4]*blocksize) < 0 then
+					riverD = math.min(riverD*glacier_factor, 1)
+				end
+			end
+
+			polygon.river_corners = {riverA, 1-riverB, 2-riverC, 1-riverD}
+
+			-- Flow directions
+			local dirA, dirB, dirC, dirD = dirs[iA], dirs[iB], dirs[iC], dirs[iD]
+			-- Determine the river flux on the edges, by testing dirs values
+			local river_west = (dirA==1 and riverA or 0) + (dirD==3 and riverD or 0)
+			local river_north = (dirA==2 and riverA or 0) + (dirB==4 and riverB or 0)
+			local river_east = 1 - (dirB==1 and riverB or 0) - (dirC==3 and riverC or 0)
+			local river_south = 1 - (dirD==2 and riverD or 0) - (dirC==4 and riverC or 0)
+
+			polygon.rivers = {river_west, river_north, river_east, river_south}
+		end
+	end
+
+	return polygons
+end
+
+return make_polygons

rivermapper.py

@@ -18,10 +18,15 @@ neighbours_dirs = np.array([
 
 neighbours_pattern = neighbours_dirs > 0
 
-def flow_dirs_lakes(dem, random=0.0625):
+def flow_dirs_lakes(dem, random=0):
+    """
+    Calculates flow direction in D4 (4 choices) for every pixel of the DEM
+    Also returns an array of lake elevation
+    """
+
     (Y, X) = dem.shape
 
-    dem_margin = np.zeros((Y+2, X+2))
+    dem_margin = np.zeros((Y+2, X+2)) # We need a margin of one pixel at every edge, to prevent crashes when scanning the neighbour pixels on the borders
     dem_margin[1:-1,1:-1] = dem
     if random > 0:
         dem_margin += np.random.random(dem_margin.shape) * random
@@ -41,10 +46,11 @@ def flow_dirs_lakes(dem, random=0.0625):
         dem_east = dem_margin[y,X]
         borders.append((dem_east, dem_east, y, X))
 
+    # Make a binary heap
     heapq.heapify(borders)
 
     dirs = np.zeros((Y+2, X+2), dtype='u1')
-    dirs[-2:,:] = 1
+    dirs[-2:,:] = 1 # Border pixels flow outside the map
     dirs[:,-2:] = 2
     dirs[ :2,:] = 3
     dirs[:, :2] = 4
@@ -56,21 +62,26 @@ def flow_dirs_lakes(dem, random=0.0625):
         heapq.heappush(borders, (alt, altmax, y, x))
 
     while len(borders) > 0:
-        (alt, altmax, y, x) = heapq.heappop(borders)
+        (alt, altmax, y, x) = heapq.heappop(borders) # Take the lowest pixel in the queue
         neighbours = dirs[y-1:y+2, x-1:x+2]
-        empty_neighbours = (neighbours == 0) * neighbours_pattern
+        empty_neighbours = (neighbours == 0) * neighbours_pattern # Find the neighbours whose flow direction is not yet defined
-        neighbours += empty_neighbours * neighbours_dirs
+        neighbours += empty_neighbours * neighbours_dirs # They flow into the pixel being studied
 
-        lake = max(alt, altmax)
+        lake = max(alt, altmax) # Set lake elevation to the maximal height of the downstream section.
         lakes[y-1,x-1] = lake
 
         coords = np.transpose(empty_neighbours.nonzero())
-        for (dy,dx) in coords-1:
+        for (dy,dx) in coords-1: # Add these neighbours into the queue
             add_point(y+dy, x+dx, lake)
 
     return dirs[1:-1,1:-1], lakes
 
 def accumulate(dirs, dem=None):
+    """
+    Calculates the quantity of water that accumulates at every pixel,
+    following flow directions.
+    """
+
     (Y, X) = dirs.shape
     dirs_margin = np.zeros((Y+2,X+2))-1
     dirs_margin[1:-1,1:-1] = dirs
@@ -79,13 +90,13 @@ def accumulate(dirs, dem=None):
     def calculate_quantity(y, x):
         if quantity[y,x] > 0:
             return quantity[y,x]
-        q = 1
+        q = 1 # Consider that every pixel contains a water quantity of 1 by default.
         neighbours = dirs_margin[y:y+3, x:x+3]
-        donors = neighbours == neighbours_dirs
+        donors = neighbours == neighbours_dirs # Identify neighbours that give their water to the pixel being studied
 
         coords = np.transpose(donors.nonzero())
         for (dy,dx) in coords-1:
-            q += calculate_quantity(y+dy, x+dx)
+            q += calculate_quantity(y+dy, x+dx) # Add water quantity of the donors pixels (this triggers calculation for these pixels, recursively)
         quantity[y, x] = q
         return q
 
@@ -96,5 +107,9 @@ def accumulate(dirs, dem=None):
     return quantity
 
 def flow(dem):
+    """
+    Calculates flow directions and water quantity
+    """
+
     dirs, lakes = flow_dirs_lakes(dem)
     return dirs, lakes, accumulate(dirs, dem)

save.py

@@ -1,8 +1,13 @@
 import numpy as np
+import zlib
 
 def save(data, fname, dtype=None):
     if dtype is not None:
         data = data.astype(dtype)
 
+    bin_data = data.tobytes()
+    bin_data_comp = zlib.compress(bin_data, 9)
+    if len(bin_data_comp) < len(bin_data):
+        bin_data = bin_data_comp
     with open(fname, 'wb') as f:
-        f.write(data.tobytes())
+        f.write(bin_data)

settings.lua

@@ -0,0 +1,53 @@
+local storage = minetest.get_mod_storage()
+local settings = minetest.settings
+
+local function get_settings(key, dtype, default)
+	if storage:contains(key) then
+		if dtype == "string" then
+			return storage:get_string(key)
+		elseif dtype == "int" then
+			return storage:get_int(key)
+		elseif dtype == "float" then
+			return storage:get_float(key)
+		elseif dtype == "bool" then
+			return storage:get_string(key) == 'true'
+		end
+	end
+
+	local conf_val = settings:get('mapgen_rivers_' .. key)
+	if conf_val then
+		if dtype == "int" then
+			conf_val = tonumber(conf_val)
+			storage:set_int(key, conf_val)
+		elseif dtype == "float" then
+			conf_val = tonumber(conf_val)
+			storage:set_float(key, conf_val)
+		elseif dtype == "string" or dtype == "bool" then
+			storage:set_string(key, conf_val)
+		end
+
+		return conf_val
+	else
+		if dtype == "int" then
+			storage:set_int(key, default)
+		elseif dtype == "float" then
+			storage:set_float(key, default)
+		elseif dtype == "string" then
+			storage:set_string(key, default)
+		elseif dtype == "bool" then
+			storage:set_string(key, tostring(default))
+		end
+
+		return default
+	end
+end
+
+mapgen_rivers.blocksize = get_settings('blocksize', 'int', 12)
+mapgen_rivers.sea_level = get_settings('sea_level', 'int', 1)
+mapgen_rivers.min_catchment = get_settings('min_catchment', 'float', 25)
+mapgen_rivers.max_catchment = get_settings('max_catchment', 'float', 40000)
+mapgen_rivers.riverbed_slope = get_settings('riverbed_slope', 'float', 0.4) * mapgen_rivers.blocksize
+--mapgen_rivers.distort = get_settings('distort', 'bool', true) To be implemented: should be possible to disable distorsion
+mapgen_rivers.glaciers = get_settings('glaciers', 'bool', true)
+mapgen_rivers.glacier_factor = get_settings('glacier_factor', 'float', 8)
+mapgen_rivers.elevation_chill = get_settings('elevation_chill', 'float', 0.25)

terrain_rivers.py

@@ -4,6 +4,7 @@ import numpy as np
 import noise
 from save import save
 from erosion import EvolutionModel
+import bounds
 import os
 import sys
 
@@ -17,30 +18,31 @@ else:
 	mapsize = 400
 
 scale = mapsize / 2
-n = np.zeros((mapsize, mapsize))
+n = np.zeros((mapsize+1, mapsize+1))
-
-#micronoise_depth = 0.05
 
+# Set noise parameters
 params = {
-    "octaves" : 8,
+    "octaves" : int(np.ceil(np.log2(mapsize)))+1,
     "persistence" : 0.5,
     "lacunarity" : 2.,
 }
 
+# Determine noise offset randomly
 xbase = np.random.randint(65536)
 ybase = np.random.randint(65536)
 
-for x in range(mapsize):
+# Generate the noise
-    for y in range(mapsize):
+for x in range(mapsize+1):
+    for y in range(mapsize+1):
         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
 
+# Initialize landscape evolution model
 print('Initializing model')
 model = EvolutionModel(nn, K=1, m=0.35, d=1, sea_level=0)
 
+# Run the model's processes: the order in which the processes are run is arbitrary and could be changed.
 print('Flow calculation 1')
 model.calculate_flow()
 
@@ -67,41 +69,32 @@ model.calculate_flow()
 
 print('Done')
 
+# Twist the grid
+bx, by = bounds.make_bounds(model.dirs, model.rivers)
+offset_x, offset_y = bounds.twist(bx, by, bounds.get_fixed(model.dirs))
+
+# Convert offset in 8-bits
+offset_x = np.clip(np.floor(offset_x * 256), -128, 127)
+offset_y = np.clip(np.floor(offset_y * 256), -128, 127)
+
+if not os.path.isdir('data'):
+    os.mkdir('data')
+os.chdir('data')
+# Save the files
 save(model.dem, 'dem', dtype='>i2')
 save(model.lakes, 'lakes', dtype='>i2')
-save(model.dirs, 'links', dtype='u1')
+save(offset_x, 'offset_x', dtype='i1')
+save(offset_y, 'offset_y', dtype='i1')
+
+save(model.dirs, 'dirs', dtype='u1')
 save(model.rivers, 'rivers', dtype='>u4')
 
 with open('size', 'w') as sfile:
-    sfile.write('{:d}\n{:d}'.format(mapsize, mapsize))
+    sfile.write('{:d}\n{:d}'.format(mapsize+1, mapsize+1))
 
+# Display the map if matplotlib is found
 try:
-    import matplotlib.pyplot as plt
+    from view_map import view_map
-
+    view_map(model.dem, model.lakes, model.rivers)
-    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

view_map.py

@@ -0,0 +1,48 @@
+#!/usr/bin/env python3
+
+import numpy as np
+import zlib
+import matplotlib.colors as mcol
+import matplotlib.pyplot as plt
+
+def view_map(dem, lakes, rivers):
+    plt.subplot(1,3,1)
+    plt.pcolormesh(dem, cmap='viridis')
+    plt.gca().set_aspect('equal', 'box')
+    plt.colorbar(orientation='horizontal')
+    plt.title('Raw elevation')
+
+    plt.subplot(1,3,2)
+    plt.pcolormesh(lakes, cmap='viridis')
+    plt.gca().set_aspect('equal', 'box')
+    plt.colorbar(orientation='horizontal')
+    plt.title('Lake surface elevation')
+
+    plt.subplot(1,3,3)
+    plt.pcolormesh(rivers, cmap='Blues', norm=mcol.LogNorm())
+    plt.gca().set_aspect('equal', 'box')
+    plt.colorbar(orientation='horizontal')
+    plt.title('Rivers flux')
+
+    plt.show()
+
+if __name__ == "__main__":
+    import sys
+    import os
+    if len(sys.argv) > 1:
+        os.chdir(sys.argv[1])
+
+    def load_map(name, dtype, shape):
+        dtype = np.dtype(dtype)
+        with open(name, 'rb') as f:
+            data = f.read()
+        if len(data) < shape[0]*shape[1]*dtype.itemsize:
+            data = zlib.decompress(data)
+        return np.frombuffer(data, dtype=dtype).reshape(shape)
+
+    shape = np.loadtxt('size', dtype='u4')
+    dem = load_map('dem', '>i2', shape)
+    lakes = load_map('lakes', '>i2', shape)
+    rivers = load_map('rivers', '>u4', shape)
+
+    view_map(dem, lakes, rivers)