46 Commits

Author SHA1 Message Date
e9fa7f9a5c Adapt to multithreading, now it is working! 2024-01-05 00:02:30 +01:00
f26fd1cccb Fixups 2024-01-04 23:54:59 +01:00
e5b8f2b3b8 Refactor to use only Settings objects for settings.
Also remove globally stored path for modpath and data.
Basically use less global variables for better adaptation to multithreading.
2024-01-04 22:47:50 +01:00
6fa1852277 Replace table settings by the mod's Settings object (to be completed) 2023-12-30 17:38:51 +01:00
11e6e72324 Started refactoring code structure to support a separate mapgen thread
Probably not working for now (untested), still some things to do
Also making the code cleaner
2023-12-28 15:27:36 +01:00
1ad8c96b8c Remove 'default' hard dependency 2022-01-17 23:20:34 +01:00
2f7098d752 Bump version (1.0.2) and add changelog 2022-01-10 12:44:33 +01:00
942a869b9f Minor fix in README 2022-01-10 12:32:38 +01:00
b3d79eaaf8 Added more comments in terrainlib_lua 2022-01-07 14:48:36 +01:00
68c19c3b94 terrainlib_lua: replaced space indents by tabs 2022-01-06 15:36:31 +01:00
417ce1bcbc Use builtin logging system and appropriate loglevels 2022-01-03 16:33:56 +01:00
c3a798933f Localize all global functions in load.lua and geometry.lua 2022-01-03 16:20:51 +01:00
0c98fc0881 Skip chunks that are fully higher than ground level 2022-01-03 16:18:48 +01:00
cb71f4400a Corrected mistake in settingtypes 2022-01-03 12:04:49 +01:00
f8f467ac3f Use local variables for math.* functions
and remove an unnecessary index calculation
2022-01-03 11:56:16 +01:00
2e29474686 Bump version (1.0.1) 2021-09-14 15:08:29 +02:00
27670addb3 Switch to singlenode mapgen if not done 2021-09-07 11:59:33 +02:00
54b94e6485 Updated screenshot 2021-08-01 18:43:44 +02:00
09de0fd298 Added reference to scientific paper 2021-07-28 19:50:27 +02:00
e3cadcdbc6 Collect garbage during pre-generation to free some space and reduce OOM problems 2021-07-26 22:28:45 +02:00
f7bc5ee0b4 Added logs during pre-generation 2021-07-25 10:14:15 +02:00
2c5e0ee5af Code style consistency: use tabs for indentation in all Lua code 2021-07-24 18:55:13 +02:00
95e87f8820 Updated README.md and added environment.yml
Switching to full Lua incredibly simplified the amount of stuff needed in README!
2021-07-24 18:49:17 +02:00
db33e58f0a Make map size a setting 2021-07-24 17:41:17 +02:00
31c5ea1025 Python map viewing: read conf file, and take world folder as input 2021-07-24 17:21:21 +02:00
c2c397c2a5 Added compatibility script
to update parameter values coming from older versions
2021-07-24 13:21:06 +02:00
9386ef51f1 Changed the way river width is determined in settings
- min_catchment now in square nodes instead of cells
- River widening power as input instead of calculating it from max_catchment
2021-07-24 13:18:58 +02:00
8ce20816e1 Some changes in default settings
- blocksize = 15 by default
- base noise is eased
- added tectonic compensation radius in settings
2021-07-24 10:23:06 +02:00
32bc9561b6 Made interactive loading optional but enabled by default 2021-06-26 16:02:00 +02:00
7e39189368 Load data on request instead of loading everything at first.
Data is cached in memory in case it is reused.
Offset values have a callback that converts them to the range ±0.5
2021-06-26 13:17:09 +02:00
636773487a Optimized load.lua by avoiding multi-argument packing 2021-06-26 11:49:26 +02:00
9cda649c93 Fixed time statistics 2021-06-26 11:26:45 +02:00
ecd1f0e08f Added time statistics and removed debug prints 2021-06-25 21:05:14 +02:00
5898354dbe Changed some default parameters
Map centered by default
Base noise higher and less wide horizontally
2021-06-25 21:05:03 +02:00
e14bc5216c Removed useless and deprecated files 2021-06-23 19:56:37 +02:00
a313244d07 Rename generate.lua -> pregenerate.lua 2021-06-23 19:53:56 +02:00
0de2f746cf Update settingtypes.txt 2021-06-06 17:17:28 +02:00
de8d685471 New settings system, use a conf file instead of datastorage
Added many missing settings, including pregeneration-related ones
TODO: update settingtypes.txt
2021-06-06 13:25:43 +02:00
51f3a2719d Generate and load map after mod loading
This has needed to globalize map tables
2021-06-05 11:24:28 +02:00
b02387944d Pre-generation: reverse X and Y directions everywhere
to make directions compatible with the mapgen code
2021-06-04 18:24:06 +02:00
74733549df Various bugfixes and workarounds
Now working in pure Lua!
Some parts of the code are very hacky (e.g. noise) and the way new and old codes have been glued together is sometimes to be rewritten.
But at least it works.
2021-06-03 23:30:04 +02:00
cb297af047 Add all code for generating a grid on world creation.
Not tested ; will likely need much testing and bugfix.
2021-06-03 20:08:57 +02:00
19efeaaff6 Globalize modpath and worldpath 2021-06-03 20:04:08 +02:00
0427b42d17 Removed Python terrainlib 2021-06-02 18:56:46 +02:00
7495d8a690 Added grid twisting (twist.lua)
Equivalent of Python terrainlib's 'bounds.py'
2021-06-02 18:42:40 +02:00
c99b8338e0 Lua Terrainlib: added first Lua files for erosion and flow routing
Tested, but not linked with the mod, yet.
2021-06-01 19:07:09 +02:00
40 changed files with 1855 additions and 1280 deletions

18
CHANGELOG.md Normal file
View File

