ant_landscape/eroder.py

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