26 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
30136bf60a Added scale (blocksize) parameter for view_map.py 2020-11-10 13:19:20 +01:00
9475b49b8d Removed duplicate calls to 2 scripts 2020-11-10 13:18:02 +01:00
36b49a7fe2 Add settings for parameters in terrain_rivers.py 2020-07-21 14:12:20 +02:00
103cd49d78 Optionally disable distorsion
by setting 'mapgen_rivers_distort = false' in minetest.conf
2020-07-21 14:01:29 +02:00
25c5cb2e1f Reverse axes order for heightmaps (iterate in Z direction first instead of X) 2020-07-21 12:46:23 +02:00
6f43430574 Added glaciers, and re-organized noise definitions 2020-05-24 12:09:21 +02:00
625768f967 Added snow and ice in function of temperature.
Uses noise parameters of builtin biomegen
2020-05-23 18:13:00 +02:00
4edd1a946e Horizontal shifting according to 3D noises:
makes slopes more irregular and natural-looking, allows overhanging.
This is done by generating an intermediate 2D elevation map and, for each node in 3D, add a 2D offset vector to the position, and seek this position on the heightmap.
2020-05-23 15:52:16 +02:00
f56857e804 Fix water not being set at lower chunk borders 2020-05-08 10:02:04 +02:00
a73a0dd80b Avoid some redundant calculation on corners
(not very significant, but why not)
2020-04-27 21:08:15 +02:00
a9ab0e53d3 Change folder structure: data files are now in a directory.
Also added a demo 400x400 map, that is overriden on pre-processing.
2020-04-26 23:29:36 +02:00
b429b302e1 Rewritten part of code to calculate river depth
Fixes bathymetry problems on turns or confluences, as well as abrupt riverbanks.
2020-04-26 22:19:05 +02:00
cd4b517585 terrain_rivers.py: mapsize is now the number of intervals
instead of the number of nodes.
2020-04-26 19:51:21 +02:00
cd90a21df4 Enhanced visualization code to display colormaps, and reuse the same code for initial and further viewing, in view_map.py 2020-04-26 18:30:29 +02:00
206c68813e Switch again to using river direction and flux instead of table of bounds 2020-04-26 18:10:23 +02:00
6af6795d90 Comment and clarify 2020-04-26 17:13:38 +02:00
49bc397718 Fix parameters for Simplex noise, to make sure the last octave has not a greater scale than 1
Also use a 401x401 grid instead of 400, so that there are 400 intervals
2020-04-26 16:52:40 +02:00
9700e948b9 Position should be strictly beyond river threshold to be a river
Prevents some wrongly placed water pixels.
2020-04-14 21:54:05 +02:00
55725ad94b Re-organized the code. All polygon-related calculations go to polygons.lua. 2020-04-14 21:11:54 +02:00
43211fc31b Removed useless functions get_point_location and geometry.area 2020-04-14 20:26:15 +02:00
14163681cc Use settings from minetest.conf 2020-04-14 18:40:51 +02:00
af7a7ce26d Compress data files (reduces size by a factor 3-4) 2020-04-13 15:59:34 +02:00
da98a538bb Draw outer corners of river turns. 2020-04-13 15:01:54 +02:00
b5db63d267 Re-added river flow data because needed for map preview 2020-04-13 15:01:44 +02:00
25 changed files with 826 additions and 364 deletions

3
.gitignore vendored
View File

@ -6,4 +6,7 @@ offset_x
offset_y
bounds_x
bounds_y
dirs
rivers
unused/
data/

View File

