Updated screenshot

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

.gitignore

@@ -9,4 +9,4 @@ bounds_y
 dirs
 rivers
 unused/
-data/
+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,36 +1,40 @@
-mapgen_rivers
-=============
+# Map Generator with Rivers
+`mapgen_rivers v1.0` by Gaël de Sailly.
 
-Procedural map generator for Minetest 5.x. Focused on river networks, and features valley erosion and lakes.
+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).
 
-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/79541028-687b3000-8089-11ea-9209-c23c15d75383.png)
+![Screenshot](https://content.minetest.net/uploads/fff09f2269.png)
+
+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 (provides biome system).
 
 # 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` and `scipy`, widely used libraries for numerical calculations
-- `noise`, doing Perlin/Simplex noises
-- optionally, `matplotlib` (for map preview)
-
-They are commonly found on `pip` or `conda` Python distributions.
+This mod should be placed in the `mods/` directory of Minetest like any other mod.
 
 # Usage
-By default, the mod contains a demo 400x400 map (so you can start the game directly), but it is recommended that you run the pre-processing script to generate a new map before world creation, if you can.
+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.
 
-## 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.
+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).
 
-If you have `matplotlib` installed, the script `view_map.py` can be used to get a map preview. Example:
-```
-./view_map.py data/
-```
+## Settings
+Settings can be found in Minetest in the `Settings` tab, `All settings` -> `Mods` -> `mapgen_rivers`.
 
-## 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.
+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
+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.
+
+It can be run from command line by passing the world folder. Example:
+```
+./view_map.py ~/.minetest/worlds/test_mg_rivers
+```

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

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

erosion.py

@@ -1,97 +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):
-    """
-    Simulate erosion by rivers.
-    This models erosion as an upstream advection of elevations ("erosion waves").
-    Advection speed depends on water flux and parameters:
-
-    v = K * flux^m
-    """
-
-    dirs = dirs.copy()
-    dirs[0,:] = 0
-    dirs[-1,:] = 0
-    dirs[:,0] = 0
-    dirs[:,-1] = 0
-
-    adv_time = 1 / (K*rivers**m) # 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 np.minimum(dem, dem_new)
-
-def diffusion(dem, time, d=1):
-    radius = d * time**.5
-    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):
-        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') # 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

heightmap.lua

@@ -1,11 +1,12 @@
-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
 
 -- Linear interpolation
 local function interp(v00, v01, v11, v10, xf, zf)
@@ -17,19 +18,18 @@ end
 local function heightmaps(minp, maxp)
 
 	local polygons = make_polygons(minp, maxp)
-	local incr = maxp.x-minp.x+1
-	local i0 = (minp.z-minp.z) * incr + (minp.x-minp.x) + 1
+	local incr = maxp.z-minp.z+1
 
 	local terrain_height_map = {}
 	local lake_height_map = {}
 
 	local i = 1
-	for x=minp.x, maxp.x do
 	for z=minp.z, maxp.z do
+		for x=minp.x, maxp.x do
 			local poly = polygons[i]
 
 			if poly then
-				local xf, zf = transform_quadri(poly.x, poly.z, x/blocksize, z/blocksize)
+				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
@@ -98,7 +98,25 @@ local function heightmaps(minp, maxp)
 					xf, zf
 				))
 
-				local lake_height = math.max(math.floor(poly.lake), terrain_height)
+				-- 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
@@ -106,8 +124,8 @@ local function heightmaps(minp, maxp)
 				terrain_height_map[i] = terrain_height
 				lake_height_map[i] = lake_height
 			else
-				terrain_height_map[i] = -31000
-				lake_height_map[i] = -31000
+				terrain_height_map[i] = MAP_BOTTOM
+				lake_height_map[i] = MAP_BOTTOM
 			end
 			i = i + 1
 		end

init.lua

@@ -1,20 +1,23 @@
 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 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
 
+if use_biomegen_mod then
+	biomegen.set_elevation_chill(elevation_chill)
+end
 dofile(modpath .. 'noises.lua')
 
-local make_polygons = dofile(modpath .. 'polygons.lua')
-
-local transform_quadri = dofile(modpath .. 'geometry.lua')
-
 local heightmaps = dofile(modpath .. 'heightmap.lua')
 
 -- Linear interpolation
