Use biomegen.generate_all

by setting 'mapgen_rivers_distort = false' in minetest.conf

heightmap.lua

@@ -7,6 +7,8 @@ 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)
@@ -17,15 +19,14 @@ end
 local function heightmaps(minp, maxp)
 
 	local polygons = make_polygons(minp, maxp)
-	local incr = maxp.x-minp.x+1
-	local i0 = (minp.z-minp.z) * incr + (minp.x-minp.x) + 1
+	local incr = maxp.z-minp.z+1
 
 	local terrain_height_map = {}
 	local lake_height_map = {}
 
 	local i = 1
-	for x=minp.x, maxp.x do
 	for z=minp.z, maxp.z do
+		for x=minp.x, maxp.x do
 			local poly = polygons[i]
 
 			if poly then
@@ -106,8 +107,8 @@ local function heightmaps(minp, maxp)
 				terrain_height_map[i] = terrain_height
 				lake_height_map[i] = lake_height
 			else
-				terrain_height_map[i] = -31000
-				lake_height_map[i] = -31000
+				terrain_height_map[i] = MAP_BOTTOM
+				lake_height_map[i] = MAP_BOTTOM
 			end
 			i = i + 1
 		end

init.lua

@@ -8,13 +8,16 @@ local blocksize = mapgen_rivers.blocksize
 local sea_level = mapgen_rivers.sea_level
 local riverbed_slope = mapgen_rivers.riverbed_slope
 local elevation_chill = mapgen_rivers.elevation_chill
+local use_distort = mapgen_rivers.distort
+local use_biomes = mapgen_rivers.biomes
+local use_biomegen_mod = use_biomes and minetest.global_exists('biomegen')
+use_biomes = use_biomes and not use_biomegen_mod
 
+if use_biomegen_mod then
+	biomegen.set_elevation_chill(elevation_chill)
+end
 dofile(modpath .. 'noises.lua')
 
-local make_polygons = dofile(modpath .. 'polygons.lua')
-
-local transform_quadri = dofile(modpath .. 'geometry.lua')
-
 local heightmaps = dofile(modpath .. 'heightmap.lua')
 
 -- Linear interpolation
@@ -48,21 +51,32 @@ local function generate(minp, maxp, seed)
 			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)
-		noise_distort_obj = minetest.get_perlin_map(mapgen_rivers.noise_params.distort_amplitude, 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
@@ -87,9 +101,12 @@ local function generate(minp, maxp, seed)
 
 		local pminp = {x=math.floor(xmin), z=math.floor(zmin)}
 		local pmaxp = {x=math.floor(xmax)+1, z=math.floor(zmax)+1}
-	local incr = pmaxp.z-pminp.z+1
-	local i_origin = 1 - pminp.x*incr - pminp.z
-	local terrain_map, lake_map = heightmaps(pminp, pmaxp)
+		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")
@@ -118,26 +135,38 @@ local function generate(minp, maxp, seed)
 		for x = minp.x, maxp.x do
 			local ivm = a:index(x, minp.y, z)
 			local ground_above = false
-			local temperature = noise_heat_map[i2d]+noise_heat_blend_map[i2d]
+			local temperature
+			if use_biomes then
+				temperature = noise_heat_map[i2d]+noise_heat_blend_map[i2d]
+			end
+			local terrain, lake
+			if not use_distort then
+				terrain = terrain_map[i2d]
+				lake = lake_map[i2d]
+			end
+
 			for y = maxp.y+1, minp.y, -1 do
+				if use_distort then
 					local xn = noise_x_map[nid]
 					local zn = noise_z_map[nid]
 					local x0 = math.floor(xn)
 					local z0 = math.floor(zn)
 
-				local i0 = i_origin + x0*incr + z0
-				local i1 = i0+incr
-				local i2 = i1+1
-				local i3 = i0+1
+					local i0 = i_origin + z0*incr + x0
+					local i1 = i0+1
+					local i2 = i1+incr
+					local i3 = i2-1
+
+					terrain = interp(terrain_map[i0], terrain_map[i1], terrain_map[i2], terrain_map[i3], xn-x0, zn-z0)
+					lake = math.min(lake_map[i0], lake_map[i1], lake_map[i2], lake_map[i3])
+				end
 
-				local terrain = interp(terrain_map[i0], terrain_map[i1], terrain_map[i2], terrain_map[i3], xn-x0, zn-z0)
 				if y <= maxp.y then
-					local lake = math.min(lake_map[i0], lake_map[i1], lake_map[i2], lake_map[i3])
 
 					local is_lake = lake > terrain
 					local ivm = a:index(x, y, z)
 					if y <= terrain then
-						if y <= terrain-1 or ground_above 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
@@ -152,8 +181,7 @@ local function generate(minp, maxp, seed)
 							end
 						end
 					elseif y <= lake and lake > sea_level then
-						local temperature_y = temperature - y*elevation_chill
-						if temperature_y >= 0 then
+						if not use_biomes or temperature - y*elevation_chill >= 0 then
 							data[ivm] = c_rwater
 						else
 							data[ivm] = c_ice
@@ -166,16 +194,29 @@ local function generate(minp, maxp, seed)
 				ground_above = y <= terrain
 
 				ivm = ivm + ystride
+				if use_distort then
 					nid = nid + incrY
 				end
-			nid = nid + incrX
-			i2d = i2d + 1
-		end
-		nid = nid + incrZ
 			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)
