Use iterative finite differences for diffusion, instead of gaussian blur

Allows variable diffusion coefficients

.gitignore

@@ -6,5 +6,7 @@ offset_x
 offset_y
 bounds_x
 bounds_y
+dirs
 rivers
 unused/
+river_data/

LICENSE

@@ -0,0 +1,165 @@
+GNU LESSER GENERAL PUBLIC LICENSE
+                       Version 3, 29 June 2007
+
+ Copyright (C) 2007 Free Software Foundation, Inc. <http://fsf.org/>
+ Everyone is permitted to copy and distribute verbatim copies
+ of this license document, but changing it is not allowed.
+
+
+  This version of the GNU Lesser General Public License incorporates
+the terms and conditions of version 3 of the GNU General Public
+License, supplemented by the additional permissions listed below.
+
+  0. Additional Definitions.
+
+  As used herein, "this License" refers to version 3 of the GNU Lesser
+General Public License, and the "GNU GPL" refers to version 3 of the GNU
+General Public License.
+
+  "The Library" refers to a covered work governed by this License,
+other than an Application or a Combined Work as defined below.
+
+  An "Application" is any work that makes use of an interface provided
+by the Library, but which is not otherwise based on the Library.
+Defining a subclass of a class defined by the Library is deemed a mode
+of using an interface provided by the Library.
+
+  A "Combined Work" is a work produced by combining or linking an
+Application with the Library.  The particular version of the Library
+with which the Combined Work was made is also called the "Linked
+Version".
+
+  The "Minimal Corresponding Source" for a Combined Work means the
+Corresponding Source for the Combined Work, excluding any source code
+for portions of the Combined Work that, considered in isolation, are
+based on the Application, and not on the Linked Version.
+
+  The "Corresponding Application Code" for a Combined Work means the
+object code and/or source code for the Application, including any data
+and utility programs needed for reproducing the Combined Work from the
+Application, but excluding the System Libraries of the Combined Work.
+
+  1. Exception to Section 3 of the GNU GPL.
+
+  You may convey a covered work under sections 3 and 4 of this License
+without being bound by section 3 of the GNU GPL.
+
+  2. Conveying Modified Versions.
+
+  If you modify a copy of the Library, and, in your modifications, a
+facility refers to a function or data to be supplied by an Application
+that uses the facility (other than as an argument passed when the
+facility is invoked), then you may convey a copy of the modified
+version:
+
+   a) under this License, provided that you make a good faith effort to
+   ensure that, in the event an Application does not supply the
+   function or data, the facility still operates, and performs
+   whatever part of its purpose remains meaningful, or
+
+   b) under the GNU GPL, with none of the additional permissions of
+   this License applicable to that copy.
+
+  3. Object Code Incorporating Material from Library Header Files.
+
+  The object code form of an Application may incorporate material from
+a header file that is part of the Library.  You may convey such object
+code under terms of your choice, provided that, if the incorporated
+material is not limited to numerical parameters, data structure
+layouts and accessors, or small macros, inline functions and templates
+(ten or fewer lines in length), you do both of the following:
+
+   a) Give prominent notice with each copy of the object code that the
+   Library is used in it and that the Library and its use are
+   covered by this License.
+
+   b) Accompany the object code with a copy of the GNU GPL and this license
+   document.
+
+  4. Combined Works.
+
+  You may convey a Combined Work under terms of your choice that,
+taken together, effectively do not restrict modification of the
+portions of the Library contained in the Combined Work and reverse
+engineering for debugging such modifications, if you also do each of
+the following:
+
+   a) Give prominent notice with each copy of the Combined Work that
+   the Library is used in it and that the Library and its use are
+   covered by this License.
+
+   b) Accompany the Combined Work with a copy of the GNU GPL and this license
+   document.
+
+   c) For a Combined Work that displays copyright notices during
+   execution, include the copyright notice for the Library among
+   these notices, as well as a reference directing the user to the
+   copies of the GNU GPL and this license document.
+
+   d) Do one of the following:
+
+       0) Convey the Minimal Corresponding Source under the terms of this
+       License, and the Corresponding Application Code in a form
+       suitable for, and under terms that permit, the user to
+       recombine or relink the Application with a modified version of
+       the Linked Version to produce a modified Combined Work, in the
+       manner specified by section 6 of the GNU GPL for conveying
+       Corresponding Source.
+
+       1) Use a suitable shared library mechanism for linking with the
+       Library.  A suitable mechanism is one that (a) uses at run time
+       a copy of the Library already present on the user's computer
+       system, and (b) will operate properly with a modified version
+       of the Library that is interface-compatible with the Linked
+       Version.
+
+   e) Provide Installation Information, but only if you would otherwise
+   be required to provide such information under section 6 of the
+   GNU GPL, and only to the extent that such information is
+   necessary to install and execute a modified version of the
+   Combined Work produced by recombining or relinking the
+   Application with a modified version of the Linked Version. (If
+   you use option 4d0, the Installation Information must accompany
+   the Minimal Corresponding Source and Corresponding Application
+   Code. If you use option 4d1, you must provide the Installation
+   Information in the manner specified by section 6 of the GNU GPL
+   for conveying Corresponding Source.)
+
+  5. Combined Libraries.
+
+  You may place library facilities that are a work based on the
+Library side by side in a single library together with other library
+facilities that are not Applications and are not covered by this
+License, and convey such a combined library under terms of your
+choice, if you do both of the following:
+
+   a) Accompany the combined library with a copy of the same work based
+   on the Library, uncombined with any other library facilities,
+   conveyed under the terms of this License.
+
+   b) Give prominent notice with the combined library that part of it
+   is a work based on the Library, and explaining where to find the
+   accompanying uncombined form of the same work.
+
+  6. Revised Versions of the GNU Lesser General Public License.
+
+  The Free Software Foundation may publish revised and/or new versions
+of the GNU Lesser General Public License from time to time. Such new
+versions will be similar in spirit to the present version, but may
+differ in detail to address new problems or concerns.
+
+  Each version is given a distinguishing version number. If the
+Library as you received it specifies that a certain numbered version
+of the GNU Lesser General Public License "or any later version"
+applies to it, you have the option of following the terms and
+conditions either of that published version or of any later version
+published by the Free Software Foundation. If the Library as you
+received it does not specify a version number of the GNU Lesser
+General Public License, you may choose any version of the GNU Lesser
+General Public License ever published by the Free Software Foundation.
+
+  If the Library as you received it specifies that a proxy can decide
+whether future versions of the GNU Lesser General Public License shall
+apply, that proxy's public statement of acceptance of any version is
+permanent authorization for you to choose that version for the
+Library.

