mtcontrib/mapgen_rivers
1
0
Fork 0
mirror of https://gitlab.com/gaelysam/mapgen_rivers.git synced 2025-05-06 19:20:35 +02:00
Code Issues Projects Releases Wiki Activity

Compare commits

base: mtcontrib:master
Branches Tags
mtcontrib:master
mtcontrib:caves
mtcontrib:dev
mtcontrib:mapgen_thread
mtcontrib:fix_rivers
mtcontrib:refactor_settings
mtcontrib:margin
mtcontrib:distinguish_sea
mtcontrib:terrainlib_lua
mtcontrib:variable_erosion
mtcontrib:old_master
mtcontrib:biomegen
mtcontrib:horznoise
mtcontrib:settings
mtcontrib:variable_width
mtcontrib:twist
mtcontrib:v1.0.2
mtcontrib:v1.0.1
mtcontrib:v1.0
..
compare: mtcontrib:v1.0.2
Branches Tags
mtcontrib:caves
mtcontrib:master
mtcontrib:dev
mtcontrib:mapgen_thread
mtcontrib:fix_rivers
mtcontrib:refactor_settings
mtcontrib:margin
mtcontrib:distinguish_sea
mtcontrib:terrainlib_lua
mtcontrib:variable_erosion
mtcontrib:old_master
mtcontrib:biomegen
mtcontrib:horznoise
mtcontrib:settings
mtcontrib:variable_width
mtcontrib:twist
mtcontrib:v1.0.2
mtcontrib:v1.0.1
mtcontrib:v1.0

No commits in common. "master" and "v1.0.2" have entirely different histories.
master ... v1.0.2

18 changed files with 859 additions and 1151 deletions
Show Stats Expand all files Collapse all files

2
README.md
View File