@ -0,0 +1,18 @@
CHANGELOG
=========
## `v1.0.2` (2022-01-10)
- Use builtin logging system and appropriate loglevels
- Skip empty chunks, when generating high above ground (~20% speedup)
- Minor optimizations (turning global variables to local...)
## `v1.0.1` (2021-09-14)
- Automatically switch to `singlenode` mapgen at init time
## `v1.0` (2021-08-01)
- Rewritten pregen code (terrainlib) in pure Lua
- Optimized grid loading
- Load grid nodes on request by default
- Changed river width settings
- Added map size in settings
- Added logs

113
README.md
View File

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

34
compatibility.lua Normal file
View File

@ -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

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

View File

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

10
environment.yml Normal file
View File

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

View File

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

View File

@ -1,3 +1,6 @@
local sqrt, abs = math.sqrt, math.abs
local unpk = unpack
local function distance_to_segment(x1, y1, x2, y2, x, y)
-- get the distance between point (x,y) and segment (x1,y1)-(x2,y2)
local a = (x1-x2)^2 + (y1-y2)^2 -- square of distance
@ -5,13 +8,13 @@ local function distance_to_segment(x1, y1, x2, y2, x, y)
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)
return sqrt(b)
elseif a + c < b then
-- The closest point of the segment is the extremity 2
return math.sqrt(c)
return 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)
return abs(x1 * (y2-y) + x2 * (y-y1) + x * (y1-y2)) / sqrt(a)
end
end
@ -19,8 +22,8 @@ local function transform_quadri(X, Y, x, y)
-- To index points in an irregular quadrilateral, giving x and y between 0 (one edge) and 1 (opposite edge)
-- X, Y 4-vectors giving the coordinates of the 4 vertices
-- x, y position to index.
local x1, x2, x3, x4 = unpack(X)
local y1, y2, y3, y4 = unpack(Y)
local x1, x2, x3, x4 = unpk(X)
local y1, y2, y3, y4 = unpk(Y)
-- Compare distance to 2 opposite edges, they give the X coordinate
local d23 = distance_to_segment(x2,y2,x3,y3,x,y)

View File

