diff --git a/generate.py b/generate.py deleted file mode 100755 index b3a1734..0000000 --- a/generate.py +++ /dev/null @@ -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!') diff --git a/terrainlib/__init__.py b/terrainlib/__init__.py deleted file mode 100644 index 8775f17..0000000 --- a/terrainlib/__init__.py +++ /dev/null @@ -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 diff --git a/terrainlib/bounds.py b/terrainlib/bounds.py deleted file mode 100644 index b179a1f..0000000 --- a/terrainlib/bounds.py +++ /dev/null @@ -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 diff --git a/terrainlib/erosion.py b/terrainlib/erosion.py deleted file mode 100644 index 9d2ed89..0000000 --- a/terrainlib/erosion.py +++ /dev/null @@ -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 diff --git a/terrainlib/rivermapper.py b/terrainlib/rivermapper.py deleted file mode 100644 index 3c93bd5..0000000 --- a/terrainlib/rivermapper.py +++ /dev/null @@ -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, 549–562, 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 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 diff --git a/terrainlib/save.py b/terrainlib/save.py deleted file mode 100644 index b21b524..0000000 --- a/terrainlib/save.py +++ /dev/null @@ -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) diff --git a/terrainlib/settings.py b/terrainlib/settings.py deleted file mode 100644 index c63b015..0000000 --- a/terrainlib/settings.py +++ /dev/null @@ -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 diff --git a/terrainlib/view.py b/view.py similarity index 100% rename from terrainlib/view.py rename to view.py diff --git a/view_map.py b/view_map.py index c3efe6d..a1216b4 100755 --- a/view_map.py +++ b/view_map.py @@ -5,7 +5,7 @@ import zlib import sys import os -from terrainlib import stats, plot +from view import stats, plot scale = 1 if len(sys.argv) > 1: