Removed Python terrainlib

This commit is contained in:
Gaël C 2021-06-02 18:56:46 +02:00
parent 7495d8a690
commit 0427b42d17
9 changed files with 1 additions and 663 deletions

View File

@ -1,179 +0,0 @@
#!/usr/bin/env python3
import numpy as np
from noise import snoise2
import os
import sys
import terrainlib
def noisemap(X, Y, scale=0.01, vscale=1.0, offset=0.0, log=False, **params):
# Determine noise offset randomly
xbase = np.random.randint(8192)-4096
ybase = np.random.randint(8192)-4096
if log:
vscale /= offset
# Generate the noise
n = np.zeros((X, Y))
for x in range(X):
for y in range(Y):
n[x,y] = snoise2(x/scale + xbase, y/scale + ybase, **params)
if log:
return np.exp(n*vscale) * offset
else:
return n*vscale + offset
### PARSE COMMAND-LINE ARGUMENTS
argc = len(sys.argv)
config_file = 'terrain_default.conf'
output_dir = 'river_data'
params_from_args = {}
i = 1 # Index of arguments
j = 1 # Number of 'orphan' arguments (the ones that are not preceded by '--something')
while i < argc:
arg = sys.argv[i]
if arg[:2] == '--':
pname = arg[2:]
v = None
split = pname.split('=', maxsplit=1)
if len(split) == 2:
pname, v = split
i += 1
elif i+1 < argc:
v = sys.argv[i+1]
i += 2
if v is not None:
if pname == 'config':
config_file = v
elif pname == 'output':
output_dir = v
else:
params_from_args[pname] = v
else:
if j == 1:
config_file = arg
elif j == 2:
output_dir = arg
i += 1
j += 1
print(config_file, output_dir)
params = terrainlib.read_config_file(config_file)
params.update(params_from_args) # Params given from args prevail against conf file
### READ SETTINGS
def get_setting(name, default):
if name in params:
return params[name]
return default
mapsize = int(get_setting('mapsize', 1000))
scale = float(get_setting('scale', 400.0))
vscale = float(get_setting('vscale', 300.0))
offset = float(get_setting('offset', 0.0))
persistence = float(get_setting('persistence', 0.6))
lacunarity = float(get_setting('lacunarity', 2.0))
K = float(get_setting('K', 0.5))
m = float(get_setting('m', 0.5))
d = float(get_setting('d', 0.5))
sea_level = float(get_setting('sea_level', 0.0))
sea_level_variations = float(get_setting('sea_level_variations', 0.0))
sea_level_variations_time = float(get_setting('sea_level_variations_time', 1.0))
flex_radius = float(get_setting('flex_radius', 20.0))
flow_method = get_setting('flow_method', 'semirandom')
time = float(get_setting('time', 10.0))
niter = int(get_setting('niter', 10))
### MAKE INITIAL TOPOGRAPHY
n = np.zeros((mapsize+1, mapsize+1))
# Set noise parameters
params = {
"offset" : offset,
"vscale" : vscale,
"scale" : scale,
"octaves" : int(np.ceil(np.log2(mapsize)))+1,
"persistence" : persistence,
"lacunarity" : lacunarity,
}
params_sealevel = {
"octaves" : 1,
"persistence" : 1,
"lacunarity" : 2,
}
if sea_level_variations != 0.0:
sea_ybase = np.random.randint(8192)-4096
sea_level_ref = snoise2(time * (1-1/niter) / sea_level_variations, sea_ybase, **params_sealevel) * sea_level_variations
params['offset'] -= (sea_level_ref + sea_level)
n = noisemap(mapsize+1, mapsize+1, **params)
### COMPUTE LANDSCAPE EVOLUTION
# Initialize landscape evolution model
print('Initializing model')
model = terrainlib.EvolutionModel(n, K=K, m=m, d=d, sea_level=sea_level, flex_radius=flex_radius, flow_method=flow_method)
terrainlib.update(model.dem, model.lakes, t=5, sea_level=model.sea_level, title='Initializing...')
dt = time/niter
# Run the model's processes: the order in which the processes are run is arbitrary and could be changed.
for i in range(niter):
disp_niter = 'Iteration {:d} of {:d}...'.format(i+1, niter)
if sea_level_variations != 0:
model.sea_level = snoise2((i*dt)/sea_level_variations_time, sea_ybase, **params_sealevel) * sea_level_variations - sea_level_ref
terrainlib.update(model.dem, model.lakes, sea_level=model.sea_level, title=disp_niter)
print(disp_niter)
print('Diffusion')
model.diffusion(dt)
print('Flow calculation')
model.calculate_flow()
terrainlib.update(model.dem, model.lakes, sea_level=model.sea_level, title=disp_niter)
print('Advection')
model.advection(dt)
print('Isostatic equilibration')
model.adjust_isostasy()
print('Last flow calculation')
model.calculate_flow()
print('Done!')
# Twist the grid
bx, by = terrainlib.make_bounds(model.dirs, model.rivers)
offset_x, offset_y = terrainlib.twist(bx, by, terrainlib.get_fixed(model.dirs))
# Convert offset in 8-bits
offset_x = np.clip(np.floor(offset_x * 256), -128, 127)
offset_y = np.clip(np.floor(offset_y * 256), -128, 127)
### SAVE OUTPUT
if not os.path.isdir(output_dir):
os.mkdir(output_dir)
os.chdir(output_dir)
# Save the files
terrainlib.save(model.dem, 'dem', dtype='>i2')
terrainlib.save(model.lakes, 'lakes', dtype='>i2')
terrainlib.save(offset_x, 'offset_x', dtype='i1')
terrainlib.save(offset_y, 'offset_y', dtype='i1')
terrainlib.save(model.dirs, 'dirs', dtype='u1')
terrainlib.save(model.rivers, 'rivers', dtype='>u4')
with open('size', 'w') as sfile:
sfile.write('{:d}\n{:d}'.format(mapsize+1, mapsize+1))
terrainlib.stats(model.dem, model.lakes)
print()
print('Grid is ready for use!')
terrainlib.plot(model.dem, model.lakes, title='Final grid, ready for use!')

