From f0dddee33c20d661c20b00f4b6779f93efc90509 Mon Sep 17 00:00:00 2001
From: Gael-de-Sailly <gael-de-sailly@netc.eu>
Date: Sat, 5 Dec 2020 14:24:50 +0100
Subject: [PATCH] Lakes map: keep initial height (reduces file size) Lake
 height is calculated for every basin, and there is a lake if lake height is
 higher than ground height. If it is lower, there is no lake. In that case, it
 was previously raised to ground level, but since this can be done in Lua, we
 can write initial lakes height in the files. This has the advantage of
 reducing file size, since there are bigger areas of equal values, that are
 more efficiently compressed.

---
 terrainlib/erosion.py     | 10 ++--------
 terrainlib/rivermapper.py |  2 +-
 terrainlib/view.py        | 10 +++++-----
 3 files changed, 8 insertions(+), 14 deletions(-)

diff --git a/terrainlib/erosion.py b/terrainlib/erosion.py
index c8b2dc2..fbed7e6 100644
--- a/terrainlib/erosion.py
+++ b/terrainlib/erosion.py
@@ -11,12 +11,6 @@ def advection(dem, dirs, rivers, time, K=1, m=0.5, sea_level=0):
     v = K * flux^m
     """
 
-    dirs = dirs.copy()
-    dirs[0,:] = 0
-    dirs[-1,:] = 0
-    dirs[:,0] = 0
-    dirs[:,-1] = 0
-
     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)
@@ -51,7 +45,7 @@ def advection(dem, dirs, rivers, time, K=1, m=0.5, sea_level=0):
             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 np.minimum(dem, dem_new)
+    return dem_new
 
 def diffusion(dem, time, d=1):
     radius = d * time**.5
@@ -80,7 +74,7 @@ class EvolutionModel:
         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)
+        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
 
diff --git a/terrainlib/rivermapper.py b/terrainlib/rivermapper.py
index 687c14d..9e09d9d 100644
--- a/terrainlib/rivermapper.py
+++ b/terrainlib/rivermapper.py
@@ -153,7 +153,7 @@ def flow(dem):
 
     waterq = accumulate_flow(dirs2)
 
-    return dirs2, np.maximum(basins[basin_id], dem), waterq
+    return dirs2, basins[basin_id], waterq
 
 def accumulate_flow(dirs):
     ndonors = np.zeros(dirs.shape, dtype=int)
diff --git a/terrainlib/view.py b/terrainlib/view.py
index d3d9c78..9b03820 100644
--- a/terrainlib/view.py
+++ b/terrainlib/view.py
@@ -22,11 +22,11 @@ if has_matplotlib:
     def view_map(dem, lakes, scale=1, title=None):
         lakes_sea = np.maximum(lakes, 0)
         water = np.maximum(lakes_sea - dem, 0)
-        max_elev = lakes_sea.max()
+        max_elev = dem.max()
         max_depth = water.max()
 
         ls = mcl.LightSource(azdeg=315, altdeg=45)
-        rgb = ls.shade(lakes_sea, cmap=cmap1, vert_exag=1/scale, blend_mode='soft', vmin=0, vmax=max_elev)
+        rgb = ls.shade(dem, cmap=cmap1, vert_exag=1/scale, blend_mode='soft', vmin=0, vmax=max_elev)
 
         (X, Y) = dem.shape
         extent = (0, Y*scale, 0, X*scale)
@@ -69,13 +69,13 @@ else:
     def plot(*args, **kwargs):
         pass
 
-def stats(dem, lake_dem, scale=1):
+def stats(dem, lakes, scale=1):
     surface = dem.size
 
-    continent = lake_dem >= 0
+    continent = np.maximum(dem, lakes) >= 0
     continent_surface = continent.sum()
 
-    lake = continent & (lake_dem>dem)
+    lake = continent & (lakes>dem)
     lake_surface = lake.sum()
 
     print('---   General    ---')