@ -1,14 +1,17 @@
local modpath = minetest.get_modpath(minetest.get_current_modname()) .. '/'
local modpath = minetest.get_modpath(minetest.get_current_modname())
local make_polygons = dofile(modpath .. 'polygons.lua')
local transform_quadri = dofile(modpath .. 'geometry.lua')
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 = tonumber(minetest.get_mapgen_setting("water_level"))
local riverbed_slope = tonumber(mapgen_rivers.settings:get("riverbed_slope")) * tonumber(mapgen_rivers.settings:get("blocksize"))
local MAP_BOTTOM = -31000
-- Localize for performance
local floor, min, max = math.floor, math.min, math.max
local unpk = unpack
-- Linear interpolation
local function interp(v00, v01, v11, v10, xf, zf)
local v0 = v01*xf + v00*(1-xf)
@ -31,11 +34,11 @@ local function heightmaps(minp, maxp)
if poly then
local xf, zf = transform_quadri(poly.x, poly.z, x, z)
local i00, i01, i11, i10 = unpack(poly.i)
local i00, i01, i11, i10 = unpk(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)
local r_west, r_north, r_east, r_south = unpk(poly.rivers)
local c_NW, c_NE, c_SE, c_SW = unpk(poly.river_corners)
-- Calculate the depth factor for each edge and corner.
-- Depth factor:
@ -65,10 +68,10 @@ local function heightmaps(minp, maxp)
-- 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)
local x0 = max(r_west, c_NW-zf, zf-c_SW)
local x1 = min(r_east, c_NE+zf, c_SE-zf)
local z0 = max(r_north, c_NW-xf, xf-c_NE)
local z1 = min(r_south, c_SW+xf, c_SE-xf)
xf = (xf-x0) / (x1-x0)
zf = (zf-z0) / (z1-z0)
elseif imax == 1 then
@ -91,7 +94,7 @@ local function heightmaps(minp, maxp)
-- Determine elevation by interpolation
local vdem = poly.dem
local terrain_height = math.floor(0.5+interp(
local terrain_height = floor(0.5+interp(
vdem[1],
vdem[2],
vdem[3],
@ -116,10 +119,10 @@ local function heightmaps(minp, maxp)
lake_id = 1
end
end
local lake_height = math.max(math.floor(poly.lake[lake_id]), terrain_height)
local lake_height = max(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)
terrain_height = min(max(lake_height, sea_level) - floor(1+depth_factor_max*riverbed_slope), terrain_height)
end
terrain_height_map[i] = terrain_height

249
init.lua
View File

@ -1,226 +1,43 @@
mapgen_rivers = {}
local modpath = minetest.get_modpath(minetest.get_current_modname()) .. '/'
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')
local heightmaps = dofile(modpath .. 'heightmap.lua')
-- Linear interpolation
local function interp(v00, v01, v11, v10, xf, zf)
local v0 = v01*xf + v00*(1-xf)
local v1 = v11*xf + v10*(1-xf)
return v1*zf + v0*(1-zf)
if minetest.get_mapgen_setting("mg_name") ~= "singlenode" then
minetest.set_mapgen_setting("mg_name", "singlenode", true)
minetest.log("warning", "[mapgen_rivers] Mapgen set to singlenode")
end
local data = {}
local modpath = minetest.get_modpath(minetest.get_current_modname())
dofile(modpath .. '/settings.lua')
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 sfile = io.open(minetest.get_worldpath() .. '/river_data/size')
if sfile then
sfile:close()
else
dofile(modpath .. '/pregenerate.lua')
collectgarbage()
end
local function generate(minp, maxp, seed)
local chulens = {
x = maxp.x-minp.x+1,
y = maxp.y-minp.y+1,
z = maxp.z-minp.z+1,
}
if not init then
mapsize = {
x = chulens.x,
y = chulens.y+1,
z = chulens.z,
}
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)
mapgen_rivers.use_mapgen_thread = minetest.settings:get_bool("mapgen_rivers_use_mapgen_thread")
mapgen_rivers.thread = 'main'
if mapgen_rivers.use_mapgen_thread then
if minetest.register_mapgen_dofile then
minetest.register_mapgen_dofile(modpath .. '/mapgen.lua')
else
terrain_map, lake_map = heightmaps(minp, maxp)
minetest.log("warning", "[mapgen_rivers] Mapgen thread not available on this Minetest version.")
mapgen_rivers.use_mapgen_thread = false
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 nid = mapsize.x*(mapsize.y-1) + 1
local incrY = -mapsize.x
local incrX = 1 - mapsize.y*incrY
local incrZ = mapsize.x*mapsize.y - mapsize.x*incrX - mapsize.x*mapsize.y*incrY
local i2d = 1
for z = minp.z, maxp.z do
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
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 + 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
if y <= maxp.y then
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
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
ground_above = y <= terrain
ivm = ivm + ystride
if use_distort then
nid = nid + incrY
end
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()
end
minetest.register_on_generated(generate)
if not mapgen_rivers.use_mapgen_thread then
dofile(modpath .. '/mapgen.lua')
end
-- Setup a metatable to load grid on request if not present
local mt = {
__index = function(_, field)
if field == 'grid' then
dofile(modpath .. '/load_grid.lua')
return mapgen_rivers.grid
end
end,
}
setmetatable(mapgen_rivers, mt)

View File

@ -1,31 +0,0 @@
local worldpath = minetest.get_worldpath() .. "/river_data/"
local function load_map(filename, bytes, signed, size)
local file = io.open(worldpath .. filename, 'rb')
local data = file:read('*all')
if #data < bytes*size then
data = minetest.decompress(data)
end
local map = {}
for i=1, size do
local i0, i1 = (i-1)*bytes+1, i*bytes
local elements = {data:byte(i0, i1)}
local n = elements[1]
if signed and n >= 128 then
n = n - 256
end
for j=2, bytes do
n = n*256 + elements[j]
end
map[i] = n
end
file:close()
return map
end
return load_map

107
load_grid.lua Normal file
View File

@ -0,0 +1,107 @@
local datapath = minetest.get_worldpath() .. "/river_data/"
local floor = math.floor
local sbyte = string.byte
local unpk = unpack
local load_map
local use_interactive_loader
if minetest.settings:has("mapgen_rivers_use_interactive_loader") then
use_interactive_loader = minetest.settings:get_bool("mapgen_rivers_use_interactive_loader")
else
use_interactive_loader = not minetest.settings:get_bool("mapgen_rivers_load_all")
end
if use_interactive_loader then
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,
}
load_map = function(filename, bytes, signed, size, converter)
local file = io.open(datapath .. 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
else
load_map = function(filename, bytes, signed, size, converter)
local file = io.open(datapath .. filename, 'rb')
local data = file:read('*all')
if #data < bytes*size then
data = minetest.decompress(data)
end
local map = {}
for i=1, size do
local i0 = (i-1)*bytes+1
local i1 = i*bytes
local elements = {data:byte(i0, i1)}
local n = sbyte(data, i0)
if signed and n >= 128 then
n = n - 256
end
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
end
local sfile = io.open(datapath..'size', 'r')
assert(sfile)
local X, Z = tonumber(sfile:read('*l')), tonumber(sfile:read('*l'))
sfile:close()
local function offset_converter(o)
return (o + 0.5) * (1/256)
end
mapgen_rivers.grid = {
size = {x=X, y=Z},
dem = load_map('dem', 2, true, X*Z),
lakes = load_map('lakes', 2, true, X*Z),
dirs = load_map('dirs', 1, false, X*Z),
rivers = load_map('rivers', 4, false, X*Z),
offset_x = load_map('offset_x', 1, true, X*Z, offset_converter),
offset_y = load_map('offset_y', 1, true, X*Z, offset_converter);
}

307
mapgen.lua Normal file
View File

@ -0,0 +1,307 @@
-- Recreate mod table if we are in a separate environment
if not minetest.global_exists("mapgen_rivers") then
mapgen_rivers = {use_mapgen_thread=true, thread='mapgen'}
mapgen_rivers.settings = Settings(minetest.get_worldpath() .. '/mapgen_rivers.conf')
end
if not mapgen_rivers.grid then
dofile(minetest.get_modpath(minetest.get_current_modname()) .. '/load_grid.lua')
end
local settings = mapgen_rivers.settings
local sea_level = tonumber(minetest.get_mapgen_setting("water_level"))
local elevation_chill = tonumber(settings:get('elevation_chill'))
local use_distort = settings:get_bool('distort')
local use_biomes = settings:get_bool('biomes')
local use_biomegen_mod = use_biomes and minetest.global_exists('biomegen')
use_biomes = use_biomes and minetest.get_modpath("default") and not use_biomegen_mod
local noiseparams = {
distort_x = settings:get_np_group('np_distort_x'),
distort_z = settings:get_np_group('np_distort_z'),
distort_amplitude = settings:get_np_group('np_distort_amplitude'),
}
if use_biomes then
noiseparams.heat = minetest.get_mapgen_setting_noiseparams("mg_biome_np_heat")
noiseparams.heat.offset = noiseparams.heat.offset + sea_level / elevation_chill
noiseparams.heat_blend = minetest.get_mapgen_setting_noiseparams("mg_biome_np_heat_blend")
end
if use_biomegen_mod then
biomegen.set_elevation_chill(elevation_chill)
end
local heightmaps = dofile(minetest.get_modpath(minetest.get_current_modname()) .. '/heightmap.lua')
-- Linear interpolation
local function interp(v00, v01, v11, v10, xf, zf)
local v0 = v01*xf + v00*(1-xf)
local v1 = v11*xf + v10*(1-xf)
return v1*zf + v0*(1-zf)
end
-- Localize for performance
local floor, min = math.floor, math.min
local data = {}
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 sumtime = 0
local sumtime2 = 0
local ngen = 0
local function init_mapgen(chulens)
mapsize = {
x = chulens.x,
y = chulens.y+1,
z = chulens.z,
}
if use_distort then
noise_x_obj = minetest.get_perlin_map(noiseparams.distort_x, mapsize)
noise_z_obj = minetest.get_perlin_map(noiseparams.distort_z, mapsize)
noise_distort_obj = minetest.get_perlin_map(noiseparams.distort_amplitude, chulens)
end
if use_biomes then
noise_heat_obj = minetest.get_perlin_map(noiseparams.heat, chulens)
noise_heat_blend_obj = minetest.get_perlin_map(noiseparams.heat_blend, chulens)
end
end
local function generate(vm, minp, maxp, seed)
minetest.log("info", ("[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,
z = maxp.z-minp.z+1,
}
if not init then
init_mapgen(chulens)
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
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=floor(xmin), z=floor(zmin)}
local pmaxp = {x=floor(xmax)+1, z=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
-- Check that there is at least one position that reaches min y
if minp.y > sea_level then
local y0 = minp.y
local is_empty = true
for i=1, #terrain_map do
if terrain_map[i] >= y0 or lake_map[i] >= y0 then
is_empty = false
break
end
end
-- If not, skip chunk
if is_empty then
local t = os.clock() - t0
ngen = ngen + 1
sumtime = sumtime + t
sumtime2 = sumtime2 + t*t
minetest.log("verbose", "[mapgen_rivers] Skipping empty chunk (fully above ground level)")
minetest.log("verbose", ("[mapgen_rivers] Done in %5.3f s"):format(t))
return
end
end
local c_stone = minetest.get_content_id("mapgen_stone")
local c_water = minetest.get_content_id("mapgen_water_source")
local c_rwater = minetest.get_content_id("mapgen_river_water_source")
local c_dirt, c_lawn, c_dirtsnow, c_snow, c_sand, c_ice
if use_biomes then
c_dirt = minetest.get_content_id("default:dirt")
c_lawn = minetest.get_content_id("default:dirt_with_grass")
c_dirtsnow = minetest.get_content_id("default:dirt_with_snow")
c_snow = minetest.get_content_id("default:snowblock")
c_sand = minetest.get_content_id("default:sand")
c_ice = minetest.get_content_id("default:ice")
end
local emin, emax = vm:get_emerged_area()
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 nid = mapsize.x*(mapsize.y-1) + 1
local incrY = -mapsize.x
local incrX = 1 - mapsize.y*incrY
local incrZ = mapsize.x*mapsize.y - mapsize.x*incrX - mapsize.x*mapsize.y*incrY
local i2d = 1
for z = minp.z, maxp.z do
for x = minp.x, maxp.x do
local ivm = a:index(x, maxp.y+1, 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
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 = floor(xn)
local z0 = floor(zn)
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 = min(lake_map[i0], lake_map[i1], lake_map[i2], lake_map[i3])
end
if y <= maxp.y then
local is_lake = lake > terrain
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
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
ground_above = y <= terrain
ivm = ivm - ystride
if use_distort then
nid = nid + incrY
end
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()
if mapgen_rivers.thread == "main" then
vm:write_to_map()
end
local t = os.clock()-t0
ngen = ngen + 1
sumtime = sumtime + t
sumtime2 = sumtime2 + t*t
minetest.log("verbose", ("[mapgen_rivers] Done in %5.3f s"):format(t))
end
if mapgen_rivers.thread == "main" then
minetest.register_on_generated(function(minp, maxp, seed)
local vm = minetest.get_mapgen_object("voxelmanip")
generate(vm, minp, maxp, seed)
end)
elseif mapgen_rivers.thread == "mapgen" then
minetest.register_on_generated(generate)
end
minetest.register_on_shutdown(function()
if ngen == 0 then
return
end
local avg = sumtime / ngen
local std = math.sqrt(sumtime2/ngen - avg*avg)
minetest.log("action", ("[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)

View File

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

View File

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

View File

@ -1,102 +1,60 @@
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 X = mapgen_rivers.grid.size.x
local Z = mapgen_rivers.grid.size.y
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 settings = mapgen_rivers.settings
local blocksize = tonumber(settings:get('blocksize'))
local min_catchment = tonumber(settings:get('min_catchment'))
local map_offset = {x=0, z=0}
if mapgen_rivers.center then
if settings:get_bool('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
-- Localize for performance
local floor, ceil, min, max, abs = math.floor, math.ceil, math.min, math.max, math.abs
min_catchment = min_catchment / (blocksize*blocksize)
local wpower = settings:get('river_widening_power')
local wfactor = 1/(2*blocksize * min_catchment^wpower)
local function river_width(flow)
flow = math.abs(flow)
flow = abs(flow)
if flow < min_catchment then
return 0
end
return math.min(wfactor * flow ^ wpower, 1)
return 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 elevation_chill = settings:get_bool('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 = settings:get_bool('glaciers')
local glacier_factor = tonumber(settings:get('glacier_factor'))
local init = false
-- On map generation, determine into which polygon every point (in 2D) will fall.
-- Also store polygon-specific data
local function make_polygons(minp, maxp)
print("Generating polygon map")
print(minp.x, maxp.x, minp.z, maxp.z)
local grid = mapgen_rivers.grid
local dem = grid.dem
local lakes = grid.lakes
local dirs = grid.dirs
local rivers = grid.rivers
local offset_x = grid.offset_x
local offset_z = grid.offset_y
if not init then
if glaciers then
@ -109,9 +67,8 @@ local function make_polygons(minp, maxp)
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)
local xpmin, xpmax = max(floor((minp.x+map_offset.x)/blocksize - 0.5), 0), min(ceil((maxp.x+map_offset.x)/blocksize + 0.5), X-2)
local zpmin, zpmax = max(floor((minp.z+map_offset.z)/blocksize - 0.5), 0), min(ceil((maxp.z+map_offset.z)/blocksize + 0.5), Z-2)
-- Iterate over the polygons
for xp = xpmin, xpmax do
@ -133,15 +90,12 @@ local function make_polygons(minp, maxp)
(offset_z[iC]+zp+1) * blocksize - map_offset.z,
(offset_z[iD]+zp+1) * blocksize - map_offset.z,
}
if xp==xpmin and zp==zpmin then
print(xp, zp, poly_x[1], poly_z[1])
end
local polygon = {x=poly_x, z=poly_z, i={iA, iB, iC, iD}}
local bounds = {} -- Will be a list of the intercepts of polygon edges for every Z position (scanline algorithm)
-- 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)
local zmin = max(floor(min(unpack(poly_z)))+1, minp.z)
local zmax = min(floor(max(unpack(poly_z))), maxp.z)
-- And initialize the arrays
for z=zmin, zmax do
bounds[z] = {}
@ -151,14 +105,14 @@ local function make_polygons(minp, maxp)
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))
local lzmin = floor(min(z1, z2))+1
local lzmax = floor(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 z=max(lzmin, minp.z), 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
@ -169,11 +123,11 @@ local function make_polygons(minp, maxp)
-- Now sort the bounds list
local zlist = bounds[z]
table.sort(zlist)
local c = math.floor(#zlist/2)
local c = 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 xmin = max(floor(zlist[l*2-1])+1, minp.x)
local xmax = min(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
@ -195,16 +149,16 @@ local function make_polygons(minp, maxp)
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)
riverA = 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)
riverB = 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)
riverC = 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)
riverD = min(riverD*glacier_factor, 1)
end
end

131
pregenerate.lua Normal file
View File

@ -0,0 +1,131 @@
local modpath = minetest.get_modpath(minetest.get_current_modname())
local EvolutionModel = dofile(modpath .. '/terrainlib_lua/erosion.lua')
local twist = dofile(modpath .. '/terrainlib_lua/twist.lua')
local blocksize = tonumber(mapgen_rivers.settings:get("blocksize"))
local tectonic_speed = tonumber(mapgen_rivers.settings:get("tectonic_speed"))
local np_base = mapgen_rivers.settings:get_np_group("np_base")
np_base.spread.x = np_base.spread.x / blocksize
np_base.spread.y = np_base.spread.y / blocksize
np_base.spread.z = np_base.spread.z / blocksize
local evol_params = {
K = tonumber(mapgen_rivers.settings:get("river_erosion_coef")),
m = tonumber(mapgen_rivers.settings:get("river_erosion_power")),
d = tonumber(mapgen_rivers.settings:get("difusive_erosion")),
}
local time = tonumber(mapgen_rivers.settings:get("evol_time"))
local time_step = tonumber(mapgen_rivers.settings:get("evol_time_step"))
local niter = math.ceil(time/time_step)
time_step = time / niter
-- Setup the model
local size = {
x = tonumber(mapgen_rivers.settings:get("grid_x_size")),
y = tonumber(mapgen_rivers.settings:get("grid_z_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()
-- Run the model
for i=1, niter do
minetest.log("info", "[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
-- Write the results in the world directory
local datapath = minetest.get_worldpath() .. "/river_data/"
minetest.mkdir(datapath)
local function write_map(filename, data, bytes)
local size = #data
local file = io.open(datapath .. filename, 'wb')
local bytelist = {}
for j=1, bytes do
bytelist[j] = 0
end
local unpk = unpack
local schar = string.char
local floor = math.floor
for i=1, size do
local n = floor(data[i])
data[i] = n
for j=bytes, 2, -1 do
bytelist[j] = n % 256
n = floor(n / 256)
end
bytelist[1] = n % 256
file:write(schar(unpk(bytelist)))
end
file:close()
end
write_map('dem', model.dem, 2)
write_map('lakes', model.lakes, 2)
write_map('dirs', model.dirs, 1)
write_map('rivers', model.rivers, 4)
write_map('offset_x', offset_x, 1)
write_map('offset_y', offset_y, 1)
local sfile = io.open(datapath .. 'size', "w")
sfile:write(size.x..'\n'..size.y)
sfile:close()
local use_interactive_loader
if minetest.settings:has("mapgen_rivers_use_interactive_loader") then
use_interactive_loader = minetest.settings:get_bool("mapgen_rivers_use_interactive_loader")
else
use_interactive_loader = not minetest.settings:get_bool("mapgen_rivers_load_all")
end
if not use_interactive_loader then
mapgen_rivers.grid = {
size = size,
dem = model.dem,
lakes = model.lakes,
dirs = model.dirs,
rivers = model.rivers,
offset_x = offset_x,
offset_y = offset_y,
}
end

42
readconfig.py Normal file
View File

@ -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'

View File

@ -1,58 +1,115 @@
local storage = minetest.get_mod_storage()
local settings = minetest.settings
local mtsettings = minetest.settings
local settings = Settings(minetest.get_worldpath() .. '/mapgen_rivers.conf')
local function get_settings(key, dtype, default)
if storage:contains(key) then
if dtype == "string" then
return storage:get_string(key)
elseif dtype == "int" then
return storage:get_int(key)
elseif dtype == "float" then
return storage:get_float(key)
elseif dtype == "bool" then
return storage:get_string(key) == 'true'
end
mapgen_rivers.version = "1.0.2"
mapgen_rivers.settings = settings
local previous_version_mt = mtsettings:get("mapgen_rivers_version") or "0.0"
local previous_version_mgr = settings:get("version") or "0.0"
if mapgen_rivers.version ~= previous_version_mt or mapgen_rivers.version ~= previous_version_mgr then
local compat_mt, compat_mgr = dofile(minetest.get_modpath(minetest.get_current_modname()) .. "/compatibility.lua")
if mapgen_rivers.version ~= previous_version_mt then
compat_mt(mtsettings)
end
local conf_val = settings:get('mapgen_rivers_' .. key)
if conf_val then
if dtype == "int" then
conf_val = tonumber(conf_val)
storage:set_int(key, conf_val)
elseif dtype == "float" then
conf_val = tonumber(conf_val)
storage:set_float(key, conf_val)
else
storage:set_string(key, conf_val)
if dtype == "bool" then
conf_val = conf_val == 'true'
end
end
return conf_val
else
if dtype == "int" then
storage:set_int(key, default)
elseif dtype == "float" then
storage:set_float(key, default)
elseif dtype == "string" then
storage:set_string(key, default)
elseif dtype == "bool" then
storage:set_string(key, tostring(default))
end
return default
if mapgen_rivers.version ~= previous_version_mgr then
compat_mgr(settings)
end
end
mapgen_rivers.center = get_settings('center', 'bool', false)
mapgen_rivers.blocksize = get_settings('blocksize', 'int', 12)
mapgen_rivers.sea_level = get_settings('sea_level', 'int', 1)
mapgen_rivers.min_catchment = get_settings('min_catchment', 'float', 25)
mapgen_rivers.max_catchment = get_settings('max_catchment', 'float', 40000)
mapgen_rivers.riverbed_slope = get_settings('riverbed_slope', 'float', 0.4) * mapgen_rivers.blocksize
mapgen_rivers.distort = get_settings('distort', 'bool', true)
mapgen_rivers.biomes = get_settings('biomes', 'bool', true)
mapgen_rivers.glaciers = get_settings('glaciers', 'bool', false)
mapgen_rivers.glacier_factor = get_settings('glacier_factor', 'float', 8)
mapgen_rivers.elevation_chill = get_settings('elevation_chill', 'float', 0.25)
mtsettings:set("mapgen_rivers_version", mapgen_rivers.version)
settings:set("version", mapgen_rivers.version)
local modified = false
local function verify_setting(name, dtype, default)
if settings:has(name) then
return
end
modified = true
local v = default
local mtname = 'mapgen_rivers_' .. name
local mthas = mtsettings:has(mtname)
if dtype == "number" or dtype == "string" then
if mthas then
v = mtsettings:get(mtname)
if dtype == "number" and tonumber(v) == nil then
v = default
end
end
settings:set(name, v)
elseif dtype == "bool" then
if mthas then
v = mtsettings:get(mtname)
end
settings:set_bool(name, v)
elseif dtype == "noise" then
if mthas then
v = mtsettings:get_np_group(mtname)
end
settings:set_np_group(name, v)
end
end
verify_setting('center', 'bool', true)
verify_setting('blocksize', 'number', 15)
verify_setting('min_catchment', 'number', 3600)
verify_setting('river_widening_power', 'number', 0.5)
verify_setting('riverbed_slope', 'number', 0.4)
verify_setting('distort', 'bool', true)
verify_setting('biomes', 'bool', true)
verify_setting('glaciers', 'bool', false)
verify_setting('glacier_factor', 'number', 8)
verify_setting('elevation_chill', 'number', 0.25)
verify_setting('grid_x_size', 'number', 1000)
verify_setting('grid_z_size', 'number', 1000)
verify_setting('river_erosion_coef', 'number', 0.5)
verify_setting('river_erosion_power', 'number', 0.4)
verify_setting('diffusive_erosion', 'number', 0.5)
verify_setting('compensation_radius', 'number', 50)
verify_setting('tectonic_speed', 'number', 70)
verify_setting('evol_time', 'number', 10)
verify_setting('evol_time_step', 'number', 1)
verify_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",
})
verify_setting('np_distort_x', 'noise', {
offset = 0,
scale = 1,
seed = -4574,
spread = {x=64, y=32, z=64},
octaves = 3,
persistence = 0.75,
lacunarity = 2,
})
verify_setting('np_distort_z', 'noise', {
offset = 0,
scale = 1,
seed = -7940,
spread = {x=64, y=32, z=64},
octaves = 3,
persistence = 0.75,
lacunarity = 2,
})
verify_setting('np_distort_amplitude', 'noise', {
offset = 0,
scale = 10,
seed = 676,
spread = {x=1024, y=1024, z=1024},
octaves = 5,
persistence = 0.5,
lacunarity = 2,
flags = "absvalue",
})
if modified then
settings:write()
end

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

@ -1,278 +0,0 @@
import numpy as np
import numpy.random as npr
from collections import defaultdict
# This file provide functions to construct the river tree from an elevation model.
# Based on a research paper:
# | Cordonnier, G., Bovy, B., and Braun, J.:
# | A versatile, linear complexity algorithm for flow routing in topographies with depressions,
# | Earth Surf. Dynam., 7, 549562, 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

View File

@ -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)

View File

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

240
terrainlib_lua/erosion.lua Normal file
View File

@ -0,0 +1,240 @@
-- erosion.lua
-- This is the main file of terrainlib_lua. It registers the EvolutionModel object and some of the
local function erode(model, time)
-- Apply river erosion on the model
-- Erosion model is based on the simplified version of the stream-power law Ey = K×A^m×S
-- where Ey is the vertical erosion speed, A catchment area of the river, S slope along the river, m and K local constants.
-- It is equivalent to considering a horizontal erosion wave travelling at Ex = K×A^m, and this latter approach allows much greather time steps so it is used here.
-- For each point, instead of moving upstream and see what point the erosion wave would reach, we move downstream and see from which point the erosion wave would reach the given point, then we can set the elevation.
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) -- Inverse of erosion speed (Ex); time needed for the erosion wave to move through the river section.
erosion_time[i] = etime
lakes[i] = mmax(lakes[i], dem[i], sea_level) -- Use lake/sea surface if higher than ground level, because rivers can not erode below.
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
-- Explore downstream until we find the point 'iw' from which the erosion wave will reach 'i'
local inext = iw
local d = dirs[iw]
-- Follow the river downstream (move 'iw')
if d == 0 then -- If no flow direction, we reach the border of the map: set elevation to the latest node's elev and abort.
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 -- We have found the node from which the erosion wave will take 'time' to arrive to 'i'.
local c = remaining / etime
new_elev = (1-c) * lakes[iw] + c * lakes[inext] -- Interpolate linearly between the two nodes
break
end
remaining = remaining - etime -- If we still don't reach the target time, decrement time and move to next point.
iw = inext
end
dem[i] = mmin(dem[i], new_elev)
end
end
local function diffuse(model, time)
-- Apply diffusion using finite differences methods
-- Adapted for small radiuses
local mmax = math.max
local dem = model.dem
local X, Y = dem.X, dem.Y
local d = model.params.d
-- 'd' is equal to 4 times the diffusion coefficient
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
-- diff should never exceed 1 per iteration.
-- If needed, we will divide the process in enough iterations so that 'ddiff' is below 1.
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
-- Laplacian Δdem × 1/4
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
end
local modpath = ""
if minetest then
modpath = minetest.get_modpath(minetest.get_current_modname()) .. "/terrainlib_lua/"
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)
-- Raises the terrain according to uplift rate (model.params.uplift)
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)
-- Adds noise to the terrain according to noise depth (model.params.noise)
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
-- Isostasy
-- This is the geological phenomenon that makes the lithosphere "float" over the underlying layers.
-- One of the key implications is that when a very large mass is removed from the ground, the lithosphere reacts by moving upward. This compensation only occurs at large scale (as the lithosphere is not flexible enough for small scale adjustments) so the implementation is using a very large-window Gaussian blur of the elevation array.
-- This implementation is quite simplistic, it does not do a mass balance of the lithosphere as this would introduce too many parameters. Instead, it defines a reference equilibrium elevation, and the ground will react toward this elevation (at the scale of the gaussian window).
-- A change in reference isostasy during the run can also be used to simulate tectonic forcing, like making a new mountain range appear.
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
-- Apply isostasy
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] -- Compute the difference between the ground level and the target level
end
-- Blur the difference map using Gaussian blur
gaussian.gaussian_blur_approx(temp, model.params.compensation_radius, 4)
for i=1, X*Y do
dem[i] = dem[i] + temp[i] -- Apply the difference
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

View File

@ -0,0 +1,88 @@
-- 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_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,
}

