ant_landscape/eroder.py

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