2 Commits

Author SHA1 Message Date
0a1c08648d Use biomegen.generate_all 2020-11-11 13:58:23 +01:00
290b998735 Added support for biomegen mod 2020-11-10 14:15:31 +01:00
29 changed files with 337 additions and 1025 deletions

2
.gitignore vendored
View File

@ -9,4 +9,4 @@ bounds_y
dirs
rivers
unused/
river_data/
data/

165
LICENSE
View File

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

122
README.md
View File

@ -1,114 +1,36 @@
# Map Generator with Rivers
`mapgen_rivers v0.0` by Gaël de Sailly.
mapgen_rivers
=============
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).
Procedural map generator for Minetest 5.x. Focused on river networks, and features valley erosion and lakes.
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.
Contains two distinct programs: Python scripts for pre-processing, and Lua scripts to generate the map on Minetest.
![Screenshot](https://user-images.githubusercontent.com/6905002/98825953-6289d980-2435-11eb-9e0b-704a95663ce0.png)
**Important to know**: Unlike most other Minetest mods, it does not contain standalone Lua code, but does part of its processing with a separate Python program (included).
- The Python part does pre-processing: it creates large-scale terrain data and applies landscape evolution algorithms, then outputs a grid of data in the mod's or world's folder. The grid is typically an array of 1000x1000 points of data, each of them representing a cell (by default 12x12 nodes). This pre-processing is long and should be run in advance.
- The Lua part does actual map generation on Minetest. It reads grid data, upscales it (by a factor 12 by default), and adds small-scale features.
# Requirements
Mod dependencies: `default` required, and [`biomegen`](https://github.com/Gael-de-Sailly/biomegen) optional.
Map pre-generation requires Python 3 with the following libraries installed:
- `numpy`, widely used library for numerical calculations
- `scipy`, a library for advanced data treatments, that is used here for Gaussian filtering
- `noise`, implementing Perlin/Simplex noises
Also, the following are optional (for map preview)
- `matplotlib`, a famous library for graphical plotting
- `colorcet` if you absolutely need better colormaps for preview :-)
They are all commonly found on `pip` or `conda` Python distributions.
![Screenshot](https://user-images.githubusercontent.com/6905002/79541028-687b3000-8089-11ea-9209-c23c15d75383.png)
# Installation
This mod should be placed in the `mods/` directory of Minetest like any other mod.
This mod should be placed in the `/mods` directory like any other Minetest mod.
The Python part relies on external libraries that you need to install:
- `numpy` and `scipy`, widely used libraries for numerical calculations
- `noise`, doing Perlin/Simplex noises
- optionally, `matplotlib` (for map preview)
They are commonly found on `pip` or `conda` Python distributions.
# 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.
By default, the mod contains a demo 400x400 map (so you can start the game directly), but it is recommended that you run the pre-processing script to generate a new map before world creation, if you can.
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.
## Pre-processing
Run the script `terrain_rivers.py` via command line. You can optionally append the map size (by default 400). Example for a 1000x1000 map:
```
./generate.py
./terrain_rivers.py 1000
```
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.
For a default 400x400 map, it should take between 1 and 2 minutes. It will generate 5 files directly in the mod folder, containing the map data.
## Parameters for `generate.py`
For a basic use you do not need to append any argument:
If you have `matplotlib` installed, the script `view_map.py` can be used to get a map preview. Example:
```
./generate.py
```
By default this will produce a 1000x1000 grid and save it in `river_data/`. Expect a computing time of about 30 minutes.
### Parameters and config files
This pre-processing takes many parameters. Instead of asking all these parameters to the end user, they are grouped in `.conf` files for usability, but the script still allows to override individual settings.
Generic usage:
```
./generate.py conf_file output_dir
./view_map.py data/
```
- `conf_file`: Path to configuration file from which parameters should be read. If omitted, attempts to read in `terrain.conf`.
- `output_dir`: Directory in which to save the grid data, defaults to `river_data/`. If it does not exist, it is created. If it already contains previous grid data, they are overwritten.
#### Config files
The mod currently includes 3 config files, providing different terrain styles:
- `terrain_default.conf` generates the standard terrain, with highest elevations around 250 with sharp peaks, and otherwise hilly terrain.
- `terrain_higher.conf` generates higher mountains (up to 400 nodes), and wider valleys.
- `terrain_original.conf` provides a terrain similar to what was generated with the first release of `mapgen_rivers`.
More work is needed to find better and more varied terrain styles.
### Complete list of parameters
Other parameters can be specified by `--parameter value`. Syntax `--parameter=value` is also supported.
| Parameter | Description | Example |
|---------------|-------------|---------|
| | **Generic parameters** |
| `mapsize` | Size of the grid, in number of cells per edge. Usually `1000`, so to have 1000x1000 cells, the grid will have 1001x1001 nodes. Note that the grid is upscaled 12x in the game (this ratio can be changed), so that a `mapsize` of 1000 will result in a 12000x12000 map by default. | `--mapsize 1000` |
| `sea_level` | Height of the sea; height below which a point is considered under water even if it is not in a closed depression. | `--sea_level 1` |
| | **Noise parameters** |
| `scale` | Horizontal variation wavlength of the largest noise octave, in grid cells (equivalent to the `spread` of a `PerlinNoise`). | `--scale 400` |
| `vscale` | Elevation coefficient, determines the approximate height difference between deepest seas and highest mountains. | `--vscale 300` |
| `offset` | Offset of the noise, will determine mean elevation. | `--offset 0` |
| `persistence` | Relative height of smaller noise octaves compared to bigger ones. | `--persistence 0.6` |
| `lacunarity` | Relative reduction of wavelength between octaves. If `lacunarity`×`persistence` is larger than 1 (usual case), smaller octaves result in higher slopes than larger ones. This case is interesting for rivers networks because slopes determine rivers position. | `--lacunarity 2` |
| | **Landscape evolution parameters**|
| `K` | Abstract erosion constant. Increasing it will increase erosive intensity. | `--K 1` |
| `m` | Parameter representing the influence of river flux on erosion. For `m=0`, small and big rivers are equal contributors to erosion. For `m=1` the erosive capability is proportional to river flux (assumed to be catchment area). Usual values: `0.25`-`0.60`. Be careful, this parameter is *highly sensitive*. | `--m 0.35` |
| `d` | Diffusion coefficient acting on sea/lake floor. Usual values `0`-`1`. | `--d 0.2` |
| `flex_radius` | Flexure radius. Wavelength over which loss/gain of mass is compensated by uplift/subsidence. This ensures that mountain ranges will not get eventually flattened by erosion, and that an equilibrium is reached. Geologically speaking, this implements [isostatic rebound](https://en.wikipedia.org/wiki/Isostasy). | `--flex_radius 20` |
| `time` | Simulated time of erosion modelling, in abstract units. | `--time 10` |
| `niter` | Number of iterations. Each iteration represents a time `time/niter`. | `--niter 10` |
| `sea_level_variations` | Amplitude of sea level variations throughout the simulation (if any). | `--sea_level_variations 10` |
| `sea_level_variations_time` | Characteristic time of variation for sea level, in the same units than `time`. Increasing it will result in slower variations between iterations. | `--sea_level_variations_time 1` |
| `flow_method` | Algorithm used for local flow calculation. Possible values are `steepest` (every node flows toward the steepest neighbour when possible), and `semirandom` (default, flow direction is determined randomly between lower neighbours, with lowest ones having greater probability). | `--flow_method semirandom` |
| | **Alternatives** |
| `config` | Another way to specify configuration file | `--config terrain_higher.conf` |
| `output` | Another way to specify output dir | `--output ~/.minetest/worlds/my_world/river_data` |
### Example
```
./generate.py terrain_higher.conf --mapsize 700 --K 0.4 --m 0.5
```
Reads parameters in `terrain_higher.conf`, and will generate a 700x700 grid using custom values for `K` and `m`.
## Map preview
If you have `matplotlib` installed, `generate.py` will automatically show the grid aspect in real time during the erosion simulation.
There is also a script to view a generated map afterwards: `view_map.py`. Its syntax is the following:
```
./view_map.py grid blocksize
```
- `grid` is the path to the grid directory to view. For example `river_data/`.
- `blocksize` is the size at which 1 grid cell will be upscaled, in order to match game coordinates. If you use default settings, use `12`.
Example:
```
./view_map.py river_data 12
```
## Map generation
Just create a Minetest world with `singlenode` mapgen, enable this mod and start the world. The data files are immediately copied in the world folder so you can re-generate them afterwards, it won't affect the old worlds.

View File

@ -6,8 +6,8 @@ def make_bounds(dirs, rivers):
"""
(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_h = np.zeros((Y, X-1), dtype='i4')
bounds_v = np.zeros((Y-1, X), dtype='i4')
bounds_v += (rivers * (dirs==1))[:-1,:]
bounds_h += (rivers * (dirs==2))[:,:-1]

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,6 +1,6 @@
import numpy as np
import scipy.ndimage as im
from .rivermapper import flow
import rivermapper as rm
def advection(dem, dirs, rivers, time, K=1, m=0.5, sea_level=0):
"""
@ -11,6 +11,12 @@ def advection(dem, dirs, rivers, time, K=1, m=0.5, sea_level=0):
v = K * flux^m
"""
dirs = dirs.copy()
dirs[0,:] = 0
dirs[-1,:] = 0
dirs[:,0] = 0
dirs[:,-1] = 0
adv_time = 1 / (K*rivers**m) # For every pixel, calculate the time an "erosion wave" will need to cross it.
dem = np.maximum(dem, sea_level)
dem_new = np.zeros(dem.shape)
@ -45,16 +51,14 @@ def advection(dem, dirs, rivers, time, K=1, m=0.5, sea_level=0):
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
return np.minimum(dem, 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'):
def __init__(self, dem, K=1, m=0.5, d=1, sea_level=0, flow=False, flex_radius=100):
self.dem = dem
#self.bedrock = dem
self.K = K
@ -63,22 +67,20 @@ class EvolutionModel:
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.dirs = np.zeros(dem.shape, dtype='u1')
self.rivers = np.zeros(dem.shape, dtype='u4')
self.flow_uptodate = False
def calculate_flow(self):
self.dirs, self.lakes, self.rivers = flow(self.dem, method=self.flow_method)
self.dirs, self.lakes, self.rivers = rm.flow(self.dem)
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)
dem = advection(self.lakes, self.dirs, self.rivers, time, K=self.K, m=self.m, sea_level=self.sea_level)
self.dem = np.minimum(dem, self.dem)
self.flow_uptodate = False

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

@ -30,7 +30,7 @@ local function heightmaps(minp, maxp)
local poly = polygons[i]
if poly then
local xf, zf = transform_quadri(poly.x, poly.z, x, z)
local xf, zf = transform_quadri(poly.x, poly.z, x/blocksize, z/blocksize)
local i00, i01, i11, i10 = unpack(poly.i)
-- Load river width on 4 edges and corners
@ -99,25 +99,7 @@ local function heightmaps(minp, maxp)
xf, zf
))
-- Spatial gradient of the interpolation
local slope_x = zf*(vdem[3]-vdem[4]) + (1-zf)*(vdem[2]-vdem[1]) < 0
local slope_z = xf*(vdem[3]-vdem[2]) + (1-xf)*(vdem[4]-vdem[1]) < 0
local lake_id = 0
if slope_x then
if slope_z then
lake_id = 3
else
lake_id = 2
end
else
if slope_z then
lake_id = 4
else
lake_id = 1
end
end
local lake_height = math.max(math.floor(poly.lake[lake_id]), terrain_height)
local lake_height = math.max(math.floor(poly.lake), terrain_height)
if imax > 0 and depth_factor_max > 0 then
terrain_height = math.min(math.max(lake_height, sea_level) - math.floor(1+depth_factor_max*riverbed_slope), terrain_height)
end

View File

@ -214,7 +214,7 @@ local function generate(minp, maxp, seed)
biomegen.generate_all(data, a, vm, minp, maxp, seed)
else
vm:set_data(data)
minetest.generate_ores(vm, minp, maxp)
mietest.generate_ores(vm, minp, maxp)
end
vm:set_lighting({day = 0, night = 0})

View File

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

View File

@ -1,5 +1,5 @@
local modpath = minetest.get_modpath(minetest.get_current_modname()) .. '/'
local mod_data_path = modpath .. 'river_data/'
local mod_data_path = modpath .. 'data/'
if not io.open(mod_data_path .. 'size', 'r') then
mod_data_path = modpath .. 'demo_data/'
end
@ -11,21 +11,21 @@ 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')
local wfile = io.open(wfilename, 'r')
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')
local mfile = io.open(mfilename, 'r')
local wfile = io.open(wfilename, 'w')
wfile:write(mfile:read("*all"))
mfile:close()
wfile:close()
end
copy_if_needed('size')
local sfile = io.open(world_data_path..'size', 'r')
local sfile = io.open(world_data_path..'size')
local X = tonumber(sfile:read('*l'))
local Z = tonumber(sfile:read('*l'))
sfile:close()
@ -60,12 +60,6 @@ local blocksize = mapgen_rivers.blocksize
local min_catchment = mapgen_rivers.min_catchment
local max_catchment = mapgen_rivers.max_catchment
local map_offset = {x=0, z=0}
if mapgen_rivers.center then
map_offset.x = blocksize*X/2
map_offset.z = blocksize*Z/2
end
-- Width coefficients: coefficients solving
-- wfactor * min_catchment ^ wpower = 1/(2*blocksize)
-- wfactor * max_catchment ^ wpower = 1
@ -95,9 +89,6 @@ local init = false
-- On map generation, determine into which polygon every point (in 2D) will fall.
-- Also store polygon-specific data
local function make_polygons(minp, maxp)
print("Generating polygon map")
print(minp.x, maxp.x, minp.z, maxp.z)
if not init then
if glaciers then
noise_heat = minetest.get_perlin(mapgen_rivers.noise_params.heat)
@ -109,9 +100,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 = math.max(math.floor(minp.x/blocksize - 0.5), 0), math.min(math.ceil(maxp.x/blocksize), X-2)
local zpmin, zpmax = math.max(math.floor(minp.z/blocksize - 0.5), 0), math.min(math.ceil(maxp.z/blocksize), Z-2)
-- Iterate over the polygons
for xp = xpmin, xpmax do
@ -121,27 +111,14 @@ local function make_polygons(minp, maxp)
local iC = index(xp+1, zp+1)
local iD = index(xp, zp+1)
-- Extract the vertices of the polygon
local poly_x = {
(offset_x[iA]+xp) * blocksize - map_offset.x,
(offset_x[iB]+xp+1) * blocksize - map_offset.x,
(offset_x[iC]+xp+1) * blocksize - map_offset.x,
(offset_x[iD]+xp) * blocksize - map_offset.x,
}
local poly_z = {
(offset_z[iA]+zp) * blocksize - map_offset.z,
(offset_z[iB]+zp) * blocksize - map_offset.z,
(offset_z[iC]+zp+1) * blocksize - map_offset.z,
(offset_z[iD]+zp+1) * blocksize - map_offset.z,
}
if xp==xpmin and zp==zpmin then
print(xp, zp, poly_x[1], poly_z[1])
end
local poly_x = {offset_x[iA]+xp, offset_x[iB]+xp+1, offset_x[iC]+xp+1, offset_x[iD]+xp}
local poly_z = {offset_z[iA]+zp, offset_z[iB]+zp, offset_z[iC]+zp+1, offset_z[iD]+zp+1}
local polygon = {x=poly_x, z=poly_z, i={iA, iB, iC, iD}}
local bounds = {} -- Will be a list of the intercepts of polygon edges for every 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 = math.max(math.floor(blocksize*math.min(unpack(poly_z)))+1, minp.z)
local zmax = math.min(math.floor(blocksize*math.max(unpack(poly_z))), maxp.z)
-- And initialize the arrays
for z=zmin, zmax do
bounds[z] = {}
@ -151,13 +128,13 @@ 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 = math.floor(blocksize*math.min(z1, z2))+1
local lzmax = math.floor(blocksize*math.max(z1, z2))
if lzmin <= lzmax then -- If there is at least one position in it
local x1, x2 = poly_x[i1], poly_x[i2]
-- Calculate coefficient of the equation defining the edge: X=aZ+b
local a = (x1-x2) / (z1-z2)
local b = (x1 - a*z1)
local b = blocksize*(x1 - a*z1)
for z=math.max(lzmin, minp.z), math.min(lzmax, maxp.z) do
-- For every Z position involved, add the intercepted X position in the table
table.insert(bounds[z], a*z+b)
@ -185,7 +162,7 @@ local function make_polygons(minp, maxp)
local poly_dem = {dem[iA], dem[iB], dem[iC], dem[iD]}
polygon.dem = poly_dem
polygon.lake = {lakes[iA], lakes[iB], lakes[iC], lakes[iD]}
polygon.lake = math.min(lakes[iA], lakes[iB], lakes[iC], lakes[iD])
-- Now, rivers.
-- Load river flux values for the 4 corners
@ -194,16 +171,16 @@ local function make_polygons(minp, maxp)
local riverC = river_width(rivers[iC])
local riverD = river_width(rivers[iD])
if glaciers then -- Widen the river
if get_temperature(poly_x[1], poly_dem[1], poly_z[1]) < 0 then
if get_temperature(poly_x[1]*blocksize, poly_dem[1], poly_z[1]*blocksize) < 0 then
riverA = math.min(riverA*glacier_factor, 1)
end
if get_temperature(poly_x[2], poly_dem[2], poly_z[2]) < 0 then
if get_temperature(poly_x[2]*blocksize, poly_dem[2], poly_z[2]*blocksize) < 0 then
riverB = math.min(riverB*glacier_factor, 1)
end
if get_temperature(poly_x[3], poly_dem[3], poly_z[3]) < 0 then
if get_temperature(poly_x[3]*blocksize, poly_dem[3], poly_z[3]*blocksize) < 0 then
riverC = math.min(riverC*glacier_factor, 1)
end
if get_temperature(poly_x[4], poly_dem[4], poly_z[4]) < 0 then
if get_temperature(poly_x[4]*blocksize, poly_dem[4], poly_z[4]*blocksize) < 0 then
riverD = math.min(riverD*glacier_factor, 1)
end
end

115
rivermapper.py Normal file
View File

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

View File

@ -22,11 +22,8 @@ local function get_settings(key, dtype, default)
elseif dtype == "float" then
conf_val = tonumber(conf_val)
storage:set_float(key, conf_val)
else
elseif dtype == "string" or dtype == "bool" then
storage:set_string(key, conf_val)
if dtype == "bool" then
conf_val = conf_val == 'true'
end
end
return conf_val
@ -45,7 +42,6 @@ local function get_settings(key, dtype, default)
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)
@ -53,6 +49,6 @@ 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.glaciers = get_settings('glaciers', 'bool', true)
mapgen_rivers.glacier_factor = get_settings('glacier_factor', 'float', 8)
mapgen_rivers.elevation_chill = get_settings('elevation_chill', 'float', 0.25)

View File

@ -1,11 +1,6 @@
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)

View File

@ -1,52 +0,0 @@
# File containing all settings for 'mapgen_rivers' mod.
# Whether the map should be centered at x=0, z=0.
mapgen_rivers_center (Center map) bool false
# Represents horizontal map scale. Every cell of the grid will be upscaled to
# a square of this size.
# For example if the grid size is 1000x1000 and block size is 12,
# the actual size of the map will be 15000.
mapgen_rivers_blocksize (Block size) float 12.0 2.0 40.0
# Sea level used by mapgen_rivers
mapgen_rivers_sea_level (Sea level) int 1
# Minimal catchment area for a river to be drawn, in grid cells
# (1 cell = blocksize x blocksize).
# Lower value means bigger river density
mapgen_rivers_min_catchment (Minimal catchment area) float 25.0 1.0 1000.0
# Catchment area in grid cells (1 grid cell = blocksize x blocksize)
# at which rivers reach their maximal width of 2*blocksize.
# Higher value means a river needs to receive more tributaries to grow in width.
mapgen_rivers_max_catchment (Maximal catchment area) float 40000.0 1000.0 10000000.0
# Lateral slope of the riverbed.
# Higher value means deeper rivers.
mapgen_rivers_riverbed_slope (Riverbed slope) float 0.4 0.0 2.0
# Enable horizontal distorsion (shearing) of landscape, to break the regularity
# of grid cells and allow overhangs.
# Distorsion uses two 3D noises and thus is intensive in terms of computing time.
mapgen_rivers_distort (Distorsion) bool true
# Enable biome generation.
# If 'biomegen' mod is installed, 'mapgen_rivers' will generate biomes from the
# native biome system. If 'biomegen' is not present, will generate only grass and
# snow.
mapgen_rivers_biomes (Biomes) bool true
# Whether to enable glaciers.
# Glaciers are widened river sections, covered by ice, that are generated in
# very cold areas.
mapgen_rivers_glaciers (Glaciers) bool false
# River channels are widened by this factor if they are a glacier.
mapgen_rivers_glacier_widening_factor (Glacier widening factor) float 8.0 1.0 20.0
# Temperature value decreases by this quantity for every node, vertically.
# This results in mountains being more covered by snow.
mapgen_rivers_elevation_chill (Elevation chill) float 0.25 0.0 5.0
# Noises: to be added. For now they are hardcoded.

View File

@ -5,13 +5,10 @@ offset = 0
persistence = 0.6
lacunarity = 2.0
K = 0.5
m = 0.5
d = 0.5
K = 1
m = 0.35
d = 0.8
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

116
terrain_rivers.py Executable file
View File

@ -0,0 +1,116 @@
#!/usr/bin/env python3
import numpy as np
import noise
from save import save
from erosion import EvolutionModel
import bounds
import os
import sys
import settings
# Always place in this script's parent directory
os.chdir(os.path.dirname(sys.argv[0]))
argc = len(sys.argv)
params = {}
if argc > 1:
if os.path.isfile(sys.argv[1]):
params = settings.read_config_file(sys.argv[1])
else:
mapsize = int(sys.argv[1])
def get_setting(name, default):
if name in params:
return params[name]
return default
mapsize = int(get_setting('mapsize', 400))
scale = float(get_setting('scale', 200.0))
vscale = float(get_setting('vscale', 200.0))
offset = float(get_setting('offset', 0.0))
persistence = float(get_setting('persistence', 0.5))
lacunarity = float(get_setting('lacunarity', 2.0))
K = float(get_setting('K', 1.0))
m = float(get_setting('m', 0.35))
d = float(get_setting('d', 1.0))
sea_level = float(get_setting('sea_level', 0.0))
flex_radius = float(get_setting('flex_radius', 20.0))
time = float(get_setting('time', 10.0))
niter = int(get_setting('niter', 10))
n = np.zeros((mapsize+1, mapsize+1))
# Set noise parameters
params = {
"octaves" : int(np.ceil(np.log2(mapsize)))+1,
"persistence" : persistence,
"lacunarity" : lacunarity,
}
# Determine noise offset randomly
xbase = np.random.randint(65536)
ybase = np.random.randint(65536)
# Generate the noise
for x in range(mapsize+1):
for y in range(mapsize+1):
n[x,y] = noise.snoise2(x/scale + xbase, y/scale + ybase, **params)
nn = n*vscale + offset
# Initialize landscape evolution model
print('Initializing model')
model = EvolutionModel(nn, K=1, m=0.35, d=1, sea_level=0, flex_radius=flex_radius)
dt = time/niter
# Run the model's processes: the order in which the processes are run is arbitrary and could be changed.
print('Initial flow calculation')
model.calculate_flow()
for i in range(niter):
print('Iteration {:d} of {:d}'.format(i+1, niter))
print('Diffusion')
model.diffusion(dt)
print('Advection')
model.advection(dt)
print('Isostatic equilibration')
model.adjust_isostasy()
print('Flow calculation')
model.calculate_flow()
print('Done')
# Twist the grid
bx, by = bounds.make_bounds(model.dirs, model.rivers)
offset_x, offset_y = bounds.twist(bx, by, bounds.get_fixed(model.dirs))
# Convert offset in 8-bits
offset_x = np.clip(np.floor(offset_x * 256), -128, 127)
offset_y = np.clip(np.floor(offset_y * 256), -128, 127)
if not os.path.isdir('data'):
os.mkdir('data')
os.chdir('data')
# Save the files
save(model.dem, 'dem', dtype='>i2')
save(model.lakes, 'lakes', dtype='>i2')
save(offset_x, 'offset_x', dtype='i1')
save(offset_y, 'offset_y', dtype='i1')
save(model.dirs, 'dirs', dtype='u1')
save(model.rivers, 'rivers', dtype='>u4')
with open('size', 'w') as sfile:
sfile.write('{:d}\n{:d}'.format(mapsize+1, mapsize+1))
# Display the map if matplotlib is found
try:
from view_map import view_map
view_map(model.dem, model.lakes, model.rivers)
except:
pass

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

View File

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