README.md

@@ -1,29 +1,114 @@
-mapgen_rivers
-=============
+# Map Generator with Rivers
+`mapgen_rivers v0.0` by Gaël de Sailly.
 
-Procedural map generator for Minetest 5.x. Still experimental and basic.
+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).
 
-Contains two distinct programs: Python scripts for pre-processing, and Lua scripts to generate the map on Minetest.
+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/79073532-7a567f00-7ce7-11ea-9791-8fb453f5175d.png)
+![Screenshot](https://user-images.githubusercontent.com/6905002/98825953-6289d980-2435-11eb-9e0b-704a95663ce0.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.
+
+# 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.
 
 # 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
-- `noise`, doing Perlin/Simplex noises
-- optionally, `matplotlib` (for map preview)
-
-They are commonly found on `pip` or `conda` Python distributions.
+This mod should be placed in the `mods/` directory of Minetest like any other mod.
 
 # Usage
-## 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).
+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.
 
-## 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.
+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.
+
+## 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.
+
+### 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`.
+
+## Map preview
+If you have `matplotlib` installed, `generate.py` will automatically show the grid aspect in real time during the erosion simulation.
+
+There is also a script to view a generated map afterwards: `view_map.py`. Its syntax is the following:
+```
+./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
+```

demo_data/size

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

erosion.py

@@ -1,86 +0,0 @@
-import numpy as np
-import scipy.ndimage as im
-import rivermapper as rm
-
-def advection(dem, dirs, rivers, time, K=1, m=0.5, sea_level=0):
-    dirs = dirs.copy()
-    dirs[0,:] = 0
-    dirs[-1,:] = 0
-    dirs[:,0] = 0
-    dirs[:,-1] = 0
-
-    adv_time = 1 / (K*rivers**m)
-    dem = np.maximum(dem, sea_level)
-    dem_new = np.zeros(dem.shape)
-
-    for y in range(dirs.shape[0]):
-        for x in range(dirs.shape[1]):
-            x0, y0 = x, y
-            x1, y1 = x, y
-            remaining = time
-            while True:
-                flow_dir = dirs[y0,x0]
-                if flow_dir == 0:
-                    remaining = 0
-                    break
-                elif flow_dir == 1:
-                    y1 += 1
-                elif flow_dir == 2:
-                    x1 += 1
-                elif flow_dir == 3:
-                    y1 -= 1
-                elif flow_dir == 4:
-                    x1 -= 1
-
-                if remaining <= adv_time[y0,x0]:
-                    break
-                remaining -= adv_time[y0,x0]
-                x0, y0 = x1, y1
-
-            c = remaining / adv_time[y0,x0]
-            dem_new[y,x] = c*dem[y1,x1] + (1-c)*dem[y0,x0]
-
-    return np.minimum(dem, dem_new)
-
-def diffusion(dem, time, d=1):
-    radius = d * time**.5
-    return im.gaussian_filter(dem, radius, mode='reflect')
-
-class EvolutionModel:
-    def __init__(self, dem, K=1, m=0.5, d=1, sea_level=0, flow=False, flex_radius=100):
-        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()
-        if flow:
-            self.calculate_flow()
-        else:
-            self.lakes = dem
-            self.dirs = np.zeros(dem.shape, dtype='u1')
-            self.rivers = np.zeros(dem.shape, dtype='u4')
-            self.flow_uptodate = False
-
-    def calculate_flow(self):
-        self.dirs, self.lakes, self.rivers = rm.flow(self.dem)
-        self.flow_uptodate = True
-
-    def advection(self, time):
-        dem = advection(self.lakes, self.dirs, self.rivers, time, K=self.K, m=self.m, sea_level=self.sea_level)
-        self.dem = np.minimum(dem, self.dem)
-        self.flow_uptodate = False
-
-    def diffusion(self, time):
-        self.dem = diffusion(self.dem, time, d=self.d)
-        self.flow_uptodate = False
-
-    def define_isostasy(self):
-        self.ref_isostasy = im.gaussian_filter(self.dem, self.flex_radius, mode='reflect')
-
-    def adjust_isostasy(self, rate=1):
-        isostasy = im.gaussian_filter(self.dem, self.flex_radius, mode='reflect')
-        correction = (self.ref_isostasy - isostasy) * rate
-        self.dem = self.dem + correction

generate.py

@@ -0,0 +1,233 @@
+#!/usr/bin/env python3
+
+import numpy as np
+from noise import snoise2, snoise3
+import os
+import sys
+
+import terrainlib
+
+class noisemap:
+    def __init__(self, X, Y, scale=0.01, vscale=1.0, tscale=1.0, offset=0.0, log=None, xbase=None, ybase=None, **params):
+        # Determine noise offset randomly
+        if xbase is None:
+            xbase = np.random.randint(8192)-4096
+        if ybase is None:
+            ybase = np.random.randint(8192)-4096
+        self.xbase = xbase
+        self.ybase = ybase
+        self.X = X
+        self.Y = Y
+        self.scale = scale
+        if log:
+            vscale /= offset
+        self.vscale = vscale
+        self.tscale = tscale
+        self.offset = offset
+        self.log = log
+        self.params = params
+
+    def get2d(self):
+        n = np.zeros((self.X, self.Y))
+        for x in range(self.X):
+            for y in range(self.Y):
+                n[x,y] = snoise2(x/self.scale + self.xbase, y/self.scale + self.ybase, **self.params)
+
+        if self.log:
+            return np.exp(n*self.vscale) * self.offset
+        else:
+            return n*self.vscale + self.offset
+
+    def get3d(self, t=0):
+        t /= self.tscale
+        n = np.zeros((self.X, self.Y))
+        for x in range(self.X):
+            for y in range(self.Y):
+                n[x,y] = snoise3(x/self.scale + self.xbase, y/self.scale + self.ybase, t, **self.params)
+
+        if self.log:
+            return np.exp(n*self.vscale) * self.offset
+        else:
+            return n*self.vscale + self.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))
+tectonics_time = float(get_setting('tectonics_time', 0.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,
+}
+
+params_K = {
+    "offset" : K,
+    "vscale" : K,
+    "scale" : 400,
+    "octaves" : 1,
+    "persistence" : 0.5,
+    "lacunarity" : 2,
+    "log" : True,
+}
+
+params_m = {
+    "offset" : m,
+    "vscale" : m*0.2,
+    "scale" : 400,
+    "octaves" : 1,
+    "persistence" : 0.5,
+    "lacunarity" : 2,
+    "log" : False,
+}
+
+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)
+
+if tectonics_time == 0.0:
+    n = noisemap(mapsize+1, mapsize+1, **params).get2d()
+else:
+    terrain_noisemap = noisemap(mapsize+1, mapsize+1, tscale=tectonics_time, **params)
+    n = terrain_noisemap.get3d()
+m_map = noisemap(mapsize+1, mapsize+1, **params_m).get2d()
+K_map = noisemap(mapsize+1, mapsize+1, **params_K).get2d()
+
+### COMPUTE LANDSCAPE EVOLUTION
+# Initialize landscape evolution model
+print('Initializing model')
+model = terrainlib.EvolutionModel(n, K=K_map, m=m_map, d=K_map, 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)
+    if tectonics_time != 0.0:
+        print('Isostasy reference redefinition')
+        model.define_isostasy(terrain_noisemap.get3d((i+1)*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!')

geometry.lua

@@ -1,45 +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
+	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)
-	-- X, Y 4-vectors giving the coordinates of the 4 nodes
+	-- 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
 
-local function area(X, Y) -- Signed area of polygon, in function of direction of rotation. Clockwise = positive.
-	local n = #X
-	local sum = X[1]*Y[n] - X[n]*Y[1]
-	for i=2, n do
-		sum = sum + X[i]*Y[i-1] - X[i-1]*Y[i]
-	end
-
-	return sum/2
-end
-
-return {
-	distance_to_segment = distance_to_segment,
-	transform_quadri = transform_quadri,
-	area = area,
-}
+return transform_quadri

heightmap.lua

@@ -0,0 +1,138 @@
+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
+
+local MAP_BOTTOM = -31000
+
+-- 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.z-minp.z+1
+
+	local terrain_height_map = {}
+	local lake_height_map = {}
+
+	local i = 1
+	for z=minp.z, maxp.z do
+		for x=minp.x, maxp.x do
+			local poly = polygons[i]
+
+			if poly then
+				local xf, zf = transform_quadri(poly.x, poly.z, x, z)
+				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
+				))
+
+				-- Spatial gradient of the interpolation
+				local slope_x = zf*(vdem[3]-vdem[4]) + (1-zf)*(vdem[2]-vdem[1]) < 0
+				local slope_z = xf*(vdem[3]-vdem[2]) + (1-xf)*(vdem[4]-vdem[1]) < 0
+				local lake_id = 0
+				if slope_x then
+					if slope_z then
+						lake_id = 3
+					else
+						lake_id = 2
+					end
+				else
+					if slope_z then
+						lake_id = 4
+					else
+						lake_id = 1
+					end
+				end
+				local lake_height = math.max(math.floor(poly.lake[lake_id]), 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] = MAP_BOTTOM
+				lake_height_map[i] = MAP_BOTTOM
+			end
+			i = i + 1
+		end
+	end
+
+	return terrain_height_map, lake_height_map
+end
+
+return heightmaps