@ -17,7 +17,7 @@ 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
No required dependency, but [`biomegen`](https://gitlab.com/gaelysam/biomegen) recommended (provides biome system).
Mod dependencies: `default` required, and [`biomegen`](https://github.com/Gael-de-Sailly/biomegen) optional (provides biome system).
# Installation
This mod should be placed in the `mods/` directory of Minetest like any other mod.

45
compatibility.lua
View File

@ -1,24 +1,3 @@
-- Fix compatibility for settings-related changes
-- Only loaded if the versions of the mod and the world mismatch
local function version_is_lower(v1, v2)
local d1, c1, d2, c2
while #v1 > 0 and #v2 > 0 do
d1, c1, v1 = v1:match("^(%d*)(%D*)(.*)$")
d2, c2, v2 = v2:match("^(%d*)(%D*)(.*)$")
d1 = tonumber(d1) or -1
d2 = tonumber(d2) or -1
if d1 ~= d2 then
return d1 < d2
end
if c1 ~= c2 then
return c1 < c2
end
end
return false
end
local function fix_min_catchment(settings, is_global)
local prefix = is_global and "mapgen_rivers_" or ""
@ -42,18 +21,6 @@ local function fix_compatibility_minetest(settings)
if previous_version == "0.0" then
fix_min_catchment(settings, true)
end
if version_is_lower(previous_version, "1.0.2-dev1") then
local blocksize = tonumber(settings:get("mapgen_rivers_blocksize") or 15)
local grid_x_size = tonumber(settings:get("mapgen_rivers_grid_x_size"))
if grid_x_size then
settings:set("mapgen_rivers_map_x_size", tostring(grid_x_size * blocksize))
end
local grid_z_size = tonumber(settings:get("mapgen_rivers_grid_z_size"))
if grid_z_size then
settings:set("mapgen_rivers_map_z_size", tostring(grid_z_size * blocksize))
end
end
end
local function fix_compatibility_mapgen_rivers(settings)
@ -62,18 +29,6 @@ local function fix_compatibility_mapgen_rivers(settings)
if previous_version == "0.0" then
fix_min_catchment(settings, false)
end
if version_is_lower(previous_version, "1.0.2-dev1") then
local blocksize = tonumber(settings:get("blocksize") or 15)
local grid_x_size = tonumber(settings:get("grid_x_size"))
if grid_x_size then
settings:set("map_x_size", tostring(grid_x_size * blocksize))
end
local grid_z_size = tonumber(settings:get("grid_z_size"))
if grid_z_size then
settings:set("map_z_size", tostring(grid_z_size * blocksize))
end
end
end
return fix_compatibility_minetest, fix_compatibility_mapgen_rivers

40
geometry.lua Normal file
View File

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

105
gridmanager.lua
View File

@ -1,105 +0,0 @@
-- Manages grid loading, writing and generation
local datapath = mapgen_rivers.world_data_path
local registered_on_grid_loaded = {}
function mapgen_rivers.register_on_grid_loaded(func)
if type(func) == "function" then
registered_on_grid_loaded[#registered_on_grid_loaded+1] = func
else
minetest.log("error", "[mapgen_rivers] register_on_grid_loaded can only register functions!")
end
end
local function on_grid_loaded_callback(grid)
for _, func in ipairs(registered_on_grid_loaded) do
func(grid)
end
end
local function offset_conv(o)
return (o + 0.5) * (1/256)
end
local floor = math.floor
local sbyte, schar = string.byte, string.char
local unpk = unpack
-- Loading files
-- Never load the full map during mapgen. Instead, create an empty lookup table
-- and read the file on-the-fly when an element is requested for the first time,
-- using __index metamethod.
local loader_mt = {
__index = function(loader, i) -- Called when accessing a missing key
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
-- Cache key for next use
loader[i] = n
return n
end,
}
local function load_file(filename, bytes, signed, size, converter)
local file = io.open(datapath .. filename, 'rb')
if file then
converter = converter or false
return setmetatable({file=file, bytes=bytes, signed=signed, size=size, conv=converter}, loader_mt)
end
end
function mapgen_rivers.load_or_generate_grid()
-- First, check whether a grid is already loaded
if mapgen_rivers.grid then
return true
end
-- If not, try to load the grid from the files
local sfile = io.open(datapath .. 'size', 'r')
if not sfile then
dofile(mapgen_rivers.modpath .. "/pregenerate.lua")
collectgarbage()
sfile = io.open(datapath .. 'size', 'r')
if not sfile then
return false
end
end
local x, z = sfile:read('*n'), sfile:read('*n')
if not x or not z then
return false
end
minetest.log("action", '[mapgen_rivers] Loading grid')
local grid = {
size = {x=x, y=z},
dem = load_file('dem', 2, true, x*z),
lakes = load_file('lakes', 2, true, x*z),
dirs = load_file('dirs', 1, false, x*z),
rivers = load_file('rivers', 4, false, x*z),
offset_x = load_file('offset_x', 1, true, x*z, offset_conv),
offset_y = load_file('offset_y', 1, true, x*z, offset_conv),
}
mapgen_rivers.grid = grid
on_grid_loaded_callback(grid)
return true
end

145
heightmap.lua
View File

@ -1,53 +1,17 @@
-- Transform polygon data into a heightmap
local modpath = mapgen_rivers.modpath
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 out_elev = mapgen_rivers.settings.margin_elev
local MAP_BOTTOM = -31000
-- Localize for performance
local floor, min, max, sqrt, abs = math.floor, math.min, math.max, math.sqrt, math.abs
local floor, min, max = math.floor, math.min, math.max
local unpk = unpack
-- Geometrical helpers
local function distance_to_segment(x1, y1, x2, y2, x, y)
-- get the distance between point (x,y) and segment (x1,y1)-(x2,y2)
local a = (x1-x2)^2 + (y1-y2)^2 -- square of distance
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 sqrt(b)
elseif a + c < b then
-- The closest point of the segment is the extremity 2
return sqrt(c)
else
-- The closest point is on the segment
return abs(x1 * (y2-y) + x2 * (y-y1) + x * (y1-y2)) / sqrt(a)
end
end
local function transform_quadri(X, Y, x, y)
-- To index points in an irregular quadrilateral, giving x and y between 0 (one edge) and 1 (opposite edge)
-- X, Y 4-vectors giving the coordinates of the 4 vertices
-- x, y position to index.
local x1, x2, x3, x4 = unpk(X)
local y1, y2, y3, y4 = unpk(Y)
-- Compare distance to 2 opposite edges, they give the X coordinate
local d23 = distance_to_segment(x2,y2,x3,y3,x,y)
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
-- Linear interpolation
local function interp(v00, v01, v11, v10, xf, zf)
local v0 = v01*xf + v00*(1-xf)
@ -55,9 +19,9 @@ local function interp(v00, v01, v11, v10, xf, zf)
return v1*zf + v0*(1-zf)
end
function mapgen_rivers.make_heightmaps(minp, maxp)
local function heightmaps(minp, maxp)
local polygons = mapgen_rivers.make_polygons(minp, maxp)
local polygons = make_polygons(minp, maxp)
local incr = maxp.z-minp.z+1
local terrain_height_map = {}
@ -82,57 +46,50 @@ function mapgen_rivers.make_heightmaps(minp, maxp)
-- = 0: on riverbank
-- > 0: inside river
local depth_factors = {
r_west - xf , -- West edge (1)
r_north - zf , -- North edge (2)
r_east - (1-xf), -- East edge (3)
r_south - (1-zf), -- South edge (4)
c_NW - xf - zf , -- North-West corner (5)
c_NE - (1-xf) - zf , -- North-East corner (6)
c_SE - (1-xf) - (1-zf), -- South-East corner (7)
c_SW - xf - (1-zf), -- South-West corner (8)
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, which determines to which of the 8 river sections (4 edges + 4 corners) the current point belongs.
-- If imax is still at 0, it means that we are not in a river.
local dpmax = 0
-- 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] > dpmax then
dpmax = depth_factors[i]
if depth_factors[i] >= depth_factor_max then
depth_factor_max = depth_factors[i]
imax = i
end
end
-- Transform the coordinates to have xfc and zfc = 0 or 1 in rivers (to avoid rivers having lateral slope and to accomodate the riverbanks smoothly)
local xfc, zfc
-- xfc:
if imax == 0 or imax == 2 or imax == 4 then -- river segment does not constrain X coordinate, so accomodate xf in function of other river sections
local x0 = max(r_west-dpmax, c_NW-zf-dpmax, c_SW-(1-zf)-dpmax, 0) -- new xf will be bounded to 0 by western riverbank
local x1 = 1-max(r_east-dpmax, c_NE-zf-dpmax, c_SE-(1-zf)-dpmax, 0) -- and bounded to 1 by eastern riverbank
if x0 >= x1 then
xfc = 0.5
else
xfc = (xf-x0) / (x1-x0)
end
elseif imax == 1 or imax == 5 or imax == 8 then -- river at the western side of the polygon
xfc = 0
else -- 3, 6, 7 : river at the eastern side of the polygon
xfc = 1
end
-- Same for zfc:
if imax == 0 or imax == 1 or imax == 3 then -- river segment does not constrain Z coordinate, so accomodate zf in function of other river sections
local z0 = max(r_north-dpmax, c_NW-xf-dpmax, c_NE-(1-xf)-dpmax, 0) -- new zf will be bounded to 0 by northern riverbank
local z1 = 1-max(r_south-dpmax, c_SW-xf-dpmax, c_SE-(1-xf)-dpmax, 0) -- and bounded to 1 by southern riverbank
if z0 >= z1 then
zfc = 0.5
else
zfc = (zf-z0) / (z1-z0)
end
elseif imax == 2 or imax == 5 or imax == 6 then -- river at the northern side of the polygon
zfc = 0
else -- 4, 7, 8 : river at the southern side of the polygon
zfc = 1
-- 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 = max(r_west, c_NW-zf, zf-c_SW)
local x1 = min(r_east, c_NE+zf, c_SE-zf)
local z0 = max(r_north, c_NW-xf, xf-c_NE)
local z1 = min(r_south, c_SW+xf, c_SE-xf)
xf = (xf-x0) / (x1-x0)
zf = (zf-z0) / (z1-z0)
elseif imax == 1 then
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
@ -142,12 +99,12 @@ function mapgen_rivers.make_heightmaps(minp, maxp)
vdem[2],
vdem[3],
vdem[4],
xfc, zfc
xf, zf
))
-- Spatial gradient of the interpolation
local slope_x = zfc*(vdem[3]-vdem[4]) + (1-zfc)*(vdem[2]-vdem[1]) < 0
local slope_z = xfc*(vdem[3]-vdem[2]) + (1-xfc)*(vdem[4]-vdem[1]) < 0
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
@ -164,15 +121,15 @@ function mapgen_rivers.make_heightmaps(minp, maxp)
end
local lake_height = max(floor(poly.lake[lake_id]), terrain_height)
if imax > 0 and dpmax > 0 then
terrain_height = min(max(lake_height, sea_level) - floor(1+dpmax*riverbed_slope), terrain_height)
if imax > 0 and depth_factor_max > 0 then
terrain_height = min(max(lake_height, sea_level) - 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] = out_elev
lake_height_map[i] = out_elev
terrain_height_map[i] = MAP_BOTTOM
lake_height_map[i] = MAP_BOTTOM
end
i = i + 1
end
@ -180,3 +137,5 @@ function mapgen_rivers.make_heightmaps(minp, maxp)
return terrain_height_map, lake_height_map
end
return heightmaps

277
init.lua
View File

@ -1,16 +1,275 @@
-- Main file, calls the other files and triggers main functions
mapgen_rivers = {}
local modpath = minetest.get_modpath(minetest.get_current_modname()) .. '/'
mapgen_rivers.modpath = modpath
mapgen_rivers.world_data_path = minetest.get_worldpath() .. '/river_data/'
dofile(modpath .. 'settings.lua')
dofile(modpath .. 'gridmanager.lua')
dofile(modpath .. 'polygons.lua')
dofile(modpath .. 'heightmap.lua')
dofile(modpath .. 'mapgen.lua')
dofile(modpath .. 'spawn.lua')
if minetest.get_mapgen_setting("mg_name") ~= "singlenode" then
minetest.set_mapgen_setting("mg_name", "singlenode", true)
minetest.log("warning", "[mapgen_rivers] Mapgen set to singlenode")
end
minetest.register_on_mods_loaded(mapgen_rivers.load_or_generate_grid)
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 use_biomegen_mod = use_biomes and minetest.global_exists('biomegen')
use_biomes = use_biomes and not use_biomegen_mod
if use_biomegen_mod then
biomegen.set_elevation_chill(elevation_chill)
end
dofile(modpath .. 'noises.lua')
local heightmaps = dofile(modpath .. 'heightmap.lua')
-- Linear interpolation
local function interp(v00, v01, v11, v10, xf, zf)
local v0 = v01*xf + v00*(1-xf)
local v1 = v11*xf + v10*(1-xf)
return v1*zf + v0*(1-zf)
end
-- Localize for performance
local floor, min = math.floor, math.min
local data = {}
local noise_x_obj, noise_z_obj, noise_distort_obj, noise_heat_obj, noise_heat_blend_obj
local noise_x_map = {}
local noise_z_map = {}
local noise_distort_map = {}
local noise_heat_map = {}
local noise_heat_blend_map = {}
local mapsize
local init = false
local sumtime = 0
local sumtime2 = 0
local ngen = 0
local function generate(minp, maxp, seed)
minetest.log("info", ("[mapgen_rivers] Generating from %s to %s"):format(minetest.pos_to_string(minp), minetest.pos_to_string(maxp)))
local chulens = {
x = maxp.x-minp.x+1,
y = maxp.y-minp.y+1,
z = maxp.z-minp.z+1,
}
if not init then
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 t0 = os.clock()
local minp2d = {x=minp.x, y=minp.z}
if use_distort then
noise_x_obj:get_3d_map_flat(minp, noise_x_map)
noise_z_obj:get_3d_map_flat(minp, noise_z_map)
noise_distort_obj:get_2d_map_flat(minp2d, noise_distort_map)
end
if use_biomes then
noise_heat_obj:get_2d_map_flat(minp2d, noise_heat_map)
noise_heat_blend_obj:get_2d_map_flat(minp2d, noise_heat_blend_map)
end
local terrain_map, lake_map, incr, i_origin
if use_distort then
local xmin, xmax, zmin, zmax = minp.x, maxp.x, minp.z, maxp.z
local i = 0
local i2d = 0
for z=minp.z, maxp.z do
for y=minp.y, maxp.y+1 do
for x=minp.x, maxp.x do
i = i+1
i2d = i2d+1
local distort = noise_distort_map[i2d]
local xv = noise_x_map[i]*distort + x
if xv < xmin then xmin = xv end
if xv > xmax then xmax = xv end
noise_x_map[i] = xv
local zv = noise_z_map[i]*distort + z
if zv < zmin then zmin = zv end
if zv > zmax then zmax = zv end
noise_z_map[i] = zv
end
i2d = i2d-chulens.x
end
end
local pminp = {x=floor(xmin), z=floor(zmin)}
local pmaxp = {x=floor(xmax)+1, z=floor(zmax)+1}
incr = pmaxp.x-pminp.x+1
i_origin = 1 - pminp.z*incr - pminp.x
terrain_map, lake_map = heightmaps(pminp, pmaxp)
else
terrain_map, lake_map = heightmaps(minp, maxp)
end
-- Check that there is at least one position that reaches min y
if minp.y > sea_level then
local y0 = minp.y
local is_empty = true
for i=1, #terrain_map do
if terrain_map[i] >= y0 or lake_map[i] >= y0 then
is_empty = false
break
end
end
-- If not, skip chunk
if is_empty then
local t = os.clock() - t0
ngen = ngen + 1
sumtime = sumtime + t
sumtime2 = sumtime2 + t*t
minetest.log("verbose", "[mapgen_rivers] Skipping empty chunk (fully above ground level)")
minetest.log("verbose", ("[mapgen_rivers] Done in %5.3f s"):format(t))
return
end
end
local c_stone = minetest.get_content_id("default:stone")
local c_dirt = minetest.get_content_id("default:dirt")
local c_lawn = minetest.get_content_id("default:dirt_with_grass")
local c_dirtsnow = minetest.get_content_id("default:dirt_with_snow")
local c_snow = minetest.get_content_id("default:snowblock")
local c_sand = minetest.get_content_id("default:sand")
local c_water = minetest.get_content_id("default:water_source")
local c_rwater = minetest.get_content_id("default:river_water_source")
local c_ice = minetest.get_content_id("default:ice")
local vm, emin, emax = minetest.get_mapgen_object("voxelmanip")
vm:get_data(data)
local a = VoxelArea:new({MinEdge = emin, MaxEdge = emax})
local ystride = a.ystride -- Tip : the ystride of a VoxelArea is the number to add to the array index to get the index of the position above. It's faster because it avoids to completely recalculate the index.
local nid = mapsize.x*(mapsize.y-1) + 1
local incrY = -mapsize.x
local incrX = 1 - mapsize.y*incrY
local incrZ = mapsize.x*mapsize.y - mapsize.x*incrX - mapsize.x*mapsize.y*incrY
local i2d = 1
for z = minp.z, maxp.z do
for x = minp.x, maxp.x do
local ivm = a:index(x, maxp.y+1, z)
local ground_above = false
local temperature
if use_biomes then
temperature = noise_heat_map[i2d]+noise_heat_blend_map[i2d]
end
local terrain, lake
if not use_distort then
terrain = terrain_map[i2d]
lake = lake_map[i2d]
end
for y = maxp.y+1, minp.y, -1 do
if use_distort then
local xn = noise_x_map[nid]
local zn = noise_z_map[nid]
local x0 = floor(xn)
local z0 = floor(zn)
local i0 = i_origin + z0*incr + x0
local i1 = i0+1
local i2 = i1+incr
local i3 = i2-1
terrain = interp(terrain_map[i0], terrain_map[i1], terrain_map[i2], terrain_map[i3], xn-x0, zn-z0)
lake = min(lake_map[i0], lake_map[i1], lake_map[i2], lake_map[i3])
end
if y <= maxp.y then
local is_lake = lake > terrain
if y <= terrain then
if not use_biomes or y <= terrain-1 or ground_above then
data[ivm] = c_stone
elseif is_lake or y < sea_level then
data[ivm] = c_sand
else
local temperature_y = temperature - y*elevation_chill
if temperature_y >= 15 then
data[ivm] = c_lawn
elseif temperature_y >= 0 then
data[ivm] = c_dirtsnow
else
data[ivm] = c_snow
end
end
elseif y <= lake and lake > sea_level then
if not use_biomes or temperature - y*elevation_chill >= 0 then
data[ivm] = c_rwater
else
data[ivm] = c_ice
end
elseif y <= sea_level then
data[ivm] = c_water
end
end
ground_above = y <= terrain
ivm = ivm - ystride
if use_distort then
nid = nid + incrY
end
end
if use_distort then
nid = nid + incrX
end
i2d = i2d + 1
end
if use_distort then
nid = nid + incrZ
end
end
if use_biomegen_mod then
biomegen.generate_all(data, a, vm, minp, maxp, seed)
else
vm:set_data(data)
minetest.generate_ores(vm, minp, maxp)
end
vm:set_lighting({day = 0, night = 0})
vm:calc_lighting()
vm:update_liquids()
vm:write_to_map()
local t = os.clock()-t0
ngen = ngen + 1
sumtime = sumtime + t
sumtime2 = sumtime2 + t*t
minetest.log("verbose", ("[mapgen_rivers] Done in %5.3f s"):format(t))
end
minetest.register_on_generated(generate)
minetest.register_on_shutdown(function()
local avg = sumtime / ngen
local std = math.sqrt(sumtime2/ngen - avg*avg)
minetest.log("action", ("[mapgen_rivers] Mapgen statistics:\n- Mapgen calls: %4d\n- Mean time: %5.3f s\n- Standard deviation: %5.3f s"):format(ngen, avg, std))
end)

98
load.lua Normal file
View File

@ -0,0 +1,98 @@
local worldpath = mapgen_rivers.world_data_path
local floor = math.floor
local sbyte, schar = string.byte, string.char
local unpk = unpack
function mapgen_rivers.load_map(filename, bytes, signed, size, converter)
local file = io.open(worldpath .. filename, 'rb')
local data = file:read('*all')
if #data < bytes*size then
data = minetest.decompress(data)
end
local map = {}
for i=1, size do
local i0 = (i-1)*bytes+1
local elements = {data:byte(i0, i1)}
local n = sbyte(data, i0)
if signed and n >= 128 then
n = n - 256
end
for j=1, bytes-1 do
n = n*256 + sbyte(data, i0+j)
end
map[i] = n
end
file:close()
if converter then
for i=1, size do
map[i] = converter(map[i])
end
end
return map
end
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 bytelist = {}
for j=1, bytes do
bytelist[j] = 0
end
for i=1, size do
local n = floor(data[i])
data[i] = n
for j=bytes, 2, -1 do
bytelist[j] = n % 256
n = floor(n / 256)
end
bytelist[1] = n % 256
file:write(schar(unpk(bytelist)))
end
file:close()
end

278
mapgen.lua
View File

@ -1,278 +0,0 @@
-- Mapgen loop and mapgen-related things
if minetest.get_mapgen_setting("mg_name") ~= "singlenode" then
minetest.set_mapgen_setting("mg_name", "singlenode", true)
minetest.log("warning", "[mapgen_rivers] Mapgen set to singlenode")
end
local sea_level = mapgen_rivers.settings.sea_level
local elevation_chill = mapgen_rivers.settings.elevation_chill
local use_distort = mapgen_rivers.settings.distort
local use_biomes = mapgen_rivers.settings.biomes
local use_biomegen_mod = use_biomes and minetest.global_exists('biomegen')
use_biomes = use_biomes and minetest.global_exists('default') and not use_biomegen_mod
if use_biomegen_mod then
biomegen.set_elevation_chill(elevation_chill)
end
-- Linear interpolation
local function interp(v00, v01, v11, v10, xf, zf)
local v0 = v01*xf + v00*(1-xf)
local v1 = v11*xf + v10*(1-xf)
return v1*zf + v0*(1-zf)
end
-- Localize for performance
local floor, min = math.floor, math.min
local data = {}
local noise_x_obj, noise_z_obj, noise_distort_obj, noise_heat_obj, noise_heat_blend_obj
local noise_x_map = {}
local noise_z_map = {}
local noise_distort_map = {}
local noise_heat_map = {}
local noise_heat_blend_map = {}
local mapsize
local init = false
local sumtime = 0
local sumtime2 = 0
local ngen = 0
function mapgen_rivers.make_chunk(minp, maxp, seed)
minetest.log("info", ("[mapgen_rivers] Generating from %s to %s"):format(minetest.pos_to_string(minp), minetest.pos_to_string(maxp)))
local chulens = {
x = maxp.x-minp.x+1,
y = maxp.y-minp.y+1,
z = maxp.z-minp.z+1,
}
if not init then
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 t0 = os.clock()
local minp2d = {x=minp.x, y=minp.z}
if use_distort then
noise_x_obj:get_3d_map_flat(minp, noise_x_map)
noise_z_obj:get_3d_map_flat(minp, noise_z_map)
noise_distort_obj:get_2d_map_flat(minp2d, noise_distort_map)
end
if use_biomes then
noise_heat_obj:get_2d_map_flat(minp2d, noise_heat_map)
noise_heat_blend_obj:get_2d_map_flat(minp2d, noise_heat_blend_map)
end
local terrain_map, lake_map, incr, i_origin
if use_distort then
local xmin, xmax, zmin, zmax = minp.x, maxp.x, minp.z, maxp.z
local i = 0
local i2d = 0
for z=minp.z, maxp.z do
for y=minp.y, maxp.y+1 do
for x=minp.x, maxp.x do
i = i+1
i2d = i2d+1
local distort = noise_distort_map[i2d]
local xv = noise_x_map[i]*distort + x
if xv < xmin then xmin = xv end
if xv > xmax then xmax = xv end
noise_x_map[i] = xv
local zv = noise_z_map[i]*distort + z
if zv < zmin then zmin = zv end
if zv > zmax then zmax = zv end
noise_z_map[i] = zv
end
i2d = i2d-chulens.x
end
i2d = i2d+chulens.x
end
local pminp = {x=floor(xmin), z=floor(zmin)}
local pmaxp = {x=floor(xmax)+1, z=floor(zmax)+1}
incr = pmaxp.x-pminp.x+1
i_origin = 1 - pminp.z*incr - pminp.x
terrain_map, lake_map = mapgen_rivers.make_heightmaps(pminp, pmaxp)
else
terrain_map, lake_map = mapgen_rivers.make_heightmaps(minp, maxp)
end
-- Check that there is at least one position that reaches min y
if minp.y > sea_level then
local y0 = minp.y
local is_empty = true
for i=1, #terrain_map do
if terrain_map[i] >= y0 or lake_map[i] >= y0 then
is_empty = false
break
end
end
-- If not, skip chunk
if is_empty then
if use_biomegen_mod and biomegen.skip_chunk then
biomegen.skip_chunk(minp, maxp)
end
local t = os.clock() - t0
ngen = ngen + 1
sumtime = sumtime + t
sumtime2 = sumtime2 + t*t
minetest.log("verbose", "[mapgen_rivers] Skipping empty chunk (fully above ground level)")
minetest.log("verbose", ("[mapgen_rivers] Done in %5.3f s"):format(t))
return
end
end
local c_stone = minetest.get_content_id("mapgen_stone")
local c_water = minetest.get_content_id("mapgen_water_source")
local c_rwater = minetest.get_content_id("mapgen_river_water_source")
local c_dirt, c_lawn, c_dirtsnow, c_snow, c_sand, c_ice
if use_biomes then
c_dirt = minetest.get_content_id("default:dirt")
c_lawn = minetest.get_content_id("default:dirt_with_grass")
c_dirtsnow = minetest.get_content_id("default:dirt_with_snow")
c_snow = minetest.get_content_id("default:snowblock")
c_sand = minetest.get_content_id("default:sand")
c_ice = minetest.get_content_id("default:ice")
end
local vm, emin, emax = minetest.get_mapgen_object("voxelmanip")
vm:get_data(data)
local a = VoxelArea:new({MinEdge = emin, MaxEdge = emax})
local ystride = a.ystride -- Tip : the ystride of a VoxelArea is the number to add to the array index to get the index of the position above. It's faster because it avoids to completely recalculate the index.
local nid = mapsize.x*(mapsize.y-1) + 1
local incrY = -mapsize.x
local incrX = 1 - mapsize.y*incrY
local incrZ = mapsize.x*mapsize.y - mapsize.x*incrX - mapsize.x*mapsize.y*incrY
local i2d = 1
for z = minp.z, maxp.z do
for x = minp.x, maxp.x do
local ivm = a:index(x, maxp.y+1, z)
local ground_above = false
local temperature
if use_biomes then
temperature = noise_heat_map[i2d]+noise_heat_blend_map[i2d]
end
local terrain, lake
if not use_distort then
terrain = terrain_map[i2d]
lake = lake_map[i2d]
end
for y = maxp.y+1, minp.y, -1 do
if use_distort then
local xn = noise_x_map[nid]
local zn = noise_z_map[nid]
local x0 = floor(xn)
local z0 = floor(zn)
local i0 = i_origin + z0*incr + x0
local i1 = i0+1
local i2 = i1+incr
local i3 = i2-1
terrain = interp(terrain_map[i0], terrain_map[i1], terrain_map[i2], terrain_map[i3], xn-x0, zn-z0)
lake = min(lake_map[i0], lake_map[i1], lake_map[i2], lake_map[i3])
end
if y <= maxp.y then
local is_lake = lake > terrain
if y <= terrain then
if not use_biomes or y <= terrain-1 or ground_above then
data[ivm] = c_stone
elseif is_lake or y < sea_level then
data[ivm] = c_sand
else
local temperature_y = temperature - y*elevation_chill
if temperature_y >= 15 then
data[ivm] = c_lawn
elseif temperature_y >= 0 then
data[ivm] = c_dirtsnow
else
data[ivm] = c_snow
end
end
elseif y <= lake and lake > sea_level then
if not use_biomes or temperature - y*elevation_chill >= 0 then
data[ivm] = c_rwater
else
data[ivm] = c_ice
end
elseif y <= sea_level then
data[ivm] = c_water
end
end
ground_above = y <= terrain
ivm = ivm - ystride
if use_distort then
nid = nid + incrY
end
end
if use_distort then
nid = nid + incrX
end
i2d = i2d + 1
end
if use_distort then
nid = nid + incrZ
end
end
if use_biomegen_mod then
biomegen.generate_all(data, a, vm, minp, maxp, seed)
else
vm:set_data(data)
minetest.generate_ores(vm, minp, maxp)
end
vm:set_lighting({day = 0, night = 0})
vm:calc_lighting()
vm:update_liquids()
vm:write_to_map()
local t = os.clock()-t0
ngen = ngen + 1
sumtime = sumtime + t
sumtime2 = sumtime2 + t*t
minetest.log("verbose", ("[mapgen_rivers] Done in %5.3f s"):format(t))
end
-- Enforce first position in mapgen callbacks
table.insert(minetest.registered_on_generateds, 1, mapgen_rivers.make_chunk)
--minetest.register_on_generated(mapgen_rivers.make_chunk)
minetest.register_on_shutdown(function()
local avg = sumtime / ngen
local std = math.sqrt(sumtime2/ngen - avg*avg)
minetest.log("action", ("[mapgen_rivers] Mapgen statistics:\n- Mapgen calls: %4d\n- Mean time: %5.3f s\n- Standard deviation: %5.3f s"):format(ngen, avg, std))
end)

3
mod.conf
View File

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

80
noises.lua Normal file
View File

@ -0,0 +1,80 @@
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', {
offset = 0,
scale = 1,
seed = -4574,
spread = {x=64, y=32, z=64},
octaves = 3,
persistence = 0.75,
lacunarity = 2,
}),
distort_z = def_setting('np_distort_z', 'noise', {
offset = 0,
scale = 1,
seed = -7940,
spread = {x=64, y=32, z=64},
octaves = 3,
persistence = 0.75,
lacunarity = 2,
}),
distort_amplitude = def_setting('np_distort_amplitude', 'noise', {
offset = 0,
scale = 10,
seed = 676,
spread = {x=1024, y=1024, z=1024},
octaves = 5,
persistence = 0.5,
lacunarity = 2,
flags = "absvalue",
}),
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
minetest.log("warning", "[mapgen_rivers] Noise " .. name .. ": 'octaves' reduced to " .. omax)
np.octaves = omax
end
end
end

