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
34 changed files with 654 additions and 1804 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.

View File

@ -1,40 +1,36 @@
# Map Generator with Rivers
`mapgen_rivers v1.0.1` by Gaël de Sailly.
mapgen_rivers
=============
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).
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://content.minetest.net/uploads/fff09f2269.png)
It used to be composed of a Python script doing pre-generation, and a Lua mod reading the pre-generation output and generating the map. The code has been rewritten in full Lua for version 1.0 (July 2021), and is now usable out-of-the-box as any other Minetest mod.
# Author and license
License: GNU LGPLv3.0
Code: Gaël de Sailly
Flow routing algorithm concept (in `terrainlib/rivermapper.lua`): Cordonnier, G., Bovy, B., & Braun, J. (2019). A versatile, linear complexity algorithm for flow routing in topographies with depressions. Earth Surface Dynamics, 7(2), 549-562.
# Requirements
Mod dependencies: `default` required, and [`biomegen`](https://github.com/Gael-de-Sailly/biomegen) optional (provides biome system).
![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
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.
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.
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).
## Settings
Settings can be found in Minetest in the `Settings` tab, `All settings` -> `Mods` -> `mapgen_rivers`.
Most settings are world-specific and a copy is made in `mapgen_rivers.conf` in the world folder, during world first use, which means that further modification of global settings will not alter existing worlds.
## Map preview
The Python script `view_map.py` can display the full map. You need to have Python 3 installed, as well as the libraries `numpy`, `matplotlib`, and optionally `colorcet`. For `conda` users, an `environment.yml` file is provided.
It can be run from command line by passing the world folder. Example:
## 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:
```
./view_map.py ~/.minetest/worlds/test_mg_rivers
./terrain_rivers.py 1000
```
For a default 400x400 map, it should take between 1 and 2 minutes. It will generate 5 files directly in the mod folder, containing the map data.
If you have `matplotlib` installed, the script `view_map.py` can be used to get a map preview. Example:
```
./view_map.py data/
```
## 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.

74
bounds.py Normal file
View File

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

View File

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

BIN
demo_data/dem Normal file

Binary file not shown.

BIN
demo_data/dirs Normal file

Binary file not shown.

BIN
demo_data/lakes Normal file

Binary file not shown.

BIN
demo_data/offset_x Normal file

Binary file not shown.

BIN
demo_data/offset_y Normal file

Binary file not shown.

BIN
demo_data/rivers Normal file

Binary file not shown.

2
demo_data/size Normal file
View File

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

View File

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

97
erosion.py Normal file
View File

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

View File

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

@ -1,20 +1,15 @@
mapgen_rivers = {}
local modpath = minetest.get_modpath(minetest.get_current_modname()) .. '/'
mapgen_rivers.modpath = modpath
mapgen_rivers.world_data_path = minetest.get_worldpath() .. '/river_data/'
if minetest.get_mapgen_setting("mg_name") ~= "singlenode" then
minetest.set_mapgen_setting("mg_name", "singlenode", true)
print("[mapgen_rivers] Mapgen set to singlenode")
end
dofile(modpath .. 'settings.lua')
local sea_level = mapgen_rivers.settings.sea_level
local elevation_chill = mapgen_rivers.settings.elevation_chill
local use_distort = mapgen_rivers.settings.distort
local use_biomes = mapgen_rivers.settings.biomes
local 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
@ -43,13 +38,7 @@ local noise_heat_blend_map = {}
local mapsize
local init = false
local sumtime = 0
local sumtime2 = 0
local ngen = 0
local function generate(minp, maxp, seed)
print(("[mapgen_rivers] Generating from %s to %s"):format(minetest.pos_to_string(minp), minetest.pos_to_string(maxp)))
local chulens = {
x = maxp.x-minp.x+1,
y = maxp.y-minp.y+1,
@ -74,7 +63,6 @@ local function generate(minp, maxp, seed)
init = true
end
local t0 = os.clock()
local minp2d = {x=minp.x, y=minp.z}
if use_distort then
noise_x_obj:get_3d_map_flat(minp, noise_x_map)
@ -142,7 +130,7 @@ local function generate(minp, maxp, seed)
local incrZ = mapsize.x*mapsize.y - mapsize.x*incrX - mapsize.x*mapsize.y*incrY
local i2d = 1
for z = minp.z, maxp.z do
for x = minp.x, maxp.x do
local ivm = a:index(x, minp.y, z)
@ -226,25 +214,13 @@ 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})
vm:calc_lighting()
vm:update_liquids()
vm:write_to_map()
local t1 = os.clock()
local t = t1-t0
ngen = ngen + 1
sumtime = sumtime + t
sumtime2 = sumtime2 + t*t
print(("[mapgen_rivers] Done in %5.3f s"):format(t))
end
minetest.register_on_generated(generate)
minetest.register_on_shutdown(function()
local avg = sumtime / ngen
local std = math.sqrt(sumtime2/ngen - avg*avg)
print(("[mapgen_rivers] Mapgen statistics:\n- Mapgen calls: %4d\n- Mean time: %5.3f s\n- Standard deviation: %5.3f s"):format(ngen, avg, std))
end)