init.lua

@@ -1,59 +1,26 @@
+mapgen_rivers = {}
+
 local modpath = minetest.get_modpath(minetest.get_current_modname()) .. '/'
-local worldpath = minetest.get_worldpath() .. '/'
-local load_map = dofile(modpath .. 'load.lua')
-local geometry = dofile(modpath .. 'geometry.lua')
 
-local function copy_if_needed(filename)
-	local wfilename = worldpath..filename
-	local wfile = io.open(wfilename, 'r')
-	if wfile then
-		wfile:close()
-		return
-	end
-	local mfilename = modpath..filename
-	local mfile = io.open(mfilename, 'r')
-	local wfile = io.open(wfilename, 'w')
-	wfile:write(mfile:read("*all"))
-	mfile:close()
-	wfile:close()
+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 use_biomegen_mod = use_biomes and minetest.global_exists('biomegen')
+use_biomes = use_biomes and not use_biomegen_mod
+
+if use_biomegen_mod then
+	biomegen.set_elevation_chill(elevation_chill)
 end
+dofile(modpath .. 'noises.lua')
 
-copy_if_needed('size')
-local sfile = io.open(worldpath..'size')
-local X = tonumber(sfile:read('*l'))
-local Z = tonumber(sfile:read('*l'))
-
-copy_if_needed('dem')
-local dem = load_map(worldpath..'dem', 2, true, X*Z)
-copy_if_needed('lakes')
-local lakes = load_map(worldpath..'lakes', 2, true, X*Z)
-copy_if_needed('bounds_x')
-local bounds_x = load_map(worldpath..'bounds_x', 4, false, (X-1)*Z)
-copy_if_needed('bounds_y')
-local bounds_z = load_map(worldpath..'bounds_y', 4, false, X*(Z-1))
-
-copy_if_needed('offset_x')
-local offset_x = load_map(worldpath..'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(worldpath..'offset_y', 1, true, X*Z)
-for k, v in ipairs(offset_z) do
-	offset_z[k] = (v+0.5)/256
-end
-
-
-local function index(x, z)
-	return z*X+x+1
-end
-
-local function get_point_location(x, z)
-	local i = index(x, z)
-	return x+offset_x[i], z+offset_z[i]
-end
+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)
@@ -62,246 +29,194 @@ end
 
 local data = {}
 