107
polygons.lua
View File

@ -1,23 +1,94 @@
-- Fetch polygons from a given areas, and compute their properties
-- and find to which polygon every point belongs
local modpath = mapgen_rivers.modpath
local mod_data_path = modpath .. 'river_data/'
if not io.open(mod_data_path .. 'size', 'r') then
mod_data_path = modpath .. 'demo_data/'
end
local blocksize = mapgen_rivers.settings.blocksize
local X = math.floor(mapgen_rivers.settings.map_x_size / blocksize)
local Z = math.floor(mapgen_rivers.settings.map_z_size / blocksize)
local world_data_path = mapgen_rivers.world_data_path
minetest.mkdir(world_data_path)
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()
minetest.log("action", '[mapgen_rivers] Generating grid, this may take a while...')
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
end
minetest.register_on_mods_loaded(function()
if load_all then
minetest.log("action", '[mapgen_rivers] Loading full grid')
else
minetest.log("action", '[mapgen_rivers] Loading grid as interactive loaders')
end
local grid = mapgen_rivers.grid
grid.dem = load_map('dem', 2, true, X*Z)
grid.lakes = load_map('lakes', 2, true, X*Z)
grid.dirs = load_map('dirs', 1, false, X*Z)
grid.rivers = load_map('rivers', 4, false, X*Z)
grid.offset_x = load_map('offset_x', 1, true, X*Z, offset_converter)
grid.offset_y = load_map('offset_y', 1, true, X*Z, offset_converter)
end)
end
mapgen_rivers.grid.size = {x=X, y=Z}
local function index(x, z)
return z*X+x+1
end
local blocksize = mapgen_rivers.settings.blocksize
local min_catchment = mapgen_rivers.settings.min_catchment
local max_catchment = mapgen_rivers.settings.max_catchment
local map_offset = {x=0, z=0}
mapgen_rivers.register_on_grid_loaded(function(grid)
X = grid.size.x
Z = grid.size.y
if mapgen_rivers.settings.center then
map_offset.x = blocksize*X/2
map_offset.z = blocksize*Z/2
end
end)
if mapgen_rivers.settings.center then
map_offset.x = blocksize*X/2
map_offset.z = blocksize*Z/2
end
-- Localize for performance
local floor, ceil, min, max, abs = math.floor, math.ceil, math.min, math.max, math.abs
@ -48,7 +119,7 @@ local init = false
-- On map generation, determine into which polygon every point (in 2D) will fall.
-- Also store polygon-specific data
function mapgen_rivers.make_polygons(minp, maxp)
local function make_polygons(minp, maxp)
local grid = mapgen_rivers.grid
local dem = grid.dem
@ -165,15 +236,15 @@ function mapgen_rivers.make_polygons(minp, maxp)
end
end
polygon.river_corners = {riverA, riverB, riverC, riverD}
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 = (dirB==1 and riverB or 0) + (dirC==3 and riverC or 0)
local river_south = (dirD==2 and riverD or 0) + (dirC==4 and riverC 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
@ -181,3 +252,5 @@ function mapgen_rivers.make_polygons(minp, maxp)
return polygons
end
return make_polygons

175
pregenerate.lua
View File

@ -1,8 +1,3 @@
-- Generate the grid using terrainlib_lua
-- Only called on first mapgen, if there is no grid yet
-- Constants
local EvolutionModel = dofile(mapgen_rivers.modpath .. '/terrainlib_lua/erosion.lua')
local twist = dofile(mapgen_rivers.modpath .. '/terrainlib_lua/twist.lua')
@ -10,7 +5,6 @@ local blocksize = mapgen_rivers.settings.blocksize
local tectonic_speed = mapgen_rivers.settings.tectonic_speed
local np_base = table.copy(mapgen_rivers.noise_params.base)
np_base.spread = vector.divide(np_base.spread, blocksize)
local evol_params = mapgen_rivers.settings.evol_params
@ -19,132 +13,69 @@ local time_step = mapgen_rivers.settings.evol_time_step
local niter = math.ceil(time/time_step)
time_step = time / niter
local use_margin = mapgen_rivers.settings.margin
local margin_width = mapgen_rivers.settings.margin_width / blocksize
local margin_elev = mapgen_rivers.settings.margin_elev
local function pregenerate(keep_loaded)
local grid = mapgen_rivers.grid
local size = grid.size
local X = math.floor(mapgen_rivers.settings.map_x_size / blocksize)
local Y = math.floor(mapgen_rivers.settings.map_z_size / blocksize)
local seed = tonumber(minetest.get_mapgen_setting("seed"))
np_base.seed = (np_base.seed or 0) + seed
local function margin(dem, width, elev)
local X, Y = dem.X, dem.Y
for i=1, width do
local c1 = ((i-1)/width) ^ 0.5
local c2 = (1-c1) * elev
local index = (i-1)*X + 1
for x=1, X do
dem[index] = dem[index] * c1 + c2
index = index + 1
end
index = i
for y=1, Y do
dem[index] = dem[index] * c1 + c2
index = index + X
end
index = X*(Y-i) + 1
for x=1, X do
dem[index] = dem[index] * c1 + c2
index = index + 1
end
index = X-i + 1
for y=1, Y do
dem[index] = dem[index] * c1 + c2
index = index + X
end
end
end
local nobj_base = PerlinNoiseMap(np_base, {x=size.x, y=1, z=size.y})
-- Generate grid
local dem = nobj_base:get_3d_map_flat({x=0, y=0, z=0})
dem.X = size.x
dem.Y = size.y
minetest.log("action", '[mapgen_rivers] Generating grid, this may take a while...')
local model = EvolutionModel(evol_params)
model.dem = dem
local ref_dem = model:define_isostasy(dem)
if X*Y > 4e6 then
minetest.log("warning", "[mapgen_rivers] You are going to generate a very large grid (>4M nodes). If you experience problems, you should increase blocksize or reduce map size.")
end
local tectonic_step = tectonic_speed * time_step
collectgarbage()
for i=1, niter do
minetest.log("info", "[mapgen_rivers] Iteration " .. i .. " of " .. niter)
local seed = tonumber(minetest.get_mapgen_setting("seed"):sub(-10))
np_base.seed = (np_base.seed or 0) + seed
local nobj_base = PerlinNoiseMap(np_base, {x=X, y=1, z=Y})
local dem = nobj_base:get_3d_map_flat({x=0, y=0, z=0})
dem.X = X
dem.Y = Y
if use_margin then
margin(dem, margin_width, margin_elev)
end
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
minetest.log("info", "[mapgen_rivers] Iteration " .. i .. " of " .. niter)
model:diffuse(time_step)
model:flow()
model:erode(time_step)
if i < niter then
if tectonic_step ~= 0 then
nobj_base:get_3d_map_flat({x=0, y=tectonic_step*i, z=0}, ref_dem)
if use_margin then
margin(ref_dem, margin_width, margin_elev)
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
model:isostasy()
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
model:flow()
local mfloor = math.floor
local mmin, mmax = math.min, math.max
local unpk, schar = unpack, string.char
local offset_x, offset_y = twist(model.dirs, model.rivers, 5)
for i=1, X*Y do
offset_x[i] = mmin(mmax(offset_x[i]*256, -128), 127)
offset_y[i] = mmin(mmax(offset_y[i]*256, -128), 127)
end
-- Write data
local datapath = mapgen_rivers.world_data_path
minetest.mkdir(datapath)
local sfile = io.open(datapath .. 'size', "w")
sfile:write(X..'\n'..Y)
sfile:close()
local function write_file(filename, data, bytes)
local file = io.open(datapath .. filename, 'wb')
local bytelist = {}
for j=1, bytes do
bytelist[j] = 0
end
for i=1, #data 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(unpk(bytelist)))
end
file:close()
end
write_file('dem', model.dem, 2)
write_file('lakes', model.lakes, 2)
write_file('dirs', model.dirs, 1)
write_file('rivers', model.rivers, 4)
write_file('offset_x', offset_x, 1)
write_file('offset_y', offset_y, 1)
return pregenerate

87
settings.lua
View File

@ -1,9 +1,7 @@
-- Read global and per-world settings
local mtsettings = minetest.settings
local mgrsettings = Settings(minetest.get_worldpath() .. '/mapgen_rivers.conf')
mapgen_rivers.version = "1.0.2-dev1"
mapgen_rivers.version = "1.0.2"
local previous_version_mt = mtsettings:get("mapgen_rivers_version") or "0.0"
local previous_version_mgr = mgrsettings:get("version") or "0.0"
@ -21,33 +19,13 @@ end
mtsettings:set("mapgen_rivers_version", mapgen_rivers.version)
mgrsettings:set("version", mapgen_rivers.version)
local defaults
do
local f = io.open(mapgen_rivers.modpath .. "/settings_default.json")
defaults = minetest.parse_json(f:read("*all"))
f:close()
end
-- Convert strings to numbers in noise params because Minetest API is not able to do it cleanly...
local function clean_np(np)
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
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 = defaults[name]
v = default
end
mgrsettings:set(name, v)
end
@ -61,7 +39,7 @@ function mapgen_rivers.define_setting(name, dtype, default)
if v == nil then
v = mtsettings:get_bool('mapgen_rivers_' .. name)
if v == nil then
v = defaults[name]
v = default
end
mgrsettings:set_bool(name, v)
end
@ -71,11 +49,10 @@ function mapgen_rivers.define_setting(name, dtype, default)
if v == nil then
v = mtsettings:get_np_group('mapgen_rivers_' .. name)
if v == nil then
v = defaults[name]
v = default
end
mgrsettings:set_np_group(name, v)
end
clean_np(v)
return v
end
end
@ -83,47 +60,33 @@ end
local def_setting = mapgen_rivers.define_setting
mapgen_rivers.settings = {
center = def_setting('center', 'bool'),
blocksize = def_setting('blocksize', 'number'),
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'),
river_widening_power = def_setting('river_widening_power', 'number'),
riverbed_slope = def_setting('riverbed_slope', 'number'),
distort = def_setting('distort', 'bool'),
biomes = def_setting('biomes', 'bool'),
glaciers = def_setting('glaciers', 'bool'),
glacier_factor = def_setting('glacier_factor', 'number'),
elevation_chill = def_setting('elevation_chill', 'number'),
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),
map_x_size = def_setting('map_x_size', 'number'),
map_z_size = def_setting('map_z_size', 'number'),
margin = def_setting('margin', 'bool'),
margin_width = def_setting('margin_width', 'number'),
margin_elev = def_setting('margin_elev', 'number'),
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'),
m = def_setting('river_erosion_power', 'number'),
d = def_setting('diffusive_erosion', 'number'),
compensation_radius = def_setting('compensation_radius', 'number'),
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'),
evol_time = def_setting('evol_time', 'number'),
evol_time_step = def_setting('evol_time_step', 'number'),
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')
}
mapgen_rivers.noise_params = {
base = def_setting("np_base", "noise"),
distort_x = def_setting("np_distort_x", "noise"),
distort_z = def_setting("np_distort_z", "noise"),
distort_amplitude = def_setting("np_distort_amplitude", "noise"),
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.settings.sea_level * mapgen_rivers.settings.elevation_chill
local function write_settings()
mgrsettings:write()
end

70
settings_default.json
View File

@ -1,70 +0,0 @@
{
"version": "1.0.2-dev1",
"center": true,
"water_level": 1,
"blocksize": 15,
"min_catchment": 3600,
"river_widening_power": 0.5,
"riverbed_slope": 0.4,
"distort": true,
"biomes": true,
"glaciers": false,
"glacier_factor": 8,
"elevation_chill": 0.25,
"map_x_size": 15000,
"map_z_size": 15000,
"margin": true,
"margin_width": 2000,
"margin_elev": -200,
"river_erosion_coef": 0.5,
"river_erosion_power": 0.4,
"diffusive_erosion": 0.5,
"compensation_radius": 50,
"tectonic_speed": 70,
"evol_time": 10,
"evol_time_step": 1,
"np_base": {
"offset": 0,
"scale": 300,
"seed": 2469,
"octaves": 8,
"spread": {"x": 2048, "y": 2048, "z": 2048},
"persist": 0.6,
"lacunarity": 2.0,
"flags": "eased"
},
"np_distort_x": {
"offset": 0,
"scale": 1,
"seed": -4574,
"octaves": 3,
"spread": {"x": 64, "y": 32, "z": 64},
"persist": 0.75,
"lacunarity": 2.0
},
"np_distort_z": {
"offset": 0,
"scale": 1,
"seed": -7940,
"octaves": 3,
"spread": {"x": 64, "y": 32, "z": 64},
"persist": 0.75,
"lacunarity": 2.0
},
"np_distort_amplitude": {
"offset": 0,
"scale": 10,
"seed": 676,
"octaves": 5,
"spread": {"x": 1024, "y": 1024, "z": 1024},
"persist": 0.5,
"lacunarity": 2.0,
"flags": "absvalue"
}
}

41
settingtypes.txt
View File

@ -3,35 +3,19 @@
# Whether the map should be centered at x=0, z=0.
mapgen_rivers_center (Center map) bool true
# Every cell of the river grid will represent a square of this size.
# A lower value will result in more detailed terrain and finer computation
# of rivers, but will be slower to generate and use more resources.
#
# WARNING: Excessively low values may cause crashes at pre-generation, due to
# memory issues
# 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 map being generated
#
# X size of the river grid will be this/blocksize
# When increasing, it is recommended to increase blocksize too
mapgen_rivers_map_x_size (Map X size) int 15000 500 66000
# 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 map being generated
#
# Z size of the river grid will be this/blocksize
# When increasing, it is recommended to increase blocksize too
mapgen_rivers_map_z_size (Map Z size) int 15000 500 66000
# If margin is enabled, elevation becomes closer to a fixed value when approaching
# the edges of the map.
mapgen_rivers_margin (Margin) bool true
# Width of the transition at map borders, in nodes
mapgen_rivers_margin_width (Margin width) float 2000.0 0.0 15000.0
# Elevation toward which to converge at map borders
mapgen_rivers_margin_elev (Margin elevation) float -200.0 -31000.0 31000.0
# 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
@ -70,6 +54,11 @@ mapgen_rivers_glacier_widening_factor (Glacier widening factor) float 8.0 1.0 20
# This results in mountains being more covered by snow.
mapgen_rivers_elevation_chill (Elevation chill) float 0.25 0.0 5.0
# 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

119
spawn.lua
View File

@ -1,119 +0,0 @@
local np_distort_x = mapgen_rivers.noise_params.distort_x
local np_distort_z = mapgen_rivers.noise_params.distort_z
local np_distort_amplitude = mapgen_rivers.noise_params.distort_amplitude
local nobj_distort_x, nobj_distort_z, nobj_distort_amplitude
local sea_level = mapgen_rivers.settings.sea_level
local distort = mapgen_rivers.settings.distort
-- 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 estimate_spawn_level(pos, use_distort)
local x, z = pos.x, pos.z
if distort and use_distort then
nobj_distort_x = nobj_distort_x or minetest.get_perlin(np_distort_x)
nobj_distort_z = nobj_distort_z or minetest.get_perlin(np_distort_z)
nobj_distort_amplitude = nobj_distort_amplitude or minetest.get_perlin(np_distort_amplitude)
local amplitude = nobj_distort_amplitude:get_2d({x=pos.x, y=pos.z})
x = x + nobj_distort_x:get_3d(pos)*amplitude
z = z + nobj_distort_z:get_3d(pos)*amplitude
end
local terrain, lakes = mapgen_rivers.make_heightmaps(
{x=math.floor(x), z=math.floor(z) },
{x=math.floor(x)+1, z=math.floor(z)+1}
)
local ex, ez = x % 1, z % 1
--local h = terrain[1]*(1-ex)*(1-ez) + terrain[2]*ex*(1-ez) + terrain[3]*(1-ex)*ez + terrain[4]*ex*ez
local h = interp(terrain[1], terrain[2], terrain[4], terrain[3], ex, ez)
local lake = math.min(lakes[1], lakes[2], lakes[3], lakes[4])
if h < lake or h <= sea_level then
return false, h
end
return true, h
end
local function get_spawn_level(x, z)
local pos = {x=x, z=z}
local suitable, y = estimate_spawn_level(pos, false)
if not suitable then
return
end
if not distort then
return math.floor(y) + 1
end
local low_bound = -math.huge
local high_bound = math.huge
local suitable_high = false
repeat
pos.y = math.max(math.min(math.floor(y+0.5), high_bound-1), low_bound+1)
suitable, y = estimate_spawn_level(pos, true)
if y > pos.y then
low_bound = pos.y
else
high_bound = pos.y
suitable_high = suitable
end
until high_bound - low_bound <= 1
if not suitable_high then
return
end
return high_bound + 1
end
minetest.get_spawn_level = get_spawn_level
local rmax = 2000
local function find_spawn_point(seed)
local level = get_spawn_level(0, 0)
if level then
return {x=0, y=level, z=0}
end
local pr = PcgRandom(seed or os.time())
local incr = 16
local r0 = 0
while r0 < rmax do
local r1 = r0 + incr
local r = pr:next(r0*r0+1, r1*r1) ^ 0.5
local a = pr:next() / 2147483648 * math.pi
local x = math.floor(math.cos(a) * r + 0.5)
local z = math.floor(math.sin(a) * r + 0.5)
level = get_spawn_level(x, z)
if level then
return {x=x, y=level, z=z}
end
r0 = r1
end
end
local function spawn_player(player)
if minetest.settings:get("static_spawnpoint") then
return
end
local spawn_point = find_spawn_point()
if spawn_point then
player:set_pos(spawn_point)
end
end
minetest.register_on_newplayer(spawn_player)
minetest.register_on_respawnplayer(spawn_player)

2
terrainlib_lua/erosion.lua
View File

@ -134,7 +134,7 @@ 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)
model.dirs, model.lakes = rivermapper.flow_routing(model.dem, model.dirs, model.lakes, 'semirandom')
model.rivers = rivermapper.accumulate(model.dirs, model.rivers)
end