View File

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

View File

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

View File

@ -1,98 +1,70 @@
local modpath = mapgen_rivers.modpath
local mod_data_path = modpath .. 'river_data/'
local modpath = minetest.get_modpath(minetest.get_current_modname()) .. '/'
local mod_data_path = modpath .. 'data/'
if not io.open(mod_data_path .. 'size', 'r') then
mod_data_path = modpath .. 'demo_data/'
end
local world_data_path = mapgen_rivers.world_data_path
local world_data_path = minetest.get_worldpath() .. '/river_data/'
minetest.mkdir(world_data_path)
dofile(modpath .. 'load.lua')
local load_map = dofile(modpath .. 'load.lua')
mapgen_rivers.grid = {}
local X = mapgen_rivers.settings.grid_x_size
local Z = mapgen_rivers.settings.grid_z_size
local function offset_converter(o)
return (o + 0.5) * (1/256)
end
local load_all = mapgen_rivers.settings.load_all
-- Try to read file 'size'
local sfile = io.open(world_data_path..'size', 'r')
local first_mapgen = true
if sfile then
X, Z = tonumber(sfile:read('*l')), tonumber(sfile:read('*l'))
sfile:close()
first_mapgen = false
end
if first_mapgen then
-- Generate a map!!
local pregenerate = dofile(mapgen_rivers.modpath .. '/pregenerate.lua')
minetest.register_on_mods_loaded(function()
print('[mapgen_rivers] Generating grid')
pregenerate(load_all)
if load_all then
local offset_x = mapgen_rivers.grid.offset_x
local offset_y = mapgen_rivers.grid.offset_y
for i=1, X*Z do
offset_x[i] = offset_converter(offset_x[i])
offset_y[i] = offset_converter(offset_y[i])
end
end
end)
end
-- if data not already loaded
if not (first_mapgen and load_all) then
local load_map
if load_all then
load_map = mapgen_rivers.load_map
else
load_map = mapgen_rivers.interactive_loader
local function copy_if_needed(filename)
local wfilename = world_data_path..filename
local wfile = io.open(wfilename, 'r')
if wfile then
wfile:close()
return
end
minetest.register_on_mods_loaded(function()
if load_all then
print('[mapgen_rivers] Loading full grid')
else
print('[mapgen_rivers] Loading grid as interactive loaders')
end
local grid = mapgen_rivers.grid
grid.dem = load_map('dem', 2, true, X*Z)
grid.lakes = load_map('lakes', 2, true, X*Z)
grid.dirs = load_map('dirs', 1, false, X*Z)
grid.rivers = load_map('rivers', 4, false, X*Z)
grid.offset_x = load_map('offset_x', 1, true, X*Z, offset_converter)
grid.offset_y = load_map('offset_y', 1, true, X*Z, offset_converter)
end)
local mfilename = mod_data_path..filename
local mfile = io.open(mfilename, 'r')
local wfile = io.open(wfilename, 'w')
wfile:write(mfile:read("*all"))
mfile:close()
wfile:close()
end
mapgen_rivers.grid.size = {x=X, y=Z}
copy_if_needed('size')
local sfile = io.open(world_data_path..'size')
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.settings.blocksize
local min_catchment = mapgen_rivers.settings.min_catchment
local max_catchment = mapgen_rivers.settings.max_catchment
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.settings.center then
map_offset.x = blocksize*X/2
map_offset.z = blocksize*Z/2
end
local min_catchment = mapgen_rivers.settings.min_catchment / (blocksize*blocksize)
local wpower = mapgen_rivers.settings.river_widening_power
local wfactor = 1/(2*blocksize * min_catchment^wpower)
-- Width coefficients: coefficients solving
-- wfactor * min_catchment ^ wpower = 1/(2*blocksize)
-- wfactor * max_catchment ^ wpower = 1
local wpower = math.log(2*blocksize)/math.log(max_catchment/min_catchment)
local wfactor = 1 / max_catchment ^ wpower
local function river_width(flow)
flow = math.abs(flow)
if flow < min_catchment then
@ -103,30 +75,20 @@ local function river_width(flow)
end
local noise_heat -- Need a large-scale noise here so no heat blend
local elevation_chill = mapgen_rivers.settings.elevation_chill
local elevation_chill = mapgen_rivers.elevation_chill
local function get_temperature(x, y, z)
local pos = {x=x, y=z}
return noise_heat:get2d(pos) - y*elevation_chill
end
local glaciers = mapgen_rivers.settings.glaciers
local glacier_factor = mapgen_rivers.settings.glacier_factor
local glaciers = mapgen_rivers.glaciers
local glacier_factor = mapgen_rivers.glacier_factor
local init = false
-- On map generation, determine into which polygon every point (in 2D) will fall.
-- Also store polygon-specific data
local function make_polygons(minp, maxp)
local grid = mapgen_rivers.grid
local dem = grid.dem
local lakes = grid.lakes
local dirs = grid.dirs
local rivers = grid.rivers
local offset_x = grid.offset_x
local offset_z = grid.offset_y
if not init then
if glaciers then
noise_heat = minetest.get_perlin(mapgen_rivers.noise_params.heat)
@ -138,8 +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)
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
@ -149,24 +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,
}
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] = {}
@ -176,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)
@ -210,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
@ -219,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

