Initial commit: working example using a basis of Simplex noise and implementing river flowing, lakes, and erosion
This commit is contained in:
Gael-de-Sailly
commit
0bf351b2f6
|
|
@ -0,0 +1,7 @@
|
|||
__pycache__/
|
||||
dem
|
||||
lakes
|
||||
links
|
||||
rivers
|
||||
size
|
||||
unused/
|
||||
|
|
@ -0,0 +1,77 @@
|
|||
import numpy as np
|
||||
import scipy.ndimage as im
|
||||
import rivermapper as rm
|
||||
|
||||
def advection(dem, dirs, rivers, time, K=1, m=0.5, sea_level=0):
|
||||
dirs = dirs.copy()
|
||||
dirs[0,:] = 0
|
||||
dirs[-1,:] = 0
|
||||
dirs[:,0] = 0
|
||||
dirs[:,-1] = 0
|
||||
|
||||
adv_time = 1 / (K*rivers**m)
|
||||
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]):
|
||||
x0, y0 = x, y
|
||||
x1, y1 = x, y
|
||||
remaining = time
|
||||
while True:
|
||||
flow_dir = dirs[y0,x0]
|
||||
if flow_dir == 0:
|
||||
remaining = 0
|
||||
break
|
||||
elif flow_dir == 1:
|
||||
y1 += 1
|
||||
elif flow_dir == 2:
|
||||
x1 += 1
|
||||
elif flow_dir == 3:
|
||||
y1 -= 1
|
||||
elif flow_dir == 4:
|
||||
x1 -= 1
|
||||
|
||||
if remaining <= adv_time[y0,x0]:
|
||||
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]
|
||||
|
||||
return np.minimum(dem, dem_new)
|
||||
|
||||
def diffusion(dem, time, d=1):
|
||||
radius = d * time**.5
|
||||
return im.gaussian_filter(dem, radius, mode='reflect')
|
||||
|
||||
class EvolutionModel:
|
||||
def __init__(self, dem, K=1, m=0.5, d=1, sea_level=0, flow=False):
|
||||
self.dem = dem
|
||||
#self.bedrock = dem
|
||||
self.K = K
|
||||
self.m = m
|
||||
self.d = d
|
||||
self.sea_level = sea_level
|
||||
#self.ref_isostasy = im.gaussian_filter(dem, radius, mode='reflect')
|
||||
if flow:
|
||||
self.calculate_flow()
|
||||
else:
|
||||
self.lakes = dem
|
||||
self.dirs = np.zeros(dem.shape, dtype='u1')
|
||||
self.rivers = np.zeros(dem.shape, dtype='u4')
|
||||
self.flow_uptodate = False
|
||||
|
||||
def calculate_flow(self):
|
||||
self.dirs, self.lakes, self.rivers = rm.flow(self.dem)
|
||||
self.flow_uptodate = True
|
||||
|
||||
def advection(self, time):
|
||||
dem = advection(self.lakes, self.dirs, self.rivers, time, K=self.K, m=self.m, sea_level=self.sea_level)
|
||||
self.dem = np.minimum(dem, self.dem)
|
||||
self.flow_uptodate = False
|
||||
|
||||
def diffusion(self, time):
|
||||
self.dem = diffusion(self.dem, time, d=self.d)
|
||||
self.flow_uptodate = False
|
||||
|
|
@ -0,0 +1,178 @@
|
|||
local modpath = minetest.get_modpath(minetest.get_current_modname()) .. '/'
|
||||
local worldpath = minetest.get_worldpath() .. '/'
|
||||
local load_map = dofile(modpath .. 'load.lua')
|
||||
|
||||
local function copy_if_needed(filename)
|
||||
local wfilename = worldpath..filename
|
||||
local wfile = io.open(wfilename, 'r')
|
||||
if wfile then
|
||||
wfile:close()
|
||||
return
|
||||
end
|
||||
local mfilename = modpath..filename
|
||||
local mfile = io.open(mfilename, 'r')
|
||||
local wfile = io.open(wfilename, 'w')
|
||||
wfile:write(mfile:read("*all"))
|
||||
mfile:close()
|
||||
wfile:close()
|
||||
end
|
||||
|
||||
copy_if_needed('size')
|
||||
local sfile = io.open(worldpath..'size')
|
||||
local X = tonumber(sfile:read('*l'))
|
||||
local Z = tonumber(sfile:read('*l'))
|
||||
|
||||
copy_if_needed('dem')
|
||||
local dem = load_map(worldpath..'dem', 2, true)
|
||||
copy_if_needed('lakes')
|
||||
local lakes = load_map(worldpath..'lakes', 2, true)
|
||||
copy_if_needed('links')
|
||||
local links = load_map(worldpath..'links', 1, false)
|
||||
copy_if_needed('rivers')
|
||||
local rivers = load_map(worldpath..'rivers', 4, false)
|
||||
|
||||
local function index(x, z)
|
||||
return z*X+x+1
|
||||
end
|
||||
|
||||
local function interp(v00, v01, v10, v11, xf, zf)
|
||||
v0 = v01*xf + v00*(1-xf)
|
||||
v1 = v11*xf + v10*(1-xf)
|
||||
return v1*zf + v0*(1-zf)
|
||||
end
|
||||
|
||||
local data = {}
|
||||
|
||||
local blocksize = 6
|
||||
local sea_level = 1
|
||||
local min_catchment = 25
|
||||
|
||||
local storage = minetest.get_mod_storage()
|
||||
if storage:contains("blocksize") then
|
||||
blocksize = storage:get_int("blocksize")
|
||||
else
|
||||
storage:set_int("blocksize", blocksize)
|
||||
end
|
||||
if storage:contains("sea_level") then
|
||||
sea_level = storage:get_int("sea_level")
|
||||
else
|
||||
storage:set_int("sea_level", sea_level)
|
||||
end
|
||||
if storage:contains("min_catchment") then
|
||||
min_catchment = storage:get_float("min_catchment")
|
||||
else
|
||||
storage:set_float("min_catchment", min_catchment)
|
||||
end
|
||||
|
||||
local function generate(minp, maxp, seed)
|
||||
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_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 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.
|
||||
|
||||
for x = minp.x, maxp.x do
|
||||
for z = minp.z, maxp.z do
|
||||
local xb = x/blocksize
|
||||
local zb = z/blocksize
|
||||
|
||||
if xb >= 0 and xb < X-1 and zb >= 0 and zb < Z-1 then
|
||||
local x0 = math.floor(xb)
|
||||
local x1 = x0+1
|
||||
local z0 = math.floor(zb)
|
||||
local z1 = z0+1
|
||||
|
||||
local xf = xb - x0
|
||||
local zf = zb - z0
|
||||
|
||||
local i00 = index(x0,z0)
|
||||
local i01 = index(x1,z0)
|
||||
local i10 = index(x0,z1)
|
||||
local i11 = index(x1,z1)
|
||||
|
||||
local terrain_height = math.floor(interp(
|
||||
dem[i00],
|
||||
dem[i01],
|
||||
dem[i10],
|
||||
dem[i11],
|
||||
xf, zf
|
||||
))
|
||||
|
||||
local lake_height = math.floor(math.min(
|
||||
lakes[i00],
|
||||
lakes[i01],
|
||||
lakes[i10],
|
||||
lakes[i11]
|
||||
))
|
||||
|
||||
local is_lake = lake_height > terrain_height
|
||||
|
||||
local is_river = false
|
||||
if xf == 0 then
|
||||
if links[i00] == 1 and rivers[i00] >= min_catchment then
|
||||
is_river = true
|
||||
elseif links[i10] == 3 and rivers[i10] >= min_catchment then
|
||||
is_river = true
|
||||
end
|
||||
end
|
||||
|
||||
if zf == 0 then
|
||||
if links[i00] == 2 and rivers[i00] >= min_catchment then
|
||||
is_river = true
|
||||
elseif links[i01] == 4 and rivers[i01] >= min_catchment then
|
||||
is_river = true
|
||||
end
|
||||
end
|
||||
|
||||
local ivm = a:index(x, minp.y-1, z)
|
||||
|
||||
if terrain_height >= minp.y then
|
||||
for y=minp.y, math.min(maxp.y, terrain_height) do
|
||||
if y == terrain_height then
|
||||
if is_lake or y <= sea_level then
|
||||
data[ivm] = c_sand
|
||||
elseif is_river then
|
||||
data[ivm] = c_rwater
|
||||
else
|
||||
data[ivm] = c_lawn
|
||||
end
|
||||
else
|
||||
data[ivm] = c_stone
|
||||
end
|
||||
ivm = ivm + ystride
|
||||
end
|
||||
end
|
||||
|
||||
if lake_height > sea_level then
|
||||
if is_lake and lake_height > minp.y then
|
||||
for y=math.max(minp.y, terrain_height+1), math.min(maxp.y, lake_height) do
|
||||
data[ivm] = c_rwater
|
||||
ivm = ivm + ystride
|
||||
end
|
||||
end
|
||||
else
|
||||
for y=math.max(minp.y, terrain_height+1), math.min(maxp.y, sea_level) do
|
||||
data[ivm] = c_water
|
||||
ivm = ivm + ystride
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
vm:set_data(data)
|
||||
minetest.generate_ores(vm, minp, maxp)
|
||||
vm:set_lighting({day = 0, night = 0})
|
||||
vm:calc_lighting()
|
||||
vm:update_liquids()
|
||||
vm:write_to_map()
|
||||
end
|
||||
|
||||
minetest.register_on_generated(generate)
|
||||
|
|
@ -0,0 +1,28 @@
|
|||
local function load_map(filename, bytes, signed)
|
||||
local file = io.open(filename, 'r')
|
||||
local data = file:read('*all')
|
||||
|
||||
local map = {}
|
||||
|
||||
local size = math.floor(#data/bytes)
|
||||
|
||||
for i=1, size do
|
||||
local i0, i1 = (i-1)*bytes+1, i*bytes
|
||||
local elements = {data:byte(i0, i1)}
|
||||
local n = elements[1]
|
||||
if signed and n >= 128 then
|
||||
n = n - 256
|
||||
end
|
||||
|
||||
for j=2, bytes do
|
||||
n = n*256 + elements[j]
|
||||
end
|
||||
|
||||
map[i] = n
|
||||
end
|
||||
file:close()
|
||||
|
||||
return map
|
||||
end
|
||||
|
||||
return load_map
|
||||
|
|
@ -0,0 +1,98 @@
|
|||
import numpy as np
|
||||
import heapq
|
||||
import sys
|
||||
|
||||
# Conventions:
|
||||
# 1 = South (+Y)
|
||||
# 2 = East (+X)
|
||||
# 3 = North (-Y)
|
||||
# 4 = West (-X)
|
||||
|
||||
sys.setrecursionlimit(65536)
|
||||
|
||||
neighbours_dirs = np.array([
|
||||
[0,1,0],
|
||||
[2,0,4],
|
||||
[0,3,0],
|
||||
], dtype='u1')
|
||||
|
||||
neighbours_pattern = neighbours_dirs > 0
|
||||
|
||||
def flow_dirs_lakes(dem):
|
||||
(Y, X) = dem.shape
|
||||
|
||||
dem_margin = np.zeros((Y+2, X+2))
|
||||
dem_margin[1:-1,1:-1] = dem
|
||||
|
||||
# Initialize: list map borders
|
||||
borders = []
|
||||
|
||||
for x in range(1,X+1):
|
||||
dem_north = dem_margin[1,x]
|
||||
borders.append((dem_north, dem_north, 1, x))
|
||||
dem_south = dem_margin[Y,x]
|
||||
borders.append((dem_south, dem_south, Y, x))
|
||||
|
||||
for y in range(2,Y):
|
||||
dem_west = dem_margin[y,1]
|
||||
borders.append((dem_west, dem_west, y, 1))
|
||||
dem_east = dem_margin[y,X]
|
||||
borders.append((dem_east, dem_east, y, X))
|
||||
|
||||
heapq.heapify(borders)
|
||||
|
||||
dirs = np.zeros((Y+2, X+2), dtype='u1')
|
||||
dirs[-2:,:] = 1
|
||||
dirs[:,-2:] = 2
|
||||
dirs[ :2,:] = 3
|
||||
dirs[:, :2] = 4
|
||||
|
||||
lakes = np.zeros((Y, X))
|
||||
|
||||
def add_point(y, x, altmax):
|
||||
alt = dem_margin[y, x]
|
||||
heapq.heappush(borders, (alt, altmax, y, x))
|
||||
|
||||
while len(borders) > 0:
|
||||
(alt, altmax, y, x) = heapq.heappop(borders)
|
||||
neighbours = dirs[y-1:y+2, x-1:x+2]
|
||||
empty_neighbours = (neighbours == 0) * neighbours_pattern
|
||||
neighbours += empty_neighbours * neighbours_dirs
|
||||
|
||||
lake = max(alt, altmax)
|
||||
lakes[y-1,x-1] = lake
|
||||
|
||||
coords = np.transpose(empty_neighbours.nonzero())
|
||||
for (dy,dx) in coords-1:
|
||||
add_point(y+dy, x+dx, lake)
|
||||
|
||||
return dirs[1:-1,1:-1], lakes
|
||||
|
||||
def accumulate(dirs, dem=None):
|
||||
(Y, X) = dirs.shape
|
||||
dirs_margin = np.zeros((Y+2,X+2))-1
|
||||
dirs_margin[1:-1,1:-1] = dirs
|
||||
quantity = np.zeros((Y, X), dtype='u4')
|
||||
|
||||
def calculate_quantity(y, x):
|
||||
if quantity[y,x] > 0:
|
||||
return quantity[y,x]
|
||||
q = 1
|
||||
neighbours = dirs_margin[y:y+3, x:x+3]
|
||||
donors = neighbours == neighbours_dirs
|
||||
|
||||
coords = np.transpose(donors.nonzero())
|
||||
for (dy,dx) in coords-1:
|
||||
q += calculate_quantity(y+dy, x+dx)
|
||||
quantity[y, x] = q
|
||||
return q
|
||||
|
||||
for x in range(X):
|
||||
for y in range(Y):
|
||||
calculate_quantity(y, x)
|
||||
|
||||
return quantity
|
||||
|
||||
def flow(dem):
|
||||
dirs, lakes = flow_dirs_lakes(dem)
|
||||
return dirs, lakes, accumulate(dirs, dem)
|
||||
|
|
@ -0,0 +1,8 @@
|
|||
import numpy as np
|
||||
|
||||
def save(data, fname, dtype=None):
|
||||
if dtype is not None:
|
||||
data = data.astype(dtype)
|
||||
|
||||
with open(fname, 'wb') as f:
|
||||
f.write(data.tobytes())
|
||||
|
|
@ -0,0 +1,101 @@
|
|||
#!/usr/bin/env python3
|
||||
|
||||
import numpy as np
|
||||
import noise
|
||||
from save import save
|
||||
from erosion import EvolutionModel
|
||||
import os
|
||||
import sys
|
||||
|
||||
# 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
|
||||
|
||||
scale = mapsize / 2
|
||||
n = np.zeros((mapsize, mapsize))
|
||||
|
||||
#micronoise_depth = 0.05
|
||||
|
||||
params = {
|
||||
"octaves" : 8,
|
||||
"persistence" : 0.5,
|
||||
"lacunarity" : 2.,
|
||||
}
|
||||
|
||||
xbase = np.random.randint(65536)
|
||||
ybase = np.random.randint(65536)
|
||||
|
||||
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
|
||||
|
||||
print('Initializing model')
|
||||
model = EvolutionModel(nn, K=1, m=0.35, d=1, sea_level=0)
|
||||
|
||||
print('Flow calculation 1')
|
||||
model.calculate_flow()
|
||||
|
||||
print('Advection 1')
|
||||
model.advection(2)
|
||||
|
||||
print('Flow calculation 2')
|
||||
model.calculate_flow()
|
||||
|
||||
print('Diffusion')
|
||||
model.diffusion(4)
|
||||
|
||||
print('Advection 2')
|
||||
model.advection(2)
|
||||
|
||||
print('Flow calculation 3')
|
||||
model.calculate_flow()
|
||||
|
||||
print('Done')
|
||||
|
||||
save(model.dem, 'dem', dtype='>i2')
|
||||
save(model.lakes, 'lakes', dtype='>i2')
|
||||
save(model.dirs, 'links', dtype='u1')
|
||||
save(model.rivers, 'rivers', dtype='>u4')
|
||||
|
||||
with open('size', 'w') as sfile:
|
||||
sfile.write('{:d}\n{:d}'.format(mapsize, mapsize))
|
||||
|
||||
try:
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
plt.subplot(2,2,1)
|
||||
plt.pcolormesh(nn, cmap='viridis')
|
||||
plt.gca().set_aspect('equal', 'box')
|
||||
#plt.colorbar(orientation='horizontal')
|
||||
plt.title('Raw elevation')
|
||||
|
||||
plt.subplot(2,2,2)
|
||||
plt.pcolormesh(model.lakes, cmap='viridis')
|
||||
plt.gca().set_aspect('equal', 'box')
|
||||
#plt.colorbar(orientation='horizontal')
|
||||
plt.title('Lake surface elevation')
|
||||
|
||||
plt.subplot(2,2,3)
|
||||
plt.pcolormesh(model.dem, cmap='viridis')
|
||||
plt.gca().set_aspect('equal', 'box')
|
||||
#plt.colorbar(orientation='horizontal')
|
||||
plt.title('Elevation after advection')
|
||||
|
||||
plt.subplot(2,2,4)
|
||||
plt.pcolormesh(model.rivers, vmin=0, vmax=mapsize**2/25, cmap='Blues')
|
||||
plt.gca().set_aspect('equal', 'box')
|
||||
#plt.colorbar(orientation='horizontal')
|
||||
plt.title('Rivers discharge')
|
||||
|
||||
plt.show()
|
||||
except:
|
||||
pass
|
||||
Loading…
Reference in New Issue