-local blocksize = 12
-local sea_level = 1
-local min_catchment = 25
-local max_catchment = 40000
-local riverbed_slope = 0.4
-
-local get_settings = dofile(modpath .. 'settings.lua')
-
-blocksize = get_settings('blocksize', 'int', blocksize)
-sea_level = get_settings('sea_level', 'int', sea_level)
-min_catchment = get_settings('min_catchment', 'float', min_catchment)
-max_catchment = get_settings('max_catchment', 'float', max_catchment)
-riverbed_slope = get_settings('riverbed_slope', 'float', riverbed_slope) * blocksize
-
--- 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_x_obj, noise_z_obj, noise_distort_obj, noise_heat_obj, noise_heat_blend_obj
+local noise_x_map = {}
+local noise_z_map = {}
+local noise_distort_map = {}
+local noise_heat_map = {}
+local noise_heat_blend_map = {}
+local mapsize
+local init = false
 
 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,
+		}
+		if use_distort then
+			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_distort_obj = minetest.get_perlin_map(mapgen_rivers.noise_params.distort_amplitude, chulens)
+		end
+		if use_biomes then
+			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)
+		end
+		init = true
+	end
+
+	local minp2d = {x=minp.x, y=minp.z}
+	if use_distort then
+		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)
+	end
+	if use_biomes then
+		noise_heat_obj:get_2d_map_flat(minp2d, noise_heat_map)
+		noise_heat_blend_obj:get_2d_map_flat(minp2d, noise_heat_blend_map)
+	end
+
+	local terrain_map, lake_map, incr, i_origin
+
+	if use_distort then
+		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}
+		incr = pmaxp.x-pminp.x+1
+		i_origin = 1 - pminp.z*incr - pminp.x
+		terrain_map, lake_map = heightmaps(pminp, pmaxp)
+	else
+		terrain_map, lake_map = heightmaps(minp, maxp)
+	end
+
 	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)
 
 	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.
-	local chulens = maxp.z - minp.z + 1
 
-	local polygons = {}
-	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)
+	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
 
-	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)
-			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 = {}
-			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)
-			for x=xmin, xmax do
-				bounds[x] = {}
+	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)
+			local ground_above = false
+			local temperature
+			if use_biomes then
+				temperature = noise_heat_map[i2d]+noise_heat_blend_map[i2d]
+			end
+			local terrain, lake
+			if not use_distort then
+				terrain = terrain_map[i2d]
+				lake = lake_map[i2d]
 			end
 