View File

@ -1,7 +0,0 @@
# Load packages and provide easy access to important functions
from .settings import read_config_file
from .erosion import EvolutionModel
from .save import save
from .bounds import make_bounds, twist, get_fixed
from .view import stats, update, plot

View File

@ -1,74 +0,0 @@
import numpy as np
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=rivers.dtype)
bounds_v = np.zeros((Y-1, X), dtype=rivers.dtype)
bounds_v += (rivers * (dirs==1))[:-1,:]
bounds_h += (rivers * (dirs==2))[:,:-1]
bounds_v -= (rivers * (dirs==3))[1:,:]
bounds_h -= (rivers * (dirs==4))[:,1:]
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
borders[0,:] |= dirs[0,:]==3
borders[:,0] |= dirs[:,0]==4
donors = np.zeros(dirs.shape, dtype='?')
donors[1:,:] |= dirs[:-1,:]==1
donors[:,1:] |= dirs[:,:-1]==2
donors[:-1,:] |= dirs[1:,:]==3
donors[:,:-1] |= dirs[:,1:]==4
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
offset_x = np.zeros((Y, X))
offset_y = np.zeros((Y, X))
for i in range(n):
force_long = np.abs(bounds_x) * (1+np.diff(offset_x, axis=1))
force_trans = np.abs(bounds_y) * np.diff(offset_x, axis=0)
force_x = np.zeros((Y, X))
force_x[:,:-1] = force_long
force_x[:,1:] -= force_long
force_x[:-1,:]+= force_trans
force_x[1:,:] -= force_trans
force_long = np.abs(bounds_y) * (1+np.diff(offset_y, axis=0))
force_trans = np.abs(bounds_x) * np.diff(offset_y, axis=1)
force_y = np.zeros((Y, X))
force_y[:-1,:] = force_long
force_y[1:,:] -= force_long
force_y[:,:-1]+= force_trans
force_y[:,1:] -= force_trans
length = np.hypot(force_x, force_y)
length[length==0] = 1
coeff = d / length * moveable # Normalize, take into account the direction only
offset_x += force_x * coeff
offset_y += force_y * coeff
return offset_x, offset_y

View File