View File

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

View File

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

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)

13
save.py Normal file
View File

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

11
settings.py Normal file
View File

@ -0,0 +1,11 @@
def read_config_file(fname):
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

View File

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

14
terrain_default.conf Normal file
View File

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

View File

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

View File

@ -1,385 +0,0 @@
-- 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 Boruvka algorithm.
-- Only flow_local and accumulate_flow are custom algorithms.
local function flow_local_semirandom(plist)
local sum = 0
for i=1, #plist do
sum = sum + plist[i]
end
--for _, p in ipairs(plist) do
--sum = sum + p
--end
if sum == 0 then
return 0
end
local r = math.random() * sum
for i=1, #plist do
local p = plist[i]
--for i, p in ipairs(plist) do
if r < p then
return i
end
r = r - p
end
return 0
end
local flow_methods = {
semirandom = flow_local_semirandom,
}
local function flow_routing(dem, dirs, lakes, method)
method = method or 'semirandom'
local flow_local = flow_methods[method] or flow_local_semirandom
dirs = dirs or {}
lakes = lakes or {}
-- Localize for performance
--local tinsert = table.insert
local tremove = table.remove
local mmax = math.max
local X, Y = dem.X, dem.Y
dirs.X = X
dirs.Y = Y
lakes.X = X
lakes.Y = Y
local i = 1
local dirs2 = {}
for i=1, X*Y do
dirs2[i] = 0
end
local singular = {}
for y=1, Y do
for x=1, X do
local zi = dem[i]
local plist = {
y<Y and mmax(zi-dem[i+X], 0) or 0, -- Southward
x<X and mmax(zi-dem[i+1], 0) or 0, -- Eastward
y>1 and mmax(zi-dem[i-X], 0) or 0, -- Northward
x>1 and mmax(zi-dem[i-1], 0) or 0, -- Westward
}
local d = flow_local(plist)
dirs[i] = d
if d == 0 then
singular[#singular+1] = i
elseif d == 1 then
dirs2[i+X] = dirs2[i+X] + 1
elseif d == 2 then
dirs2[i+1] = dirs2[i+1] + 2
elseif d == 3 then
dirs2[i-X] = dirs2[i-X] + 4
elseif d == 4 then
dirs2[i-1] = dirs2[i-1] + 8
end
i = i + 1
end
end
-- Compute basins and links
local nbasins = #singular
local basin_id = {}
local links = {}
local basin_links
local function add_link(i1, i2, b1, isY)
local b2
if i2 == 0 then
b2 = 0
else
b2 = basin_id[i2]
if b2 == 0 then
return
end
end
if b2 ~= b1 then
local elev = i2 == 0 and dem[i1] or mmax(dem[i1], dem[i2])
local l2 = basin_links[b2]
if not l2 then
l2 = {}
basin_links[b2] = l2
end
if not l2.elev or l2.elev > elev then
l2.elev = elev
l2.i = mmax(i1,i2)
l2.is_y = isY
l2[1] = b2
l2[2] = b1
end
end
end
for i=1, X*Y do
basin_id[i] = 0
end
--for ib, s in ipairs(singular) do
for ib=1, nbasins do
--local s = singular[ib]
local queue = {singular[ib]}
basin_links = {}
links[#links+1] = basin_links
--tinsert(links, basin_links)
while #queue > 0 do
local i = tremove(queue)
basin_id[i] = ib
local d = dirs2[i]
if d >= 8 then -- River coming from East
d = d - 8
queue[#queue+1] = i+1
--tinsert(queue, i+X)
elseif i%X > 0 then
add_link(i, i+1, ib, false)
else
add_link(i, 0, ib, false)
end
if d >= 4 then -- River coming from South
d = d - 4
queue[#queue+1] = i+X
--tinsert(queue, i+1)
elseif i <= X*(Y-1) then
add_link(i, i+X, ib, true)
else
add_link(i, 0, ib, true)
end
if d >= 2 then -- River coming from West
d = d - 2
queue[#queue+1] = i-1
--tinsert(queue, i-X)
elseif i%X ~= 1 then
add_link(i, i-1, ib, false)
else
add_link(i, 0, ib, false)
end
if d >= 1 then -- River coming from North
queue[#queue+1] = i-X
--tinsert(queue, i-1)
elseif i > X then
add_link(i, i-X, ib, true)
else
add_link(i, 0, ib, true)
end
end
end
dirs2 = nil
links[0] = {}
local nlinks = {}
for i=0, nbasins do
nlinks[i] = 0
end
--for ib1, blinks in ipairs(links) do
for ib1=1, #links do
for ib2, link in pairs(links[ib1]) do
if ib2 < ib1 then
links[ib2][ib1] = link
nlinks[ib1] = nlinks[ib1] + 1
nlinks[ib2] = nlinks[ib2] + 1
end
end
end
local lowlevel = {}
for i, n in pairs(nlinks) do
if n <= 8 then
lowlevel[i] = links[i]
end
end
local basin_graph = {}
for n=1, nbasins do
local b1, lnk1 = next(lowlevel)
lowlevel[b1] = nil
local b2
local lowest = math.huge
local lnk1 = links[b1]
local i = 0
for bn, bdata in pairs(lnk1) do
i = i + 1
if bdata.elev < lowest then
lowest = bdata.elev
b2 = bn
end
end
-- Add link to the graph
local bound = lnk1[b2]
local bb1, bb2 = bound[1], bound[2]
if not basin_graph[bb1] then
basin_graph[bb1] = {}
end
if not basin_graph[bb2] then
basin_graph[bb2] = {}
end
basin_graph[bb1][bb2] = bound
basin_graph[bb2][bb1] = bound
-- Merge basin b1 into b2
local lnk2 = links[b2]
-- First, remove the link between b1 and b2
lnk1[b2] = nil
lnk2[b1] = nil
nlinks[b2] = nlinks[b2] - 1
if nlinks[b2] == 8 then
lowlevel[b2] = lnk2
end
-- Look for basin 1's neighbours, and add them to basin 2 if they have a lower pass
for bn, bdata in pairs(lnk1) do
local lnkn = links[bn]
lnkn[b1] = nil
if lnkn[b2] then
nlinks[bn] = nlinks[bn] - 1
if nlinks[bn] == 8 then
lowlevel[bn] = lnkn
end
else
nlinks[b2] = nlinks[b2] + 1
if nlinks[b2] == 9 then
lowlevel[b2] = nil
end
end
if not lnkn[b2] or lnkn[b2].elev > bdata.elev then
lnkn[b2] = bdata
lnk2[bn] = bdata
end
end
end
local queue = {[0] = -math.huge}
local basin_lake = {}
for n=1, nbasins do
basin_lake[n] = 0
end
local reverse = {3, 4, 1, 2, [0]=0}
for n=1, nbasins do
local b1, elev1 = next(queue)
queue[b1] = nil
basin_lake[b1] = elev1
for b2, bound in pairs(basin_graph[b1]) do
-- Make b2 flow into b1
local i = bound.i
local dir = bound.is_y and 3 or 4
if basin_id[i] ~= b2 then
dir = dir - 2
if bound.is_y then
i = i - X
else
i = i - 1
end
elseif b1 == 0 then
dir = 0
end
repeat
dir, dirs[i] = dirs[i], dir
if dir == 1 then
i = i + X
elseif dir == 2 then
i = i + 1
elseif dir == 3 then
i = i - X
elseif dir == 4 then
i = i - 1
end
dir = reverse[dir]
until dir == 0
-- Add b2 into the queue
queue[b2] = mmax(elev1, bound.elev)
basin_graph[b2][b1] = nil
end
basin_graph[b1] = nil
end
for i=1, X*Y do
lakes[i] = basin_lake[basin_id[i]]
end
return dirs, lakes
end
local function accumulate(dirs, waterq)
waterq = waterq or {}
local X, Y = dirs.X, dirs.Y
--local tinsert = table.insert
local ndonors = {}
local waterq = {X=X, Y=Y}
for i=1, X*Y do
ndonors[i] = 0
waterq[i] = 1
end
--for i1, dir in ipairs(dirs) do
for i1=1, X*Y do
local i2
local dir = dirs[i1]
if dir == 1 then
i2 = i1+X
elseif dir == 2 then
i2 = i1+1
elseif dir == 3 then
i2 = i1-X
elseif dir == 4 then
i2 = i1-1
end
if i2 then
ndonors[i2] = ndonors[i2] + 1
end
end
for i1=1, X*Y do
if ndonors[i1] == 0 then
local i2 = i1
local dir = dirs[i2]
local w = waterq[i2]
while dir > 0 do
if dir == 1 then
i2 = i2 + X
elseif dir == 2 then
i2 = i2 + 1
elseif dir == 3 then
i2 = i2 - X
elseif dir == 4 then
i2 = i2 - 1
end
w = w + waterq[i2]
waterq[i2] = w
if ndonors[i2] > 1 then
ndonors[i2] = ndonors[i2] - 1
break
end
dir = dirs[i2]
end
end
end
return waterq
end
return {
flow_routing = flow_routing,
accumulate = accumulate,
flow_methods = flow_methods,
}

View File

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

102
view.py
View File

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

View File

@ -2,36 +2,51 @@
import numpy as np
import zlib
import sys
import os
import matplotlib.colors as mcol
import matplotlib.pyplot as plt
from view import stats, plot
from readconfig import read_conf_file
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')
os.chdir(sys.argv[1])
conf = read_conf_file('mapgen_rivers.conf')
if 'center' in conf:
center = conf['center'] == 'true'
else:
center = True
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')
if 'blocksize' in conf:
blocksize = float(conf['blocksize'])
else:
blocksize = 15.0
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')
def load_map(name, dtype, shape):
dtype = np.dtype(dtype)
with open(name, 'rb') as f:
data = f.read()
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('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)
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=blocksize)
plot(dem, lakes, scale=blocksize, center=center)
view_map(dem, lakes, rivers, scale)