Added optional sea level variations for the simulation.

This results in more varied coastline morphologies.

README.md

@@ -76,6 +76,8 @@ Other parameters can be specified by `--parameter value`. Syntax `--parameter=va
 | `flex_radius` | Flexure radius. Wavelength over which loss/gain of mass is compensated by uplift/subsidence. This ensures that mountain ranges will not get eventually flattened by erosion, and that an equilibrium is reached. Geologically speaking, this implements [isostatic rebound](https://en.wikipedia.org/wiki/Isostasy). | `--flex_radius 20` |
 | `time`        | Simulated time of erosion modelling, in abstract units. | `--time 10` |
 | `niter`       | Number of iterations. Each iteration represents a time `time/niter`. | `--niter 10` |
+| `sea_level_variations` | Amplitude of sea level variations throughout the simulation (if any). | `--sea_level_variations 10` |
+| `sea_level_variations_time` | Characteristic time of variation for sea level, in the same units than `time`. Increasing it will result in slower variations between iterations. | `--sea_level_variations_time 1` |
 | | **Alternatives** |
 | `config`      | Another way to specify configuration file | `--config terrain_higher.conf` |
 | `output`      | Another way to specify output dir | `--output ~/.minetest/worlds/my_world/river_data` |

generate.py

@@ -65,6 +65,8 @@ K = float(get_setting('K', 1.0))
 m = float(get_setting('m', 0.35))
 d = float(get_setting('d', 0.2))
 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))
 
 time = float(get_setting('time', 10.0))
@@ -80,9 +82,19 @@ params = {
     "lacunarity" : lacunarity,
 }
 
+params_sealevel = {
+    "octaves" : 1,
+    "persistence" : 1,
+    "lacunarity" : 2,
+}
+
 # Determine noise offset randomly
 xbase = np.random.randint(8192)-4096
 ybase = np.random.randint(8192)-4096
+if sea_level_variations != 0.0:
+    sea_ybase = np.random.randint(8192)-4096
+    sea_level_ref = noise.snoise2(time * (1-1/niter) / sea_level_variations, sea_ybase, **params_sealevel) * sea_level_variations
+    offset -= (sea_level_ref + sea_level)
 
 # Generate the noise
 for x in range(mapsize+1):
@@ -95,7 +107,7 @@ nn = n*vscale + offset
 # Initialize landscape evolution model
 print('Initializing model')
 model = terrainlib.EvolutionModel(nn, K=K, m=m, d=d, sea_level=sea_level, flex_radius=flex_radius)
-terrainlib.update(model.dem, model.lakes, t=5, title='Initializing...')
+terrainlib.update(model.dem, model.lakes, t=5, sea_level=model.sea_level, title='Initializing...')
 
 dt = time/niter
 
@@ -103,13 +115,15 @@ dt = time/niter
 
 for i in range(niter):
     disp_niter = 'Iteration {:d} of {:d}...'.format(i+1, niter)
-    terrainlib.update(model.dem, model.lakes, title=disp_niter)
+    if sea_level_variations != 0:
+        model.sea_level = noise.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, title=disp_niter)
+    terrainlib.update(model.dem, model.lakes, sea_level=model.sea_level, title=disp_niter)
     print('Advection')
     model.advection(dt)
     print('Isostatic equilibration')

terrainlib/erosion.py

@@ -49,6 +49,8 @@ def advection(dem, dirs, rivers, time, K=1, m=0.5, sea_level=0):
 
 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:

terrainlib/view.py

@@ -19,14 +19,16 @@ except ImportError: # No module matplotlib
     has_matplotlib = False
 
 if has_matplotlib:
-    def view_map(dem, lakes, scale=1, title=None):
-        lakes_sea = np.maximum(lakes, 0)
+    def view_map(dem, lakes, scale=1, sea_level=0.0, title=None):
+        lakes_sea = np.maximum(lakes, sea_level)
         water = np.maximum(lakes_sea - dem, 0)
         max_elev = dem.max()
         max_depth = water.max()
 
         ls = mcl.LightSource(azdeg=315, altdeg=45)
-        rgb = ls.shade(dem, cmap=cmap1, vert_exag=1/scale, blend_mode='soft', vmin=0, vmax=max_elev)
+        norm_ground = plt.Normalize(vmin=sea_level, vmax=max_elev)
+        norm_sea = plt.Normalize(vmin=0, vmax=max_depth)
+        rgb = ls.shade(dem, cmap=cmap1, vert_exag=1/scale, blend_mode='soft', norm=norm_ground)
 
         (X, Y) = dem.shape
         extent = (0, Y*scale, 0, X*scale)
@@ -34,10 +36,10 @@ if has_matplotlib:
         alpha = (water > 0).astype('u1')
         plt.imshow(np.flipud(water), alpha=np.flipud(alpha), cmap=cmap2, extent=extent, vmin=0, vmax=max_depth, interpolation='antialiased')
 
-        sm1 = plt.cm.ScalarMappable(cmap=cmap1, norm=plt.Normalize(vmin=0, vmax=max_elev))
+        sm1 = plt.cm.ScalarMappable(cmap=cmap1, norm=norm_ground)
         plt.colorbar(sm1).set_label('Elevation')
 
-        sm2 = plt.cm.ScalarMappable(cmap=cmap2, norm=plt.Normalize(vmin=0, vmax=max_depth))
+        sm2 = plt.cm.ScalarMappable(cmap=cmap2, norm=norm_sea)
         plt.colorbar(sm2).set_label('Water depth')
 
         plt.xlabel('X')