ant_landscape/eroder.py

   1 # SPDX-License-Identifier: GPL-2.0-or-later
   2
   3 from time import time
   4 import unittest
   5 import sys
   6 import os
   7 from random import random as rand, shuffle
   8 import numpy as np
   9
  10 numexpr_available = False
  11
  12
  13 def getmemsize():
  14     return 0.0
  15
  16
  17 def getptime():
  18     return time()
  19
  20
  21 class Grid:
  22
  23     def __init__(self, size=10, dtype=np.single):
  24         self.center = np.zeros([size, size], dtype)
  25         self.water = None
  26         self.sediment = None
  27         self.scour = None
  28         self.flowrate = None
  29         self.sedimentpct = None
  30         self.sedimentpct = None
  31         self.capacity = None
  32         self.avalanced = None
  33         self.minx = None
  34         self.miny = None
  35         self.maxx = None
  36         self.maxy = None
  37         self.zscale = 1
  38         self.maxrss = 0.0
  39         self.sequence = [0, 1, 2, 3]
  40         self.watermax = 1.0
  41         self.flowratemax = 1.0
  42         self.scourmax = 1.0
  43         self.sedmax = 1.0
  44         self.scourmin = 1.0
  45
  46     def init_water_and_sediment(self):
  47         if self.water is None:
  48             self.water = np.zeros(self.center.shape, dtype=np.single)
  49         if self.sediment is None:
  50             self.sediment = np.zeros(self.center.shape, dtype=np.single)
  51         if self.scour is None:
  52             self.scour = np.zeros(self.center.shape, dtype=np.single)
  53         if self.flowrate is None:
  54             self.flowrate = np.zeros(self.center.shape, dtype=np.single)
  55         if self.sedimentpct is None:
  56             self.sedimentpct = np.zeros(self.center.shape, dtype=np.single)
  57         if self.capacity is None:
  58             self.capacity = np.zeros(self.center.shape, dtype=np.single)
  59         if self.avalanced is None:
  60             self.avalanced = np.zeros(self.center.shape, dtype=np.single)
  61
  62     def __str__(self):
  63         return ''.join(self.__str_iter__(fmt="%.3f"))
  64
  65     def __str_iter__(self, fmt):
  66         for row in self.center[::]:
  67             values = []
  68             for v in row:
  69                 values.append(fmt % v)
  70             yield ' '.join(values) + '\n'
  71
  72     @staticmethod
  73     def fromFile(filename):
  74         if filename == '-':
  75             filename = sys.stdin
  76         g = Grid()
  77         g.center = np.loadtxt(filename, np.single)
  78         return g
  79
  80     def toFile(self, filename, fmt="%.3f"):
  81         if filename == '-':
  82             filename = sys.stdout.fileno()
  83         with open(filename, "w") as f:
  84             for line in self.__str_iter__(fmt):
  85                 f.write(line)
  86
  87     def raw(self, format="%.3f"):
  88         fstr = format + " " + format + " " + format + " "
  89         a = self.center / self.zscale
  90         minx = 0.0 if self.minx is None else self.minx
  91         miny = 0.0 if self.miny is None else self.miny
  92         maxx = 1.0 if self.maxx is None else self.maxx
  93         maxy = 1.0 if self.maxy is None else self.maxy
  94         dx = (maxx - minx) / (a.shape[0] - 1)
  95         dy = (maxy - miny) / (a.shape[1] - 1)
  96         for row in range(a.shape[0] - 1):
  97             row0 = miny + row * dy
  98             row1 = row0 + dy
  99             for col in range(a.shape[1] - 1):
 100                 col0 = minx + col * dx
 101                 col1 = col0 + dx
 102                 yield (fstr % (row0, col0, a[row][col]) +
 103                        fstr % (row0, col1, a[row][col + 1]) +
 104                        fstr % (row1, col0, a[row + 1][col]) + "\n")
 105                 yield (fstr % (row0, col1, a[row][col + 1]) +
 106                        fstr % (row1, col0, a[row + 1][col]) +
 107                        fstr % (row1, col1, a[row + 1][col + 1]) + "\n")
 108
 109     def toRaw(self, filename, infomap=None):
 110         with open(filename if type(filename) == str else sys.stdout.fileno(), "w") as f:
 111             f.writelines(self.raw())
 112         if infomap:
 113             with open(os.path.splitext(filename)[0] + ".inf" if type(filename) == str else sys.stdout.fileno(), "w") as f:
 114                 f.writelines("\n".join("%-15s: %s" % t for t in sorted(infomap.items())))
 115
 116     @staticmethod
 117     def fromRaw(filename):
 118         """initialize a grid from a Blender .raw file.
 119         currently supports just rectangular grids of all triangles
 120         """
 121         g = Grid.fromFile(filename)
 122         # we assume tris and an axis aligned grid
 123         g.center = np.reshape(g.center, (-1, 3))
 124         g._sort()
 125         return g
 126
 127     def _sort(self, expfact):
 128         # keep unique vertices only by creating a set and sort first on x then on y coordinate
 129         # using rather slow python sort but couldn't wrap my head around np.lexsort
 130         verts = sorted(list({tuple(t) for t in self.center[::]}))
 131         x = set(c[0] for c in verts)
 132         y = set(c[1] for c in verts)
 133         nx = len(x)
 134         ny = len(y)
 135         self.minx = min(x)
 136         self.maxx = max(x)
 137         self.miny = min(y)
 138         self.maxy = max(y)
 139         xscale = (self.maxx - self.minx) / (nx - 1)
 140         yscale = (self.maxy - self.miny) / (ny - 1)
 141         # note: a purely flat plane cannot be scaled
 142         if (yscale != 0.0) and (abs(xscale / yscale) - 1.0 > 1e-3):
 143             raise ValueError("Mesh spacing not square %d x %d  %.4f x %4.f" % (nx, ny, xscale, yscale))
 144         self.zscale = 1.0
 145         if abs(yscale) > 1e-6:
 146             self.zscale = 1.0 / yscale
 147
 148         # keep just the z-values and null any offset
 149         # we might catch a reshape error that will occur if nx*ny != # of vertices
 150         # (if we are not dealing with a heightfield but with a mesh with duplicate
 151         # x,y coords, like an axis aligned cube
 152         self.center = np.array([c[2] for c in verts], dtype=np.single).reshape(nx, ny)
 153         self.center = (self.center - np.amin(self.center)) * self.zscale
 154         if self.rainmap is not None:
 155             rmscale = np.max(self.center)
 156             self.rainmap = expfact + (1 - expfact) * (self.center / rmscale)
 157
 158     @staticmethod
 159     def fromBlenderMesh(me, vg, expfact):
 160         g = Grid()
 161         g.center = np.asarray(list(tuple(v.co) for v in me.vertices), dtype=np.single)
 162         g.rainmap = None
 163         if vg is not None:
 164             for v in me.vertices:
 165                 vg.add([v.index], 0.0, 'ADD')
 166             g.rainmap = np.asarray(list((v.co[0], v.co[1], vg.weight(v.index)) for v in me.vertices), dtype=np.single)
 167         g._sort(expfact)
 168         return g
 169
 170     def setrainmap(self, rainmap):
 171         self.rainmap = rainmap
 172
 173     def _verts(self, surface):
 174         a = surface / self.zscale
 175         minx = 0.0 if self.minx is None else self.minx
 176         miny = 0.0 if self.miny is None else self.miny
 177         maxx = 1.0 if self.maxx is None else self.maxx
 178         maxy = 1.0 if self.maxy is None else self.maxy
 179         dx = (maxx - minx) / (a.shape[0] - 1)
 180         dy = (maxy - miny) / (a.shape[1] - 1)
 181         for row in range(a.shape[0]):
 182             row0 = miny + row * dy
 183             for col in range(a.shape[1]):
 184                 col0 = minx + col * dx
 185                 yield (row0, col0, a[row][col])
 186
 187     def _faces(self):
 188         nrow, ncol = self.center.shape
 189         for row in range(nrow - 1):
 190             for col in range(ncol - 1):
 191                 vi = row * ncol + col
 192                 yield (vi, vi + ncol, vi + 1)
 193                 yield (vi + 1, vi + ncol, vi + ncol + 1)
 194
 195     def toBlenderMesh(self, me):
 196         # pass me as argument so that we don't need to import bpy and create a dependency
 197         # the docs state that from_pydata takes iterators as arguments but it will
 198         # fail with generators because it does len(arg)
 199         me.from_pydata(list(self._verts(self.center)), [], list(self._faces()))
 200
 201     def toWaterMesh(self, me):
 202         # pass me as argument so that we don't need to import bpy and create a dependency
 203         # the docs state that from_pydata takes iterators as arguments but it will
 204         # fail with generators because it does len(arg)
 205         me.from_pydata(list(self._verts(self.water)), [], list(self._faces()))
 206
 207     def peak(self, value=1):
 208         nx, ny = self.center.shape
 209         self.center[int(nx / 2), int(ny / 2)] += value
 210
 211     def shelf(self, value=1):
 212         nx, ny = self.center.shape
 213         self.center[:nx / 2] += value
 214
 215     def mesa(self, value=1):
 216         nx, ny = self.center.shape
 217         self.center[nx / 4:3 * nx / 4, ny / 4:3 * ny / 4] += value
 218
 219     def random(self, value=1):
 220         self.center += np.random.random_sample(self.center.shape) * value
 221
 222     def neighborgrid(self):
 223         self.up = np.roll(self.center, -1, 0)
 224         self.down = np.roll(self.center, 1, 0)
 225         self.left = np.roll(self.center, -1, 1)
 226         self.right = np.roll(self.center, 1, 1)
 227
 228     def zeroedge(self, quantity=None):
 229         c = self.center if quantity is None else quantity
 230         c[0, :] = 0
 231         c[-1, :] = 0
 232         c[:, 0] = 0
 233         c[:, -1] = 0
 234
 235     def diffuse(self, Kd, IterDiffuse, numexpr):
 236         self.zeroedge()
 237         c = self.center[1:-1, 1:-1]
 238         up = self.center[:-2, 1:-1]
 239         down = self.center[2:, 1:-1]
 240         left = self.center[1:-1, :-2]
 241         right = self.center[1:-1, 2:]
 242         if(numexpr and numexpr_available):
 243             self.center[1:-1, 1:-1] = ne.evaluate('c + Kd * (up + down + left + right - 4.0 * c)')
 244         else:
 245             self.center[1:-1, 1:-1] = c + (Kd / IterDiffuse) * (up + down + left + right - 4.0 * c)
 246         self.maxrss = max(getmemsize(), self.maxrss)
 247         return self.center
 248
 249     def avalanche(self, delta, iterava, prob, numexpr):
 250         self.zeroedge()
 251         c = self.center[1:-1, 1:-1]
 252         up = self.center[:-2, 1:-1]
 253         down = self.center[2:, 1:-1]
 254         left = self.center[1:-1, :-2]
 255         right = self.center[1:-1, 2:]
 256         where = np.where
 257
 258         if(numexpr and numexpr_available):
 259             self.center[1:-1, 1:-1] = ne.evaluate('c + where((up   -c) > delta ,(up   -c -delta)/2, 0) \
 260                  + where((down -c) > delta ,(down -c -delta)/2, 0)  \
 261                  + where((left -c) > delta ,(left -c -delta)/2, 0)  \
 262                  + where((right-c) > delta ,(right-c -delta)/2, 0)  \
 263                  + where((up   -c) < -delta,(up   -c +delta)/2, 0)  \
 264                  + where((down -c) < -delta,(down -c +delta)/2, 0)  \
 265                  + where((left -c) < -delta,(left -c +delta)/2, 0)  \
 266                  + where((right-c) < -delta,(right-c +delta)/2, 0)')
 267         else:
 268             sa = (
 269                 # incoming
 270                 where((up - c) > delta, (up - c - delta) / 2, 0)
 271                 + where((down - c) > delta, (down - c - delta) / 2, 0)
 272                 + where((left - c) > delta, (left - c - delta) / 2, 0)
 273                 + where((right - c) > delta, (right - c - delta) / 2, 0)
 274                 # outgoing
 275                 + where((up - c) < -delta, (up - c + delta) / 2, 0)
 276                 + where((down - c) < -delta, (down - c + delta) / 2, 0)
 277                 + where((left - c) < -delta, (left - c + delta) / 2, 0)
 278                 + where((right - c) < -delta, (right - c + delta) / 2, 0)
 279             )
 280             randarray = np.random.randint(0, 100, sa.shape) * 0.01
 281             sa = where(randarray < prob, sa, 0)
 282             self.avalanced[1:-1, 1:-1] = self.avalanced[1:-1, 1:-1] + sa / iterava
 283             self.center[1:-1, 1:-1] = c + sa / iterava
 284
 285         self.maxrss = max(getmemsize(), self.maxrss)
 286         return self.center
 287
 288     def rain(self, amount=1, variance=0, userainmap=False):
 289         self.water += (1.0 - np.random.random(self.water.shape) * variance) * \
 290             (amount if ((self.rainmap is None) or (not userainmap)) else self.rainmap * amount)
 291
 292     def spring(self, amount, px, py, radius):
 293         # px, py and radius are all fractions
 294         nx, ny = self.center.shape
 295         rx = max(int(nx * radius), 1)
 296         ry = max(int(ny * radius), 1)
 297         px = int(nx * px)
 298         py = int(ny * py)
 299         self.water[px - rx:px + rx + 1, py - ry:py + ry + 1] += amount
 300
 301     def river(self, Kc, Ks, Kdep, Ka, Kev, numexpr):
 302         zeros = np.zeros
 303         where = np.where
 304         min = np.minimum
 305         max = np.maximum
 306         abs = np.absolute
 307         arctan = np.arctan
 308         sin = np.sin
 309
 310         center = (slice(1, -1, None), slice(1, -1, None))
 311         up = (slice(None, -2, None), slice(1, -1, None))
 312         down = (slice(2, None, None), slice(1, -1, None))
 313         left = (slice(1, -1, None), slice(None, -2, None))
 314         right = (slice(1, -1, None), slice(2, None, None))
 315
 316         water = self.water
 317         rock = self.center
 318         sediment = self.sediment
 319         height = rock + water
 320
 321         # !! this gives a runtime warning for division by zero
 322         verysmallnumber = 0.0000000001
 323         water += verysmallnumber
 324         sc = where(water > verysmallnumber, sediment / water, 0)
 325
 326         sdw = zeros(water[center].shape)
 327         svdw = zeros(water[center].shape)
 328         sds = zeros(water[center].shape)
 329         angle = zeros(water[center].shape)
 330         for d in (up, down, left, right):
 331             if(numexpr and numexpr_available):
 332                 hdd = height[d]
 333                 hcc = height[center]
 334                 dw = ne.evaluate('hdd-hcc')
 335                 inflow = ne.evaluate('dw > 0')
 336                 wdd = water[d]
 337                 wcc = water[center]
 338                 # nested where() represent min() and max()
 339                 dw = ne.evaluate('where(inflow, where(wdd<dw, wdd, dw), where(-wcc>dw, -wcc, dw))/4.0')
 340                 sdw = ne.evaluate('sdw + dw')
 341                 scd = sc[d]
 342                 scc = sc[center]
 343                 rockd = rock[d]
 344                 rockc = rock[center]
 345                 sds = ne.evaluate('sds + dw * where(inflow, scd, scc)')
 346                 svdw = ne.evaluate('svdw + abs(dw)')
 347                 angle = ne.evaluate('angle + arctan(abs(rockd-rockc))')
 348             else:
 349                 dw = (height[d] - height[center])
 350                 inflow = dw > 0
 351                 dw = where(inflow, min(water[d], dw), max(-water[center], dw)) / 4.0
 352                 sdw = sdw + dw
 353                 sds = sds + dw * where(inflow, sc[d], sc[center])
 354                 svdw = svdw + abs(dw)
 355                 angle = angle + np.arctan(abs(rock[d] - rock[center]))
 356
 357         if numexpr and numexpr_available:
 358             wcc = water[center]
 359             scc = sediment[center]
 360             rcc = rock[center]
 361             water[center] = ne.evaluate('wcc + sdw')
 362             sediment[center] = ne.evaluate('scc + sds')
 363             sc = ne.evaluate('where(wcc>0, scc/wcc, 2000*Kc)')
 364             fKc = ne.evaluate('Kc*sin(Ka*angle)*svdw')
 365             ds = ne.evaluate('where(sc > fKc, -Kd * scc, Ks * svdw)')
 366             rock[center] = ne.evaluate('rcc - ds')
 367             # there isn't really a bottom to the rock but negative values look ugly
 368             rock[center] = ne.evaluate('where(rcc<0,0,rcc)')
 369             sediment[center] = ne.evaluate('scc + ds')
 370         else:
 371             wcc = water[center]
 372             scc = sediment[center]
 373             rcc = rock[center]
 374             water[center] = wcc * (1 - Kev) + sdw
 375             sediment[center] = scc + sds
 376             sc = where(wcc > 0, scc / wcc, 2 * Kc)
 377             fKc = Kc * svdw
 378             ds = where(fKc > sc, (fKc - sc) * Ks, (fKc - sc) * Kdep) * wcc
 379             self.flowrate[center] = svdw
 380             self.scour[center] = ds
 381             self.sedimentpct[center] = sc
 382             self.capacity[center] = fKc
 383             sediment[center] = scc + ds + sds
 384
 385     def flow(self, Kc, Ks, Kz, Ka, numexpr):
 386         zeros = np.zeros
 387         where = np.where
 388         min = np.minimum
 389         max = np.maximum
 390         abs = np.absolute
 391         arctan = np.arctan
 392         sin = np.sin
 393
 394         center = (slice(1, -1, None), slice(1, -1, None))
 395         rock = self.center
 396         ds = self.scour[center]
 397         rcc = rock[center]
 398         rock[center] = rcc - ds * Kz
 399         # there isn't really a bottom to the rock but negative values look ugly
 400         rock[center] = where(rcc < 0, 0, rcc)
 401
 402     def rivergeneration(
 403             self,
 404             rainamount,
 405             rainvariance,
 406             userainmap,
 407             Kc,
 408             Ks,
 409             Kdep,
 410             Ka,
 411             Kev,
 412             Kspring,
 413             Kspringx,
 414             Kspringy,
 415             Kspringr,
 416             numexpr,
 417     ):
 418         self.init_water_and_sediment()
 419         self.rain(rainamount, rainvariance, userainmap)
 420         self.zeroedge(self.water)
 421         self.zeroedge(self.sediment)
 422         self.river(Kc, Ks, Kdep, Ka, Kev, numexpr)
 423         self.watermax = np.max(self.water)
 424
 425     def fluvial_erosion(
 426             self,
 427             rainamount,
 428             rainvariance,
 429             userainmap,
 430             Kc,
 431             Ks,
 432             Kdep,
 433             Ka,
 434             Kspring,
 435             Kspringx,
 436             Kspringy,
 437             Kspringr,
 438             numexpr,
 439     ):
 440         self.flow(Kc, Ks, Kdep, Ka, numexpr)
 441         self.flowratemax = np.max(self.flowrate)
 442         self.scourmax = np.max(self.scour)
 443         self.scourmin = np.min(self.scour)
 444         self.sedmax = np.max(self.sediment)
 445
 446     def analyze(self):
 447         self.neighborgrid()
 448         # just looking at up and left to avoid needless double calculations
 449         slopes = np.concatenate((np.abs(self.left - self.center), np.abs(self.up - self.center)))
 450         return '\n'.join(["%-15s: %.3f" % t for t in [
 451             ('height average', np.average(self.center)),
 452             ('height median', np.median(self.center)),
 453             ('height max', np.max(self.center)),
 454             ('height min', np.min(self.center)),
 455             ('height std', np.std(self.center)),
 456             ('slope average', np.average(slopes)),
 457             ('slope median', np.median(slopes)),
 458             ('slope max', np.max(slopes)),
 459             ('slope min', np.min(slopes)),
 460             ('slope std', np.std(slopes))
 461         ]])
 462
 463
 464 class TestGrid(unittest.TestCase):
 465
 466     def test_diffuse(self):
 467         g = Grid(5)
 468         g.peak(1)
 469         self.assertEqual(g.center[2, 2], 1.0)
 470         g.diffuse(0.1, numexpr=False)
 471         for n in [(2, 1), (2, 3), (1, 2), (3, 2)]:
 472             self.assertAlmostEqual(g.center[n], 0.1)
 473         self.assertAlmostEqual(g.center[2, 2], 0.6)
 474
 475     def test_diffuse_numexpr(self):
 476         g = Grid(5)
 477         g.peak(1)
 478         g.diffuse(0.1, numexpr=False)
 479         h = Grid(5)
 480         h.peak(1)
 481         h.diffuse(0.1, numexpr=True)
 482         self.assertEqual(list(g.center.flat), list(h.center.flat))
 483
 484     def test_avalanche_numexpr(self):
 485         g = Grid(5)
 486         g.peak(1)
 487         g.avalanche(0.1, numexpr=False)
 488         h = Grid(5)
 489         h.peak(1)
 490         h.avalanche(0.1, numexpr=True)
 491         print(g)
 492         print(h)
 493         np.testing.assert_almost_equal(g.center, h.center)
 494
 495
 496 if __name__ == "__main__":
 497
 498     import argparse
 499
 500     parser = argparse.ArgumentParser(description='Erode a terrain while assuming zero boundary conditions.')
 501     parser.add_argument('-I', dest='iterations', type=int, default=1, help='the number of iterations')
 502     parser.add_argument('-Kd', dest='Kd', type=float, default=0.01, help='Diffusion constant')
 503     parser.add_argument('-Kh', dest='Kh', type=float, default=6, help='Maximum stable cliff height')
 504     parser.add_argument('-Kp', dest='Kp', type=float, default=0.1, help='Avalanche probability for unstable cliffs')
 505     parser.add_argument('-Kr', dest='Kr', type=float, default=0.1, help='Average amount of rain per iteration')
 506     parser.add_argument('-Kspring', dest='Kspring', type=float, default=0.0,
 507                         help='Average amount of wellwater per iteration')
 508     parser.add_argument('-Kspringx', dest='Kspringx', type=float, default=0.5, help='relative x position of spring')
 509     parser.add_argument('-Kspringy', dest='Kspringy', type=float, default=0.5, help='relative y position of spring')
 510     parser.add_argument('-Kspringr', dest='Kspringr', type=float, default=0.02, help='radius of spring')
 511     parser.add_argument('-Kdep', dest='Kdep', type=float, default=0.1, help='Sediment deposition constant')
 512     parser.add_argument('-Ks', dest='Ks', type=float, default=0.1, help='Soil softness constant')
 513     parser.add_argument('-Kc', dest='Kc', type=float, default=1.0, help='Sediment capacity')
 514     parser.add_argument('-Ka', dest='Ka', type=float, default=2.0, help='Slope dependency of erosion')
 515     parser.add_argument(
 516         '-ri',
 517         action='store_true',
 518         dest='rawin',
 519         default=False,
 520         help='use Blender raw format for input')
 521     parser.add_argument(
 522         '-ro',
 523         action='store_true',
 524         dest='rawout',
 525         default=False,
 526         help='use Blender raw format for output')
 527     parser.add_argument('-i', action='store_true', dest='useinputfile', default=False,
 528                         help='use an inputfile (instead of just a synthesized grid)')
 529     parser.add_argument(
 530         '-t',
 531         action='store_true',
 532         dest='timingonly',
 533         default=False,
 534         help='do not write anything to an output file')
 535     parser.add_argument('-infile', type=str, default="-", help='input filename')
 536     parser.add_argument('-outfile', type=str, default="-", help='output filename')
 537     parser.add_argument('-Gn', dest='gridsize', type=int, default=20, help='Gridsize (always square)')
 538     parser.add_argument('-Gp', dest='gridpeak', type=float, default=0, help='Add peak with given height')
 539     parser.add_argument('-Gs', dest='gridshelf', type=float, default=0, help='Add shelve with given height')
 540     parser.add_argument('-Gm', dest='gridmesa', type=float, default=0, help='Add mesa with given height')
 541     parser.add_argument('-Gr', dest='gridrandom', type=float, default=0,
 542                         help='Add random values between 0 and given value')
 543     parser.add_argument('-m', dest='threads', type=int, default=1, help='number of threads to use')
 544     parser.add_argument('-u', action='store_true', dest='unittest', default=False, help='perform unittests')
 545     parser.add_argument('-a', action='store_true', dest='analyze', default=False,
 546                         help='show some statistics of input and output meshes')
 547     parser.add_argument('-d', action='store_true', dest='dump', default=False,
 548                         help='show sediment and water meshes at end of run')
 549     parser.add_argument('-n', action='store_true', dest='usenumexpr', default=False, help='use numexpr optimizations')
 550
 551     args = parser.parse_args()
 552     print("\nInput arguments:")
 553     print("\n".join("%-15s: %s" % t for t in sorted(vars(args).items())), file=sys.stderr)
 554
 555     if args.unittest:
 556         unittest.main(argv=[sys.argv[0]])
 557         sys.exit(0)
 558
 559     if args.useinputfile:
 560         if args.rawin:
 561             grid = Grid.fromRaw(args.infile)
 562         else:
 563             grid = Grid.fromFile(args.infile)
 564     else:
 565         grid = Grid(args.gridsize)
 566
 567     if args.gridpeak > 0:
 568         grid.peak(args.gridpeak)
 569     if args.gridmesa > 0:
 570         grid.mesa(args.gridmesa)
 571     if args.gridshelf > 0:
 572         grid.shelf(args.gridshelf)
 573     if args.gridrandom > 0:
 574         grid.random(args.gridrandom)
 575
 576     if args.analyze:
 577         print('\nstatistics of the input grid:\n\n', grid.analyze(), file=sys.stderr, sep='')
 578     t = getptime()
 579     for g in range(args.iterations):
 580         if args.Kd > 0:
 581             grid.diffuse(args.Kd, args.usenumexpr)
 582         if args.Kh > 0 and args.Kp > rand():
 583             grid.avalanche(args.Kh, args.usenumexpr)
 584         if args.Kr > 0 or args.Kspring > 0:
 585             grid.fluvial_erosion(
 586                 args.Kr,
 587                 args.Kc,
 588                 args.Ks,
 589                 args.Kdep,
 590                 args.Ka,
 591                 args.Kspring,
 592                 args.Kspringx,
 593                 args.Kspringy,
 594                 args.Kspringr,
 595                 args.usenumexpr,
 596             )
 597     t = getptime() - t
 598     print("\nElapsed time: %.1f seconds, max memory %.1f Mb.\n" % (t, grid.maxrss), file=sys.stderr)
 599     if args.analyze:
 600         print('\nstatistics of the output grid:\n\n', grid.analyze(), file=sys.stderr, sep='')
 601
 602     if not args.timingonly:
 603         if args.rawout:
 604             grid.toRaw(args.outfile, vars(args))
 605         else:
 606             grid.toFile(args.outfile)
 607
 608     if args.dump:
 609         print("sediment\n", np.array_str(grid.sediment, precision=3), file=sys.stderr)
 610         print("water\n", np.array_str(grid.water, precision=3), file=sys.stderr)
 611         print("sediment concentration\n", np.array_str(grid.sediment / grid.water, precision=3), file=sys.stderr)