-			local i1 = 4
-			for i2=1, 4 do -- Loop on 4 edges
-				local x1, x2 = poly_x[i1], poly_x[i2]
-				local lxmin = math.floor(blocksize*math.min(x1, x2))+1
-				local lxmax = math.floor(blocksize*math.max(x1, x2))
-				if lxmin <= lxmax then
-					local z1, z2 = poly_z[i1], poly_z[i2]
-					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
-						table.insert(bounds[x], a*x+b)
-					end
-				end
-				i1 = i2
-			end
-			for x=xmin, xmax do
-				local xlist = bounds[x]
-				table.sort(xlist)
-				local c = math.floor(#xlist/2)
-				for l=1, c do
-					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
-						polygons[i] = polygon
-						i = i + 1
-					end
-				end
-			end
+			for y = maxp.y+1, minp.y, -1 do
+				if use_distort then
+					local xn = noise_x_map[nid]
+					local zn = noise_z_map[nid]
+					local x0 = math.floor(xn)
+					local z0 = math.floor(zn)
 
-			polygon.dem = {dem[iA], dem[iB], dem[iC], dem[iD]}
-			polygon.lake = math.min(lakes[iA], lakes[iB], lakes[iC], lakes[iD])
+					local i0 = i_origin + z0*incr + x0
+					local i1 = i0+1
+					local i2 = i1+incr
+					local i3 = i2-1
 
-			local river_west = river_width(bounds_z[iA])
-			local river_north = river_width(bounds_x[iA-zp])
-			local river_east = 1-river_width(bounds_z[iB])
-			local river_south = 1-river_width(bounds_x[iD-zp-1])
-			if river_west > river_east then
-				local mean = (river_west + river_east) / 2
-				river_west = mean
-				river_east = mean
-			end
-			if river_north > river_south then
-				local mean = (river_north + river_south) / 2
-				river_north = mean
-				river_south = mean
-			end
-			polygon.rivers = {river_west, river_north, river_east, river_south}
-
-			local around = {0,0,0,0,0,0,0,0}
-			if zp > 0 then
-				around[1] = river_width(bounds_z[iA-X])
-				around[2] = river_width(bounds_z[iB-X])
-			end
-			if xp < X-2 then
-				around[3] = river_width(bounds_x[iB-zp])
-				around[4] = river_width(bounds_x[iC-zp-1])
-			end
-			if zp < Z-2 then
-				around[5] = river_width(bounds_z[iC])
-				around[6] = river_width(bounds_z[iD])
-			end
-			if xp > 0 then
-				around[7] = river_width(bounds_x[iD-zp-2])
-				around[8] = river_width(bounds_x[iA-zp-1])
-			end
-			polygon.river_corners = {math.max(around[8], around[1]), math.max(around[2], around[3]), math.max(around[4], around[5]), math.max(around[6], around[7])}
-		end
-	end
-
-	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 = geometry.transform_quadri(poly.x, poly.z, x/blocksize, z/blocksize)
-				local i00, i01, i11, i10 = unpack(poly.i)
-
-				local is_river = false
-				local depth_factor = 0
-				local r_west, r_north, r_east, r_south = unpack(poly.rivers)
-				if xf >= r_east then
-					is_river = true
-					depth_factor = xf-r_east
-					xf = 1
-				elseif xf <= r_west then
-					is_river = true
-					depth_factor = r_west-xf
-					xf = 0
-				end
-				if zf >= r_south then
-					is_river = true
-					depth_factor = zf-r_south
-					zf = 1
-				elseif zf <= r_north then
-					is_river = true
-					depth_factor = r_north-zf
-					zf = 0
+					terrain = interp(terrain_map[i0], terrain_map[i1], terrain_map[i2], terrain_map[i3], xn-x0, zn-z0)
+					lake = math.min(lake_map[i0], lake_map[i1], lake_map[i2], lake_map[i3])
 				end
 
-				if not is_river then
-					local c_NW, c_NE, c_SE, c_SW = unpack(poly.river_corners)
-					if xf+zf <= c_NW then
-						is_river = true
-						depth_factor = c_NW-xf-zf
-						xf, zf = 0, 0
-					elseif 1-xf+zf <= c_NE then
-						is_river = true
-						depth_factor = c_NE-1+xf-zf
-						xf, zf = 1, 0
-					elseif 2-xf-zf <= c_SE then
-						is_river = true
-						depth_factor = c_SE-2+xf+zf
-						xf, zf = 1, 1
-					elseif xf+1-zf <= c_SW then
-						is_river = true
-						depth_factor = c_SW-xf-1+zf
-						xf, zf = 0, 1
-					end
-				end
+				if y <= maxp.y then
 
-				if not is_river then
-					xf = (xf-r_west) / (r_east-r_west)
-					zf = (zf-r_north) / (r_south-r_north)
-				end
-
-				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 is_river then
-					terrain_height = math.min(math.max(lake_height, sea_level) - math.floor(1+depth_factor*riverbed_slope), terrain_height)
-				end
-				local is_lake = lake_height > terrain_height
-				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
-							else
-								data[ivm] = c_lawn
-							end
-						else
+					local is_lake = lake > terrain
+					local ivm = a:index(x, y, z)
+					if y <= terrain then
+						if not use_biomes or y <= terrain-1 or ground_above then
 							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
 						end
-						ivm = ivm + ystride
+					elseif y <= lake and lake > sea_level then
+						if not use_biomes or temperature - y*elevation_chill >= 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
-					if is_lake and lake_height > minp.y then
-						for y=math.max(minp.y, terrain_height+1), math.min(maxp.y, lake_height) do
-							data[ivm] = c_rwater
-							ivm = ivm + ystride
-						end
-					end
-				else
-					for y=math.max(minp.y, terrain_height+1), math.min(maxp.y, sea_level) do
-						data[ivm] = c_water
-						ivm = ivm + ystride
-					end
+				ground_above = y <= terrain
+
+				ivm = ivm + ystride
+				if use_distort then
+					nid = nid + incrY
 				end
 			end
-			i = i + 1
+
+			if use_distort then
+				nid = nid + incrX
+			end
+			i2d = i2d + 1
+		end
+
+		if use_distort then
+			nid = nid + incrZ
 		end
 	end
 
-	vm:set_data(data)
-	minetest.generate_ores(vm, minp, maxp)
+	if use_biomegen_mod then
+		biomegen.generate_all(data, a, vm, minp, maxp, seed)
+	else
+		vm:set_data(data)
+		minetest.generate_ores(vm, minp, maxp)
+	end
+
 	vm:set_lighting({day = 0, night = 0})
 	vm:calc_lighting()
 	vm:update_liquids()

load.lua

@@ -1,5 +1,7 @@
+local worldpath = minetest.get_worldpath() .. "/river_data/"
+
 local function load_map(filename, bytes, signed, size)
-	local file = io.open(filename, 'r')
+	local file = io.open(worldpath .. filename, 'rb')
 	local data = file:read('*all')
 	if #data < bytes*size then
 		data = minetest.decompress(data)

mod.conf

@@ -0,0 +1,4 @@
+name = mapgen_rivers
+title = Map generator with realistic rivers
+depends = default
+optional_depends = biomegen

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,228 @@
+local modpath = minetest.get_modpath(minetest.get_current_modname()) .. '/'
+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/'
+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, '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()
+end
+
+copy_if_needed('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
+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
+
+local map_offset = {x=0, z=0}
+if mapgen_rivers.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 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)
+	print("Generating polygon map")
+	print(minp.x, maxp.x, minp.z, maxp.z)
+
+	if not init then
+		if glaciers then
+			noise_heat = minetest.get_perlin(mapgen_rivers.noise_params.heat)
+		end
+		init = true
+	end
+
+	local chulens = maxp.x - minp.x + 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+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
+		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)   * blocksize - map_offset.x,
+				(offset_x[iB]+xp+1) * blocksize - map_offset.x,
+				(offset_x[iC]+xp+1) * blocksize - map_offset.x,
+				(offset_x[iD]+xp)   * blocksize - map_offset.x,
+			}
+			local poly_z = {
+				(offset_z[iA]+zp)   * blocksize - map_offset.z,
+				(offset_z[iB]+zp)   * blocksize - map_offset.z,
+				(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)
+			-- Calculate the min and max Z positions 
+			local zmin = math.max(math.floor(math.min(unpack(poly_z)))+1, minp.z)
+			local zmax = math.min(math.floor(math.max(unpack(poly_z))), maxp.z)
+			-- And initialize the arrays
+			for z=zmin, zmax do
+				bounds[z] = {}
+			end
+
+			local i1 = 4
+			for i2=1, 4 do -- Loop on 4 edges
+				local z1, z2 = poly_z[i1], poly_z[i2]
+				-- Calculate the integer Z positions over which this edge spans
+				local lzmin = math.floor(math.min(z1, z2))+1
+				local lzmax = math.floor(math.max(z1, z2))
+				if lzmin <= lzmax then -- If there is at least one position in it
+					local x1, x2 = poly_x[i1], poly_x[i2]
+					-- Calculate coefficient of the equation defining the edge: X=aZ+b
+					local a = (x1-x2) / (z1-z2)
+					local b = (x1 - a*z1)
+					for z=math.max(lzmin, minp.z), math.min(lzmax, maxp.z) do
+						-- For every Z position involved, add the intercepted X position in the table
+						table.insert(bounds[z], a*z+b)
+					end
+				end
+				i1 = i2
+			end
+			for z=zmin, zmax do
+				-- Now sort the bounds list
+				local zlist = bounds[z]
+				table.sort(zlist)
+				local c = math.floor(#zlist/2)
+				for l=1, c do
+					-- Take pairs of X coordinates: all positions between them belong to the polygon.
+					local xmin = math.max(math.floor(zlist[l*2-1])+1, minp.x)
+					local xmax = math.min(math.floor(zlist[l*2]), maxp.x)
+					local i = (z-minp.z) * chulens + (xmin-minp.x) + 1
+					for x=xmin, xmax 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 = {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], poly_dem[1], poly_z[1]) < 0 then
+					riverA = math.min(riverA*glacier_factor, 1)
+				end
+				if get_temperature(poly_x[2], poly_dem[2], poly_z[2]) < 0 then
+					riverB = math.min(riverB*glacier_factor, 1)
+				end
+				if get_temperature(poly_x[3], poly_dem[3], poly_z[3]) < 0 then
+					riverC = math.min(riverC*glacier_factor, 1)
+				end
+				if get_temperature(poly_x[4], poly_dem[4], poly_z[4]) < 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

@@ -1,100 +0,0 @@
-import numpy as np
-import heapq
-import sys
-
-# Conventions:
-# 1 = South (+Y)
-# 2 = East  (+X)
-# 3 = North (-Y)
-# 4 = West  (-X)
-
-sys.setrecursionlimit(65536)
-
-neighbours_dirs = np.array([
-    [0,1,0],
-    [2,0,4],
-    [0,3,0],
-], dtype='u1')
-
-neighbours_pattern = neighbours_dirs > 0
-
-def flow_dirs_lakes(dem, random=0):
-    (Y, X) = dem.shape
-
-    dem_margin = np.zeros((Y+2, X+2))
-    dem_margin[1:-1,1:-1] = dem
-    if random > 0:
-        dem_margin += np.random.random(dem_margin.shape) * random
-
-    # Initialize: list map borders
-    borders = []
-
-    for x in range(1,X+1):
-        dem_north = dem_margin[1,x]
-        borders.append((dem_north, dem_north, 1, x))
-        dem_south = dem_margin[Y,x]
-        borders.append((dem_south, dem_south, Y, x))
-
-    for y in range(2,Y):
-        dem_west = dem_margin[y,1]
-        borders.append((dem_west, dem_west, y, 1))
-        dem_east = dem_margin[y,X]
-        borders.append((dem_east, dem_east, y, X))
-
-    heapq.heapify(borders)
-
-    dirs = np.zeros((Y+2, X+2), dtype='u1')
-    dirs[-2:,:] = 1
-    dirs[:,-2:] = 2
-    dirs[ :2,:] = 3
-    dirs[:, :2] = 4
-
-    lakes = np.zeros((Y, X))
-
-    def add_point(y, x, altmax):
-        alt = dem_margin[y, x]
-        heapq.heappush(borders, (alt, altmax, y, x))
-
-    while len(borders) > 0:
-        (alt, altmax, y, x) = heapq.heappop(borders)
-        neighbours = dirs[y-1:y+2, x-1:x+2]
-        empty_neighbours = (neighbours == 0) * neighbours_pattern
-        neighbours += empty_neighbours * neighbours_dirs
-
-        lake = max(alt, altmax)
-        lakes[y-1,x-1] = lake
-
-        coords = np.transpose(empty_neighbours.nonzero())
-        for (dy,dx) in coords-1:
-            add_point(y+dy, x+dx, lake)
-
-    return dirs[1:-1,1:-1], lakes
-
-def accumulate(dirs, dem=None):
-    (Y, X) = dirs.shape
-    dirs_margin = np.zeros((Y+2,X+2))-1
-    dirs_margin[1:-1,1:-1] = dirs
-    quantity = np.zeros((Y, X), dtype='i4')
-
-    def calculate_quantity(y, x):
-        if quantity[y,x] > 0:
-            return quantity[y,x]
-        q = 1
-        neighbours = dirs_margin[y:y+3, x:x+3]
-        donors = neighbours == neighbours_dirs
-
-        coords = np.transpose(donors.nonzero())
-        for (dy,dx) in coords-1:
-            q += calculate_quantity(y+dy, x+dx)
-        quantity[y, x] = q
-        return q
-
-    for x in range(X):
-        for y in range(Y):
-            calculate_quantity(y, x)
-
-    return quantity
-
-def flow(dem):
-    dirs, lakes = flow_dirs_lakes(dem)
-    return dirs, lakes, accumulate(dirs, dem)

settings.lua

@@ -9,6 +9,8 @@ local function get_settings(key, dtype, default)
 			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
 
@@ -20,8 +22,11 @@ local function get_settings(key, dtype, default)
 		elseif dtype == "float" then
 			conf_val = tonumber(conf_val)
 			storage:set_float(key, conf_val)
-		elseif dtype == "string" then
+		else
 			storage:set_string(key, conf_val)
+			if dtype == "bool" then
+				conf_val = conf_val == 'true'
+			end
 		end
 
 		return conf_val
@@ -32,10 +37,22 @@ local function get_settings(key, dtype, default)
 			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
 
-return get_settings
+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)