View File

@ -0,0 +1,442 @@
-- 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, 549562, 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 Borůvka algorithm.
-- Only flow_local and accumulate_flow are custom algorithms.
local function flow_local_semirandom(plist)
-- Determines how water should flow at 1 node scale.
-- The straightforward approach would be "Water will flow to the lowest of the 4 neighbours", but here water flows to one of the lower neighbours, chosen randomly, but probability depends on height difference.
-- This makes rivers better follow the curvature of the topography at large scale, and be less biased by pure N/E/S/W directions.
-- 'plist': array of downward height differences (0 if upward)
local sum = 0
for i=1, #plist do
sum = sum + plist[i] -- Sum of probabilities
end
if sum == 0 then
return 0
end
local r = math.random() * sum
for i=1, #plist do
local p = plist[i]
if r < p then
return i
end
r = r - p
end
return 0
end
-- Maybe implement more flow methods in the future?
local flow_methods = {
semirandom = flow_local_semirandom,
}
-- Applies all steps of the flow routing, to calculate flow direction for every node, and lake surface elevation.
-- It's quite a hard piece of code, but we will go step by step and explain what's going on, so stay with me and... let's goooooooo!
local function flow_routing(dem, dirs, lakes, method) -- 'dirs' and 'lakes' are optional tables to reuse for memory optimization, they may contain any data.
method = method or 'semirandom'
local flow_local = flow_methods[method] or flow_local_semirandom
dirs = dirs or {}
lakes = lakes or {}
-- Localize for performance
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
----------------------------------------
-- STEP 1: Find local flow directions --
----------------------------------------
-- Use the local flow function and fill the flow direction tables
local singular = {}
for y=1, Y do
for x=1, X do
local zi = dem[i]
local plist = { -- Get the height difference of the 4 neighbours (and 0 if uphill)
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': Direction toward which water flow
-- 'dirs2': Directions from which water comes
dirs[i] = d
if d == 0 then -- If water can't flow from this node, add it to the list of singular nodes that will be resolved later
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
--------------------------------------
-- STEP 2: Compute basins and links --
--------------------------------------
-- Now water can flow until it reaches a singular node (which is in most cases the bottom of a depression)
-- We will calculate the drainage basin of every singular node (all the nodes from which the water will flow in this singular node, directly or indirectly), make an adjacency list of basins, and find the lowest pass between each pair of adjacent basins (they are potential lake outlets)
local nbasins = #singular
local basin_id = {}
local links = {}
local basin_links
-- Function to analyse a link between two nodes
local function add_link(i1, i2, b1, isY)
-- i1, i2: coordinates of two nodes
-- b1: basin that contains i1
-- isY: whether the link is in Y direction
local b2
-- Note that basin number #0 represents the outside of the map; or if the coordinate is inside the map, means that the basin number is uninitialized.
if i2 == 0 then -- If outside the map
b2 = 0
else
b2 = basin_id[i2]
if b2 == 0 then -- If basin of i2 is not already computed, skip
return
end
end
if b2 ~= b1 then -- If these two nodes don't belong to the same basin, we have found a link between two adjacent basins
local elev = i2 == 0 and dem[i1] or mmax(dem[i1], dem[i2]) -- Elevation of the highest of the two sides of the link (or only i1 if b2 is map outside)
local l2 = basin_links[b2]
if not l2 then
l2 = {}
basin_links[b2] = l2
end
if not l2.elev or l2.elev > elev then -- If this link is lower than the lowest registered link between these two basins, register it as the new lowest pass
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=1, nbasins do
-- Here we will recursively search upstream from the singular node to determine its drainage basin
local queue = {singular[ib]} -- Start with the singular node, then this queue will be filled with water donors neighbours
basin_links = {}
links[#links+1] = basin_links
while #queue > 0 do
local i = tremove(queue)
basin_id[i] = ib
local d = dirs2[i] -- Get the directions water is coming from
-- Iterate through the 4 directions
if d >= 8 then -- River coming from the East
d = d - 8
queue[#queue+1] = i+1
-- If no river is coming from the East, we might be at the limit of two basins, thus we need to test adjacency.
elseif i%X > 0 then
add_link(i, i+1, ib, false)
else -- If the eastern neighbour is outside the map
add_link(i, 0, ib, false)
end
if d >= 4 then -- River coming from the South
d = d - 4
queue[#queue+1] = i+X
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 the West
d = d - 2
queue[#queue+1] = i-1
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 the North
queue[#queue+1] = i-X
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
-- Iterate through pairs of adjacent basins, and make the links reciprocal
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
-----------------------------------------------------
-- STEP 3: Compute minimal spanning tree of basins --
-----------------------------------------------------
-- We've got an adjacency list of basins with the elevation of their links.
-- We will build a minimal spanning tree of the basins (where costs are the elevation of the links). As demonstrated by Cordonnier et al., this finds the outlets of the basins, where water would naturally flow. This does not tell in which direction water is flowing, however.
-- We will use a version of Borůvka's algorithm, with Mareš' optimizations to approach linear complexity (see paper).
-- The concept of Borůvka's algorithm is to take elements and merge them with their lowest neighbour, until all elements are merged.
-- Mareš' optimizations mainly consist in skipping elements that have over 8 links, until extra links are removed when other elements are merged.
-- Note that for this step we are only working on basins, not grid nodes.
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
-- Iterate in lowlevel but its contents may change during the loop
-- 'next' called with only one argument always returns an element if table is not empty
local b1, lnk1 = next(lowlevel)
lowlevel[b1] = nil
local b2
local lowest = math.huge
local lnk1 = links[b1]
local i = 0
-- Look for lowest link
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, in both directions
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
-- When the number of links is changing, we need to check whether the basin can be added to / removed from 'lowlevel'
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 -- If bassin bn is also linked to b2
nlinks[bn] = nlinks[bn] - 1 -- Then bassin bn is losing a link because it keeps only one link toward b1/b2 after the merge
if nlinks[bn] == 8 then
lowlevel[bn] = lnkn
end
else -- If bn was linked to b1 but not to b2
nlinks[b2] = nlinks[b2] + 1 -- Then b2 is gaining a link to bn because of the merge
if nlinks[b2] == 9 then
lowlevel[b2] = nil
end
end
if not lnkn[b2] or lnkn[b2].elev > bdata.elev then -- If the link b1-bn will become the new lowest link between b2 and bn, redirect the link to b2
lnkn[b2] = bdata
lnk2[bn] = bdata
end
end
end
--------------------------------------------------------------
-- STEP 4: Orient basin graph, and grid nodes inside basins --
--------------------------------------------------------------
-- We will finally solve those freaking singular nodes.
-- To orient the basin graph, we will consider that the ultimate basin water should flow into is the map outside (basin #0). We will start from it and recursively walk upstream to the neighbouring basins, using only links that are in the minimal spanning tree. This gives the flow direction of the links, and thus, the outlet of every basin.
-- This will also give lake elevation, which is the highest link encountered between map outside and the given basin on the spanning tree.
-- And within each basin, we need to modify flow directions to connect the singular node to the outlet.
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) -- Pop from queue
queue[b1] = nil
basin_lake[b1] = elev1
-- Iterate through b1's neighbours (according to the spanning tree)
for b2, bound in pairs(basin_graph[b1]) do
-- Make b2 flow into b1
local i = bound.i -- Get the coordinate of the link (which is the basin's outlet)
local dir = bound.is_y and 3 or 4 -- And get the direction (S/E/N/W)
if basin_id[i] ~= b2 then
dir = dir - 2
-- Coordinate 'i' refers to the side of the link with the highest X/Y position. In case it is in the wrong basin, take the other side by decrementing by one row/column.
if bound.is_y then
i = i - X
else
i = i - 1
end
elseif b1 == 0 then
dir = 0
end
-- Use the flow directions computed in STEP 2 to find the route from the outlet position to the singular node, and reverse this route to make the singular node flow into the outlet
-- This can make the river flow uphill, which may seem unnatural, but it can only happen below a lake (because outlet elevation defines lake surface elevation)
repeat
-- Assign i's direction to 'dir', and get i's former direction
dir, dirs[i] = dirs[i], dir
-- Move i by following its former flow direction (downstream)
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
-- Reverse the flow direction for the next node, which will flow into i
dir = reverse[dir]
until dir == 0 -- Stop when reaching the singular node
-- Add basin b2 into the queue, and keep the highest link elevation, that will define the elevation of the lake in b2
queue[b2] = mmax(elev1, bound.elev)
-- Remove b1 from b2's neighbours to avoid coming back to b1
basin_graph[b2][b1] = nil
end
basin_graph[b1] = nil
end
-- Every node will be assigned the lake elevation of the basin it belongs to.
-- If lake elevation is lower than ground elevation, it simply means that there is no lake here.
for i=1, X*Y do
lakes[i] = basin_lake[basin_id[i]]
end
-- That's it!
return dirs, lakes
end
local function accumulate(dirs, waterq)
-- Calculates the river flow by determining the surface of the catchment area for every node
-- This means: how many nodes will give their water to that given node, directly or indirectly?
-- This is obtained by following rivers downstream and summing up the flow of every tributary, starting with a value of 1 at the sources.
-- This will give non-zero values for every node but only large values will be considered to be rivers.
waterq = waterq or {}
local X, Y = dirs.X, dirs.Y
local ndonors = {}
local waterq = {X=X, Y=Y}
for i=1, X*Y do
ndonors[i] = 0
waterq[i] = 1
end
-- Calculate the number of direct donors
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
-- Find sources (nodes that have no donor)
if ndonors[i1] == 0 then
local i2 = i1
local dir = dirs[i2]
local w = waterq[i2]
-- Follow the water flow downstream: move 'i2' to the next node according to its flow direction
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
-- Increment the water quantity of i2
w = w + waterq[i2]
waterq[i2] = w
-- Stop on an unresolved confluence (node with >1 donors) and decrease the number of remaining donors
-- When the ndonors of a confluence has decreased to 1, it means that its water quantity has already been incremented by its tributaries, so it can be resolved like a standard river section. However, do not decrease ndonors to zero to avoid considering it as a source.
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,
}

102
terrainlib_lua/twist.lua Normal file
View File

@ -0,0 +1,102 @@
-- twist.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 File

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

View File

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