@@ -35,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,
@@ -48,21 +57,33 @@ local function generate(minp, maxp, seed)
 			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)
-		noise_distort_obj = minetest.get_perlin_map(mapgen_rivers.noise_params.distort_amplitude, chulens)
+		end
 		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)
 		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
@@ -87,9 +108,12 @@ local function generate(minp, maxp, seed)
 
 		local pminp = {x=math.floor(xmin), z=math.floor(zmin)}
 		local pmaxp = {x=math.floor(xmax)+1, z=math.floor(zmax)+1}
-	local incr = pmaxp.z-pminp.z+1
-	local i_origin = 1 - pminp.x*incr - pminp.z
-	local terrain_map, lake_map = heightmaps(pminp, pmaxp)
+		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")
@@ -118,26 +142,38 @@ local function generate(minp, maxp, seed)
 		for x = minp.x, maxp.x do
 			local ivm = a:index(x, minp.y, z)
 			local ground_above = false
-			local temperature = noise_heat_map[i2d]+noise_heat_blend_map[i2d]
+			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
+
 			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)
 
-				local i0 = i_origin + x0*incr + z0
-				local i1 = i0+incr
-				local i2 = i1+1
-				local i3 = i0+1
+					local i0 = i_origin + z0*incr + x0
+					local i1 = i0+1
+					local i2 = i1+incr
+					local i3 = i2-1
+
+					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
 
-				local terrain = interp(terrain_map[i0], terrain_map[i1], terrain_map[i2], terrain_map[i3], xn-x0, zn-z0)
 				if y <= maxp.y then
-					local lake = math.min(lake_map[i0], lake_map[i1], lake_map[i2], lake_map[i3])
 
 					local is_lake = lake > terrain
 					local ivm = a:index(x, y, z)
 					if y <= terrain then
-						if y <= terrain-1 or ground_above 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
@@ -152,8 +188,7 @@ local function generate(minp, maxp, seed)
 							end
 						end
 					elseif y <= lake and lake > sea_level then
-						local temperature_y = temperature - y*elevation_chill
-						if temperature_y >= 0 then
+						if not use_biomes or temperature - y*elevation_chill >= 0 then
 							data[ivm] = c_rwater
 						else
 							data[ivm] = c_ice
@@ -166,20 +201,45 @@ local function generate(minp, maxp, seed)
 				ground_above = y <= terrain
 
 				ivm = ivm + ystride
+				if use_distort then
 					nid = nid + incrY
 				end
-			nid = nid + incrX
-			i2d = i2d + 1
-		end
-		nid = nid + incrZ
 			end
 
+			if use_distort then
+				nid = nid + incrX
+			end
+			i2d = i2d + 1
+		end
+
+		if use_distort then
+			nid = nid + incrZ
+		end
+	end
+
+	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()
 	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)
-	local file = io.open(worldpath .. filename, 'r')
+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

mod.conf

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

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,70 +1,98 @@
-local modpath = minetest.get_modpath(minetest.get_current_modname()) .. '/'
-local mod_data_path = modpath .. 'data/'
+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, 'r')
-	if wfile then
-		wfile:close()
-		return
-	end
-	local mfilename = mod_data_path..filename
-	local mfile = io.open(mfilename, 'r')
-	local wfile = io.open(wfilename, 'w')
-	wfile:write(mfile:read("*all"))
-	mfile:close()
-	wfile:close()
+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 sfile = io.open(world_data_path..'size')
-local X = tonumber(sfile:read('*l'))
-local Z = tonumber(sfile:read('*l'))
+local load_all = mapgen_rivers.settings.load_all
+
+-- Try to read file 'size'
+local sfile = io.open(world_data_path..'size', 'r')
+local first_mapgen = true
+if sfile then
+	X, Z = tonumber(sfile:read('*l')), 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
+	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
 
--- 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 map_offset = {x=0, z=0}
+if mapgen_rivers.settings.center then
+	map_offset.x = blocksize*X/2
+	map_offset.z = blocksize*Z/2
+end
+
+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
@@ -75,20 +103,30 @@ 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)
+
+	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
 			noise_heat = minetest.get_perlin(mapgen_rivers.noise_params.heat)