settingtypes.txt

@@ -0,0 +1,52 @@
+# 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
+
+#    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
+
+#    Sea level used by mapgen_rivers
+mapgen_rivers_sea_level (Sea level) int 1
+
+#    Minimal catchment area for a river to be drawn, in grid cells
+#    (1 cell = blocksize x blocksize).
+#    Lower value means bigger river density
+mapgen_rivers_min_catchment (Minimal catchment area) float 25.0 1.0 1000.0
+
+#    Catchment area in grid cells (1 grid cell = blocksize x blocksize)
+#    at which rivers reach their maximal width of 2*blocksize.
+#    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
+
+#    Lateral slope of the riverbed.
+#    Higher value means deeper rivers.
+mapgen_rivers_riverbed_slope (Riverbed slope) float 0.4 0.0 2.0
+
+#    Enable horizontal distorsion (shearing) of landscape, to break the regularity
+#    of grid cells and allow overhangs.
+#    Distorsion uses two 3D noises and thus is intensive in terms of computing time.
+mapgen_rivers_distort (Distorsion) bool true
+
+#    Enable biome generation.
+#    If 'biomegen' mod is installed, 'mapgen_rivers' will generate biomes from the
+#    native biome system. If 'biomegen' is not present, will generate only grass and
+#    snow.
+mapgen_rivers_biomes (Biomes) bool true
+
+#    Whether to enable glaciers.
+#    Glaciers are widened river sections, covered by ice, that are generated in
+#    very cold areas.
+mapgen_rivers_glaciers (Glaciers) bool false
+
+#    River channels are widened by this factor if they are a glacier.
+mapgen_rivers_glacier_widening_factor (Glacier widening factor) float 8.0 1.0 20.0
+
+#    Temperature value decreases by this quantity for every node, vertically.
+#    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.