@ -1,95 +0,0 @@
import numpy as np
import scipy.ndimage as im
from .rivermapper import flow
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
"""
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
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]: # 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] # If between 2 pixels, perform linear interpolation.
return dem_new
def diffusion(dem, time, d=1):
radius = d * time**.5
if radius == 0:
return dem
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, flow_method='semirandom'):
self.dem = dem
#self.bedrock = dem
self.K = K
self.m = m
self.d = d
self.sea_level = sea_level
self.flex_radius = flex_radius
self.define_isostasy()
self.flow_method = flow_method
#set_flow_method(flow_method)
if flow:
self.calculate_flow()
else:
self.lakes = dem
self.dirs = np.zeros(dem.shape, dtype=int)
self.rivers = np.zeros(dem.shape, dtype=int)
self.flow_uptodate = False
def calculate_flow(self):
self.dirs, self.lakes, self.rivers = flow(self.dem, method=self.flow_method)
self.flow_uptodate = True
def advection(self, time):
dem = advection(np.maximum(self.dem, 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
def define_isostasy(self):
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') # Calculate blurred DEM
correction = (self.ref_isostasy - isostasy) * rate # Compare it with the reference isostasy
self.dem = self.dem + correction # Adjust

View File

@ -1,278 +0,0 @@
import numpy as np
import numpy.random as npr
from collections import defaultdict
# This file provide functions to construct the river tree from an elevation model.
# Based on a research paper:
# | Cordonnier, G., Bovy, B., and Braun, J.:
# | A versatile, linear complexity algorithm for flow routing in topographies with depressions,
# | Earth Surf. Dynam., 7, 549562, https://doi.org/10.5194/esurf-7-549-2019, 2019.
# 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 Boruvka algorithm.
# Only flow_local and accumulate_flow are custom algorithms.
# Define two different method for local flow routing
def flow_local_steepest(plist):
vmax = 0.0
imax = 0.0
for i, p in enumerate(plist):
if p > vmax:
vmax = p
imax = i
if vmax > 0.0:
return imax+1
return 0
def flow_local_semirandom(plist):
"""
Determines a flow direction based on denivellation for every neighbouring node.
Denivellation must be positive for downward and zero for flat or upward:
dz = max(zref-z, 0)
"""
psum = sum(plist)
if psum == 0:
return 0
r = npr.random() * psum
for i, p in enumerate(plist):
if r < p:
return i+1
r -= p
flow_local_methods = {
'steepest' : flow_local_steepest,
'semirandom' : flow_local_semirandom,
}
def flow(dem, method='semirandom'):
if method in flow_local_methods:
flow_local = flow_local_methods[method]
else:
raise KeyError('Flow method \'{}\' does not exist'.format(method))
# Flow locally
dirs1 = np.zeros(dem.shape, dtype=int)
dirs2 = np.zeros(dem.shape, dtype=int)
(X, Y) = dem.shape
Xmax, Ymax = X-1, Y-1
singular = []
for x in range(X):
z0 = z1 = z2 = dem[x,0]
for y in range(Y):
z0 = z1
z1 = z2
if y < Ymax:
z2 = dem[x, y+1]
plist = [
max(z1-dem[x+1,y],0) if x<Xmax else 0, # 1: x -> x+1
max(z1-z2,0), # 2: y -> y+1
max(z1-dem[x-1,y],0) if x>0 else 0, # 3: x -> x-1
max(z1-z0,0), # 4: y -> y-1
]
pdir = flow_local(plist)
dirs2[x,y] = pdir
if pdir == 0:
singular.append((x,y))
elif pdir == 1:
dirs1[x+1,y] += 1
elif pdir == 2:
dirs1[x,y+1] += 2
elif pdir == 3:
dirs1[x-1,y] += 4
elif pdir == 4:
dirs1[x,y-1] += 8
# Compute basins
basin_id = np.zeros(dem.shape, dtype=int)
stack = []
for i, s in enumerate(singular):
queue = [s]
while queue:
x, y = queue.pop()
basin_id[x,y] = i
d = int(dirs1[x,y])
if d & 1:
queue.append((x-1,y))
if d & 2:
queue.append((x,y-1))
if d & 4:
queue.append((x+1,y))
if d & 8:
queue.append((x,y+1))
del dirs1
# Link basins
nsing = len(singular)
links = {}
def add_link(b0, b1, elev, bound):
b = (min(b0,b1),max(b0,b1))
if b not in links or links[b][0] > elev:
links[b] = (elev, bound)
for x in range(X):
b0 = basin_id[x,0]
add_link(-1, b0, dem[x,0], (True, x, 0))
for y in range(1,Y):
b1 = basin_id[x,y]
if b0 != b1:
add_link(b0, b1, max(dem[x,y-1],dem[x,y]), (True, x, y))
b0 = b1
add_link(-1, b1, dem[x,Ymax], (True, x, Y))
for y in range(Y):
b0 = basin_id[0,y]
add_link(-1, b0, dem[0,y], (False, 0, y))
for x in range(1,X):
b1 = basin_id[x,y]
if b0 != b1:
add_link(b0, b1, max(dem[x-1,y],dem[x,y]), (False, x, y))
b0 = b1
add_link(-1, b1, dem[Xmax,y], (False, X, y))
# Computing basin tree
graph = planar_boruvka(links)
basin_links = defaultdict(dict)
for elev, b1, b2, bound in graph:
basin_links[b1][b2] = basin_links[b2][b1] = (elev, bound)
basins = np.zeros(nsing+1)
stack = [(-1, float('-inf'))]
# Applying basin flowing
dir_reverse = (0, 3, 4, 1, 2)
while stack:
b1, elev1 = stack.pop()
basins[b1] = elev1
for b2, (elev2, bound) in basin_links[b1].items():
stack.append((b2, max(elev1, elev2)))
# Reverse flow direction in b2 (TODO)
isY, x, y = bound
backward = True # Whether water will escape the basin in +X/+Y direction
if not (x < X and y < Y and basin_id[x,y] == b2):
if isY:
y -= 1
else:
x -= 1
backward = False
d = 2*backward + isY + 1
while d > 0:
d, dirs2[x,y] = dirs2[x,y], d
if d == 1:
x += 1
elif d == 2:
y += 1
elif d == 3:
x -= 1
elif d == 4:
y -= 1
d = dir_reverse[d]
del basin_links[b2][b1]
del basin_links[b1]
# Calculating water quantity
dirs2[-1,:][dirs2[-1,:]==1] = 0
dirs2[:,-1][dirs2[:,-1]==2] = 0
dirs2[0,:][dirs2[0,:]==3] = 0
dirs2[:,0][dirs2[:,0]==4] = 0
waterq = accumulate_flow(dirs2)
return dirs2, basins[basin_id], waterq
def accumulate_flow(dirs):
ndonors = np.zeros(dirs.shape, dtype=int)
ndonors[1:,:] += dirs[:-1,:] == 1
ndonors[:,1:] += dirs[:,:-1] == 2
ndonors[:-1,:] += dirs[1:,:] == 3
ndonors[:,:-1] += dirs[:,1:] == 4
waterq = np.ones(dirs.shape, dtype=int)
(X, Y) = dirs.shape
rangeX = range(X)
rangeY = range(Y)
for x in rangeX:
for y in rangeY:
if ndonors[x,y] > 0:
continue
xw, yw = x, y
w = waterq[xw, yw]
while 1:
d = dirs[xw, yw]
if d <= 0:
break
elif d == 1:
xw += 1
elif d == 2:
yw += 1
elif d == 3:
xw -= 1
elif d == 4:
yw -= 1
w += waterq[xw, yw]
waterq[xw, yw] = w
if ndonors[xw, yw] > 1:
ndonors[xw, yw] -= 1
break
return waterq
def planar_boruvka(links):
# Compute basin tree
basin_list = defaultdict(dict)
for (b1, b2), (elev, bound) in links.items():
basin_list[b1][b2] = basin_list[b2][b1] = (elev, b1, b2, bound)
threshold = 8
lowlevel = {}
for k, v in basin_list.items():
if len(v) <= threshold:
lowlevel[k] = v
basin_graph = []
n = len(basin_list)
while n > 1:
(b1, lnk1) = lowlevel.popitem()
b2 = min(lnk1, key=lnk1.get)
lnk2 = basin_list[b2]
# Add link to the graph
basin_graph.append(lnk1[b2])
# Union : merge basin 1 into basin 2
# First, delete the direct link
del lnk1[b2]
del lnk2[b1]
# Look for basin 1's neighbours, and add them to basin 2 if they have a lower pass
for k, v in lnk1.items():
bk = basin_list[k]
if k in lnk2 and lnk2[k] < v:
del bk[b1]
else:
lnk2[k] = v
bk[b2] = bk.pop(b1)
if k not in lowlevel and len(bk) <= threshold:
lowlevel[k] = bk
if b2 in lowlevel:
if len(lnk2) > threshold:
del lowlevel[b2]
elif len(lnk2) <= threshold:
lowlevel[b2] = lnk2
del lnk1
n -= 1
return basin_graph

View File

@ -1,13 +0,0 @@
import numpy as np
import zlib
def save(data, fname, dtype=None):
if dtype is not None:
data = data.astype(dtype)
bin_data = data.tobytes()
bin_data_comp = zlib.compress(bin_data, 9)
if len(bin_data_comp) < len(bin_data):
bin_data = bin_data_comp
with open(fname, 'wb') as f:
f.write(bin_data)

View File

@ -1,16 +0,0 @@
import os.path
def read_config_file(fname):
settings = {}
if not os.path.isfile(fname):
return 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

View File

@ -5,7 +5,7 @@ import zlib
import sys import sys
import os import os
from terrainlib import stats, plot from view import stats, plot
scale = 1 scale = 1
if len(sys.argv) > 1: if len(sys.argv) > 1: