From 6af6795d90f26b2f23eff0fea567f3b9f152651a Mon Sep 17 00:00:00 2001 From: Gael-de-Sailly Date: Sun, 26 Apr 2020 17:13:38 +0200 Subject: [PATCH] Comment and clarify --- README.md | 8 ++++---- bounds.py | 16 ++++++++++++++-- erosion.py | 27 +++++++++++++++++++-------- geometry.lua | 10 ++++++++-- init.lua | 6 ++++-- polygons.lua | 22 ++++++++++++++++++++-- rivermapper.py | 35 +++++++++++++++++++++++++---------- terrain_rivers.py | 15 ++++++++++----- view_map.py | 2 +- 9 files changed, 105 insertions(+), 36 deletions(-) diff --git a/README.md b/README.md index 5f25a03..c0c5217 100644 --- a/README.md +++ b/README.md @@ -1,17 +1,17 @@ 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) @@ -23,7 +23,7 @@ Run the script `terrain_rivers.py` via command line. You can optionally append t ``` ./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 401x401 map, it should take between 1 and 2 minutes. It will generate 5 files directly in the mod folder, containing the map 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. diff --git a/bounds.py b/bounds.py index 5a43532..830a87d 100644 --- a/bounds.py +++ b/bounds.py @@ -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 diff --git a/erosion.py b/erosion.py index bd5879e..7f830cf 100644 --- a/erosion.py +++ b/erosion.py @@ -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 diff --git a/geometry.lua b/geometry.lua index 773a39f..fcf3ce4 100644 --- a/geometry.lua +++ b/geometry.lua @@ -1,27 +1,33 @@ 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) diff --git a/init.lua b/init.lua index 4575cd5..957fc64 100644 --- a/init.lua +++ b/init.lua @@ -12,6 +12,7 @@ local make_polygons = dofile(modpath .. 'polygons.lua') local transform_quadri = dofile(modpath .. 'geometry.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) @@ -44,6 +45,7 @@ local function generate(minp, maxp, seed) local xf, zf = transform_quadri(poly.x, poly.z, x/blocksize, z/blocksize) local i00, i01, i11, i10 = unpack(poly.i) + -- Test the 4 edges to see whether we are in a river or not local is_river = false local depth_factor = 0 local r_west, r_north, r_east, r_south = unpack(poly.rivers) @@ -66,7 +68,7 @@ local function generate(minp, maxp, seed) zf = 0 end - if not is_river then + if not is_river then -- Test corners also local c_NW, c_NE, c_SE, c_SW = unpack(poly.river_corners) if xf+zf < c_NW then is_river = true @@ -87,7 +89,7 @@ local function generate(minp, maxp, seed) end end - if not is_river then + if not is_river then -- Redefine indicesto have 0/1 on the riverbanks (avoids ugly edge cuts, at least for small rivers) xf = (xf-r_west) / (r_east-r_west) zf = (zf-r_north) / (r_south-r_north) end diff --git a/polygons.lua b/polygons.lua index 22eae57..443f3ac 100644 --- a/polygons.lua +++ b/polygons.lua @@ -44,6 +44,7 @@ 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 @@ -66,26 +67,33 @@ local function river_width(flow) return math.min(wfactor * flow ^ wpower, 1) end +-- On map generation, determine into which polygon every point (in 2D) will fall. +-- Also store polygon-specific data local function make_polygons(minp, maxp) local chulens = maxp.z - minp.z + 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 = {} + 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) + -- And initialize the arrays for x=xmin, xmax do bounds[x] = {} end @@ -93,27 +101,33 @@ local function make_polygons(minp, maxp) 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 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) 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 + -- Fill the map at these places polygons[i] = polygon i = i + 1 end @@ -123,10 +137,13 @@ local function make_polygons(minp, maxp) polygon.dem = {dem[iA], dem[iB], dem[iC], dem[iD]} polygon.lake = math.min(lakes[iA], lakes[iB], lakes[iC], lakes[iD]) + -- Now, rivers. + -- Start by finding the river width (if any) for the polygon's 4 edges. 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]) + -- Only if opposite rivers overlap (should be rare) if river_west > river_east then local mean = (river_west + river_east) / 2 river_west = mean @@ -139,6 +156,7 @@ local function make_polygons(minp, maxp) end polygon.rivers = {river_west, river_north, river_east, river_south} + -- Look for river corners local around = {0,0,0,0,0,0,0,0} if zp > 0 then around[1] = river_width(bounds_z[iA-X]) diff --git a/rivermapper.py b/rivermapper.py index b494748..db1c563 100644 --- a/rivermapper.py +++ b/rivermapper.py @@ -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) diff --git a/terrain_rivers.py b/terrain_rivers.py index 2df2f7f..108476b 100755 --- a/terrain_rivers.py +++ b/terrain_rivers.py @@ -20,28 +20,29 @@ else: scale = (mapsize-1) / 2 n = np.zeros((mapsize, mapsize)) -#micronoise_depth = 0.05 - +# Set noise parameters params = { "octaves" : int(np.ceil(np.log2(mapsize-1)))+1, "persistence" : 0.5, "lacunarity" : 2., } +# Determine noise offset randomly xbase = np.random.randint(65536) ybase = np.random.randint(65536) +# Generate the noise for x in range(mapsize): for y in range(mapsize): 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 +# Initialize landscape evolution model print('Initializing model') model = EvolutionModel(nn, K=1, m=0.35, d=1, sea_level=0) +# 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() @@ -68,12 +69,15 @@ 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)) +# Convert offset in 8-bits offset_x = np.clip(np.floor(ox * 256), -128, 127) offset_y = np.clip(np.floor(oy * 256), -128, 127) +# Save the files save(model.dem, 'dem', dtype='>i2') save(model.lakes, 'lakes', dtype='>i2') save(np.abs(bx), 'bounds_x', dtype='>i4') @@ -86,6 +90,7 @@ save(model.rivers, 'rivers', dtype='>u4') with open('size', 'w') as sfile: sfile.write('{:d}\n{:d}'.format(mapsize, mapsize)) +# Display the map if matplotlib is found try: import matplotlib.pyplot as plt @@ -111,7 +116,7 @@ try: 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.title('Rivers flux') plt.show() except: diff --git a/view_map.py b/view_map.py index 475dffa..cbf8aab 100755 --- a/view_map.py +++ b/view_map.py @@ -34,6 +34,6 @@ 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.title('Rivers flux') plt.show()