336
terrainlib_lua/rivermapper.lua
View File

@ -10,18 +10,50 @@
-- Big thanks to them for releasing this paper under a free license ! :)
-- The algorithm here makes use of most of the paper's concepts, including the Planar Borůvka algorithm.
-- Only flow_local and accumulate_flow are custom algorithms.
local function flow_local_semirandom(plist)
-- Determines how water should flow at 1 node scale.
-- The straightforward approach would be "Water will flow to the lowest of the 4 neighbours", but here water flows to one of the lower neighbours, chosen randomly, but probability depends on height difference.
-- This makes rivers better follow the curvature of the topography at large scale, and be less biased by pure N/E/S/W directions.
-- 'plist': array of downward height differences (0 if upward)
local sum = 0
for i=1, #plist do
sum = sum + plist[i] -- Sum of probabilities
end
if sum == 0 then
return 0
end
local r = math.random() * sum
for i=1, #plist do
local p = plist[i]
if r < p then
return i
end
r = r - p
end
return 0
end
-- Maybe implement more flow methods in the future?
local flow_methods = {
semirandom = flow_local_semirandom,
}
-- Applies all steps of the flow routing, to calculate flow direction for every node, and lake surface elevation.
-- It's quite a hard piece of code, but we will go step by step and explain what's going on, so stay with me and... let's goooooooo!
local function flow_routing(dem, dirs, lakes) -- 'dirs' and 'lakes' are optional tables to reuse for memory optimization, they may contain any data.
local function flow_routing(dem, dirs, lakes, method) -- 'dirs' and 'lakes' are optional tables to reuse for memory optimization, they may contain any data.
method = method or 'semirandom'
local flow_local = flow_methods[method] or flow_local_semirandom
dirs = dirs or {}
lakes = lakes or {}
-- Localize for performance
local tremove = table.remove
local mmax = math.max
local mrand = math.random
local X, Y = dem.X, dem.Y
dirs.X = X
@ -42,29 +74,14 @@ local function flow_routing(dem, dirs, lakes) -- 'dirs' and 'lakes' are optional
for y=1, Y do
for x=1, X do
local zi = dem[i]
-- Determine how water should flow at 1 node scale.
-- The straightforward approach would be "Water will flow to the lowest of the 4 neighbours", but here water flows to one of the lower neighbours, chosen randomly, with probability depending on height difference.
-- This makes rivers better follow the curvature of the topography at large scale, and be less biased by pure N/E/S/W directions.
local pSouth = y<Y and mmax(zi-dem[i+X], 0) or 0
local pEast = x<X and mmax(zi-dem[i+1], 0) or 0
local pNorth = y>1 and mmax(zi-dem[i-X], 0) or 0
local pWest = x>1 and mmax(zi-dem[i-1], 0) or 0
local d = 0
local sum = pSouth + pEast + pNorth + pWest
local r = mrand() * sum
if sum > 0 then
if r < pSouth then
d = 1
elseif r-pSouth < pEast then
d = 2
elseif r-pSouth-pEast < pNorth then
d = 3
else
d = 4
end
end
local plist = { -- Get the height difference of the 4 neighbours (and 0 if uphill)
y<Y and mmax(zi-dem[i+X], 0) or 0, -- Southward
x<X and mmax(zi-dem[i+1], 0) or 0, -- Eastward
y>1 and mmax(zi-dem[i-X], 0) or 0, -- Northward
x>1 and mmax(zi-dem[i-1], 0) or 0, -- Westward
}
local d = flow_local(plist)
-- 'dirs': Direction toward which water flow
-- 'dirs2': Directions from which water comes
dirs[i] = d
@ -93,172 +110,89 @@ local function flow_routing(dem, dirs, lakes) -- 'dirs' and 'lakes' are optional
local links = {}
local basin_links
-- Function to analyse a link between two nodes
local function add_link(i1, i2, b1, isY)
-- i1, i2: coordinates of two nodes
-- b1: basin that contains i1
-- isY: whether the link is in Y direction
local b2
-- Note that basin number #0 represents the outside of the map; or if the coordinate is inside the map, means that the basin number is uninitialized.
if i2 == 0 then -- If outside the map
b2 = 0
else
b2 = basin_id[i2]
if b2 == 0 then -- If basin of i2 is not already computed, skip
return
end
end
if b2 ~= b1 then -- If these two nodes don't belong to the same basin, we have found a link between two adjacent basins
local elev = i2 == 0 and dem[i1] or mmax(dem[i1], dem[i2]) -- Elevation of the highest of the two sides of the link (or only i1 if b2 is map outside)
local l2 = basin_links[b2]
if not l2 then
l2 = {}
basin_links[b2] = l2
end
if not l2.elev or l2.elev > elev then -- If this link is lower than the lowest registered link between these two basins, register it as the new lowest pass
l2.elev = elev
l2.i = mmax(i1,i2)
l2.is_y = isY
l2[1] = b2
l2[2] = b1
end
end
end
for i=1, X*Y do
basin_id[i] = 0
end
local cur = nbasins
local ib = 0
while cur > 0 do
local i = singular[cur]
cur = cur - 1
if dirs[i] == 0 then
basin_links = {}
links[#links+1] = basin_links
ib = ib + 1
end
for ib=1, nbasins do
-- Here we will recursively search upstream from the singular node to determine its drainage basin
local queue = {singular[ib]} -- Start with the singular node, then this queue will be filled with water donors neighbours
basin_links = {}
links[#links+1] = basin_links
while #queue > 0 do
local i = tremove(queue)
basin_id[i] = ib
local d = dirs2[i] -- Get the directions water is coming from
-- Iterate through the 4 directions
-- Loop is unrolled on purpose, for performance (critical part!)
----------
-- EAST --
----------
if d >= 8 then -- River coming from the East
d = d - 8
cur = cur + 1
singular[cur] = i+1
queue[#queue+1] = i+1
-- If no river is coming from the East, we might be at the limit of two basins, thus we need to test adjacency.
elseif i%X > 0 then
if basin_id[i+1] ~= ib and basin_id[i+1] ~= 0 then
local b2 = basin_id[i+1]
local elev = mmax(dem[i], dem[i+1]) -- Elevation of the highest of the two sides of the link (or only i1 if b2 is map outside)
local l2 = basin_links[b2]
if not l2 then
l2 = {b2, ib, elev=elev, i=i+1, is_y=false}
basin_links[b2] = l2 -- Potential non-linear complexity here
elseif l2.elev > elev then -- If this link is lower than the lowest registered link between these two basins, register it as the new lowest pass
l2.elev = elev
l2.i = i+1
l2.is_y = false
l2[1] = b2
l2[2] = ib
end
end
add_link(i, i+1, ib, false)
else -- If the eastern neighbour is outside the map
local l2 = basin_links[0]
if not l2 then
l2 = {0, ib, elev=dem[i], i=i, is_y=false}
basin_links[0] = l2
elseif l2.elev > dem[i] then
l2.elev = dem[i]
l2.i = i
l2.is_y = false
l2[1] = 0
l2[2] = ib
end
add_link(i, 0, ib, false)
end
-----------
-- SOUTH --
-----------
if d >= 4 then -- River coming from the South
d = d - 4
cur = cur + 1
singular[cur] = i+X
queue[#queue+1] = i+X
elseif i <= X*(Y-1) then
if basin_id[i+X] ~= ib and basin_id[i+X] ~= 0 then
local b2 = basin_id[i+X]
local elev = mmax(dem[i], dem[i+X])
local l2 = basin_links[b2]
if not l2 then
l2 = {b2, ib, elev=elev, i=i+X, is_y=true}
basin_links[b2] = l2
elseif l2.elev > elev then
l2.elev = elev
l2.i = i+X
l2.is_y = true
l2[1] = b2
l2[2] = ib
end
end
add_link(i, i+X, ib, true)
else
local l2 = basin_links[0]
if not l2 then
l2 = {0, ib, elev=dem[i], i=i, is_y=true}
basin_links[0] = l2
elseif l2.elev > dem[i] then
l2.elev = dem[i]
l2.i = i
l2.is_y = true
l2[1] = 0
l2[2] = ib
end
add_link(i, 0, ib, true)
end
----------
-- WEST --
----------
if d >= 2 then -- River coming from the West
d = d - 2
cur = cur + 1
singular[cur] = i-1
queue[#queue+1] = i-1
elseif i%X ~= 1 then
if basin_id[i-1] ~= ib and basin_id[i-1] ~= 0 then
local b2 = basin_id[i-1]
local elev = mmax(dem[i], dem[i-1])
local l2 = basin_links[b2]
if not l2 then
l2 = {b2, ib, elev=elev, i=i, is_y=false}
basin_links[b2] = l2
elseif l2.elev > elev then
l2.elev = elev
l2.i = i
l2.is_y = false
l2[1] = b2
l2[2] = ib
end
end
add_link(i, i-1, ib, false)
else
local l2 = basin_links[0]
if not l2 then
l2 = {0, ib, elev=dem[i], i=i, is_y=false}
basin_links[0] = l2
elseif l2.elev > dem[i] then
l2.elev = dem[i]
l2.i = i
l2.is_y = false
l2[1] = 0
l2[2] = ib
end
add_link(i, 0, ib, false)
end
-----------
-- NORTH --
-----------
if d >= 1 then -- River coming from the North
cur = cur + 1
singular[cur] = i-X
queue[#queue+1] = i-X
elseif i > X then
if basin_id[i-X] ~= ib and basin_id[i-X] ~= 0 then
local b2 = basin_id[i-X]
local elev = mmax(dem[i], dem[i-X])
local l2 = basin_links[b2]
if not l2 then
l2 = {b2, ib, elev=elev, i=i, is_y=true}
basin_links[b2] = l2
elseif l2.elev > elev then
l2.elev = elev
l2.i = i
l2.is_y = true
l2[1] = b2
l2[2] = ib
end
end
add_link(i, i-X, ib, true)
else
local l2 = basin_links[0]
if not l2 then
l2 = {0, ib, elev=dem[i], i=i, is_y=true}
basin_links[0] = l2
elseif l2.elev > dem[i] then
l2.elev = dem[i]
l2.i = i
l2.is_y = true
l2[1] = 0
l2[2] = ib
end
add_link(i, 0, ib, true)
end
end
end
dirs2 = nil
@ -289,31 +223,26 @@ local function flow_routing(dem, dirs, lakes) -- 'dirs' and 'lakes' are optional
-- Mareš' optimizations mainly consist in skipping elements that have over 8 links, until extra links are removed when other elements are merged.
-- Note that for this step we are only working on basins, not grid nodes.
local lowlevel = {}
cur = 0
local ref = singular -- Reuse table
for i=0, nbasins do
if nlinks[i] <= 8 then
cur = cur + 1
lowlevel[cur] = i
ref[i] = cur
for i, n in pairs(nlinks) do
if n <= 8 then
lowlevel[i] = links[i]
end
end
local basin_graph = {}
for i=0, nbasins do
basin_graph[i] = {} -- Initialize (to ensure subtables don't go in the hash part)
end
for i=1, nbasins do
for n=1, nbasins do
-- Iterate in lowlevel but its contents may change during the loop
local b1 = lowlevel[cur]
cur = cur - 1
local lnk1 = links[b1]
-- 'next' called with only one argument always returns an element if table is not empty
local b1, lnk1 = next(lowlevel)
lowlevel[b1] = nil
local b2
local lowest = math.huge
local lnk1 = links[b1]
local i = 0
-- Look for lowest link
for bn, bdata in pairs(lnk1) do
i = i + 1
if bdata.elev < lowest then
lowest = bdata.elev
b2 = bn
@ -323,7 +252,13 @@ local function flow_routing(dem, dirs, lakes) -- 'dirs' and 'lakes' are optional
-- Add link to the graph, in both directions
local bound = lnk1[b2]
local bb1, bb2 = bound[1], bound[2]
basin_graph[bb1][bb2] = bound -- Potential non-linear complexity here
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
@ -334,9 +269,7 @@ local function flow_routing(dem, dirs, lakes) -- 'dirs' and 'lakes' are optional
nlinks[b2] = nlinks[b2] - 1
-- When the number of links is changing, we need to check whether the basin can be added to / removed from 'lowlevel'
if nlinks[b2] == 8 then
cur = cur + 1
lowlevel[cur] = b2
ref[b2] = cur
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
@ -346,16 +279,12 @@ local function flow_routing(dem, dirs, lakes) -- 'dirs' and 'lakes' are optional
if lnkn[b2] then -- If bassin bn is also linked to b2
nlinks[bn] = nlinks[bn] - 1 -- Then bassin bn is losing a link because it keeps only one link toward b1/b2 after the merge
if nlinks[bn] == 8 then
cur = cur + 1
lowlevel[cur] = bn
ref[bn] = cur
lowlevel[bn] = lnkn
end
else -- If bn was linked to b1 but not to b2
nlinks[b2] = nlinks[b2] + 1 -- Then b2 is gaining a link to bn because of the merge
if nlinks[b2] == 9 then
lowlevel[ref[b2]] = lowlevel[cur]
ref[lowlevel[cur]] = ref[b2]
cur = cur - 1
lowlevel[b2] = nil
end
end
@ -373,17 +302,15 @@ local function flow_routing(dem, dirs, lakes) -- 'dirs' and 'lakes' are optional
-- To orient the basin graph, we will consider that the ultimate basin water should flow into is the map outside (basin #0). We will start from it and recursively walk upstream to the neighbouring basins, using only links that are in the minimal spanning tree. This gives the flow direction of the links, and thus, the outlet of every basin.
-- This will also give lake elevation, which is the highest link encountered between map outside and the given basin on the spanning tree.
-- And within each basin, we need to modify flow directions to connect the singular node to the outlet.
local queue = {0}
local queuevalues = {-math.huge}
cur = 1
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}
while cur > 0 do
local b1, elev1 = queue[cur], queuevalues[cur] -- Pop from queue
cur = cur - 1
for n=1, nbasins do
local b1, elev1 = next(queue) -- Pop from queue
queue[b1] = nil
basin_lake[b1] = elev1
-- Iterate through b1's neighbours (according to the spanning tree)
for b2, bound in pairs(basin_graph[b1]) do
@ -422,9 +349,7 @@ local function flow_routing(dem, dirs, lakes) -- 'dirs' and 'lakes' are optional
until dir == 0 -- Stop when reaching the singular node
-- Add basin b2 into the queue, and keep the highest link elevation, that will define the elevation of the lake in b2
cur = cur + 1
queue[cur] = b2
queuevalues[cur] = mmax(elev1, bound.elev)
queue[b2] = mmax(elev1, bound.elev)
-- Remove b1 from b2's neighbours to avoid coming back to b1
basin_graph[b2][b1] = nil
end
@ -447,25 +372,31 @@ local function accumulate(dirs, waterq)
-- This means: how many nodes will give their water to that given node, directly or indirectly?
-- This is obtained by following rivers downstream and summing up the flow of every tributary, starting with a value of 1 at the sources.
-- This will give non-zero values for every node but only large values will be considered to be rivers.
waterq = waterq or {}
local X, Y = dirs.X, dirs.Y
waterq = waterq or {X=X, Y=Y}
local ndonors = {}
local waterq = {X=X, Y=Y}
for i=1, X*Y do
ndonors[i] = 0
waterq[i] = 1
end
-- Calculate the number of direct donors
for i=1, X*Y do
if dirs[i] == 1 then
ndonors[i+X] = ndonors[i+X] + 1
elseif dirs[i] == 2 then
ndonors[i+1] = ndonors[i+1] + 1
elseif dirs[i] == 3 then
ndonors[i-X] = ndonors[i-X] + 1
elseif dirs[i] == 4 then
ndonors[i-1] = ndonors[i-1] + 1
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
@ -507,4 +438,5 @@ end
return {
flow_routing = flow_routing,
accumulate = accumulate,
flow_methods = flow_methods,
}