Implement simple tectonics.

Uplift and subsidence are determined with a noise, at every iteration.
There is no distinctive pattern like tectonic plates, just vertical movements
disturbing rivers from their equilibrium state, and thus creating more diversity.
More lakes and waterfalls especially.

generate.py

@@ -111,6 +111,7 @@ 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))
+tectonics_time = float(get_setting('tectonics_time', 0.0))
 flow_method = get_setting('flow_method', 'semirandom')
 
 time = float(get_setting('time', 10.0))
@@ -160,7 +161,11 @@ if sea_level_variations != 0.0:
     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).get2d()
+if tectonics_time == 0.0:
+    n = noisemap(mapsize+1, mapsize+1, **params).get2d()
+else:
+    terrain_noisemap = noisemap(mapsize+1, mapsize+1, tscale=tectonics_time, **params)
+    n = terrain_noisemap.get3d()
 m_map = noisemap(mapsize+1, mapsize+1, **params_m).get2d()
 K_map = noisemap(mapsize+1, mapsize+1, **params_K).get2d()
 
@@ -187,6 +192,9 @@ for i in range(niter):
     terrainlib.update(model.dem, model.lakes, sea_level=model.sea_level, title=disp_niter)
     print('Advection')
     model.advection(dt)
+    if tectonics_time != 0.0:
+        print('Isostasy reference redefinition')
+        model.define_isostasy(terrain_noisemap.get3d((i+1)*dt))
     print('Isostatic equilibration')
     model.adjust_isostasy()
 

terrainlib/erosion.py

@@ -86,8 +86,10 @@ class EvolutionModel:
         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 define_isostasy(self, dem=None):
+        if dem is None:
+            dem = self.dem
+        self.ref_isostasy = im.gaussian_filter(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