Updated screenshot

Switching to full Lua incredibly simplified the amount of stuff needed in README!

README.md

@@ -1,114 +1,40 @@
 # Map Generator with Rivers
-`mapgen_rivers v0.0` by Gaël de Sailly.
+`mapgen_rivers v1.0` by Gaël de Sailly.
 
-Procedural map generator for Minetest 5.x. It aims to create realistic and nice-looking landscapes for the game, focused on river networks. It is based on algorithms modelling water flow and river erosion at a broad scale, similar to some used by researchers in Earth Sciences. It is taking some inspiration from [Fastscape](https://github.com/fastscape-lem/fastscape).
+Semi-procedural map generator for Minetest 5.x. It aims to create realistic and nice-looking landscapes for the game, focused on river networks. It is based on algorithms modelling water flow and river erosion at a broad scale, similar to some used by researchers in Earth Sciences. It is taking some inspiration from [Fastscape](https://github.com/fastscape-lem/fastscape).
 
 Its main particularity compared to conventional Minetest mapgens is that rivers that flow strictly downhill, and combine together to form wider rivers, until they reach the sea. Another notable feature is the possibility of large lakes above sea level.
 
-![Screenshot](https://user-images.githubusercontent.com/6905002/98825953-6289d980-2435-11eb-9e0b-704a95663ce0.png)
+![Screenshot](https://content.minetest.net/uploads/fff09f2269.png)
 
-**Important to know**: Unlike most other Minetest mods, it does not contain standalone Lua code, but does part of its processing with a separate Python program (included).
-- The Python part does pre-processing: it creates large-scale terrain data and applies landscape evolution algorithms, then outputs a grid of data in the mod's or world's folder. The grid is typically an array of 1000x1000 points of data, each of them representing a cell (by default 12x12 nodes). This pre-processing is long and should be run in advance.
-- The Lua part does actual map generation on Minetest. It reads grid data, upscales it (by a factor 12 by default), and adds small-scale features.
+It used to be composed of a Python script doing pre-generation, and a Lua mod reading the pre-generation output and generating the map. The code has been rewritten in full Lua for version 1.0 (July 2021), and is now usable out-of-the-box as any other Minetest mod.
+
+# Author and license
+License: GNU LGPLv3.0
+
+Code: Gaël de Sailly
+Flow routing algorithm concept (in `terrainlib/rivermapper.lua`): Cordonnier, G., Bovy, B., & Braun, J. (2019). A versatile, linear complexity algorithm for flow routing in topographies with depressions. Earth Surface Dynamics, 7(2), 549-562.
 
 # Requirements
-Mod dependencies: `default` required, and [`biomegen`](https://github.com/Gael-de-Sailly/biomegen) optional.
-
-Map pre-generation requires Python 3 with the following libraries installed:
-- `numpy`, widely used library for numerical calculations
-- `scipy`, a library for advanced data treatments, that is used here for Gaussian filtering
-- `noise`, implementing Perlin/Simplex noises
-
-Also, the following are optional (for map preview)
-- `matplotlib`, a famous library for graphical plotting
-- `colorcet` if you absolutely need better colormaps for preview :-)
-
-They are all commonly found on `pip` or `conda` Python distributions.
+Mod dependencies: `default` required, and [`biomegen`](https://github.com/Gael-de-Sailly/biomegen) optional (provides biome system).
 
 # Installation
 This mod should be placed in the `mods/` directory of Minetest like any other mod.
 
 # Usage
-By default, the mod contains a demo 400x400 grid (so you can start the game directly), but it is recommended to run the pre-processing script to generate a new grid before world creation, if you can.
+It is recommended to use it **only in new worlds, with `singlenode` mapgen**. On first start, it runs pre-generation to produce a grid, from which the map will be generated. This usually takes a few seconds, but depending on custom settings this can grow considerably longer.
 
-1. Run the script `generate.py` to generate a grid, preferentially from inside the mod's directory, but you can also run it directly in a Minetest world. See next paragraph for details about parameters.
-```
-./generate.py
-```
-2. Start Minetest, create a world with `singlenode` mapgen, enable `mapgen_rivers` mod, and launch the game. If you generated a grid in the world directory, it will copy it. If not, it will use the demo grid.
+By default, it only generates a 15k x 15k map, centered around the origin. To obtain a bigger map, you can increase grid size and/or block size in settings, but this can be more ressource-intensive (as the map has to be loaded in full at pre-generation).
 
-## Parameters for `generate.py`
-For a basic use you do not need to append any argument:
-```
-./generate.py
-```
-By default this will produce a 1000x1000 grid and save it in `river_data/`. Expect a computing time of about 30 minutes.
+## Settings
+Settings can be found in Minetest in the `Settings` tab, `All settings` -> `Mods` -> `mapgen_rivers`.
 
-### Parameters and config files
-This pre-processing takes many parameters. Instead of asking all these parameters to the end user, they are grouped in `.conf` files for usability, but the script still allows to override individual settings.
-
-Generic usage:
-```
-./generate.py conf_file output_dir
-```
-
-- `conf_file`: Path to configuration file from which parameters should be read. If omitted, attempts to read in `terrain.conf`.
-- `output_dir`: Directory in which to save the grid data, defaults to `river_data/`. If it does not exist, it is created. If it already contains previous grid data, they are overwritten.
-
-#### Config files
-The mod currently includes 3 config files, providing different terrain styles:
-- `terrain_default.conf` generates the standard terrain, with highest elevations around 250 with sharp peaks, and otherwise hilly terrain.
-- `terrain_higher.conf` generates higher mountains (up to 400 nodes), and wider valleys.
-- `terrain_original.conf` provides a terrain similar to what was generated with the first release of `mapgen_rivers`.
-
-More work is needed to find better and more varied terrain styles.
-
-### Complete list of parameters
-Other parameters can be specified by `--parameter value`. Syntax `--parameter=value` is also supported.
-
-| Parameter     | Description | Example |
-|---------------|-------------|---------|
-| | **Generic parameters** |
-| `mapsize`     | Size of the grid, in number of cells per edge. Usually `1000`, so to have 1000x1000 cells, the grid will have 1001x1001 nodes. Note that the grid is upscaled 12x in the game (this ratio can be changed), so that a `mapsize` of 1000 will result in a 12000x12000 map by default. | `--mapsize 1000` |
-| `sea_level`   | Height of the sea; height below which a point is considered under water even if it is not in a closed depression. | `--sea_level 1` |
-| | **Noise parameters** |
-| `scale`       | Horizontal variation wavlength of the largest noise octave, in grid cells (equivalent to the `spread` of a `PerlinNoise`). | `--scale 400` |
-| `vscale`      | Elevation coefficient, determines the approximate height difference between deepest seas and highest mountains. | `--vscale 300` |
-| `offset`      | Offset of the noise, will determine mean elevation. | `--offset 0` |
-| `persistence` | Relative height of smaller noise octaves compared to bigger ones. | `--persistence 0.6` |
-| `lacunarity`  | Relative reduction of wavelength between octaves. If `lacunarity`×`persistence` is larger than 1 (usual case), smaller octaves result in higher slopes than larger ones. This case is interesting for rivers networks because slopes determine rivers position. | `--lacunarity 2` |
-| | **Landscape evolution parameters**|
-| `K`           | Abstract erosion constant. Increasing it will increase erosive intensity. | `--K 1` |
-| `m`           | Parameter representing the influence of river flux on erosion. For `m=0`, small and big rivers are equal contributors to erosion. For `m=1` the erosive capability is proportional to river flux (assumed to be catchment area). Usual values: `0.25`-`0.60`. Be careful, this parameter is *highly sensitive*. | `--m 0.35` |
-| `d`           | Diffusion coefficient acting on sea/lake floor. Usual values `0`-`1`. | `--d 0.2` |
-| `flex_radius` | Flexure radius. Wavelength over which loss/gain of mass is compensated by uplift/subsidence. This ensures that mountain ranges will not get eventually flattened by erosion, and that an equilibrium is reached. Geologically speaking, this implements [isostatic rebound](https://en.wikipedia.org/wiki/Isostasy). | `--flex_radius 20` |
-| `time`        | Simulated time of erosion modelling, in abstract units. | `--time 10` |
-| `niter`       | Number of iterations. Each iteration represents a time `time/niter`. | `--niter 10` |
-| `sea_level_variations` | Amplitude of sea level variations throughout the simulation (if any). | `--sea_level_variations 10` |
-| `sea_level_variations_time` | Characteristic time of variation for sea level, in the same units than `time`. Increasing it will result in slower variations between iterations. | `--sea_level_variations_time 1` |
-| `flow_method` | Algorithm used for local flow calculation. Possible values are `steepest` (every node flows toward the steepest neighbour when possible), and `semirandom` (default, flow direction is determined randomly between lower neighbours, with lowest ones having greater probability). | `--flow_method semirandom` |
-| | **Alternatives** |
-| `config`      | Another way to specify configuration file | `--config terrain_higher.conf` |
-| `output`      | Another way to specify output dir | `--output ~/.minetest/worlds/my_world/river_data` |
-
-### Example
-```
-./generate.py terrain_higher.conf --mapsize 700 --K 0.4 --m 0.5
-```
-Reads parameters in `terrain_higher.conf`, and will generate a 700x700 grid using custom values for `K` and `m`.
+Most settings are world-specific and a copy is made in `mapgen_rivers.conf` in the world folder, during world first use, which means that further modification of global settings will not alter existing worlds.
 
 ## Map preview
-If you have `matplotlib` installed, `generate.py` will automatically show the grid aspect in real time during the erosion simulation.
+The Python script `view_map.py` can display the full map. You need to have Python 3 installed, as well as the libraries `numpy`, `matplotlib`, and optionally `colorcet`. For `conda` users, an `environment.yml` file is provided.
 
-There is also a script to view a generated map afterwards: `view_map.py`. Its syntax is the following:
+It can be run from command line by passing the world folder. Example:
 ```
-./view_map.py grid blocksize
-```
-
-- `grid` is the path to the grid directory to view. For example `river_data/`.
-- `blocksize` is the size at which 1 grid cell will be upscaled, in order to match game coordinates. If you use default settings, use `12`.
-
-Example:
-```
-./view_map.py river_data 12
+./view_map.py ~/.minetest/worlds/test_mg_rivers
 ```

compatibility.lua

@@ -0,0 +1,34 @@
+local function fix_min_catchment(settings, is_global)
+	local prefix = is_global and "mapgen_rivers_" or ""
+
+	local min_catchment = settings:get(prefix.."min_catchment")
+	if min_catchment then
+		min_catchment = tonumber(min_catchment)
+		local blocksize = tonumber(settings:get(prefix.."blocksize") or 15)
+		settings:set(prefix.."min_catchment", tonumber(min_catchment) * blocksize*blocksize)
+		local max_catchment = settings:get(prefix.."max_catchment")
+		if max_catchment then
+			max_catchment = tonumber(max_catchment)
+			local wpower = math.log(2*blocksize)/math.log(max_catchment/min_catchment)
+			settings:set(prefix.."river_widening_power", wpower)
+		end
+	end
+end
+
+local function fix_compatibility_minetest(settings)
+	local previous_version = settings:get("mapgen_rivers_version") or "0.0"
+
+	if previous_version == "0.0" then
+		fix_min_catchment(settings, true)
+	end
+end
+
+local function fix_compatibility_mapgen_rivers(settings)
+	local previous_version = settings:get("version") or "0.0"
+
+	if previous_version == "0.0" then
+		fix_min_catchment(settings, false)
+	end
+end
+
+return fix_compatibility_minetest, fix_compatibility_mapgen_rivers

demo_data/size

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

environment.yml

@@ -0,0 +1,10 @@
+name: mapgen_rivers
+
+channels:
+    - conda-forge
+
+dependencies:
+    - python
+    - matplotlib
+    - numpy
+    - colorcet

generate.py

@@ -1,179 +0,0 @@
-#!/usr/bin/env python3
-
-import numpy as np
-from noise import snoise2
-import os
-import sys
-
-import terrainlib
-
-def noisemap(X, Y, scale=0.01, vscale=1.0, offset=0.0, log=False, **params):
-    # Determine noise offset randomly
-    xbase = np.random.randint(8192)-4096
-    ybase = np.random.randint(8192)-4096
-
-    if log:
-        vscale /= offset
-
-    # Generate the noise
-    n = np.zeros((X, Y))
-    for x in range(X):
-        for y in range(Y):
-            n[x,y] = snoise2(x/scale + xbase, y/scale + ybase, **params)
-
-    if log:
-        return np.exp(n*vscale) * offset
-    else:
-        return n*vscale + offset
-
-### PARSE COMMAND-LINE ARGUMENTS
-argc = len(sys.argv)
-
-config_file = 'terrain_default.conf'
-output_dir = 'river_data'
-params_from_args = {}
-i = 1 # Index of arguments
-j = 1 # Number of 'orphan' arguments (the ones that are not preceded by '--something')
-while i < argc:
-    arg = sys.argv[i]
-    if arg[:2] == '--':
-        pname = arg[2:]
-        v = None
-        split = pname.split('=', maxsplit=1)
-        if len(split) == 2:
-            pname, v = split
-            i += 1
-        elif i+1 < argc:
-            v = sys.argv[i+1]
-            i += 2
-
-        if v is not None:
-            if pname == 'config':
-                config_file = v
-            elif pname == 'output':
-                output_dir = v
-            else:
-                params_from_args[pname] = v
-    else:
-        if j == 1:
-            config_file = arg
-        elif j == 2:
-            output_dir = arg
-        i += 1
-        j += 1
-
-print(config_file, output_dir)
-
-params = terrainlib.read_config_file(config_file)
-params.update(params_from_args) # Params given from args prevail against conf file
-
-### READ SETTINGS
-def get_setting(name, default):
-    if name in params:
-        return params[name]
-    return default
-
-mapsize = int(get_setting('mapsize', 1000))
-scale = float(get_setting('scale', 400.0))
-vscale = float(get_setting('vscale', 300.0))
-offset = float(get_setting('offset', 0.0))
-persistence = float(get_setting('persistence', 0.6))
-lacunarity = float(get_setting('lacunarity', 2.0))
-
-K = float(get_setting('K', 0.5))
-m = float(get_setting('m', 0.5))
-d = float(get_setting('d', 0.5))
-sea_level = float(get_setting('sea_level', 0.0))
-sea_level_variations = float(get_setting('sea_level_variations', 0.0))
-sea_level_variations_time = float(get_setting('sea_level_variations_time', 1.0))
-flex_radius = float(get_setting('flex_radius', 20.0))
-flow_method = get_setting('flow_method', 'semirandom')
-
-time = float(get_setting('time', 10.0))
-niter = int(get_setting('niter', 10))
-
-### MAKE INITIAL TOPOGRAPHY
-n = np.zeros((mapsize+1, mapsize+1))
-
-# Set noise parameters
-params = {
-    "offset" : offset,
-    "vscale" : vscale,
-    "scale" : scale,
-    "octaves" : int(np.ceil(np.log2(mapsize)))+1,
-    "persistence" : persistence,
-    "lacunarity" : lacunarity,
-}
-
-params_sealevel = {
-    "octaves" : 1,
-    "persistence" : 1,
-    "lacunarity" : 2,
-}
-
-if sea_level_variations != 0.0:
-    sea_ybase = np.random.randint(8192)-4096
-    sea_level_ref = snoise2(time * (1-1/niter) / sea_level_variations, sea_ybase, **params_sealevel) * sea_level_variations
-    params['offset'] -= (sea_level_ref + sea_level)
-
-n = noisemap(mapsize+1, mapsize+1, **params)
-
-### COMPUTE LANDSCAPE EVOLUTION
-# Initialize landscape evolution model
-print('Initializing model')
-model = terrainlib.EvolutionModel(n, K=K, m=m, d=d, sea_level=sea_level, flex_radius=flex_radius, flow_method=flow_method)
-terrainlib.update(model.dem, model.lakes, t=5, sea_level=model.sea_level, title='Initializing...')
-
-dt = time/niter
-
-# Run the model's processes: the order in which the processes are run is arbitrary and could be changed.
-
-for i in range(niter):
-    disp_niter = 'Iteration {:d} of {:d}...'.format(i+1, niter)
-    if sea_level_variations != 0:
-        model.sea_level = snoise2((i*dt)/sea_level_variations_time, sea_ybase, **params_sealevel) * sea_level_variations - sea_level_ref
-    terrainlib.update(model.dem, model.lakes, sea_level=model.sea_level, title=disp_niter)
-    print(disp_niter)
-    print('Diffusion')
-    model.diffusion(dt)
-    print('Flow calculation')
-    model.calculate_flow()
-    terrainlib.update(model.dem, model.lakes, sea_level=model.sea_level, title=disp_niter)
-    print('Advection')
-    model.advection(dt)
-    print('Isostatic equilibration')
-    model.adjust_isostasy()
-
-print('Last flow calculation')
-model.calculate_flow()
-
-print('Done!')
-
-# Twist the grid
-bx, by = terrainlib.make_bounds(model.dirs, model.rivers)
-offset_x, offset_y = terrainlib.twist(bx, by, terrainlib.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)
-
-### SAVE OUTPUT
-if not os.path.isdir(output_dir):
-    os.mkdir(output_dir)
-os.chdir(output_dir)
-# Save the files
-terrainlib.save(model.dem, 'dem', dtype='>i2')
-terrainlib.save(model.lakes, 'lakes', dtype='>i2')
-terrainlib.save(offset_x, 'offset_x', dtype='i1')
-terrainlib.save(offset_y, 'offset_y', dtype='i1')
-
-terrainlib.save(model.dirs, 'dirs', dtype='u1')
-terrainlib.save(model.rivers, 'rivers', dtype='>u4')
-
-with open('size', 'w') as sfile:
-    sfile.write('{:d}\n{:d}'.format(mapsize+1, mapsize+1))
-
-terrainlib.stats(model.dem, model.lakes)
-print()
-print('Grid is ready for use!')
-terrainlib.plot(model.dem, model.lakes, title='Final grid, ready for use!')

heightmap.lua

@@ -1,11 +1,10 @@
-local modpath = minetest.get_modpath(minetest.get_current_modname()) .. '/'
+local modpath = mapgen_rivers.modpath
 
 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
+local sea_level = mapgen_rivers.settings.sea_level
+local riverbed_slope = mapgen_rivers.settings.riverbed_slope * mapgen_rivers.settings.blocksize
 
 local MAP_BOTTOM = -31000
 

init.lua

@@ -1,15 +1,15 @@
 mapgen_rivers = {}
 
 local modpath = minetest.get_modpath(minetest.get_current_modname()) .. '/'
+mapgen_rivers.modpath = modpath
+mapgen_rivers.world_data_path = minetest.get_worldpath() .. '/river_data/'
 
 dofile(modpath .. 'settings.lua')
 
-local blocksize = mapgen_rivers.blocksize
-local sea_level = mapgen_rivers.sea_level
-local riverbed_slope = mapgen_rivers.riverbed_slope
-local elevation_chill = mapgen_rivers.elevation_chill
-local use_distort = mapgen_rivers.distort
-local use_biomes = mapgen_rivers.biomes
+local sea_level = mapgen_rivers.settings.sea_level
+local elevation_chill = mapgen_rivers.settings.elevation_chill
+local use_distort = mapgen_rivers.settings.distort
+local use_biomes = mapgen_rivers.settings.biomes
 local use_biomegen_mod = use_biomes and minetest.global_exists('biomegen')
 use_biomes = use_biomes and not use_biomegen_mod
 
@@ -38,7 +38,13 @@ local noise_heat_blend_map = {}
 local mapsize
 local init = false
 
+local sumtime = 0
+local sumtime2 = 0
+local ngen = 0
+
 local function generate(minp, maxp, seed)
+	print(("[mapgen_rivers] Generating from %s to %s"):format(minetest.pos_to_string(minp), minetest.pos_to_string(maxp)))
+
 	local chulens = {
 		x = maxp.x-minp.x+1,
 		y = maxp.y-minp.y+1,
@@ -63,6 +69,7 @@ local function generate(minp, maxp, seed)
 		init = true
 	end
 
+	local t0 = os.clock()
 	local minp2d = {x=minp.x, y=minp.z}
 	if use_distort then
 		noise_x_obj:get_3d_map_flat(minp, noise_x_map)
@@ -130,7 +137,7 @@ local function generate(minp, maxp, seed)
 	local incrZ = mapsize.x*mapsize.y - mapsize.x*incrX - mapsize.x*mapsize.y*incrY
 
 	local i2d = 1
-	
+
 	for z = minp.z, maxp.z do
 		for x = minp.x, maxp.x do
 			local ivm = a:index(x, minp.y, z)
@@ -221,6 +228,18 @@ local function generate(minp, maxp, seed)
 	vm:calc_lighting()
 	vm:update_liquids()
 	vm:write_to_map()
+	local t1 = os.clock()
+
+	local t = t1-t0
+	ngen = ngen + 1
+	sumtime = sumtime + t
+	sumtime2 = sumtime2 + t*t
+	print(("[mapgen_rivers] Done in %5.3f s"):format(t))
 end
 
 minetest.register_on_generated(generate)
+minetest.register_on_shutdown(function()
+	local avg = sumtime / ngen
+	local std = math.sqrt(sumtime2/ngen - avg*avg)
+	print(("[mapgen_rivers] Mapgen statistics:\n- Mapgen calls: %4d\n- Mean time: %5.3f s\n- Standard deviation: %5.3f s"):format(ngen, avg, std))
+end)

load.lua

@@ -1,31 +1,100 @@
-local worldpath = minetest.get_worldpath() .. "/river_data/"
+local worldpath = mapgen_rivers.world_data_path
 
-local function load_map(filename, bytes, signed, size)
+function mapgen_rivers.load_map(filename, bytes, signed, size, converter)
 	local file = io.open(worldpath .. filename, 'rb')
 	local data = file:read('*all')
 	if #data < bytes*size then
 		data = minetest.decompress(data)
 	end
+	local sbyte = string.byte
 
 	local map = {}
 
 	for i=1, size do
-		local i0, i1 = (i-1)*bytes+1, i*bytes
+		local i0 = (i-1)*bytes+1
 		local elements = {data:byte(i0, i1)}
-		local n = elements[1]
+		local n = sbyte(data, i0)
 		if signed and n >= 128 then
 			n = n - 256
 		end
 
-		for j=2, bytes do
-			n = n*256 + elements[j]
+		for j=1, bytes-1 do
+			n = n*256 + sbyte(data, i0+j)
 		end
 
 		map[i] = n
 	end
 	file:close()
 
+	if converter then
+		for i=1, size do
+			map[i] = converter(map[i])
+		end
+	end
+
 	return map
 end
 
-return load_map
+local sbyte = string.byte
+
+local loader_mt = {
+	__index = function(loader, i)
+		local file = loader.file
+		local bytes = loader.bytes
+		file:seek('set', (i-1)*bytes)
+		local strnum = file:read(bytes)
+
+		local n = sbyte(strnum, 1)
+		if loader.signed and n >= 128 then
+			n = n - 256
+		end
+
+		for j=2, bytes do
+			n = n*256 + sbyte(strnum, j)
+		end
+
+		if loader.conv then
+			n = loader.conv(n)
+		end
+		loader[i] = n
+		return n
+	end,
+}
+
+function mapgen_rivers.interactive_loader(filename, bytes, signed, size, converter)
+	local file = io.open(worldpath .. filename, 'rb')
+	if file then
+		minetest.register_on_shutdown(function()
+			file:close()
+		end)
+		converter = converter or false
+		return setmetatable({file=file, bytes=bytes, signed=signed, size=size, conv=converter}, loader_mt)
+	end
+end
+
+function mapgen_rivers.write_map(filename, data, bytes)
+	local size = #data
+	local file = io.open(worldpath .. filename, 'wb')
+	local mfloor = math.floor
+	local schar = string.char
+	local upack = unpack
+
+	local bytelist = {}
+	for j=1, bytes do
+		bytelist[j] = 0
+	end
+
+	for i=1, size do
+		local n = mfloor(data[i])
+		data[i] = n
+		for j=bytes, 2, -1 do
+			bytelist[j] = n % 256
+			n = mfloor(n / 256)
+		end
+		bytelist[1] = n % 256
+
+		file:write(schar(upack(bytelist)))
+	end
+
+	file:close()
+end

noises.lua

@@ -1,5 +1,18 @@
+local def_setting = mapgen_rivers.define_setting
+
 mapgen_rivers.noise_params = {
-	distort_x = {
+	base = def_setting('np_base', 'noise', {
+		offset = 0,
+		scale = 300,
+		seed = 2469,
+		octaves = 8,
+		spread = {x=2048, y=2048, z=2048},
+		persist = 0.6,
+		lacunarity = 2,
+		flags = "eased",
+	}),
+
+	distort_x = def_setting('np_distort_x', 'noise', {
 		offset = 0,
 		scale = 1,
 		seed = -4574,
@@ -7,9 +20,9 @@ mapgen_rivers.noise_params = {
 		octaves = 3,
 		persistence = 0.75,
 		lacunarity = 2,
-	},
+	}),
 
-	distort_z = {
+	distort_z = def_setting('np_distort_z', 'noise', {
 		offset = 0,
 		scale = 1,
 		seed = -7940,
@@ -17,9 +30,9 @@ mapgen_rivers.noise_params = {
 		octaves = 3,
 		persistence = 0.75,
 		lacunarity = 2,
-	},
+	}),
 
-	distort_amplitude = {
+	distort_amplitude = def_setting('np_distort_amplitude', 'noise', {
 		offset = 0,
 		scale = 10,
 		seed = 676,
@@ -28,10 +41,40 @@ mapgen_rivers.noise_params = {
 		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
+-- Convert to number because Minetest API is not able to do it cleanly...
+for name, np in pairs(mapgen_rivers.noise_params) do
+	for field, value in pairs(np) do
+		if field ~= 'flags' and type(value) == 'string' then
+			np[field] = tonumber(value) or value
+		elseif field == 'spread' then
+			for dir, v in pairs(value) do
+				value[dir] = tonumber(v) or v
+			end
+		end
+	end
+end
+
+local heat = mapgen_rivers.noise_params.heat
+local base = mapgen_rivers.noise_params.base
+local settings = mapgen_rivers.settings
+heat.offset = heat.offset + settings.sea_level * settings.elevation_chill
+base.spread.x = base.spread.x / settings.blocksize
+base.spread.y = base.spread.y / settings.blocksize
+base.spread.z = base.spread.z / settings.blocksize
+
+for name, np in pairs(mapgen_rivers.noise_params) do
+	local lac = np.lacunarity or 2
+	if lac > 1 then
+		local omax = math.floor(math.log(math.min(np.spread.x, np.spread.y, np.spread.z)) / math.log(lac))+1
+		if np.octaves > omax then
+			print("[mapgen_rivers] Noise " .. name .. ": 'octaves' reduced to " .. omax)
+			np.octaves = omax
+		end
+	end
+end

polygons.lua

@@ -1,76 +1,98 @@
-local modpath = minetest.get_modpath(minetest.get_current_modname()) .. '/'
+local modpath = mapgen_rivers.modpath
 local mod_data_path = modpath .. 'river_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/'
+local world_data_path = mapgen_rivers.world_data_path
 minetest.mkdir(world_data_path)
 
-local load_map = dofile(modpath .. 'load.lua')
+dofile(modpath .. 'load.lua')
 
-local function copy_if_needed(filename)
-	local wfilename = world_data_path..filename
-	local wfile = io.open(wfilename, 'rb')
-	if wfile then
-		wfile:close()
-		return
-	end
-	local mfilename = mod_data_path..filename
-	local mfile = io.open(mfilename, 'rb')
-	local wfile = io.open(wfilename, 'wb')
-	wfile:write(mfile:read("*all"))
-	mfile:close()
-	wfile:close()
+mapgen_rivers.grid = {}
+
+local X = mapgen_rivers.settings.grid_x_size
+local Z = mapgen_rivers.settings.grid_z_size
+
+local function offset_converter(o)
+	return (o + 0.5) * (1/256)
 end
 
-copy_if_needed('size')
+local load_all = mapgen_rivers.settings.load_all
+
+-- Try to read file 'size'
 local sfile = io.open(world_data_path..'size', 'r')
-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
+local first_mapgen = true
+if sfile then
+	X, Z = tonumber(sfile:read('*l')), tonumber(sfile:read('*l'))
+	sfile:close()
+	first_mapgen = false
 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
+if first_mapgen then
+	-- Generate a map!!
+	local pregenerate = dofile(mapgen_rivers.modpath .. '/pregenerate.lua')
+	minetest.register_on_mods_loaded(function()
+		print('[mapgen_rivers] Generating grid')
+		pregenerate(load_all)
+
+		if load_all then
+			local offset_x = mapgen_rivers.grid.offset_x
+			local offset_y = mapgen_rivers.grid.offset_y
+			for i=1, X*Z do
+				offset_x[i] = offset_converter(offset_x[i])
+				offset_y[i] = offset_converter(offset_y[i])
+			end
+		end
+	end)
 end
 
--- To index a flat array representing a 2D map
+-- if data not already loaded
+if not (first_mapgen and load_all) then
+	local load_map
+	if load_all then
+		load_map = mapgen_rivers.load_map
+	else
+		load_map = mapgen_rivers.interactive_loader
+	end
+
+	minetest.register_on_mods_loaded(function()
+		if load_all then
+			print('[mapgen_rivers] Loading full grid')
+		else
+			print('[mapgen_rivers] Loading grid as interactive loaders')
+		end
+		local grid = mapgen_rivers.grid
+
+		grid.dem = load_map('dem', 2, true, X*Z)
+		grid.lakes = load_map('lakes', 2, true, X*Z)
+		grid.dirs = load_map('dirs', 1, false, X*Z)
+		grid.rivers = load_map('rivers', 4, false, X*Z)
+
+		grid.offset_x = load_map('offset_x', 1, true, X*Z, offset_converter)
+		grid.offset_y = load_map('offset_y', 1, true, X*Z, offset_converter)
+	end)
+end
+
+mapgen_rivers.grid.size = {x=X, y=Z}
+
 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
+local blocksize = mapgen_rivers.settings.blocksize
+local min_catchment = mapgen_rivers.settings.min_catchment
+local max_catchment = mapgen_rivers.settings.max_catchment
 
 local map_offset = {x=0, z=0}
-if mapgen_rivers.center then
+if mapgen_rivers.settings.center then
 	map_offset.x = blocksize*X/2
 	map_offset.z = blocksize*Z/2
 end
 
--- 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 min_catchment = mapgen_rivers.settings.min_catchment / (blocksize*blocksize)
+local wpower = mapgen_rivers.settings.river_widening_power
+local wfactor = 1/(2*blocksize * min_catchment^wpower)
 local function river_width(flow)
 	flow = math.abs(flow)
 	if flow < min_catchment then
@@ -81,22 +103,29 @@ local function river_width(flow)
 end
 
 local noise_heat -- Need a large-scale noise here so no heat blend
-local elevation_chill = mapgen_rivers.elevation_chill
+local elevation_chill = mapgen_rivers.settings.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 glaciers = mapgen_rivers.settings.glaciers
+local glacier_factor = mapgen_rivers.settings.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)
-	print("Generating polygon map")
-	print(minp.x, maxp.x, minp.z, maxp.z)
+
+	local grid = mapgen_rivers.grid
+	local dem = grid.dem
+	local lakes = grid.lakes
+	local dirs = grid.dirs
+	local rivers = grid.rivers
+
+	local offset_x = grid.offset_x
+	local offset_z = grid.offset_y
 
 	if not init then
 		if glaciers then
@@ -111,7 +140,6 @@ local function make_polygons(minp, maxp)
 	-- 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+map_offset.x)/blocksize - 0.5), 0), math.min(math.ceil((maxp.x+map_offset.x)/blocksize + 0.5), X-2)
 	local zpmin, zpmax = math.max(math.floor((minp.z+map_offset.z)/blocksize - 0.5), 0), math.min(math.ceil((maxp.z+map_offset.z)/blocksize + 0.5), Z-2)
-	print(xpmin, xpmax, zpmin, zpmax)
 
 	-- Iterate over the polygons
 	for xp = xpmin, xpmax do
@@ -133,9 +161,6 @@ local function make_polygons(minp, maxp)
 				(offset_z[iC]+zp+1) * blocksize - map_offset.z,
 				(offset_z[iD]+zp+1) * blocksize - map_offset.z,
 			}
-			if xp==xpmin and zp==zpmin then
-				print(xp, zp, poly_x[1], poly_z[1])
-			end
 			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 Z position (scanline algorithm)

pregenerate.lua

@@ -0,0 +1,81 @@
+local EvolutionModel = dofile(mapgen_rivers.modpath .. '/terrainlib_lua/erosion.lua')
+local twist = dofile(mapgen_rivers.modpath .. '/terrainlib_lua/twist.lua')
+
+local blocksize = mapgen_rivers.settings.blocksize
+local tectonic_speed = mapgen_rivers.settings.tectonic_speed
+
+local np_base = table.copy(mapgen_rivers.noise_params.base)
+
+local evol_params = mapgen_rivers.settings.evol_params
+
+local time = mapgen_rivers.settings.evol_time
+local time_step = mapgen_rivers.settings.evol_time_step
+local niter = math.ceil(time/time_step)
+time_step = time / niter
+
+local function pregenerate(keep_loaded)
+	local grid = mapgen_rivers.grid
+	local size = grid.size
+
+	local seed = tonumber(minetest.get_mapgen_setting("seed"))
+	np_base.seed = (np_base.seed or 0) + seed
+
+	local nobj_base = PerlinNoiseMap(np_base, {x=size.x, y=1, z=size.y})
+
+	local dem = nobj_base:get_3d_map_flat({x=0, y=0, z=0})
+	dem.X = size.x
+	dem.Y = size.y
+
+	local model = EvolutionModel(evol_params)
+	model.dem = dem
+	local ref_dem = model:define_isostasy(dem)
+
+	local tectonic_step = tectonic_speed * time_step
+	collectgarbage()
+	for i=1, niter do
+		print("[mapgen_rivers] Iteration " .. i .. " of " .. niter)
+
+		model:diffuse(time_step)
+		model:flow()
+		model:erode(time_step)
+		if i < niter then
+			if tectonic_step ~= 0 then
+				nobj_base:get_3d_map_flat({x=0, y=tectonic_step*i, z=0}, ref_dem)
+			end
+			model:isostasy()
+		end
+
+		collectgarbage()
+	end
+	model:flow()
+
+	local mfloor = math.floor
+	local mmin, mmax = math.min, math.max
+	local offset_x, offset_y = twist(model.dirs, model.rivers, 5)
+	for i=1, size.x*size.y do
+		offset_x[i] = mmin(mmax(offset_x[i]*256, -128), 127)
+		offset_y[i] = mmin(mmax(offset_y[i]*256, -128), 127)
+	end
+
+	mapgen_rivers.write_map('dem', model.dem, 2)
+	mapgen_rivers.write_map('lakes', model.lakes, 2)
+	mapgen_rivers.write_map('dirs', model.dirs, 1)
+	mapgen_rivers.write_map('rivers', model.rivers, 4)
+	mapgen_rivers.write_map('offset_x', offset_x, 1)
+	mapgen_rivers.write_map('offset_y', offset_y, 1)
+	local sfile = io.open(mapgen_rivers.world_data_path .. 'size', "w")
+	sfile:write(size.x..'\n'..size.y)
+	sfile:close()
+
+	if keep_loaded then
+		grid.dem = model.dem
+		grid.lakes = model.lakes
+		grid.dirs = model.dirs
+		grid.rivers = model.rivers
+		grid.offset_x = offset_x
+		grid.offset_y = offset_y
+	end
+	collectgarbage()
+end
+
+return pregenerate

readconfig.py

@@ -0,0 +1,42 @@
+def read_conf_file(filename):
+    f = open(filename, 'r')
+    return read_conf(f)
+
+def read_conf(f, end_tag=None):
+    conf = {}
+    while True:
+        line = f.readline()
+        if len(line) == 0:
+            return conf
+        line = line.strip()
+        if line == end_tag:
+            return conf
+        if len(line) == 0 or line[0] == '#':
+            continue
+
+        eqpos = line.find('=')
+        if eqpos < 0:
+            continue
+            
+        name, value = line[:eqpos].rstrip(), line[eqpos+1:].lstrip()
+        if value == '{':
+            # Group
+            conf[name] = read_conf(f, end_tag='}')
+
+        elif value == '"""':
+            # Multiline
+            conf[value] = read_multiline(f)
+
+        else:
+            conf[name] = value
+
+def read_multiline(f):
+    mline = ''
+    while True:
+        line = f.readline()
+        if len(line) == 0:
+            return mline
+        line = line.strip()
+        if line == '"""':
+            return mline
+        mline += line + '\n'

settings.lua

@@ -1,58 +1,95 @@
-local storage = minetest.get_mod_storage()
-local settings = minetest.settings
+local mtsettings = minetest.settings
+local mgrsettings = Settings(minetest.get_worldpath() .. '/mapgen_rivers.conf')
 
-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
+mapgen_rivers.version = "1.0"
+
+local previous_version_mt = mtsettings:get("mapgen_rivers_version") or "0.0"
+local previous_version_mgr = mgrsettings:get("version") or "0.0"
+
+if mapgen_rivers.version ~= previous_version_mt or mapgen_rivers.version ~= previous_version_mgr then
+	local compat_mt, compat_mgr = dofile(minetest.get_modpath(minetest.get_current_modname()) .. "/compatibility.lua")
+	if mapgen_rivers.version ~= previous_version_mt then
+		compat_mt(mtsettings)
 	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)
-		else
-			storage:set_string(key, conf_val)
-			if dtype == "bool" then
-				conf_val = conf_val == 'true'
-			end
-		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
+	if mapgen_rivers.version ~= previous_version_mgr then
+		compat_mgr(mgrsettings)
 	end
 end
 
-mapgen_rivers.center = get_settings('center', 'bool', false)
-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)
-mapgen_rivers.biomes = get_settings('biomes', 'bool', true)
-mapgen_rivers.glaciers = get_settings('glaciers', 'bool', false)
-mapgen_rivers.glacier_factor = get_settings('glacier_factor', 'float', 8)
-mapgen_rivers.elevation_chill = get_settings('elevation_chill', 'float', 0.25)
+mtsettings:set("mapgen_rivers_version", mapgen_rivers.version)
+mgrsettings:set("version", mapgen_rivers.version)
+
+function mapgen_rivers.define_setting(name, dtype, default)
+	if dtype == "number" or dtype == "string" then
+		local v = mgrsettings:get(name)
+		if v == nil then
+			v = mtsettings:get('mapgen_rivers_' .. name)
+			if v == nil then
+				v = default
+			end
+			mgrsettings:set(name, v)
+		end
+		if dtype == "number" then
+			return tonumber(v)
+		else
+			return v
+		end
+	elseif dtype == "bool" then
+		local v = mgrsettings:get_bool(name)
+		if v == nil then
+			v = mtsettings:get_bool('mapgen_rivers_' .. name)
+			if v == nil then
+				v = default
+			end
+			mgrsettings:set_bool(name, v)
+		end
+		return v
+	elseif dtype == "noise" then
+		local v = mgrsettings:get_np_group(name)
+		if v == nil then
+			v = mtsettings:get_np_group('mapgen_rivers_' .. name)
+			if v == nil then
+				v = default
+			end
+			mgrsettings:set_np_group(name, v)
+		end
+		return v
+	end
+end
+
+local def_setting = mapgen_rivers.define_setting
+
+mapgen_rivers.settings = {
+	center = def_setting('center', 'bool', true),
+	blocksize = def_setting('blocksize', 'number', 15),
+	sea_level = tonumber(minetest.get_mapgen_setting('water_level')),
+	min_catchment = def_setting('min_catchment', 'number', 3600),
+	river_widening_power = def_setting('river_widening_power', 'number', 0.5),
+	riverbed_slope = def_setting('riverbed_slope', 'number', 0.4),
+	distort = def_setting('distort', 'bool', true),
+	biomes = def_setting('biomes', 'bool', true),
+	glaciers = def_setting('glaciers', 'bool', false),
+	glacier_factor = def_setting('glacier_factor', 'number', 8),
+	elevation_chill = def_setting('elevation_chill', 'number', 0.25),
+
+	grid_x_size = def_setting('grid_x_size', 'number', 1000),
+	grid_z_size = def_setting('grid_z_size', 'number', 1000),
+	evol_params = {
+		K = def_setting('river_erosion_coef', 'number', 0.5),
+		m = def_setting('river_erosion_power', 'number', 0.4),
+		d = def_setting('diffusive_erosion', 'number', 0.5),
+		compensation_radius = def_setting('compensation_radius', 'number', 50),
+	},
+	tectonic_speed = def_setting('tectonic_speed', 'number', 70),
+	evol_time = def_setting('evol_time', 'number', 10),
+	evol_time_step = def_setting('evol_time_step', 'number', 1),
+
+	load_all = mtsettings:get_bool('mapgen_rivers_load_all')
+}
+
+local function write_settings()
+	mgrsettings:write()
+end
+
+minetest.register_on_mods_loaded(write_settings)
+minetest.register_on_shutdown(write_settings)

settingtypes.txt

@@ -1,26 +1,31 @@
 # File containing all settings for 'mapgen_rivers' mod.
 
 #    Whether the map should be centered at x=0, z=0.
-mapgen_rivers_center (Center map) bool false
+mapgen_rivers_center (Center map) bool true
 
 #    Represents horizontal map scale. Every cell of the grid will be upscaled to
 #    a square of this size.
 #    For example if the grid size is 1000x1000 and block size is 12,
-#    the actual size of the map will be 15000.
-mapgen_rivers_blocksize (Block size) float 12.0 2.0 40.0
+#    the actual size of the map will be 12000.
+mapgen_rivers_blocksize (Block size) float 15.0 2.0 100.0
 
-#    Sea level used by mapgen_rivers
-mapgen_rivers_sea_level (Sea level) int 1
+#    X size of the grid being generated
+#    Actual size of the map is grid_x_size * blocksize
+mapgen_rivers_grid_x_size (Grid X size) int 1000 50 5000
 
-#    Minimal catchment area for a river to be drawn, in grid cells
-#    (1 cell = blocksize x blocksize).
+#    Z size of the grid being generated
+#    Actual size of the map is grid_z_size * blocksize
+mapgen_rivers_grid_z_size (Grid Z size) int 1000 50 5000
+
+#    Minimal catchment area for a river to be drawn, in square nodes
 #    Lower value means bigger river density
-mapgen_rivers_min_catchment (Minimal catchment area) float 25.0 1.0 1000.0
+mapgen_rivers_min_catchment (Minimal catchment area) float 3600.0 100.0 1000000.0
 
-#    Catchment area in grid cells (1 grid cell = blocksize x blocksize)
-#    at which rivers reach their maximal width of 2*blocksize.
+#    Coefficient describing how rivers widen when merging.
+#    Riwer width is a power law W = a*D^p. D is river flow and p is this parameter.
 #    Higher value means a river needs to receive more tributaries to grow in width.
-mapgen_rivers_max_catchment (Maximal catchment area) float 40000.0 1000.0 10000000.0
+#    Note that a river can never exceed 2*blocksize.
+mapgen_rivers_river_widening_power (River widening power) float 0.5 0.0 1.0
 
 #    Lateral slope of the riverbed.
 #    Higher value means deeper rivers.
@@ -49,4 +54,64 @@ mapgen_rivers_glacier_widening_factor (Glacier widening factor) float 8.0 1.0 20
 #    This results in mountains being more covered by snow.
 mapgen_rivers_elevation_chill (Elevation chill) float 0.25 0.0 5.0
 
-# Noises: to be added. For now they are hardcoded.
+#    If enabled, loads all grid data in memory at init time.
+#    If disabled, data will be loaded on request and cached in memory.
+#    It's recommended to disable it for very large maps (> 2000 grid nodes or so)
+mapgen_rivers_load_all (Load all data in memory) bool false
+
+[Landscape evolution parameters]
+
+#    Modelled landscape evolution time, in arbitrary units
+mapgen_rivers_evol_time (Landscape evolution time) float 10.0 0.0 100.0
+
+#    Model time steps in arbitrary units
+#    Smaller values will result in more time steps to be necessary to
+#    complete the simulation, taking more time.
+mapgen_rivers_evol_time_step (Landscape evolution time step) float 1.0 0.0 50.0
+
+#    To adjust river erosion proportionnally.
+#    This type of erosion acts by deepening the valleys.
+mapgen_rivers_river_erosion_coef (River erosion coefficient) float 0.5 0.0 10.0
+
+#    Represents how much river erosion depends on river flow (catchment area).
+#    Catchment area is elevated to this power.
+#    Extreme cases: 0.0 -> All rivers have the same erosive capabilities
+#    1.0 -> Erosion is proportional to river flow
+#    Reasonable values are generally between 0.4 and 0.7.
+#
+#    This parameter is extremely sensitive, and changes may require to adjust
+#    'river_erosion_coef' as well.
+mapgen_rivers_river_erosion_power (River erosion power) float 0.4 0.0 1.0
+
+#    Intensity of diffusive erosion.
+#    Smoothes peaks and valleys, and tends to prevent sharp cliffs from forming.
+mapgen_rivers_diffusive_erosion (Diffusive erosion) float 0.5 0.0 10.0
+
+#    Radius of compensation for isostatic/tectonic processes
+#    Tectonic uplift forces will have a diffuse effect over this radius
+mapgen_rivers_compensation_radius (Tectonic compensation radius) float 50 1.0 1000.0
+
+#    Speed of evolution of tectonic conditions between steps
+#    Higher values means tectonics will be very different from one step to the other,
+#    resulting in geologically unstable and more varied landforms (plateau, gorge, lake...)
+mapgen_rivers_tectonic_speed (Tectonic speed) float 70 0 10000
+
+[Noises]
+
+#    Y level of terrain at a very large scale. Only used during pre-generation.
+#    X and Z axes correspond to map's X and Z directions, and Y axis is time.
+#    Successive XZ slices of this noise represent successive tectonic states.
+mapgen_rivers_np_base (Terrain base noise) noise_params_3d 0, 300, (2048, 2048, 2048), 2469, 8, 0.6, 2.0, eased
+
+#    This noise will shear the terrain on the X axis,
+#    to break the regularity of the river grid.
+mapgen_rivers_np_distort_x (X-axis distorsion noise) noise_params_3d 0, 1, (64, 32, 64), -4574, 3, 0.75, 2.0
+
+#    This noise will shear the terrain on the Z axis,
+#    to break the regularity of the river grid.
+mapgen_rivers_np_distort_z (Z-axis distorsion noise) noise_params_3d 0, 1, (64, 32, 64), -7940, 3, 0.75, 2.0
+
+#    Amplitude of the distorsion.
+#    Too small values may leave the grid pattern apparent,
+#    and too high values could make the terrain insanely twisted.
+mapgen_rivers_np_distort_amplitude (Distorsion amplitude noise) noise_params_2d 0, 10, (1024, 1024, 1024), 676, 5, 0.5, 2.0, absvalue

terrain_default.conf

@@ -1,17 +0,0 @@
-mapsize = 1000
-scale = 400
-vscale = 300
-offset = 0
-persistence = 0.6
-lacunarity = 2.0
-
-K = 0.5
-m = 0.5
-d = 0.5
-sea_level = 0
-sea_level_variations = 8
-sea_level_variations_time = 2
-flex_radius = 20
-
-time = 10
-niter = 10

terrain_higher.conf

@@ -1,17 +0,0 @@
-mapsize = 1000
-scale = 400
-vscale = 600
-offset = 0
-persistence = 0.65
-lacunarity = 2.0
-
-K = 0.5
-m = 0.45
-d = 0.55
-sea_level = 0
-sea_level_variations = 12
-sea_level_variations_time = 2
-flex_radius = 50
-
-time = 10
-niter = 10

terrain_original.conf

@@ -1,16 +0,0 @@
-mapsize = 1000
-scale = 400
-vscale = 300
-offset = 0
-persistence = 0.6
-lacunarity = 2.0
-flow_method = steepest
-
-K = 1
-m = 0.35
-d = 0
-sea_level = 0
-flex_radius = 20
-
-time = 10
-niter = 10

terrainlib/__init__.py

@@ -1,7 +0,0 @@
-# Load packages and provide easy access to important functions
-
-from .settings import read_config_file
-from .erosion import EvolutionModel
-from .save import save
-from .bounds import make_bounds, twist, get_fixed
-from .view import stats, update, plot

terrainlib/bounds.py

@@ -1,74 +0,0 @@
-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=rivers.dtype)
-    bounds_v = np.zeros((Y-1, X), dtype=rivers.dtype)
-
-    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

terrainlib/erosion.py

@@ -1,95 +0,0 @@
-import numpy as np
-import scipy.ndimage as im
-from .rivermapper import flow
-
-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
-    """
-
-    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
-                    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]: # 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] # If between 2 pixels, perform linear interpolation.
-
-    return dem_new
-
-def diffusion(dem, time, d=1):
-    radius = d * time**.5
-    if radius == 0:
-        return dem
-    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, flow_method='semirandom'):
-        self.dem = dem
-        #self.bedrock = dem
-        self.K = K
-        self.m = m
-        self.d = d
-        self.sea_level = sea_level
-        self.flex_radius = flex_radius
-        self.define_isostasy()
-        self.flow_method = flow_method
-        #set_flow_method(flow_method)
-        if flow:
-            self.calculate_flow()
-        else:
-            self.lakes = dem
-            self.dirs = np.zeros(dem.shape, dtype=int)
-            self.rivers = np.zeros(dem.shape, dtype=int)
-            self.flow_uptodate = False
-
-    def calculate_flow(self):
-        self.dirs, self.lakes, self.rivers = flow(self.dem, method=self.flow_method)
-        self.flow_uptodate = True
-
-    def advection(self, time):
-        dem = advection(np.maximum(self.dem, 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
-
-    def define_isostasy(self):
-        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') # Calculate blurred DEM
-        correction = (self.ref_isostasy - isostasy) * rate # Compare it with the reference isostasy
-        self.dem = self.dem + correction # Adjust

terrainlib/rivermapper.py

@@ -1,278 +0,0 @@
-import numpy as np
-import numpy.random as npr
-from collections import defaultdict
-
-# This file provide functions to construct the river tree from an elevation model.
-# Based on a research paper:
-#   | Cordonnier, G., Bovy, B., and Braun, J.:
-#   | A versatile, linear complexity algorithm for flow routing in topographies with depressions,
-#   | Earth Surf. Dynam., 7, 549–562, https://doi.org/10.5194/esurf-7-549-2019, 2019.
-# Big thanks to them for releasing this paper under a free license ! :)
-
-# The algorithm here makes use of most of the paper's concepts, including the Planar Boruvka algorithm.
-# Only flow_local and accumulate_flow are custom algorithms.
-
-# Define two different method for local flow routing
-def flow_local_steepest(plist):
-    vmax = 0.0
-    imax = 0.0
-    for i, p in enumerate(plist):
-        if p > vmax:
-            vmax = p
-            imax = i
-    if vmax > 0.0:
-        return imax+1
-    return 0
-
-def flow_local_semirandom(plist):
-    """
-    Determines a flow direction based on denivellation for every neighbouring node.
-    Denivellation must be positive for downward and zero for flat or upward:
-    dz = max(zref-z, 0)
-    """
-    psum = sum(plist)
-    if psum == 0:
-        return 0
-    r = npr.random() * psum
-    for i, p in enumerate(plist):
-        if r < p:
-            return i+1
-        r -= p
-
-flow_local_methods = {
-    'steepest' : flow_local_steepest,
-    'semirandom' : flow_local_semirandom,
-}
-
-def flow(dem, method='semirandom'):
-    if method in flow_local_methods:
-        flow_local = flow_local_methods[method]
-    else:
-        raise KeyError('Flow method \'{}\' does not exist'.format(method))
-
-    # Flow locally
-    dirs1 = np.zeros(dem.shape, dtype=int)
-    dirs2 = np.zeros(dem.shape, dtype=int)
-    (X, Y) = dem.shape
-    Xmax, Ymax = X-1, Y-1
-    singular = []
-    for x in range(X):
-        z0 = z1 = z2 = dem[x,0]
-        for y in range(Y):
-            z0 = z1
-            z1 = z2
-            if y < Ymax:
-                z2 = dem[x, y+1]
-
-            plist = [
-                max(z1-dem[x+1,y],0) if x<Xmax else 0, # 1: x -> x+1
-                max(z1-z2,0),                          # 2: y -> y+1
-                max(z1-dem[x-1,y],0) if x>0    else 0, # 3: x -> x-1
-                max(z1-z0,0),                          # 4: y -> y-1
-            ]
-            
-            pdir = flow_local(plist)
-            dirs2[x,y] = pdir
-            if pdir == 0:
-                singular.append((x,y))
-            elif pdir == 1:
-                dirs1[x+1,y] += 1
-            elif pdir == 2:
-                dirs1[x,y+1] += 2
-            elif pdir == 3:
-                dirs1[x-1,y] += 4
-            elif pdir == 4:
-                dirs1[x,y-1] += 8
-
-    # Compute basins
-    basin_id = np.zeros(dem.shape, dtype=int)
-    stack = []
-
-    for i, s in enumerate(singular):
-        queue = [s]
-        while queue:
-            x, y = queue.pop()
-            basin_id[x,y] = i
-            d = int(dirs1[x,y])
-
-            if d & 1:
-                queue.append((x-1,y))
-            if d & 2:
-                queue.append((x,y-1))
-            if d & 4:
-                queue.append((x+1,y))
-            if d & 8:
-                queue.append((x,y+1))
-
-    del dirs1
-
-    # Link basins
-    nsing = len(singular)
-    links = {}
-    def add_link(b0, b1, elev, bound):
-        b = (min(b0,b1),max(b0,b1))
-        if b not in links or links[b][0] > elev:
-            links[b] = (elev, bound)
-
-    for x in range(X):
-        b0 = basin_id[x,0]
-        add_link(-1, b0, dem[x,0], (True, x, 0))
-        for y in range(1,Y):
-            b1 = basin_id[x,y]
-            if b0 != b1:
-                add_link(b0, b1, max(dem[x,y-1],dem[x,y]), (True, x, y))
-            b0 = b1
-        add_link(-1, b1, dem[x,Ymax], (True, x, Y))
-    for y in range(Y):
-        b0 = basin_id[0,y]
-        add_link(-1, b0, dem[0,y], (False, 0, y))
-        for x in range(1,X):
-            b1 = basin_id[x,y]
-            if b0 != b1:
-                add_link(b0, b1, max(dem[x-1,y],dem[x,y]), (False, x, y))
-            b0 = b1
-        add_link(-1, b1, dem[Xmax,y], (False, X, y))
-
-    # Computing basin tree
-    graph = planar_boruvka(links)
-
-    basin_links = defaultdict(dict)
-    for elev, b1, b2, bound in graph:
-        basin_links[b1][b2] = basin_links[b2][b1] = (elev, bound)
-    basins = np.zeros(nsing+1)
-    stack = [(-1, float('-inf'))]
-
-    # Applying basin flowing
-    dir_reverse = (0, 3, 4, 1, 2)
-    while stack:
-        b1, elev1 = stack.pop()
-        basins[b1] = elev1
-
-        for b2, (elev2, bound) in basin_links[b1].items():
-            stack.append((b2, max(elev1, elev2)))
-
-            # Reverse flow direction in b2 (TODO)
-            isY, x, y = bound
-            backward = True # Whether water will escape the basin in +X/+Y direction
-            if not (x < X and y < Y and basin_id[x,y] == b2):
-                if isY:
-                    y -= 1
-                else:
-                    x -= 1
-                backward = False
-            d = 2*backward + isY + 1
-            while d > 0:
-                d, dirs2[x,y] = dirs2[x,y], d
-                if d == 1:
-                    x += 1
-                elif d == 2:
-                    y += 1
-                elif d == 3:
-                    x -= 1
-                elif d == 4:
-                    y -= 1
-                d = dir_reverse[d]
-
-            del basin_links[b2][b1]
-        del basin_links[b1]
-
-    # Calculating water quantity
-    dirs2[-1,:][dirs2[-1,:]==1] = 0
-    dirs2[:,-1][dirs2[:,-1]==2] = 0
-    dirs2[0,:][dirs2[0,:]==3] = 0
-    dirs2[:,0][dirs2[:,0]==4] = 0
-
-    waterq = accumulate_flow(dirs2)
-
-    return dirs2, basins[basin_id], waterq
-
-def accumulate_flow(dirs):
-    ndonors = np.zeros(dirs.shape, dtype=int)
-    ndonors[1:,:] += dirs[:-1,:] == 1
-    ndonors[:,1:] += dirs[:,:-1] == 2
-    ndonors[:-1,:] += dirs[1:,:] == 3
-    ndonors[:,:-1] += dirs[:,1:] == 4
-    waterq = np.ones(dirs.shape, dtype=int)
-
-    (X, Y) = dirs.shape
-    rangeX = range(X)
-    rangeY = range(Y)
-    for x in rangeX:
-        for y in rangeY:
-            if ndonors[x,y] > 0:
-                continue
-            xw, yw = x, y
-            w = waterq[xw, yw]
-            while 1:
-                d = dirs[xw, yw]
-                if d <= 0:
-                    break
-                elif d == 1:
-                    xw += 1
-                elif d == 2:
-                    yw += 1
-                elif d == 3:
-                    xw -= 1
-                elif d == 4:
-                    yw -= 1
-
-                w += waterq[xw, yw]
-                waterq[xw, yw] = w
-
-                if ndonors[xw, yw] > 1:
-                    ndonors[xw, yw] -= 1
-                    break
-
-    return waterq
-
-def planar_boruvka(links):
-    # Compute basin tree
-
-    basin_list = defaultdict(dict)
-
-    for (b1, b2), (elev, bound) in links.items():
-        basin_list[b1][b2] = basin_list[b2][b1] = (elev, b1, b2, bound)
-
-    threshold = 8
-    lowlevel = {}
-    for k, v in basin_list.items():
-        if len(v) <= threshold:
-            lowlevel[k] = v
-
-    basin_graph = []
-    n = len(basin_list)
-    while n > 1:
-        (b1, lnk1) = lowlevel.popitem()
-        b2 = min(lnk1, key=lnk1.get)
-        lnk2 = basin_list[b2]
-
-        # Add link to the graph
-        basin_graph.append(lnk1[b2])
-
-        # Union : merge basin 1 into basin 2
-        # First, delete the direct link
-        del lnk1[b2]
-        del lnk2[b1]
-
-        # Look for basin 1's neighbours, and add them to basin 2 if they have a lower pass
-        for k, v in lnk1.items():
-            bk = basin_list[k]
-            if k in lnk2 and lnk2[k] < v:
-                del bk[b1]
-            else:
-                lnk2[k] = v
-                bk[b2] = bk.pop(b1)
-
-            if k not in lowlevel and len(bk) <= threshold:
-                lowlevel[k] = bk
-
-        if b2 in lowlevel:
-            if len(lnk2) > threshold:
-                del lowlevel[b2]
-        elif len(lnk2) <= threshold:
-            lowlevel[b2] = lnk2
-        del lnk1
-
-        n -= 1
-
-    return basin_graph

terrainlib/save.py

@@ -1,13 +0,0 @@
-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(bin_data)

terrainlib/settings.py

@@ -1,16 +0,0 @@
-import os.path
-
-def read_config_file(fname):
-    settings = {}
-
-    if not os.path.isfile(fname):
-        return settings
-
-    with open(fname, 'r') as f:
-        for line in f:
-            slist = line.split('=', 1)
-            if len(slist) >= 2:
-                prefix, suffix = slist
-                settings[prefix.strip()] = suffix.strip()
-
-    return settings

terrainlib_lua/erosion.lua

@@ -0,0 +1,221 @@
+-- erosion.lua
+
+local function erode(model, time)
+	--local tinsert = table.insert
+	local mmin, mmax = math.min, math.max
+	local dem = model.dem
+	local dirs = model.dirs
+	local lakes = model.lakes
+	local rivers = model.rivers
+	local sea_level = model.params.sea_level
+	local K = model.params.K
+	local m = model.params.m
+	local X, Y = dem.X, dem.Y
+	local scalars = type(K) == "number" and type(m) == "number"
+
+	local erosion_time
+	if model.params.variable_erosion then
+		erosion_time = {}
+	else
+		erosion_time = model.erosion_time or {}
+	end
+
+	if scalars then
+		for i=1, X*Y do
+			local etime = 1 / (K*rivers[i]^m)
+			erosion_time[i] = etime
+			lakes[i] = mmax(lakes[i], dem[i], sea_level)
+		end
+	else
+		for i=1, X*Y do
+			local etime = 1 / (K[i]*rivers[i]^m[i])
+			erosion_time[i] = etime
+			lakes[i] = mmax(lakes[i], dem[i], sea_level)
+		end
+	end
+
+	for i=1, X*Y do
+		local iw = i
+		local remaining = time
+		local new_elev
+		while true do
+			local inext = iw
+			local d = dirs[iw]
+
+			if d == 0 then
+				new_elev = lakes[iw]
+				break
+			elseif d == 1 then
+				inext = iw+X
+			elseif d == 2 then
+				inext = iw+1
+			elseif d == 3 then
+				inext = iw-X
+			elseif d == 4 then
+				inext = iw-1
+			end
+
+			local etime = erosion_time[iw]
+			if remaining <= etime then
+				local c = remaining / etime
+				new_elev = (1-c) * lakes[iw] + c * lakes[inext]
+				break
+			end
+
+			remaining = remaining - etime
+			iw = inext
+		end
+
+		dem[i] = mmin(dem[i], new_elev)
+	end
+end
+
+local function diffuse(model, time)
+    local mmax = math.max
+    local dem = model.dem
+    local X, Y = dem.X, dem.Y
+    local d = model.params.d
+    local dmax = d
+    if type(d) == "table" then
+        dmax = -math.huge
+        for i=1, X*Y do
+            dmax = mmax(dmax, d[i])
+        end
+    end
+
+    local diff = dmax * time
+    local niter = math.floor(diff) + 1
+    local ddiff = diff / niter
+
+    local temp = {}
+    for n=1, niter do
+        local i = 1
+        for y=1, Y do
+            local iN = (y==1) and 0 or -X
+            local iS = (y==Y) and 0 or X
+            for x=1, X do
+                local iW = (x==1) and 0 or -1
+                local iE = (x==X) and 0 or 1
+                temp[i] = (dem[i+iN]+dem[i+iE]+dem[i+iS]+dem[i+iW])*0.25 - dem[i]
+                i = i + 1
+            end
+        end
+
+        for i=1, X*Y do
+            dem[i] = dem[i] + temp[i]*ddiff
+        end
+    end
+    -- TODO Test this
+end
+
+local modpath = ""
+if minetest then
+    if minetest.global_exists('mapgen_rivers') then
+        modpath = mapgen_rivers.modpath .. "terrainlib_lua/"
+    else
+        modpath = minetest.get_modpath(minetest.get_current_modname()) .. "terrainlib_lua/"
+    end
+end
+
+local rivermapper = dofile(modpath .. "rivermapper.lua")
+local gaussian = dofile(modpath .. "gaussian.lua")
+
+local function flow(model)
+	model.dirs, model.lakes = rivermapper.flow_routing(model.dem, model.dirs, model.lakes, 'semirandom')
+	model.rivers = rivermapper.accumulate(model.dirs, model.rivers)
+end
+
+local function uplift(model, time)
+	local dem = model.dem
+	local X, Y = dem.X, dem.Y
+	local uplift_rate = model.params.uplift
+	if type(uplift_rate) == "number" then
+		local uplift_total = uplift_rate * time
+		for i=1, X*Y do
+			dem[i] = dem[i] + uplift_total
+		end
+	else
+		for i=1, X*Y do
+			dem[i] = dem[i] + uplift_rate[i]*time
+		end
+	end
+end
+
+local function noise(model, time)
+	local random = math.random
+	local dem = model.dem
+	local noise_depth = model.params.noise * 2 * time
+	local X, Y = dem.X, dem.Y
+	for i=1, X*Y do
+		dem[i] = dem[i] + (random()-0.5) * noise_depth
+	end
+end
+
+local function define_isostasy(model, ref, link)
+    ref = ref or model.dem
+    if link then
+        model.isostasy_ref = ref
+        return
+    end
+
+    local X, Y = ref.X, ref.Y
+    local ref2 = model.isostasy_ref or {X=X, Y=Y}
+    model.isostasy_ref = ref2
+    for i=1, X*Y do
+        ref2[i] = ref[i]
+    end
+
+    return ref2
+end
+
+local function isostasy(model)
+	local dem = model.dem
+	local X, Y = dem.X, dem.Y
+	local temp = {X=X, Y=Y}
+	local ref = model.isostasy_ref
+	for i=1, X*Y do
+		temp[i] = ref[i] - dem[i]
+	end
+
+	gaussian.gaussian_blur_approx(temp, model.params.compensation_radius, 4)
+	for i=1, X*Y do
+		dem[i] = dem[i] + temp[i]
+	end
+end
+
+local evol_model_mt = {
+	erode = erode,
+	diffuse = diffuse,
+	flow = flow,
+	uplift = uplift,
+	noise = noise,
+	isostasy = isostasy,
+	define_isostasy = define_isostasy,
+}
+
+evol_model_mt.__index = evol_model_mt
+
+local defaults = {
+	K = 1,
+	m = 0.5,
+	d = 1,
+	variable_erosion = false,
+	sea_level = 0,
+	uplift = 10,
+	noise = 0.001,
+	compensation_radius = 50,
+}
+
+local function EvolutionModel(params)
+	params = params or {}
+	local o = {params = params}
+	for k, v in pairs(defaults) do
+		if params[k] == nil then
+			params[k] = v
+		end
+	end
+	o.dem = params.dem
+	return setmetatable(o, evol_model_mt)
+end
+
+return EvolutionModel

terrainlib_lua/gaussian.lua

@@ -0,0 +1,102 @@
+-- gaussian.lua
+
+local function get_box_size(sigma, n)
+	local v = sigma^2 / n
+	local r_ideal = ((12*v + 1) ^ 0.5 - 1) / 2
+	local r_down = math.floor(r_ideal)
+	local r_up = math.ceil(r_ideal)
+	local v_down = ((2*r_down+1)^2 - 1) / 12
+	local v_up = ((2*r_up+1)^2 - 1) / 12
+	local m_ideal = (v - v_down) / (v_up - v_down) * n
+	local m = math.floor(m_ideal+0.5)
+
+	local sizes = {}
+	for i=1, n do
+		sizes[i] = i<=m and 2*r_up+1 or 2*r_down+1
+	end
+
+	return sizes
+end
+
+local function box_blur_1d(map, size, first, incr, len, map2)
+	local n = math.ceil(size/2)
+	first = first or 1
+	incr = incr or 1
+	len = len or math.floor((#map-first)/incr)+1
+	local last = first + (len-1)*incr
+
+	local nth = first+(n-1)*incr
+	local sum = 0
+	for i=first, nth, incr do
+		if i == first then
+			sum = sum + map[i]
+		else
+			sum = sum + 2*map[i]
+		end
+	end
+
+	local i_left = nth
+	local incr_left = -incr
+	local i_right = nth
+	local incr_right = incr
+
+	map2 = map2 or {}
+	for i=first, last, incr do
+		map2[i] = sum / size
+		i_right = i_right + incr_right
+		sum = sum - map[i_left] + map[i_right]
+		i_left = i_left + incr_left
+
+		if i_left == first then
+			incr_left = incr
+		end
+		if i_right == last then
+			incr_right = -incr
+		end
+	end
+
+	return map2
+end
+
+local function box_blur_2d(map1, size, map2)
+	local X, Y = map1.X, map1.Y
+	map2 = map2 or {}
+	for y=1, Y do
+		box_blur_1d(map1, size, (y-1)*X+1, 1, X, map2)
+	end
+	for x=1, X do
+		box_blur_1d(map2, size, x, X, Y, map1)
+	end
+
+	return map1
+end
+
+--[[local function gaussian_blur(map, std, tail)
+	local exp = math.exp
+
+	local kernel_mid = math.ceil(std*tail) + 1
+	local kernel_size = kernel_mid * 2 - 1
+	local kernel = {}
+	local cst1 = 1/(std*(2*math.pi)^0.5)
+	local cst2 = -1/(2*std^2)
+	for i=1, kernel_size do
+		kernel[i] = cst1 * exp((i-kernel_mid)^2 * cst2)
+	end
+
+	]]
+
+local function gaussian_blur_approx(map, sigma, n, map2)
+	map2 = map2 or {}
+	local sizes = get_box_size(sigma, n)
+	for i=1, n do
+		box_blur_2d(map, sizes[i], map2)
+	end
+	return map
+end
+
+return {
+	get_box_size = get_box_size,
+	box_blur_1d = box_blur_1d,
+	box_blur_2d = box_blur_2d,
+	gaussian_blur_approx = gaussian_blur_approx,
+}

terrainlib_lua/rivermapper.lua

@@ -0,0 +1,385 @@
+-- rivermapper.lua
+
+-- This file provide functions to construct the river tree from an elevation model.
+-- Based on a research paper:
+--
+--     Cordonnier, G., Bovy, B., and Braun, J.:
+--     A versatile, linear complexity algorithm for flow routing in topographies with depressions,
+--     Earth Surf. Dynam., 7, 549–562, https://doi.org/10.5194/esurf-7-549-2019, 2019.
+--
+-- Big thanks to them for releasing this paper under a free license ! :)
+
+-- The algorithm here makes use of most of the paper's concepts, including the Planar Boruvka algorithm.
+-- Only flow_local and accumulate_flow are custom algorithms.
+
+local function flow_local_semirandom(plist)
+	local sum = 0
+	for i=1, #plist do
+		sum = sum + plist[i]
+	end
+	--for _, p in ipairs(plist) do
+		--sum = sum + p
+	--end
+	if sum == 0 then
+		return 0
+	end
+	local r = math.random() * sum
+	for i=1, #plist do
+		local p = plist[i]
+	--for i, p in ipairs(plist) do
+		if r < p then
+			return i
+		end
+		r = r - p
+	end
+	return 0
+end
+
+local flow_methods = {
+	semirandom = flow_local_semirandom,
+}
+
+local function flow_routing(dem, dirs, lakes, method)
+	method = method or 'semirandom'
+	local flow_local = flow_methods[method] or flow_local_semirandom
+
+	dirs = dirs or {}
+	lakes = lakes or {}
+
+	-- Localize for performance
+	--local tinsert = table.insert
+	local tremove = table.remove
+	local mmax = math.max
+
+	local X, Y = dem.X, dem.Y
+	dirs.X = X
+	dirs.Y = Y
+	lakes.X = X
+	lakes.Y = Y
+	local i = 1
+	local dirs2 = {}
+	for i=1, X*Y do
+		dirs2[i] = 0
+	end
+
+	local singular = {}
+	for y=1, Y do
+		for x=1, X do
+			local zi = dem[i]
+			local plist = {
+				y<Y and mmax(zi-dem[i+X], 0) or 0, -- Southward
+				x<X and mmax(zi-dem[i+1], 0) or 0, -- Eastward
+				y>1 and mmax(zi-dem[i-X], 0) or 0, -- Northward
+				x>1 and mmax(zi-dem[i-1], 0) or 0, -- Westward
+			}
+
+			local d = flow_local(plist)
+			dirs[i] = d
+			if d == 0 then
+				singular[#singular+1] = i
+			elseif d == 1 then
+				dirs2[i+X] = dirs2[i+X] + 1
+			elseif d == 2 then
+				dirs2[i+1] = dirs2[i+1] + 2
+			elseif d == 3 then
+				dirs2[i-X] = dirs2[i-X] + 4
+			elseif d == 4 then
+				dirs2[i-1] = dirs2[i-1] + 8
+			end
+			i = i + 1
+		end
+	end
+
+	-- Compute basins and links
+	local nbasins = #singular
+	local basin_id = {}
+	local links = {}
+	local basin_links
+
+	local function add_link(i1, i2, b1, isY)
+		local b2
+		if i2 == 0 then
+			b2 = 0
+		else
+			b2 = basin_id[i2]
+			if b2 == 0 then
+				return
+			end
+		end
+		if b2 ~= b1 then
+			local elev = i2 == 0 and dem[i1] or mmax(dem[i1], dem[i2])
+			local l2 = basin_links[b2]
+			if not l2 then
+				l2 = {}
+				basin_links[b2] = l2
+			end
+			if not l2.elev or l2.elev > elev then
+				l2.elev = elev
+				l2.i = mmax(i1,i2)
+				l2.is_y = isY
+				l2[1] = b2
+				l2[2] = b1
+			end
+		end
+	end
+
+	for i=1, X*Y do
+		basin_id[i] = 0
+	end
+	--for ib, s in ipairs(singular) do
+	for ib=1, nbasins do
+		--local s = singular[ib]
+		local queue = {singular[ib]}
+		basin_links = {}
+		links[#links+1] = basin_links
+		--tinsert(links, basin_links)
+		while #queue > 0 do
+			local i = tremove(queue)
+			basin_id[i] = ib
+			local d = dirs2[i]
+
+			if d >= 8 then -- River coming from East
+				d = d - 8
+				queue[#queue+1] = i+1
+				--tinsert(queue, i+X)
+			elseif i%X > 0 then
+				add_link(i, i+1, ib, false)
+			else
+				add_link(i, 0, ib, false)
+			end
+
+			if d >= 4 then -- River coming from South
+				d = d - 4
+				queue[#queue+1] = i+X
+				--tinsert(queue, i+1)
+			elseif i <= X*(Y-1) then
+				add_link(i, i+X, ib, true)
+			else
+				add_link(i, 0, ib, true)
+			end
+
+			if d >= 2 then -- River coming from West
+				d = d - 2
+				queue[#queue+1] = i-1
+				--tinsert(queue, i-X)
+			elseif i%X ~= 1 then
+				add_link(i, i-1, ib, false)
+			else
+				add_link(i, 0, ib, false)
+			end
+
+			if d >= 1 then -- River coming from North
+				queue[#queue+1] = i-X
+				--tinsert(queue, i-1)
+			elseif i > X then
+				add_link(i, i-X, ib, true)
+			else
+				add_link(i, 0, ib, true)
+			end
+		end
+	end
+	dirs2 = nil
+
+	links[0] = {}
+	local nlinks = {}
+	for i=0, nbasins do
+		nlinks[i] = 0
+	end
+
+	--for ib1, blinks in ipairs(links) do
+	for ib1=1, #links do
+		for ib2, link in pairs(links[ib1]) do
+			if ib2 < ib1 then
+				links[ib2][ib1] = link
+				nlinks[ib1] = nlinks[ib1] + 1
+				nlinks[ib2] = nlinks[ib2] + 1
+			end
+		end
+	end
+
+	local lowlevel = {}
+	for i, n in pairs(nlinks) do
+		if n <= 8 then
+			lowlevel[i] = links[i]
+		end
+	end
+
+	local basin_graph = {}
+	for n=1, nbasins do
+		local b1, lnk1 = next(lowlevel)
+		lowlevel[b1] = nil
+
+		local b2
+		local lowest = math.huge
+		local lnk1 = links[b1]
+		local i = 0
+		for bn, bdata in pairs(lnk1) do
+			i = i + 1
+			if bdata.elev < lowest then
+				lowest = bdata.elev
+				b2 = bn
+			end
+		end
+
+		-- Add link to the graph
+		local bound = lnk1[b2]
+		local bb1, bb2 = bound[1], bound[2]
+		if not basin_graph[bb1] then
+			basin_graph[bb1] = {}
+		end
+		if not basin_graph[bb2] then
+			basin_graph[bb2] = {}
+		end
+		basin_graph[bb1][bb2] = bound
+		basin_graph[bb2][bb1] = bound
+
+		-- Merge basin b1 into b2
+		local lnk2 = links[b2]
+		-- First, remove the link between b1 and b2
+		lnk1[b2] = nil
+		lnk2[b1] = nil
+		nlinks[b2] = nlinks[b2] - 1
+		if nlinks[b2] == 8 then
+			lowlevel[b2] = lnk2
+		end
+		-- Look for basin 1's neighbours, and add them to basin 2 if they have a lower pass
+		for bn, bdata in pairs(lnk1) do
+			local lnkn = links[bn]
+			lnkn[b1] = nil
+
+			if lnkn[b2] then
+				nlinks[bn] = nlinks[bn] - 1
+				if nlinks[bn] == 8 then
+					lowlevel[bn] = lnkn
+				end
+			else
+				nlinks[b2] = nlinks[b2] + 1
+				if nlinks[b2] == 9 then
+					lowlevel[b2] = nil
+				end
+			end
+
+			if not lnkn[b2] or lnkn[b2].elev > bdata.elev then
+				lnkn[b2] = bdata
+				lnk2[bn] = bdata
+			end
+		end
+	end
+
+	local queue = {[0] = -math.huge}
+	local basin_lake = {}
+	for n=1, nbasins do
+		basin_lake[n] = 0
+	end
+	local reverse = {3, 4, 1, 2, [0]=0}
+	for n=1, nbasins do
+		local b1, elev1 = next(queue)
+		queue[b1] = nil
+		basin_lake[b1] = elev1
+		for b2, bound in pairs(basin_graph[b1]) do
+			-- Make b2 flow into b1
+			local i = bound.i
+			local dir = bound.is_y and 3 or 4
+			if basin_id[i] ~= b2 then
+				dir = dir - 2
+				if bound.is_y then
+					i = i - X
+				else
+					i = i - 1
+				end
+			elseif b1 == 0 then
+				dir = 0
+			end
+
+			repeat
+				dir, dirs[i] = dirs[i], dir
+				if dir == 1 then
+					i = i + X
+				elseif dir == 2 then
+					i = i + 1
+				elseif dir == 3 then
+					i = i - X
+				elseif dir == 4 then
+					i = i - 1
+				end
+				dir = reverse[dir]
+			until dir == 0
+			-- Add b2 into the queue
+			queue[b2] = mmax(elev1, bound.elev)
+			basin_graph[b2][b1] = nil
+		end
+		basin_graph[b1] = nil
+	end
+
+	for i=1, X*Y do
+		lakes[i] = basin_lake[basin_id[i]]
+	end
+
+	return dirs, lakes
+end
+
+local function accumulate(dirs, waterq)
+	waterq = waterq or {}
+	local X, Y = dirs.X, dirs.Y
+	--local tinsert = table.insert
+
+	local ndonors = {}
+	local waterq = {X=X, Y=Y}
+	for i=1, X*Y do
+		ndonors[i] = 0
+		waterq[i] = 1
+	end
+
+	--for i1, dir in ipairs(dirs) do
+	for i1=1, X*Y do
+		local i2
+		local dir = dirs[i1]
+		if dir == 1 then
+			i2 = i1+X
+		elseif dir == 2 then
+			i2 = i1+1
+		elseif dir == 3 then
+			i2 = i1-X
+		elseif dir == 4 then
+			i2 = i1-1
+		end
+		if i2 then
+			ndonors[i2] = ndonors[i2] + 1
+		end
+	end
+
+	for i1=1, X*Y do
+		if ndonors[i1] == 0 then
+			local i2 = i1
+			local dir = dirs[i2]
+			local w = waterq[i2]
+			while dir > 0 do
+				if dir == 1 then
+					i2 = i2 + X
+				elseif dir == 2 then
+					i2 = i2 + 1
+				elseif dir == 3 then
+					i2 = i2 - X
+				elseif dir == 4 then
+					i2 = i2 - 1
+				end
+				w = w + waterq[i2]
+				waterq[i2] = w
+
+				if ndonors[i2] > 1 then
+					ndonors[i2] = ndonors[i2] - 1
+					break
+				end
+				dir = dirs[i2]
+			end
+		end
+	end
+
+	return waterq
+end
+
+return {
+	flow_routing = flow_routing,
+	accumulate = accumulate,
+	flow_methods = flow_methods,
+}

terrainlib_lua/twist.lua

@@ -0,0 +1,102 @@
+-- bounds.lua
+
+local function get_bounds(dirs, rivers)
+	local X, Y = dirs.X, dirs.Y
+	local bounds_x = {X=X, Y=Y}
+	local bounds_y = {X=X, Y=Y}
+	for i=1, X*Y do
+		bounds_x[i] = 0
+		bounds_y[i] = 0
+	end
+
+	for i=1, X*Y do
+		local dir = dirs[i]
+		local river = rivers[i]
+		if dir == 1 then -- South (+Y)
+			bounds_y[i] = river
+		elseif dir == 2 then -- East (+X)
+			bounds_x[i] = river
+		elseif dir == 3 then -- North (-Y)
+			bounds_y[i-X] = river
+		elseif dir == 4 then -- West (-X)
+			bounds_x[i-1] = river
+		end
+	end
+
+	return bounds_x, bounds_y
+end
+
+local function twist(dirs, rivers, n)
+	n = n or 5
+	local X, Y = dirs.X, dirs.Y
+	local bounds_x, bounds_y = get_bounds(dirs, rivers)
+	local dn = 0.5 / n
+
+	local offset_x = {X=X, Y=Y}
+	local offset_y = {X=X, Y=Y}
+	local offset_x_alt = {X=X, Y=Y}
+	local offset_y_alt = {X=X, Y=Y}
+	for i=1, X*Y do
+		offset_x[i] = 0
+		offset_y[i] = 0
+	end
+
+	for nn=1, n do
+		local i = 1
+		for y=1, Y do
+			for x=1, X do
+				local ox, oy = offset_x[i], offset_y[i]
+				if dirs[i] ~= 0 and rivers[i] > 1 then
+					local sum_fx = 0
+					local sum_fy = 0
+					local sum_w = 0
+					local b
+					if x < X then
+						b = bounds_x[i]
+						sum_fx = sum_fx + b*(offset_x[i+1]+1)
+						sum_fy = sum_fy + b*offset_y[i+1]
+						sum_w = sum_w + b
+					end
+					if y < Y then
+						b = bounds_y[i]
+						sum_fx = sum_fx + b*offset_x[i+X]
+						sum_fy = sum_fy + b*(offset_y[i+X]+1)
+						sum_w = sum_w + b
+					end
+					if x > 1 then
+						b = bounds_x[i-1]
+						sum_fx = sum_fx + b*(offset_x[i-1]-1)
+						sum_fy = sum_fy + b*offset_y[i-1]
+						sum_w = sum_w + b
+					end
+					if y > 1 then
+						b = bounds_y[i-X]
+						sum_fx = sum_fx + b*offset_x[i-X]
+						sum_fy = sum_fy + b*(offset_y[i-X]-1)
+						sum_w = sum_w + b
+					end
+					local fx, fy = sum_fx/sum_w - ox, sum_fy/sum_w - oy
+					local fd = (fx*fx+fy*fy) ^ 0.5
+					if fd > dn then
+						local c = dn/fd
+						fx, fy = fx*c, fy*c
+					end
+
+					offset_x_alt[i] = ox+fx
+					offset_y_alt[i] = oy+fy
+				else
+					offset_x_alt[i] = ox
+					offset_y_alt[i] = oy
+				end
+
+				i = i + 1
+			end
+		end
+		offset_x, offset_x_alt = offset_x_alt, offset_x
+		offset_y, offset_y_alt = offset_y_alt, offset_y
+	end
+
+	return offset_x, offset_y
+end
+
+return twist

view.py

@@ -19,7 +19,7 @@ except ImportError: # No module matplotlib
     has_matplotlib = False
 
 if has_matplotlib:
-    def view_map(dem, lakes, scale=1, sea_level=0.0, title=None):
+    def view_map(dem, lakes, scale=1, center=False, sea_level=0.0, title=None):
         lakes_sea = np.maximum(lakes, sea_level)
         water = np.maximum(lakes_sea - dem, 0)
         max_elev = dem.max()
@@ -31,7 +31,10 @@ if has_matplotlib:
         rgb = ls.shade(dem, cmap=cmap1, vert_exag=1/scale, blend_mode='soft', norm=norm_ground)
 
         (X, Y) = dem.shape
-        extent = (0, Y*scale, 0, X*scale)
+        if center:
+            extent = (-(Y+1)*scale/2, (Y-1)*scale/2, -(X+1)*scale/2, (X-1)*scale/2)
+        else:
+            extent = (-0.5*scale, (Y-0.5)*scale, -0.5*scale, (X-0.5)*scale)
         plt.imshow(np.flipud(rgb), extent=extent, interpolation='antialiased')
         alpha = (water > 0).astype('u1')
         plt.imshow(np.flipud(water), alpha=np.flipud(alpha), cmap=cmap2, extent=extent, vmin=0, vmax=max_depth, interpolation='antialiased')

view_map.py

@@ -5,13 +5,20 @@ import zlib
 import sys
 import os
 
-from terrainlib import stats, plot
+from view import stats, plot
+from readconfig import read_conf_file
 
-scale = 1
-if len(sys.argv) > 1:
-    os.chdir(sys.argv[1])
-if len(sys.argv) > 2:
-    scale = int(sys.argv[2])
+os.chdir(sys.argv[1])
+conf = read_conf_file('mapgen_rivers.conf')
+if 'center' in conf:
+    center = conf['center'] == 'true'
+else:
+    center = True
+
+if 'blocksize' in conf:
+    blocksize = float(conf['blocksize'])
+else:
+    blocksize = 15.0
 
 def load_map(name, dtype, shape):
     dtype = np.dtype(dtype)
@@ -21,9 +28,10 @@ def load_map(name, dtype, shape):
             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)
+shape = np.loadtxt('river_data/size', dtype='u4')
+shape = (shape[1], shape[0])
+dem = load_map('river_data/dem', '>i2', shape)
+lakes = load_map('river_data/lakes', '>i2', shape)
 
-stats(dem, lakes, scale=scale)
-plot(dem, lakes, scale=scale)
+stats(dem, lakes, scale=blocksize)
+plot(dem, lakes, scale=blocksize, center=center)