@ -1,29 +1,36 @@
mapgen_rivers
=============
Procedural map generator for Minetest 5.x. Still experimental and basic.
Procedural map generator for Minetest 5.x. Focused on river networks, and features valley erosion and lakes.
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/79073532-7a567f00-7ce7-11ea-9791-8fb453f5175d.png)
![Screenshot](https://user-images.githubusercontent.com/6905002/79541028-687b3000-8089-11ea-9209-c23c15d75383.png)
# Installation
This mod should be placed in the `/mods` directory like any other Minetest mod.
The Python part relies on external libraries that you need to install:
- `numpy`, a widely used library for numerical calculations
- `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 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.
## Pre-processing
Run the script `terrain_rivers.py` via command line. You can optionally append the map size (by default 400). Example for a 1000x1000 map:
```
./terrain_rivers.py 1000
```
For a default 400x400 map, it should take between 1 and 2 minutes. It will generate 5 files directly in the mod folder, containing the map data (1.4 MB for the default size).
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.

View File

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

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

@ -3,22 +3,33 @@ 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)
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
@ -32,19 +43,19 @@ def advection(dem, dirs, rivers, time, K=1, m=0.5, sea_level=0):
elif flow_dir == 4:
x1 -= 1
if remaining <= adv_time[y0,x0]:
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]
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')
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):
@ -78,9 +89,9 @@ class EvolutionModel:
self.flow_uptodate = False
def define_isostasy(self):
self.ref_isostasy = im.gaussian_filter(self.dem, self.flex_radius, mode='reflect')
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')
correction = (self.ref_isostasy - isostasy) * rate
self.dem = self.dem + correction
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,45 +1,37 @@
local function distance_to_segment(x1, y1, x2, y2, x, y)
-- get the distance between point (x,y) and segment (x1,y1)-(x2,y2)
local a = (x1-x2)^2 + (y1-y2)^2
local a = (x1-x2)^2 + (y1-y2)^2 -- square of distance
local b = (x1-x)^2 + (y1-y)^2
local c = (x2-x)^2 + (y2-y)^2
if a + b < c then
-- The closest point of the segment is the extremity 1
return math.sqrt(b)
elseif a + c < b then
-- The closest point of the segment is the extremity 2
return math.sqrt(c)
else
-- The closest point is on the segment
return math.abs(x1 * (y2-y) + x2 * (y-y1) + x * (y1-y2)) / math.sqrt(a)
end
end
local function transform_quadri(X, Y, x, y)
-- X, Y 4-vectors giving the coordinates of the 4 nodes
-- To index points in an irregular quadrilateral, giving x and y between 0 (one edge) and 1 (opposite edge)
-- X, Y 4-vectors giving the coordinates of the 4 vertices
-- x, y position to index.
local x1, x2, x3, x4 = unpack(X)
local y1, y2, y3, y4 = unpack(Y)
-- Compare distance to 2 opposite edges, they give the X coordinate
local d23 = distance_to_segment(x2,y2,x3,y3,x,y)
local d41 = distance_to_segment(x4,y4,x1,y1,x,y)
local xc = d41 / (d23+d41)
-- Same for the 2 other edges, they give the Y coordinate
local d12 = distance_to_segment(x1,y1,x2,y2,x,y)
local d34 = distance_to_segment(x3,y3,x4,y4,x,y)
local yc = d12 / (d12+d34)
return xc, yc
end
local function area(X, Y) -- Signed area of polygon, in function of direction of rotation. Clockwise = positive.
local n = #X
local sum = X[1]*Y[n] - X[n]*Y[1]
for i=2, n do
sum = sum + X[i]*Y[i-1] - X[i-1]*Y[i]
end
return sum/2
end
return {
distance_to_segment = distance_to_segment,
transform_quadri = transform_quadri,
area = area,
}
return transform_quadri

120
heightmap.lua Normal file
View File

@ -0,0 +1,120 @@
local modpath = minetest.get_modpath(minetest.get_current_modname()) .. '/'
local make_polygons = dofile(modpath .. 'polygons.lua')
local transform_quadri = dofile(modpath .. 'geometry.lua')
local blocksize = mapgen_rivers.blocksize
local sea_level = mapgen_rivers.sea_level
local riverbed_slope = mapgen_rivers.riverbed_slope
local MAP_BOTTOM = -31000
-- Linear interpolation
local function interp(v00, v01, v11, v10, xf, zf)
local v0 = v01*xf + v00*(1-xf)
local v1 = v11*xf + v10*(1-xf)
return v1*zf + v0*(1-zf)
end
local function heightmaps(minp, maxp)
local polygons = make_polygons(minp, maxp)
local incr = maxp.z-minp.z+1
local terrain_height_map = {}
local lake_height_map = {}
local i = 1
for z=minp.z, maxp.z do
for x=minp.x, maxp.x do
local poly = polygons[i]
if poly then
local xf, zf = transform_quadri(poly.x, poly.z, x/blocksize, z/blocksize)
local i00, i01, i11, i10 = unpack(poly.i)
-- Load river width on 4 edges and corners
local r_west, r_north, r_east, r_south = unpack(poly.rivers)
local c_NW, c_NE, c_SE, c_SW = unpack(poly.river_corners)
-- Calculate the depth factor for each edge and corner.
-- Depth factor:
-- < 0: outside river
-- = 0: on riverbank
-- > 0: inside river
local depth_factors = {
r_west - xf,
r_north - zf,
xf - r_east,
zf - r_south,
c_NW-xf-zf,
xf-zf-c_NE,
xf+zf-c_SE,
zf-xf-c_SW,
}
-- Find the maximal depth factor and determine to which river it belongs
local depth_factor_max = 0
local imax = 0
for i=1, 8 do
if depth_factors[i] >= depth_factor_max then
depth_factor_max = depth_factors[i]
imax = i
end
end
-- Transform the coordinates to have xf and zf = 0 or 1 in rivers (to avoid rivers having lateral slope and to accomodate the surrounding smoothly)
if imax == 0 then
local x0 = math.max(r_west, c_NW-zf, zf-c_SW)
local x1 = math.min(r_east, c_NE+zf, c_SE-zf)
local z0 = math.max(r_north, c_NW-xf, xf-c_NE)
local z1 = math.min(r_south, c_SW+xf, c_SE-xf)
xf = (xf-x0) / (x1-x0)
zf = (zf-z0) / (z1-z0)
elseif imax == 1 then
xf = 0
elseif imax == 2 then
zf = 0
elseif imax == 3 then
xf = 1
elseif imax == 4 then
zf = 1
elseif imax == 5 then
xf, zf = 0, 0
elseif imax == 6 then
xf, zf = 1, 0
elseif imax == 7 then
xf, zf = 1, 1
elseif imax == 8 then
xf, zf = 0, 1
end
-- Determine elevation by interpolation
local vdem = poly.dem
local terrain_height = math.floor(0.5+interp(
vdem[1],
vdem[2],
vdem[3],
vdem[4],
xf, zf
))
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
terrain_height_map[i] = terrain_height
lake_height_map[i] = lake_height
else
terrain_height_map[i] = MAP_BOTTOM
lake_height_map[i] = MAP_BOTTOM
end
i = i + 1
end
end
return terrain_height_map, lake_height_map
end
return heightmaps

401
init.lua
View File

@ -1,59 +1,26 @@
mapgen_rivers = {}
local modpath = minetest.get_modpath(minetest.get_current_modname()) .. '/'
local worldpath = minetest.get_worldpath() .. '/'
local load_map = dofile(modpath .. 'load.lua')
local geometry = dofile(modpath .. 'geometry.lua')
local function copy_if_needed(filename)
local wfilename = worldpath..filename
local wfile = io.open(wfilename, 'r')
if wfile then
wfile:close()
return
end
local mfilename = modpath..filename
local mfile = io.open(mfilename, 'r')
local wfile = io.open(wfilename, 'w')
wfile:write(mfile:read("*all"))
mfile:close()
wfile:close()
dofile(modpath .. 'settings.lua')
local blocksize = mapgen_rivers.blocksize
local sea_level = mapgen_rivers.sea_level
local riverbed_slope = mapgen_rivers.riverbed_slope
local elevation_chill = mapgen_rivers.elevation_chill
local use_distort = mapgen_rivers.distort
local use_biomes = mapgen_rivers.biomes
local use_biomegen_mod = use_biomes and minetest.global_exists('biomegen')
use_biomes = use_biomes and not use_biomegen_mod
if use_biomegen_mod then
biomegen.set_elevation_chill(elevation_chill)
end
dofile(modpath .. 'noises.lua')
copy_if_needed('size')
local sfile = io.open(worldpath..'size')
local X = tonumber(sfile:read('*l'))
local Z = tonumber(sfile:read('*l'))
copy_if_needed('dem')
local dem = load_map(worldpath..'dem', 2, true)
copy_if_needed('lakes')
local lakes = load_map(worldpath..'lakes', 2, true)
copy_if_needed('bounds_x')
local bounds_x = load_map(worldpath..'bounds_x', 4, false)
copy_if_needed('bounds_y')
local bounds_z = load_map(worldpath..'bounds_y', 4, false)
copy_if_needed('offset_x')
local offset_x = load_map(worldpath..'offset_x', 1, true)
for k, v in ipairs(offset_x) do
offset_x[k] = (v+0.5)/256
end
copy_if_needed('offset_y')
local offset_z = load_map(worldpath..'offset_y', 1, true)
for k, v in ipairs(offset_z) do
offset_z[k] = (v+0.5)/256
end
local function index(x, z)
return z*X+x+1
end
local function get_point_location(x, z)
local i = index(x, z)
return x+offset_x[i], z+offset_z[i]
end
local heightmaps = dofile(modpath .. 'heightmap.lua')
-- Linear interpolation
local function interp(v00, v01, v11, v10, xf, zf)
local v0 = v01*xf + v00*(1-xf)
local v1 = v11*xf + v10*(1-xf)
@ -62,216 +29,194 @@ end
local data = {}
local blocksize = 12
local sea_level = 1
local min_catchment = 25
local max_catchment = 40000
local storage = minetest.get_mod_storage()
if storage:contains("blocksize") then
blocksize = storage:get_int("blocksize")
else
storage:set_int("blocksize", blocksize)
end
if storage:contains("sea_level") then
sea_level = storage:get_int("sea_level")
else
storage:set_int("sea_level", sea_level)
end
if storage:contains("min_catchment") then
min_catchment = storage:get_float("min_catchment")
else
storage:set_float("min_catchment", min_catchment)
end
if storage:contains("max_catchment") then
max_catchment = storage:get_float("max_catchment")
else
storage:set_float("max_catchment", max_catchment)
end
-- Width coefficients: coefficients solving
-- wfactor * min_catchment ^ wpower = 1/(2*blocksize)
-- wfactor * max_catchment ^ wpower = 1
local wpower = math.log(2*blocksize)/math.log(max_catchment/min_catchment)
local wfactor = 1 / max_catchment ^ wpower
local function river_width(flow)
flow = math.abs(flow)
if flow < min_catchment then
return 0
end
return math.min(wfactor * flow ^ wpower, 1)
end
local noise_x_obj, noise_z_obj, noise_distort_obj, noise_heat_obj, noise_heat_blend_obj
local noise_x_map = {}
local noise_z_map = {}
local noise_distort_map = {}
local noise_heat_map = {}
local noise_heat_blend_map = {}
local mapsize
local init = false
local function generate(minp, maxp, seed)
local chulens = {
x = maxp.x-minp.x+1,
y = maxp.y-minp.y+1,
z = maxp.z-minp.z+1,
}
if not init then
mapsize = {
x = chulens.x,
y = chulens.y+1,
z = chulens.z,
}
if use_distort then
noise_x_obj = minetest.get_perlin_map(mapgen_rivers.noise_params.distort_x, mapsize)
noise_z_obj = minetest.get_perlin_map(mapgen_rivers.noise_params.distort_z, mapsize)
noise_distort_obj = minetest.get_perlin_map(mapgen_rivers.noise_params.distort_amplitude, chulens)
end
if use_biomes then
noise_heat_obj = minetest.get_perlin_map(mapgen_rivers.noise_params.heat, chulens)
noise_heat_blend_obj = minetest.get_perlin_map(mapgen_rivers.noise_params.heat_blend, chulens)
end
init = true
end
local minp2d = {x=minp.x, y=minp.z}
if use_distort then
noise_x_obj:get_3d_map_flat(minp, noise_x_map)
noise_z_obj:get_3d_map_flat(minp, noise_z_map)
noise_distort_obj:get_2d_map_flat(minp2d, noise_distort_map)
end
if use_biomes then
noise_heat_obj:get_2d_map_flat(minp2d, noise_heat_map)
noise_heat_blend_obj:get_2d_map_flat(minp2d, noise_heat_blend_map)
end
local terrain_map, lake_map, incr, i_origin
if use_distort then
local xmin, xmax, zmin, zmax = minp.x, maxp.x, minp.z, maxp.z
local i = 0
local i2d = 0
for z=minp.z, maxp.z do
for y=minp.y, maxp.y+1 do
for x=minp.x, maxp.x do
i = i+1
i2d = i2d+1
local distort = noise_distort_map[i2d]
local xv = noise_x_map[i]*distort + x
if xv < xmin then xmin = xv end
if xv > xmax then xmax = xv end
noise_x_map[i] = xv
local zv = noise_z_map[i]*distort + z
if zv < zmin then zmin = zv end
if zv > zmax then zmax = zv end
noise_z_map[i] = zv
end
i2d = i2d-chulens.x
end
end
local pminp = {x=math.floor(xmin), z=math.floor(zmin)}
local pmaxp = {x=math.floor(xmax)+1, z=math.floor(zmax)+1}
incr = pmaxp.x-pminp.x+1
i_origin = 1 - pminp.z*incr - pminp.x
terrain_map, lake_map = heightmaps(pminp, pmaxp)
else
terrain_map, lake_map = heightmaps(minp, maxp)
end
local c_stone = minetest.get_content_id("default:stone")
local c_dirt = minetest.get_content_id("default:dirt")
local c_lawn = minetest.get_content_id("default:dirt_with_grass")
local c_dirtsnow = minetest.get_content_id("default:dirt_with_snow")
local c_snow = minetest.get_content_id("default:snowblock")
local c_sand = minetest.get_content_id("default:sand")
local c_water = minetest.get_content_id("default:water_source")
local c_rwater = minetest.get_content_id("default:river_water_source")
local c_ice = minetest.get_content_id("default:ice")
local vm, emin, emax = minetest.get_mapgen_object("voxelmanip")
vm:get_data(data)
local a = VoxelArea:new({MinEdge = emin, MaxEdge = emax})
local ystride = a.ystride -- Tip : the ystride of a VoxelArea is the number to add to the array index to get the index of the position above. It's faster because it avoids to completely recalculate the index.
local chulens = maxp.z - minp.z + 1
local polygons = {}
local xpmin, xpmax = math.max(math.floor(minp.x/blocksize - 0.5), 0), math.min(math.ceil(maxp.x/blocksize), X-2)
local zpmin, zpmax = math.max(math.floor(minp.z/blocksize - 0.5), 0), math.min(math.ceil(maxp.z/blocksize), Z-2)
local nid = mapsize.x*(mapsize.y-1) + 1
local incrY = -mapsize.x
local incrX = 1 - mapsize.y*incrY
local incrZ = mapsize.x*mapsize.y - mapsize.x*incrX - mapsize.x*mapsize.y*incrY
for xp = xpmin, xpmax do
for zp=zpmin, zpmax do
local iA = index(xp, zp)
local iB = index(xp+1, zp)
local iC = index(xp+1, zp+1)
local iD = index(xp, zp+1)
local poly_x = {offset_x[iA]+xp, offset_x[iB]+xp+1, offset_x[iC]+xp+1, offset_x[iD]+xp}
local poly_z = {offset_z[iA]+zp, offset_z[iB]+zp, offset_z[iC]+zp+1, offset_z[iD]+zp+1}
local polygon = {x=poly_x, z=poly_z, i={iA, iB, iC, iD}}
local bounds = {}
local xmin = math.max(math.floor(blocksize*math.min(unpack(poly_x)))+1, minp.x)
local xmax = math.min(math.floor(blocksize*math.max(unpack(poly_x))), maxp.x)
for x=xmin, xmax do
bounds[x] = {}
local i2d = 1
for z = minp.z, maxp.z do
for x = minp.x, maxp.x do
local ivm = a:index(x, minp.y, z)
local ground_above = false
local temperature
if use_biomes then
temperature = noise_heat_map[i2d]+noise_heat_blend_map[i2d]
end
local terrain, lake
if not use_distort then
terrain = terrain_map[i2d]
lake = lake_map[i2d]
end
local i1 = 4
for i2=1, 4 do -- Loop on 4 edges
local x1, x2 = poly_x[i1], poly_x[i2]
local lxmin = math.floor(blocksize*math.min(x1, x2))+1
local lxmax = math.floor(blocksize*math.max(x1, x2))
if lxmin <= lxmax then
local z1, z2 = poly_z[i1], poly_z[i2]
local a = (z1-z2) / (x1-x2)
local b = blocksize*(z1 - a*x1)
for x=math.max(lxmin, minp.x), math.min(lxmax, maxp.x) do
table.insert(bounds[x], a*x+b)
end
end
i1 = i2
end
for x=xmin, xmax do
local xlist = bounds[x]
table.sort(xlist)
local c = math.floor(#xlist/2)
for l=1, c do
local zmin = math.max(math.floor(xlist[l*2-1])+1, minp.z)
local zmax = math.min(math.floor(xlist[l*2]), maxp.z)
local i = (x-minp.x) * chulens + (zmin-minp.z) + 1
for z=zmin, zmax do
polygons[i] = polygon
i = i + 1
end
end
end
for y = maxp.y+1, minp.y, -1 do
if use_distort then
local xn = noise_x_map[nid]
local zn = noise_z_map[nid]
local x0 = math.floor(xn)
local z0 = math.floor(zn)
polygon.dem = {dem[iA], dem[iB], dem[iC], dem[iD]}
polygon.lake = math.min(lakes[iA], lakes[iB], lakes[iC], lakes[iD])
local i0 = i_origin + z0*incr + x0
local i1 = i0+1
local i2 = i1+incr
local i3 = i2-1
local river_west = river_width(bounds_z[iA])
local river_north = river_width(bounds_x[iA-zp])
local river_east = 1-river_width(bounds_z[iB])
local river_south = 1-river_width(bounds_x[iD-zp-1])
if river_west > river_east then
local mean = (river_west + river_east) / 2
river_west = mean
river_east = mean
end
if river_north > river_south then
local mean = (river_north + river_south) / 2
river_north = mean
river_south = mean
end
polygon.rivers = {river_west, river_north, river_east, river_south}
end
end
local i = 1
for x = minp.x, maxp.x do
for z = minp.z, maxp.z do
local poly = polygons[i]
if poly then
local xf, zf = geometry.transform_quadri(poly.x, poly.z, x/blocksize, z/blocksize)
local i00, i01, i11, i10 = unpack(poly.i)
local is_river = false
local r_west, r_north, r_east, r_south = unpack(poly.rivers)
if xf >= r_east then
is_river = true
xf = 1
elseif xf <= r_west then
is_river = true
xf = 0
end
if zf >= r_south then
is_river = true
zf = 1
elseif zf <= r_north then
is_river = true
zf = 0
terrain = interp(terrain_map[i0], terrain_map[i1], terrain_map[i2], terrain_map[i3], xn-x0, zn-z0)
lake = math.min(lake_map[i0], lake_map[i1], lake_map[i2], lake_map[i3])
end
if not is_river then
xf = (xf-r_west) / (r_east-r_west)
zf = (zf-r_north) / (r_south-r_north)
end
if y <= maxp.y then
local vdem = poly.dem
local terrain_height = math.floor(0.5+interp(
vdem[1],
vdem[2],
vdem[3],
vdem[4],
xf, zf
))
local lake_height = math.floor(poly.lake)
local is_lake = lake_height > terrain_height
local ivm = a:index(x, minp.y-1, z)
if terrain_height >= minp.y then
for y=minp.y, math.min(maxp.y, terrain_height) do
if y == terrain_height then
if is_lake or y <= sea_level then
data[ivm] = c_sand
elseif is_river then
data[ivm] = c_rwater
else
data[ivm] = c_lawn
end
else
local is_lake = lake > terrain
local ivm = a:index(x, y, z)
if y <= terrain then
if not use_biomes or y <= terrain-1 or ground_above then
data[ivm] = c_stone
elseif is_lake or y < sea_level then
data[ivm] = c_sand
else
local temperature_y = temperature - y*elevation_chill
if temperature_y >= 15 then
data[ivm] = c_lawn
elseif temperature_y >= 0 then
data[ivm] = c_dirtsnow
else
data[ivm] = c_snow
end
end
ivm = ivm + ystride
elseif y <= lake and lake > sea_level then
if not use_biomes or temperature - y*elevation_chill >= 0 then
data[ivm] = c_rwater
else
data[ivm] = c_ice
end
elseif y <= sea_level then
data[ivm] = c_water
end
end
if lake_height > sea_level then
if is_lake and lake_height > minp.y then
for y=math.max(minp.y, terrain_height+1), math.min(maxp.y, lake_height) do
data[ivm] = c_rwater
ivm = ivm + ystride
end
end
else
for y=math.max(minp.y, terrain_height+1), math.min(maxp.y, sea_level) do
data[ivm] = c_water
ivm = ivm + ystride
end
ground_above = y <= terrain
ivm = ivm + ystride
if use_distort then
nid = nid + incrY
end
end
i = i + 1
if use_distort then
nid = nid + incrX
end
i2d = i2d + 1
end
if use_distort then
nid = nid + incrZ
end
end
vm:set_data(data)
minetest.generate_ores(vm, minp, maxp)
if use_biomegen_mod then
biomegen.generate_all(data, a, vm, minp, maxp, seed)
else
vm:set_data(data)
mietest.generate_ores(vm, minp, maxp)
end
vm:set_lighting({day = 0, night = 0})
vm:calc_lighting()
vm:update_liquids()

View File

@ -1,11 +1,14 @@
local function load_map(filename, bytes, signed)
local file = io.open(filename, 'r')
local worldpath = minetest.get_worldpath() .. "/river_data/"
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 map = {}
local size = math.floor(#data/bytes)
for i=1, size do
local i0, i1 = (i-1)*bytes+1, i*bytes
local elements = {data:byte(i0, i1)}

4
mod.conf Normal file
View File

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

37
noises.lua Normal file
View File

@ -0,0 +1,37 @@
mapgen_rivers.noise_params = {
distort_x = {
offset = 0,
scale = 1,
seed = -4574,
spread = {x=64, y=32, z=64},
octaves = 3,
persistence = 0.75,
lacunarity = 2,
},
distort_z = {
offset = 0,
scale = 1,
seed = -7940,
spread = {x=64, y=32, z=64},
octaves = 3,
persistence = 0.75,
lacunarity = 2,
},
distort_amplitude = {
offset = 0,
scale = 10,
seed = 676,
spread = {x=1024, y=1024, z=1024},
octaves = 5,
persistence = 0.5,
lacunarity = 2,
flags = "absvalue",
},
heat = minetest.get_mapgen_setting_noiseparams('mg_biome_np_heat'),
heat_blend = minetest.get_mapgen_setting_noiseparams('mg_biome_np_heat_blend'),
}
mapgen_rivers.noise_params.heat.offset = mapgen_rivers.noise_params.heat.offset + mapgen_rivers.sea_level*mapgen_rivers.elevation_chill

205
polygons.lua Normal file
View File

@ -0,0 +1,205 @@
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 = minetest.get_worldpath() .. '/river_data/'
minetest.mkdir(world_data_path)
local load_map = dofile(modpath .. 'load.lua')
local function copy_if_needed(filename)
local wfilename = world_data_path..filename
local wfile = io.open(wfilename, 'r')
if wfile then
wfile:close()
return
end
local mfilename = mod_data_path..filename
local mfile = io.open(mfilename, 'r')
local wfile = io.open(wfilename, 'w')
wfile:write(mfile:read("*all"))
mfile:close()
wfile:close()
end
copy_if_needed('size')
local sfile = io.open(world_data_path..'size')
local X = tonumber(sfile:read('*l'))
local Z = tonumber(sfile:read('*l'))
sfile:close()
copy_if_needed('dem')
local dem = load_map('dem', 2, true, X*Z)
copy_if_needed('lakes')
local lakes = load_map('lakes', 2, true, X*Z)
copy_if_needed('dirs')
local dirs = load_map('dirs', 1, false, X*Z)
copy_if_needed('rivers')
local rivers = load_map('rivers', 4, false, X*Z)
copy_if_needed('offset_x')
local offset_x = load_map('offset_x', 1, true, X*Z)
for k, v in ipairs(offset_x) do
offset_x[k] = (v+0.5)/256
end
copy_if_needed('offset_y')
local offset_z = load_map('offset_y', 1, true, X*Z)
for k, v in ipairs(offset_z) do
offset_z[k] = (v+0.5)/256
end
-- To index a flat array representing a 2D map
local function index(x, z)
return z*X+x+1
end
local blocksize = mapgen_rivers.blocksize
local min_catchment = mapgen_rivers.min_catchment
local max_catchment = mapgen_rivers.max_catchment
-- Width coefficients: coefficients solving
-- wfactor * min_catchment ^ wpower = 1/(2*blocksize)
-- wfactor * max_catchment ^ wpower = 1
local wpower = math.log(2*blocksize)/math.log(max_catchment/min_catchment)
local wfactor = 1 / max_catchment ^ wpower
local function river_width(flow)
flow = math.abs(flow)
if flow < min_catchment then
return 0
end
return math.min(wfactor * flow ^ wpower, 1)
end
local noise_heat -- Need a large-scale noise here so no heat blend
local elevation_chill = mapgen_rivers.elevation_chill
local function get_temperature(x, y, z)
local pos = {x=x, y=z}
return noise_heat:get2d(pos) - y*elevation_chill
end
local glaciers = mapgen_rivers.glaciers
local glacier_factor = mapgen_rivers.glacier_factor
local init = false
-- On map generation, determine into which polygon every point (in 2D) will fall.
-- Also store polygon-specific data
local function make_polygons(minp, maxp)
if not init then
if glaciers then
noise_heat = minetest.get_perlin(mapgen_rivers.noise_params.heat)
end
init = true
end
local chulens = maxp.x - minp.x + 1
local polygons = {}
-- Determine the minimum and maximum coordinates of the polygons that could be on the chunk, knowing that they have an average size of 'blocksize' and a maximal offset of 0.5 blocksize.
local xpmin, xpmax = math.max(math.floor(minp.x/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
for zp=zpmin, zpmax do
local iA = index(xp, zp)
local iB = index(xp+1, zp)
local iC = index(xp+1, zp+1)
local iD = index(xp, zp+1)
-- Extract the vertices of the polygon
local poly_x = {offset_x[iA]+xp, 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(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] = {}
end
local i1 = 4
for i2=1, 4 do -- Loop on 4 edges
local z1, z2 = poly_z[i1], poly_z[i2]
-- Calculate the integer Z positions over which this edge spans
local lzmin = math.floor(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 = 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)
end
end
i1 = i2
end
for z=zmin, zmax do
-- Now sort the bounds list
local zlist = bounds[z]
table.sort(zlist)
local c = math.floor(#zlist/2)
for l=1, c do
-- Take pairs of X coordinates: all positions between them belong to the polygon.
local xmin = math.max(math.floor(zlist[l*2-1])+1, minp.x)
local xmax = math.min(math.floor(zlist[l*2]), maxp.x)
local i = (z-minp.z) * chulens + (xmin-minp.x) + 1
for x=xmin, xmax do
-- Fill the map at these places
polygons[i] = polygon
i = i + 1
end
end
end
local poly_dem = {dem[iA], dem[iB], dem[iC], dem[iD]}
polygon.dem = poly_dem
polygon.lake = math.min(lakes[iA], lakes[iB], lakes[iC], lakes[iD])
-- Now, rivers.
-- Load river flux values for the 4 corners
local riverA = river_width(rivers[iA])
local riverB = river_width(rivers[iB])
local riverC = river_width(rivers[iC])
local riverD = river_width(rivers[iD])
if glaciers then -- Widen the river
if get_temperature(poly_x[1]*blocksize, poly_dem[1], poly_z[1]*blocksize) < 0 then
riverA = math.min(riverA*glacier_factor, 1)
end
if get_temperature(poly_x[2]*blocksize, poly_dem[2], poly_z[2]*blocksize) < 0 then
riverB = math.min(riverB*glacier_factor, 1)
end
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]*blocksize, poly_dem[4], poly_z[4]*blocksize) < 0 then
riverD = math.min(riverD*glacier_factor, 1)
end
end
polygon.river_corners = {riverA, 1-riverB, 2-riverC, 1-riverD}
-- Flow directions
local dirA, dirB, dirC, dirD = dirs[iA], dirs[iB], dirs[iC], dirs[iD]
-- Determine the river flux on the edges, by testing dirs values
local river_west = (dirA==1 and riverA or 0) + (dirD==3 and riverD or 0)
local river_north = (dirA==2 and riverA or 0) + (dirB==4 and riverB or 0)
local river_east = 1 - (dirB==1 and riverB or 0) - (dirC==3 and riverC or 0)
local river_south = 1 - (dirD==2 and riverD or 0) - (dirC==4 and riverC or 0)
polygon.rivers = {river_west, river_north, river_east, river_south}
end
end
return polygons
end
return make_polygons

View File

@ -19,9 +19,14 @@ neighbours_dirs = np.array([
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))
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
@ -41,10 +46,11 @@ def flow_dirs_lakes(dem, random=0):
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
dirs[-2:,:] = 1 # Border pixels flow outside the map
dirs[:,-2:] = 2
dirs[ :2,:] = 3
dirs[:, :2] = 4
@ -56,21 +62,26 @@ def flow_dirs_lakes(dem, random=0):
heapq.heappush(borders, (alt, altmax, y, x))
while len(borders) > 0:
(alt, altmax, y, x) = heapq.heappop(borders)
(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
neighbours += empty_neighbours * neighbours_dirs
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)
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:
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
@ -79,13 +90,13 @@ def accumulate(dirs, dem=None):
def calculate_quantity(y, x):
if quantity[y,x] > 0:
return quantity[y,x]
q = 1
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
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)
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
@ -96,5 +107,9 @@ def accumulate(dirs, dem=None):
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

@ -1,8 +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(data.tobytes())
f.write(bin_data)

54
settings.lua Normal file
View File

@ -0,0 +1,54 @@
local storage = minetest.get_mod_storage()
local settings = minetest.settings
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
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
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

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

View File

@ -7,110 +7,110 @@ 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:
mapsize = int(sys.argv[1])
else:
mapsize = 400
if os.path.isfile(sys.argv[1]):
params = settings.read_config_file(sys.argv[1])
else:
mapsize = int(sys.argv[1])
scale = mapsize / 2
n = np.zeros((mapsize, mapsize))
def get_setting(name, default):
if name in params:
return params[name]
return default
#micronoise_depth = 0.05
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.log2(mapsize)),
"persistence" : 0.5,
"lacunarity" : 2.,
"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)
for x in range(mapsize):
for y in range(mapsize):
# 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)
#micronoise = np.random.rand(mapsize, mapsize)
#nn = np.exp(n*2) + micronoise*micronoise_depth
nn = n*mapsize/5 + mapsize/20
nn = n*vscale + offset
# Initialize landscape evolution model
print('Initializing model')
model = EvolutionModel(nn, K=1, m=0.35, d=1, sea_level=0)
model = EvolutionModel(nn, K=1, m=0.35, d=1, sea_level=0, flex_radius=flex_radius)
print('Flow calculation 1')
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()
print('Advection 1')
model.advection(2)
print('Isostatic equilibration 1')
model.adjust_isostasy()
print('Flow calculation 2')
model.calculate_flow()
print('Diffusion')
model.diffusion(4)
print('Advection 2')
model.advection(2)
print('Isostatic equilibration 2')
model.adjust_isostasy()
print('Flow calculation 3')
model.calculate_flow()
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)
ox, oy = bounds.twist(bx, by, bounds.get_fixed(model.dirs))
offset_x, offset_y = bounds.twist(bx, by, bounds.get_fixed(model.dirs))
offset_x = np.clip(np.floor(ox * 256), -128, 127)
offset_y = np.clip(np.floor(oy * 256), -128, 127)
# 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(np.abs(bx), 'bounds_x', dtype='>i4')
save(np.abs(by), 'bounds_y', dtype='>i4')
save(offset_x, 'offset_x', dtype='i1')
save(offset_y, 'offset_y', dtype='i1')
save(model.dirs, 'dirs', dtype='u1')
save(model.rivers, 'rivers', dtype='>u4')
with open('size', 'w') as sfile:
sfile.write('{:d}\n{:d}'.format(mapsize, mapsize))
sfile.write('{:d}\n{:d}'.format(mapsize+1, mapsize+1))
# Display the map if matplotlib is found
try:
import matplotlib.pyplot as plt
plt.subplot(2,2,1)
plt.pcolormesh(nn, cmap='viridis')
plt.gca().set_aspect('equal', 'box')
#plt.colorbar(orientation='horizontal')
plt.title('Raw elevation')
plt.subplot(2,2,2)
plt.pcolormesh(model.lakes, cmap='viridis')
plt.gca().set_aspect('equal', 'box')
#plt.colorbar(orientation='horizontal')
plt.title('Lake surface elevation')
plt.subplot(2,2,3)
plt.pcolormesh(model.dem, cmap='viridis')
plt.gca().set_aspect('equal', 'box')
#plt.colorbar(orientation='horizontal')
plt.title('Elevation after advection')
plt.subplot(2,2,4)
plt.pcolormesh(model.rivers, vmin=0, vmax=mapsize**2/25, cmap='Blues')
plt.gca().set_aspect('equal', 'box')
#plt.colorbar(orientation='horizontal')
plt.title('Rivers discharge')
plt.show()
from view_map import view_map
view_map(model.dem, model.lakes, model.rivers)
except:
pass

View File

@ -1,30 +1,52 @@
#!/usr/bin/env python3
import numpy as np
import zlib
import matplotlib.colors as mcol
import matplotlib.pyplot as plt
shape = np.loadtxt('size', dtype='u4')
n = shape[0] * shape[1]
dem = np.fromfile('dem', dtype='>i2').reshape(shape)
lakes = np.fromfile('lakes', dtype='>i2').reshape(shape)
rivers = np.fromfile('rivers', dtype='>u4').reshape(shape)
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')
plt.subplot(1,3,1)
plt.pcolormesh(dem, cmap='viridis')
plt.gca().set_aspect('equal', 'box')
#plt.colorbar(orientation='horizontal')
plt.title('Raw elevation')
plt.subplot(1,3,2)
plt.pcolormesh(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')
plt.subplot(1,3,2)
plt.pcolormesh(lakes, cmap='viridis')
plt.gca().set_aspect('equal', 'box')
#plt.colorbar(orientation='horizontal')
plt.title('Lake surface elevation')
plt.subplot(1,3,3)
plt.pcolormesh(np.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.subplot(1,3,3)
plt.pcolormesh(np.log(rivers), vmin=0, vmax=np.log(n/25), cmap='Blues')
plt.gca().set_aspect('equal', 'box')
#plt.colorbar(orientation='horizontal')
plt.title('Rivers discharge')
plt.show()
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)
rivers = load_map('rivers', '>u4', shape)
view_map(dem, lakes, rivers, scale)