terrain_default.conf

@@ -0,0 +1,17 @@
+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

@@ -0,0 +1,17 @@
+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

@@ -0,0 +1,16 @@
+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

terrain_rivers.py

@@ -1,118 +0,0 @@
-#!/usr/bin/env python3
-
-import numpy as np
-import noise
-from save import save
-from erosion import EvolutionModel
-import bounds
-import os
-import sys
-
-# Always place in this script's parent directory
-os.chdir(os.path.dirname(sys.argv[0]))
-argc = len(sys.argv)
-
-if argc > 1:
-	mapsize = int(sys.argv[1])
-else:
-	mapsize = 400
-
-scale = mapsize / 2
-n = np.zeros((mapsize, mapsize))
-
-#micronoise_depth = 0.05
-
-params = {
-    "octaves" : int(np.log2(mapsize)),
-    "persistence" : 0.5,
-    "lacunarity" : 2.,
-}
-
-xbase = np.random.randint(65536)
-ybase = np.random.randint(65536)
-
-for x in range(mapsize):
-    for y in range(mapsize):
-        n[x,y] = noise.snoise2(x/scale + xbase, y/scale + ybase, **params)
-
-#micronoise = np.random.rand(mapsize, mapsize)
-#nn = np.exp(n*2) + micronoise*micronoise_depth
-nn = n*mapsize/5 + mapsize/20
-
-print('Initializing model')
-model = EvolutionModel(nn, K=1, m=0.35, d=1, sea_level=0)
-
-print('Flow calculation 1')
-model.calculate_flow()
-
-print('Advection 1')
-model.advection(2)
-
-print('Isostatic equilibration 1')
-model.adjust_isostasy()
-
-print('Flow calculation 2')
-model.calculate_flow()
-
-print('Diffusion')
-model.diffusion(4)
-
-print('Advection 2')
-model.advection(2)
-
-print('Isostatic equilibration 2')
-model.adjust_isostasy()
-
-print('Flow calculation 3')
-model.calculate_flow()
-
-print('Done')
-
-bx, by = bounds.make_bounds(model.dirs, model.rivers)
-ox, oy = bounds.twist(bx, by, bounds.get_fixed(model.dirs))
-
-offset_x = np.clip(np.floor(ox * 256), -128, 127)
-offset_y = np.clip(np.floor(oy * 256), -128, 127)
-
-save(model.dem, 'dem', dtype='>i2')
-save(model.lakes, 'lakes', dtype='>i2')
-save(np.abs(bx), 'bounds_x', dtype='>i4')
-save(np.abs(by), 'bounds_y', dtype='>i4')
-save(offset_x, 'offset_x', dtype='i1')
-save(offset_y, 'offset_y', dtype='i1')
-
-save(model.rivers, 'rivers', dtype='>u4')
-
-with open('size', 'w') as sfile:
-    sfile.write('{:d}\n{:d}'.format(mapsize, mapsize))
-
-try:
-    import matplotlib.pyplot as plt
-
-    plt.subplot(2,2,1)
-    plt.pcolormesh(nn, cmap='viridis')
-    plt.gca().set_aspect('equal', 'box')
-    #plt.colorbar(orientation='horizontal')
-    plt.title('Raw elevation')
-
-    plt.subplot(2,2,2)
-    plt.pcolormesh(model.lakes, cmap='viridis')
-    plt.gca().set_aspect('equal', 'box')
-    #plt.colorbar(orientation='horizontal')
-    plt.title('Lake surface elevation')
-
-    plt.subplot(2,2,3)
-    plt.pcolormesh(model.dem, cmap='viridis')
-    plt.gca().set_aspect('equal', 'box')
-    #plt.colorbar(orientation='horizontal')
-    plt.title('Elevation after advection')
-
-    plt.subplot(2,2,4)
-    plt.pcolormesh(model.rivers, vmin=0, vmax=mapsize**2/25, cmap='Blues')
-    plt.gca().set_aspect('equal', 'box')
-    #plt.colorbar(orientation='horizontal')
-    plt.title('Rivers discharge')
-
-    plt.show()
-except:
-    pass

terrainlib/__init__.py

@@ -0,0 +1,7 @@
+# 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,10 +1,13 @@
 import numpy as np
-import matplotlib.pyplot as plt
 
 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_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]
@@ -14,6 +17,10 @@ def make_bounds(dirs, rivers):
     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
@@ -28,6 +35,11 @@ def get_fixed(dirs):
     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
@@ -55,7 +67,7 @@ def twist(bounds_x, bounds_y, fixed, d=0.1, n=5):
 
         length = np.hypot(force_x, force_y)
         length[length==0] = 1
-        coeff = d / length * moveable
+        coeff = d / length * moveable # Normalize, take into account the direction only
         offset_x += force_x * coeff
         offset_y += force_y * coeff
 