@@ -96,12 +134,12 @@ local function make_polygons(minp, maxp)
 		init = true
 	end
 
-	local chulens = maxp.z - minp.z + 1
+	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/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 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)
 
 	-- Iterate over the polygons
 	for xp = xpmin, xpmax do
@@ -111,48 +149,58 @@ local function make_polygons(minp, maxp)
 			local iC = index(xp+1, zp+1)
 			local iD = index(xp, zp+1)
 			-- Extract the vertices of the polygon
-			local poly_x = {offset_x[iA]+xp, offset_x[iB]+xp+1, offset_x[iC]+xp+1, offset_x[iD]+xp}
-			local poly_z = {offset_z[iA]+zp, offset_z[iB]+zp, offset_z[iC]+zp+1, offset_z[iD]+zp+1}
+			local 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,
+			}
 			local polygon = {x=poly_x, z=poly_z, i={iA, iB, iC, iD}}
 
-			local bounds = {} -- Will be a list of the intercepts of polygon edges for every X position (scanline algorithm)
-			-- Calculate the min and max X positions 
-			local xmin = math.max(math.floor(blocksize*math.min(unpack(poly_x)))+1, minp.x)
-			local xmax = math.min(math.floor(blocksize*math.max(unpack(poly_x))), maxp.x)
+			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 x=xmin, xmax do
-				bounds[x] = {}
+			for z=zmin, zmax do
+				bounds[z] = {}
 			end
 
 			local i1 = 4
 			for i2=1, 4 do -- Loop on 4 edges
-				local x1, x2 = poly_x[i1], poly_x[i2]
-				-- Calculate the integer X positions over which this edge spans
-				local lxmin = math.floor(blocksize*math.min(x1, x2))+1
-				local lxmax = math.floor(blocksize*math.max(x1, x2))
-				if lxmin <= lxmax then -- If there is at least one position in it
 				local z1, z2 = poly_z[i1], poly_z[i2]
-					-- Calculate coefficient of the equation defining the edge: Z=aX+b
-					local a = (z1-z2) / (x1-x2)
-					local b = blocksize*(z1 - a*x1)
-					for x=math.max(lxmin, minp.x), math.min(lxmax, maxp.x) do
-						-- For every X position involved, add the intercepted Z position in the table
-						table.insert(bounds[x], a*x+b)
+				-- 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 x=xmin, xmax do
-				-- Now sort the bounds list
-				local xlist = bounds[x]
-				table.sort(xlist)
-				local c = math.floor(#xlist/2)
-				for l=1, c do
-					-- Take pairs of Z coordinates: all positions between them belong to the polygon.
-					local zmin = math.max(math.floor(xlist[l*2-1])+1, minp.z)
-					local zmax = math.min(math.floor(xlist[l*2]), maxp.z)
-					local i = (x-minp.x) * chulens + (zmin-minp.z) + 1
 			for z=zmin, zmax do
+				-- 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
@@ -162,7 +210,7 @@ local function make_polygons(minp, maxp)
 
 			local poly_dem = {dem[iA], dem[iB], dem[iC], dem[iD]}
 			polygon.dem = poly_dem
-			polygon.lake = math.min(lakes[iA], lakes[iB], lakes[iC], lakes[iD])
+			polygon.lake = {lakes[iA], lakes[iB], lakes[iC], lakes[iD]}
 
 			-- Now, rivers.
 			-- Load river flux values for the 4 corners
@@ -171,16 +219,16 @@ local function make_polygons(minp, maxp)
 			local riverC = river_width(rivers[iC])
 			local riverD = river_width(rivers[iD])
 			if glaciers then -- Widen the river
-				if get_temperature(poly_x[1]*blocksize, poly_dem[1], poly_z[1]*blocksize) < 0 then
+				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]*blocksize, poly_dem[2], poly_z[2]*blocksize) < 0 then
+				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]*blocksize, poly_dem[3], poly_z[3]*blocksize) < 0 then
+				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]*blocksize, poly_dem[4], poly_z[4]*blocksize) < 0 then
+				if get_temperature(poly_x[4], poly_dem[4], poly_z[4]) < 0 then
 					riverD = math.min(riverD*glacier_factor, 1)
 				end
 			end

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'