+		mietest.generate_ores(vm, minp, maxp)
+	end
+
 	vm:set_lighting({day = 0, night = 0})
 	vm:calc_lighting()
 	vm:update_liquids()

mod.conf

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

polygons.lua

@@ -96,7 +96,7 @@ local function make_polygons(minp, maxp)
 		init = true
 	end
 
-	local chulens = maxp.z - minp.z + 1
+	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.
@@ -115,44 +115,44 @@ local function make_polygons(minp, maxp)
 			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 X position (scanline algorithm)
-			-- Calculate the min and max X positions 
-			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)
+			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 x=xmin, xmax do
-				bounds[x] = {}
+			for z=zmin, zmax do
+				bounds[z] = {}
 			end
 
 			local i1 = 4
 			for i2=1, 4 do -- Loop on 4 edges
-				local x1, x2 = poly_x[i1], poly_x[i2]
-				-- Calculate the integer X positions over which this edge spans
-				local lxmin = math.floor(blocksize*math.min(x1, x2))+1
-				local lxmax = math.floor(blocksize*math.max(x1, x2))
-				if lxmin <= lxmax then -- If there is at least one position in it
 				local z1, z2 = poly_z[i1], poly_z[i2]
-					-- Calculate coefficient of the equation defining the edge: Z=aX+b
-					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
-						-- For every X position involved, add the intercepted Z position in the table
-						table.insert(bounds[x], a*x+b)
+				-- 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 x=xmin, xmax do
-				-- Now sort the bounds list
-				local xlist = bounds[x]
-				table.sort(xlist)
-				local c = math.floor(#xlist/2)
-				for l=1, c do
-					-- Take pairs of Z coordinates: all positions between them belong to the polygon.
-					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
+				-- 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

settings.lua

@@ -47,7 +47,8 @@ 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) To be implemented: should be possible to disable distorsion
+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)

settings.py

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

terrain_default.conf

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

terrain_rivers.py

@@ -7,24 +7,48 @@ 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)
 
-if argc > 1:
-	mapsize = int(sys.argv[1])
-else:
-	mapsize = 400
+params = {}
+
+if argc > 1:
+	if os.path.isfile(sys.argv[1]):
+		params = settings.read_config_file(sys.argv[1])
+	else:
+		mapsize = int(sys.argv[1])
+
+def get_setting(name, default):
+	if name in params:
+		return params[name]
+	return default
+
+mapsize = int(get_setting('mapsize', 400))
+scale = float(get_setting('scale', 200.0))
+vscale = float(get_setting('vscale', 200.0))
+offset = float(get_setting('offset', 0.0))
+persistence = float(get_setting('persistence', 0.5))
+lacunarity = float(get_setting('lacunarity', 2.0))
+
+K = float(get_setting('K', 1.0))
+m = float(get_setting('m', 0.35))
+d = float(get_setting('d', 1.0))
+sea_level = float(get_setting('sea_level', 0.0))
+flex_radius = float(get_setting('flex_radius', 20.0))
+
+time = float(get_setting('time', 10.0))
+niter = int(get_setting('niter', 10))
 
-scale = mapsize / 2
 n = np.zeros((mapsize+1, mapsize+1))
 
 # Set noise parameters
 params = {
     "octaves" : int(np.ceil(np.log2(mapsize)))+1,
-    "persistence" : 0.5,
-    "lacunarity" : 2.,
+    "persistence" : persistence,
+    "lacunarity" : lacunarity,
 }
 
 # Determine noise offset randomly
@@ -36,35 +60,27 @@ for x in range(mapsize+1):
     for y in range(mapsize+1):
         n[x,y] = noise.snoise2(x/scale + xbase, y/scale + ybase, **params)
 
-nn = n*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)
+
+dt = time/niter
 
 # Run the model's processes: the order in which the processes are run is arbitrary and could be changed.
-print('Flow calculation 1')
-model.calculate_flow()
-
-print('Advection 1')
-model.advection(2)
-
-print('Isostatic equilibration 1')
-model.adjust_isostasy()
-
-print('Flow calculation 2')
+print('Initial flow calculation')
 model.calculate_flow()
 
+for i in range(niter):
+    print('Iteration {:d} of {:d}'.format(i+1, niter))
     print('Diffusion')
-model.diffusion(4)
-
-print('Advection 2')
-model.advection(2)
-
-print('Isostatic equilibration 2')
+    model.diffusion(dt)
+    print('Advection')
+    model.advection(dt)
+    print('Isostatic equilibration')
     model.adjust_isostasy()
-
-print('Flow calculation 3')
+    print('Flow calculation')
     model.calculate_flow()
 
 print('Done')

view_map.py

@@ -5,21 +5,21 @@ import zlib
 import matplotlib.colors as mcol
 import matplotlib.pyplot as plt
 
-def view_map(dem, lakes, rivers):
+def view_map(dem, lakes, rivers, scale):
     plt.subplot(1,3,1)
-    plt.pcolormesh(dem, cmap='viridis')
+    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,2)
-    plt.pcolormesh(lakes, cmap='viridis')
+    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,3)
-    plt.pcolormesh(rivers, cmap='Blues', norm=mcol.LogNorm())
+    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')
@@ -29,8 +29,12 @@ def view_map(dem, lakes, rivers):
 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)
@@ -45,4 +49,4 @@ if __name__ == "__main__":
     lakes = load_map('lakes', '>i2', shape)
     rivers = load_map('rivers', '>u4', shape)
 
-    view_map(dem, lakes, rivers)
+    view_map(dem, lakes, rivers, scale)