terrainlib/erosion.py

@@ -0,0 +1,121 @@
+import numpy as np
+import scipy.ndimage as im
+import scipy.signal as si
+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
+
+second_derivative_matrix = np.array([
+    [0., 0.25, 0.],
+    [0.25,-1., 0.25],
+    [0., 0.25, 0.],
+])
+
+diff_max = 1.0
+
+def diffusion(dem, time, d=1.0):
+    if isinstance(d, np.ndarray):
+        dmax = d.max()
+    else:
+        dmax = d
+    diff = time * dmax
+    print(diff)
+    niter = int(diff//diff_max) + 1
+    ddiff = d * (time / niter)
+
+    #print('{:d} iterations'.format(niter))
+    for i in range(niter):
+        dem[1:-1,1:-1] += si.convolve2d(dem, second_derivative_matrix, mode='valid') * ddiff
+        #print('iteration {:d}'.format(i+1))
+
+    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
+        if isinstance(d, np.ndarray):
+            self.d = d[1:-1,1:-1]
+        else:
+            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, dem=None):
+        if dem is None:
+            dem = self.dem
+        self.ref_isostasy = im.gaussian_filter(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

@@ -0,0 +1,278 @@
+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/settings.py

@@ -0,0 +1,16 @@
+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/view.py

@@ -0,0 +1,99 @@
+#!/usr/bin/env python3
+
+import numpy as np
+import sys, traceback
+
+has_matplotlib = True
+try:
+    import matplotlib.colors as mcl
+    import matplotlib.pyplot as plt
+    try:
+        import colorcet as cc
+        cmap1 = cc.cm.CET_L11
+        cmap2 = cc.cm.CET_L12
+    except ImportError: # No module colorcet
+        import matplotlib.cm as cm
+        cmap1 = cm.summer
+        cmap2 = cm.Blues
+except ImportError: # No module matplotlib
+    has_matplotlib = False
+
+if has_matplotlib:
+    def view_map(dem, lakes, scale=1, sea_level=0.0, title=None):
+        lakes_sea = np.maximum(lakes, sea_level)
+        water = np.maximum(lakes_sea - dem, 0)
+        max_elev = dem.max()
+        max_depth = water.max()
+
+        ls = mcl.LightSource(azdeg=315, altdeg=45)
+        norm_ground = plt.Normalize(vmin=sea_level, vmax=max_elev)
+        norm_sea = plt.Normalize(vmin=0, vmax=max_depth)
+        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)
+        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')
+
+        sm1 = plt.cm.ScalarMappable(cmap=cmap1, norm=norm_ground)
+        plt.colorbar(sm1).set_label('Elevation')
+
+        sm2 = plt.cm.ScalarMappable(cmap=cmap2, norm=norm_sea)
+        plt.colorbar(sm2).set_label('Water depth')
+
+        plt.xlabel('X')
+        plt.ylabel('Z')
+
+        if title is not None:
+            plt.title(title, fontweight='bold')
+
+    def update(*args, t=0.01, **kwargs):
+        try:
+            plt.clf()
+            view_map(*args, **kwargs)
+            plt.pause(t)
+        except:
+            traceback.print_exception(*sys.exc_info())
+
+    def plot(*args, **kwargs):
+        try:
+            plt.clf()
+            view_map(*args, **kwargs)
+            plt.pause(0.01)
+            plt.show()
+        except Exception as e:
+            traceback.print_exception(*sys.exc_info())
+
+else:
+    def update(*args, **kwargs):
+        pass
+    def plot(*args, **kwargs):
+        pass
+
+def stats(dem, lakes, scale=1):
+    surface = dem.size
+
+    continent = np.maximum(dem, lakes) >= 0
+    continent_surface = continent.sum()
+
+    lake = continent & (lakes>dem)
+    lake_surface = lake.sum()
+
+    print('---   General    ---')
+    print('Grid size:    {:5d}x{:5d}'.format(dem.shape[0], dem.shape[1]))
+    if scale > 1:
+        print('Map size:     {:5d}x{:5d}'.format(int(dem.shape[0]*scale), int(dem.shape[1]*scale)))
+    print()
+    print('---   Surfaces   ---')
+    print('Continents:        {:6.2%}'.format(continent_surface/surface))
+    print('-> Ground:         {:6.2%}'.format((continent_surface-lake_surface)/surface))
+    print('-> Lakes:          {:6.2%}'.format(lake_surface/surface))
+    print('Oceans:            {:6.2%}'.format(1-continent_surface/surface))
+    print()
+    print('---  Elevations  ---')
+    print('Mean elevation:      {:4.0f}'.format(dem.mean()))
+    print('Mean ocean depth:    {:4.0f}'.format((dem*~continent).sum()/(surface-continent_surface)))
+    print('Mean continent elev: {:4.0f}'.format((dem*continent).sum()/continent_surface))
+    print('Lowest elevation:    {:4.0f}'.format(dem.min()))
+    print('Highest elevation:   {:4.0f}'.format(dem.max()))

view_map.py

@@ -2,38 +2,28 @@
 
 import numpy as np
 import zlib
-import matplotlib.pyplot as plt
+import sys
+import os
+
+from terrainlib import stats, plot
+
+scale = 1
+if len(sys.argv) > 1:
+    os.chdir(sys.argv[1])
+if len(sys.argv) > 2:
+    scale = int(sys.argv[2])
 
 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)
+        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')
-n = shape[0] * shape[1]
 dem = load_map('dem', '>i2', shape)
 lakes = load_map('lakes', '>i2', shape)
-rivers = load_map('rivers', '>u4', shape)
 
-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(np.log(rivers), vmin=0, vmax=np.log(n/25), cmap='Blues')
-plt.gca().set_aspect('equal', 'box')
-#plt.colorbar(orientation='horizontal')
-plt.title('Rivers discharge')
-
-plt.show()
+stats(dem, lakes, scale=scale)
+plot(dem, lakes, scale=scale)