rivermapper.py

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

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)

settings.lua

@@ -1,53 +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'
+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
+	if mapgen_rivers.version ~= previous_version_mgr then
+		compat_mgr(mgrsettings)
 	end
 end
 
-	local conf_val = settings:get('mapgen_rivers_' .. key)
-	if conf_val then
-		if dtype == "int" then
-			conf_val = tonumber(conf_val)
-			storage:set_int(key, conf_val)
-		elseif dtype == "float" then
-			conf_val = tonumber(conf_val)
-			storage:set_float(key, conf_val)
-		elseif dtype == "string" or dtype == "bool" then
-			storage:set_string(key, conf_val)
-		end
+mtsettings:set("mapgen_rivers_version", mapgen_rivers.version)
+mgrsettings:set("version", mapgen_rivers.version)
 
-		return conf_val
+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
-		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)
+			return v
+		end
 	elseif dtype == "bool" then
-			storage:set_string(key, tostring(default))
+		local v = mgrsettings:get_bool(name)
+		if v == nil then
+			v = mtsettings:get_bool('mapgen_rivers_' .. name)
+			if v == nil then
+				v = default
 			end
-
-		return default
+			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
 
-mapgen_rivers.blocksize = get_settings('blocksize', 'int', 12)
-mapgen_rivers.sea_level = get_settings('sea_level', 'int', 1)
-mapgen_rivers.min_catchment = get_settings('min_catchment', 'float', 25)
-mapgen_rivers.max_catchment = get_settings('max_catchment', 'float', 40000)
-mapgen_rivers.riverbed_slope = get_settings('riverbed_slope', 'float', 0.4) * mapgen_rivers.blocksize
---mapgen_rivers.distort = get_settings('distort', 'bool', true) To be implemented: should be possible to disable distorsion
-mapgen_rivers.glaciers = get_settings('glaciers', 'bool', true)
-mapgen_rivers.glacier_factor = get_settings('glacier_factor', 'float', 8)
-mapgen_rivers.elevation_chill = get_settings('elevation_chill', 'float', 0.25)
+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

@@ -0,0 +1,117 @@
+# 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 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 12000.
+mapgen_rivers_blocksize (Block size) float 15.0 2.0 100.0
+
+#    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
+
+#    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 3600.0 100.0 1000000.0
+
+#    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.
+#    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.
+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
+
+#    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_rivers.py

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

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

@@ -0,0 +1,102 @@
+#!/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, 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()
+        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
+        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')
+
+        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,35 +2,23 @@
 
 import numpy as np
 import zlib
-import matplotlib.colors as mcol
-import matplotlib.pyplot as plt
-
-def view_map(dem, lakes, rivers):
-    plt.subplot(1,3,1)
-    plt.pcolormesh(dem, cmap='viridis')
-    plt.gca().set_aspect('equal', 'box')
-    plt.colorbar(orientation='horizontal')
-    plt.title('Raw elevation')
-
-    plt.subplot(1,3,2)
-    plt.pcolormesh(lakes, cmap='viridis')
-    plt.gca().set_aspect('equal', 'box')
-    plt.colorbar(orientation='horizontal')
-    plt.title('Lake surface elevation')
-
-    plt.subplot(1,3,3)
-    plt.pcolormesh(rivers, cmap='Blues', norm=mcol.LogNorm())
-    plt.gca().set_aspect('equal', 'box')
-    plt.colorbar(orientation='horizontal')
-    plt.title('Rivers flux')
-
-    plt.show()
-
-if __name__ == "__main__":
 import sys
 import os
-    if len(sys.argv) > 1:
+
+from view import stats, plot
+from readconfig import read_conf_file
+
 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)
@@ -40,9 +28,10 @@ if __name__ == "__main__":
             data = zlib.decompress(data)
         return np.frombuffer(data, dtype=dtype).reshape(shape)
 
-    shape = np.loadtxt('size', dtype='u4')
-    dem = load_map('dem', '>i2', shape)
-    lakes = load_map('lakes', '>i2', shape)
-    rivers = load_map('rivers', '>u4', shape)
+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)
 
-    view_map(dem, lakes, rivers)
+stats(dem, lakes, scale=blocksize)
+plot(dem, lakes, scale=blocksize, center=center)