diff options
Diffstat (limited to 'ant_landscape/eroder.py')
| -rw-r--r-- | ant_landscape/eroder.py | 704 |
1 files changed, 324 insertions, 380 deletions
diff --git a/ant_landscape/eroder.py b/ant_landscape/eroder.py index d36aef67..04e04bd7 100644 --- a/ant_landscape/eroder.py +++ b/ant_landscape/eroder.py @@ -2,7 +2,7 @@ # # erode.py -- a script to simulate erosion of height fields # (c) 2014 Michel J. Anders (varkenvarken) -# now with some modifications by Ian Huish (nerk) +# with some modifications by Ian Huish (nerk) # # This program is free software; you can redistribute it and/or # modify it under the terms of the GNU General Public License @@ -20,37 +20,29 @@ # # ##### END GPL LICENSE BLOCK ##### + from time import time import unittest import sys import os -# import resource # so much for platform independence, this only works on unix :-( from random import random as rand, shuffle import numpy as np -#from .perlin import pnoise numexpr_available = False -# Sorry, nerk can't handle numexpr at this time! -#try: -# import numexpr as ne -# numexpr_available = True -#except ImportError: -# pass + def getmemsize(): - return 0.0 - #return resource.getrusage(resource.RUSAGE_SELF).ru_maxrss*resource.getpagesize()/(1024.0*1024.0) - + return 0.0 + + def getptime(): - #r = resource.getrusage(resource.RUSAGE_SELF) - #return r.ru_utime + r.ru_stime - return time() - + return time() + + class Grid: def __init__(self, size=10, dtype=np.single): - self.center = np.zeros([size,size], dtype) - #print("Centre\n", np.array_str(self.center,precision=3), file=sys.stderr) + self.center = np.zeros([size, size], dtype) self.water = None self.sediment = None self.scour = None @@ -59,19 +51,20 @@ class Grid: self.sedimentpct = None self.capacity = None self.avalanced = None - self.minx=None - self.miny=None - self.maxx=None - self.maxy=None - self.zscale=1 - self.maxrss=0.0 - self.sequence=[0,1,2,3] + self.minx = None + self.miny = None + self.maxx = None + self.maxy = None + self.zscale = 1 + self.maxrss = 0.0 + self.sequence = [0, 1, 2, 3] self.watermax = 1.0 self.flowratemax = 1.0 self.scourmax = 1.0 self.sedmax = 1.0 self.scourmin = 1.0 - + + def init_water_and_sediment(self): if self.water is None: self.water = np.zeros(self.center.shape, dtype=np.single) @@ -88,9 +81,11 @@ class Grid: if self.avalanced is None: self.avalanced = np.zeros(self.center.shape, dtype=np.single) + def __str__(self): return ''.join(self.__str_iter__(fmt="%.3f")) - + + def __str_iter__(self, fmt): for row in self.center[::]: values=[] @@ -98,20 +93,24 @@ class Grid: values.append(fmt%v) yield ' '.join(values) + '\n' + @staticmethod def fromFile(filename): - if filename == '-' : filename = sys.stdin + if filename == '-': + filename = sys.stdin g=Grid() g.center=np.loadtxt(filename,np.single) return g + def toFile(self, filename, fmt="%.3f"): if filename == '-' : filename = sys.stdout.fileno() with open(filename,"w") as f: for line in self.__str_iter__(fmt): f.write(line) - + + def raw(self,format="%.3f"): fstr=format+" "+ format+" "+ format+" " a=self.center / self.zscale @@ -134,6 +133,7 @@ class Grid: fstr%(row1 ,col0 ,a[row+1][col ])+ fstr%(row1 ,col1 ,a[row+1][col+1])+"\n") + def toRaw(self, filename, infomap=None): with open(filename if type(filename) == str else sys.stdout.fileno() , "w") as f: f.writelines(self.raw()) @@ -141,85 +141,81 @@ class Grid: with open(os.path.splitext(filename)[0]+".inf" if type(filename) == str else sys.stdout.fileno() , "w") as f: f.writelines("\n".join("%-15s: %s"%t for t in sorted(infomap.items()))) + @staticmethod def fromRaw(filename): """initialize a grid from a Blender .raw file. currenly suports just rectangular grids of all triangles """ - g=Grid.fromFile(filename) - # we assume tris and an axis aligned grid - g.center=np.reshape(g.center,(-1,3)) + g = Grid.fromFile(filename) + # we assume tris and an axis aligned grid + g.center = np.reshape(g.center,(-1,3)) g._sort() return g + def _sort(self, expfact): - # keep unique vertices only by creating a set and sort first on x then on y coordinate - # using rather slow python sort but couldn;t wrap my head around np.lexsort + # keep unique vertices only by creating a set and sort first on x then on y coordinate + # using rather slow python sort but couldn;t wrap my head around np.lexsort verts = sorted(list({ tuple(t) for t in self.center[::] })) - x=set(c[0] for c in verts) - y=set(c[1] for c in verts) - nx=len(x) - ny=len(y) - self.minx=min(x) - self.maxx=max(x) - self.miny=min(y) - self.maxy=max(y) - xscale=(self.maxx-self.minx)/(nx-1) - yscale=(self.maxy-self.miny)/(ny-1) + x = set(c[0] for c in verts) + y = set(c[1] for c in verts) + nx = len(x) + ny = len(y) + self.minx = min(x) + self.maxx = max(x) + self.miny = min(y) + self.maxy = max(y) + xscale = (self.maxx-self.minx)/(nx-1) + yscale = (self.maxy-self.miny)/(ny-1) # note: a purely flat plane cannot be scaled - if (yscale != 0.0) and (abs(xscale/yscale) - 1.0 > 1e-3) : raise ValueError("Mesh spacing not square %d x %d %.4f x %4.f"%(nx,ny,xscale,yscale)) - self.zscale=1.0 + if (yscale != 0.0) and (abs(xscale/yscale) - 1.0 > 1e-3): + raise ValueError("Mesh spacing not square %d x %d %.4f x %4.f"%(nx,ny,xscale,yscale)) + self.zscale = 1.0 if abs(yscale) > 1e-6 : - self.zscale=1.0/yscale + self.zscale = 1.0/yscale - # keep just the z-values and null any ofsset - # 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 - self.center=np.array([c[2] for c in verts],dtype=np.single).reshape(nx,ny) - self.center=(self.center-np.amin(self.center))*self.zscale + # keep just the z-values and null any ofsset + # 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 + self.center = np.array([c[2] for c in verts],dtype=np.single).reshape(nx,ny) + self.center = (self.center-np.amin(self.center))*self.zscale if self.rainmap is not None: - #rainmap = sorted(list({ tuple(t) for t in self.rainmap[::] })) - #self.rainmap=np.array([c[2] for c in rainmap],dtype=np.single).reshape(nx,ny) rmscale = np.max(self.center) - #self.rainmap = (self.center/rmscale) * np.exp(expfact*((self.center/rmscale)-1)) self.rainmap = expfact + (1-expfact)*(self.center/rmscale) + @staticmethod def fromBlenderMesh(me, vg, expfact): - g=Grid() - g.center=np.asarray(list(tuple(v.co) for v in me.vertices), dtype=np.single ) - g.rainmap=None - print("VertexGroup\n",vg, file=sys.stderr) + g = Grid() + g.center = np.asarray(list(tuple(v.co) for v in me.vertices), dtype=np.single ) + g.rainmap = None if vg is not None: for v in me.vertices: vg.add([v.index],0.0,'ADD') g.rainmap=np.asarray(list( (v.co[0], v.co[1], vg.weight(v.index)) for v in me.vertices), dtype=np.single ) g._sort(expfact) - #print("CentreMesh\n", np.array_str(g.center,precision=3), file=sys.stderr) - #print('rainmap',np.max(g.rainmap),np.min(g.rainmap)) return g -# def rainmapcolor(me, vg): -# if vg is not None: -# for v in me.vertices: - def setrainmap(self, rainmap): self.rainmap = rainmap + def _verts(self, surface): - a=surface / self.zscale - minx=0.0 if self.minx is None else self.minx - miny=0.0 if self.miny is None else self.miny - maxx=1.0 if self.maxx is None else self.maxx - maxy=1.0 if self.maxy is None else self.maxy - dx=(maxx-minx)/(a.shape[0]-1) - dy=(maxy-miny)/(a.shape[1]-1) + a = surface / self.zscale + minx = 0.0 if self.minx is None else self.minx + miny = 0.0 if self.miny is None else self.miny + maxx = 1.0 if self.maxx is None else self.maxx + maxy = 1.0 if self.maxy is None else self.maxy + dx = (maxx - minx) / (a.shape[0] - 1) + dy = (maxy - miny) / (a.shape[1] - 1) for row in range(a.shape[0]): - row0=miny+row*dy + row0 = miny + row * dy for col in range(a.shape[1]): - col0=minx+col*dx + col0 = minx + col * dx yield (row0 ,col0 ,a[row ][col ]) + def _faces(self): nrow, ncol = self.center.shape for row in range(nrow-1): @@ -228,61 +224,70 @@ class Grid: yield (vi, vi+ncol, vi+1) yield (vi+1, vi+ncol, vi+ncol+1) - def toBlenderMesh(self, me): # pass me as argument so that we don't need to import bpy and create a dependency + + def toBlenderMesh(self, me): + # pass me as argument so that we don't need to import bpy and create a dependency # the docs state that from_pydata takes iterators as arguments but it will fail with generators because it does len(arg) me.from_pydata(list(self._verts(self.center)),[],list(self._faces())) - def toWaterMesh(self, me): # pass me as argument so that we don't need to import bpy and create a dependency + + def toWaterMesh(self, me): + # pass me as argument so that we don't need to import bpy and create a dependency # the docs state that from_pydata takes iterators as arguments but it will fail with generators because it does len(arg) me.from_pydata(list(self._verts(self.water)),[],list(self._faces())) + def peak(self, value=1): nx,ny = self.center.shape self.center[int(nx/2),int(ny/2)] += value + def shelf(self, value=1): nx,ny = self.center.shape self.center[:nx/2] += value + def mesa(self, value=1): nx,ny = self.center.shape self.center[nx/4:3*nx/4,ny/4:3*ny/4] += value + def random(self, value=1): self.center += np.random.random_sample(self.center.shape)*value + def neighborgrid(self): - self.up=np.roll(self.center,-1,0) - self.down=np.roll(self.center,1,0) - self.left=np.roll(self.center,-1,1) - self.right=np.roll(self.center,1,1) + self.up = np.roll(self.center,-1,0) + self.down = np.roll(self.center,1,0) + self.left = np.roll(self.center,-1,1) + self.right = np.roll(self.center,1,1) + def zeroedge(self, quantity=None): c = self.center if quantity is None else quantity - c[0,:]=0 - c[-1,:]=0 - c[:,0]=0 - c[:,-1]=0 + c[0,:] = 0 + c[-1,:] = 0 + c[:,0] = 0 + c[:,-1] = 0 + def diffuse(self, Kd, IterDiffuse, numexpr): self.zeroedge() - c = self.center[1:-1,1:-1] - up = self.center[ :-2,1:-1] - down = self.center[2: ,1:-1] - left = self.center[1:-1, :-2] + c = self.center[1:-1,1:-1] + up = self.center[ :-2,1:-1] + down = self.center[2: ,1:-1] + left = self.center[1:-1, :-2] right = self.center[1:-1,2: ] if(numexpr and numexpr_available): self.center[1:-1,1:-1] = ne.evaluate('c + Kd * (up + down + left + right - 4.0 * c)') else: self.center[1:-1,1:-1] = c + (Kd/IterDiffuse) * (up + down + left + right - 4.0 * c) - print("diffuse: ", Kd) self.maxrss = max(getmemsize(), self.maxrss) return self.center + def avalanche(self, delta, iterava, prob, numexpr): self.zeroedge() - #print(self.center) - c = self.center[1:-1,1:-1] up = self.center[ :-2,1:-1] down = self.center[2: ,1:-1] @@ -301,12 +306,12 @@ class Grid: + where((right-c) < -delta,(right-c +delta)/2, 0)') else: sa = ( - # incoming + # incoming where((up -c) > delta ,(up -c -delta)/2, 0) + where((down -c) > delta ,(down -c -delta)/2, 0) + where((left -c) > delta ,(left -c -delta)/2, 0) + where((right-c) > delta ,(right-c -delta)/2, 0) - # outgoing + # outgoing + where((up -c) < -delta,(up -c +delta)/2, 0) + where((down -c) < -delta,(down -c +delta)/2, 0) + where((left -c) < -delta,(left -c +delta)/2, 0) @@ -317,14 +322,16 @@ class Grid: self.avalanced[1:-1,1:-1] = self.avalanced[1:-1,1:-1] + sa/iterava self.center[1:-1,1:-1] = c + sa/iterava - #print(self.center) self.maxrss = max(getmemsize(), self.maxrss) return self.center + def rain(self, amount=1, variance=0, userainmap=False): self.water += (1.0 - np.random.random(self.water.shape) * variance) * (amount if ((self.rainmap is None) or (not userainmap)) else self.rainmap * amount) - def spring(self, amount, px, py, radius): # px, py and radius are all fractions + + def spring(self, amount, px, py, radius): + # px, py and radius are all fractions nx, ny = self.center.shape rx = max(int(nx*radius),1) ry = max(int(ny*radius),1) @@ -332,320 +339,257 @@ class Grid: py = int(ny*py) self.water[px-rx:px+rx+1,py-ry:py+ry+1] += amount + def river(self, Kc, Ks, Kdep, Ka, Kev, numexpr): + zeros = np.zeros + where = np.where + min = np.minimum + max = np.maximum + abs = np.absolute + arctan = np.arctan + sin = np.sin + + center = (slice( 1, -1,None),slice( 1, -1,None)) + up = (slice(None, -2,None),slice( 1, -1,None)) + down = (slice( 2, None,None),slice( 1, -1,None)) + left = (slice( 1, -1,None),slice(None, -2,None)) + right = (slice( 1, -1,None),slice( 2,None,None)) + + water = self.water + rock = self.center + sediment = self.sediment + height = rock + water + + # !! this gives a runtime warning for division by zero + verysmallnumber = 0.0000000001 + water += verysmallnumber + sc = where(water > verysmallnumber, sediment / water, 0) + + sdw = zeros(water[center].shape) + svdw = zeros(water[center].shape) + sds = zeros(water[center].shape) + angle = zeros(water[center].shape) + for d in (up,down,left,right): + if(numexpr and numexpr_available): + hdd = height[d] + hcc = height[center] + dw = ne.evaluate('hdd-hcc') + inflow = ne.evaluate('dw > 0') + wdd = water[d] + wcc = water[center] + dw = ne.evaluate('where(inflow, where(wdd<dw, wdd, dw), where(-wcc>dw, -wcc, dw))/4.0') # nested where() represent min() and max() + sdw = ne.evaluate('sdw + dw') + scd = sc[d] + scc = sc[center] + rockd= rock[d] + rockc= rock[center] + sds = ne.evaluate('sds + dw * where(inflow, scd, scc)') + svdw = ne.evaluate('svdw + abs(dw)') + angle= ne.evaluate('angle + arctan(abs(rockd-rockc))') + else: + dw = (height[d]-height[center]) + inflow = dw > 0 + dw = where(inflow, min(water[d], dw), max(-water[center], dw))/4.0 + sdw = sdw + dw + sds = sds + dw * where(inflow, sc[d], sc[center]) + svdw = svdw + abs(dw) + angle= angle + np.arctan(abs(rock[d]-rock[center])) - zeros = np.zeros - where = np.where - min = np.minimum - max = np.maximum - abs = np.absolute - arctan = np.arctan - sin = np.sin - - center = (slice( 1, -1,None),slice( 1, -1,None)) - #print("CentreSlice\n", np.array_str(center,precision=3), file=sys.stderr) - up = (slice(None, -2,None),slice( 1, -1,None)) - down = (slice( 2, None,None),slice( 1, -1,None)) - left = (slice( 1, -1,None),slice(None, -2,None)) - right = (slice( 1, -1,None),slice( 2,None,None)) - - water = self.water - rock = self.center - sediment = self.sediment - height = rock + water - sc = where(water>0, sediment/water, 0) ##!! this gives a runtime warning for division by zero - sdw = zeros(water[center].shape) - svdw = zeros(water[center].shape) - sds = zeros(water[center].shape) - angle = zeros(water[center].shape) - #print(water[center]) - for d in (up,down,left,right): - if(numexpr and numexpr_available): - hdd = height[d] - hcc = height[center] - dw = ne.evaluate('hdd-hcc') - inflow = ne.evaluate('dw > 0') - wdd = water[d] + if(numexpr and numexpr_available): + wcc = water[center] + scc = sediment[center] + rcc = rock[center] + water[center] = ne.evaluate('wcc + sdw') + sediment[center] = ne.evaluate('scc + sds') + sc = ne.evaluate('where(wcc>0, scc/wcc, 2000*Kc)') + fKc = ne.evaluate('Kc*sin(Ka*angle)*svdw') + ds = ne.evaluate('where(sc > fKc, -Kd * scc, Ks * svdw)') + rock[center] = ne.evaluate('rcc - ds') + # there isn't really a bottom to the rock but negative values look ugly + rock[center] = ne.evaluate('where(rcc<0,0,rcc)') + sediment[center] = ne.evaluate('scc + ds') + else: wcc = water[center] - dw = ne.evaluate('where(inflow, where(wdd<dw, wdd, dw), where(-wcc>dw, -wcc, dw))/4.0') # nested where() represent min() and max() - sdw = ne.evaluate('sdw + dw') - scd = sc[d] - scc = sc[center] - rockd= rock[d] - rockc= rock[center] - sds = ne.evaluate('sds + dw * where(inflow, scd, scc)') - svdw = ne.evaluate('svdw + abs(dw)') - angle= ne.evaluate('angle + arctan(abs(rockd-rockc))') - else: - dw = (height[d]-height[center]) - inflow = dw > 0 - dw = where(inflow, min(water[d], dw), max(-water[center], dw))/4.0 - sdw = sdw + dw - sds = sds + dw * where(inflow, sc[d], sc[center]) - svdw = svdw + abs(dw) - angle= angle + np.arctan(abs(rock[d]-rock[center])) - - if(numexpr and numexpr_available): - wcc = water[center] - scc = sediment[center] - rcc = rock[center] - water[center] = ne.evaluate('wcc + sdw') - sediment[center] = ne.evaluate('scc + sds') - sc = ne.evaluate('where(wcc>0, scc/wcc, 2000*Kc)') - fKc = ne.evaluate('Kc*sin(Ka*angle)*svdw') - ds = ne.evaluate('where(sc > fKc, -Kd * scc, Ks * svdw)') - rock[center] = ne.evaluate('rcc - ds') - rock[center] = ne.evaluate('where(rcc<0,0,rcc)') # there isn't really a bottom to the rock but negative values look ugly - sediment[center] = ne.evaluate('scc + ds') - else: - wcc = water[center] - scc = sediment[center] - rcc = rock[center] - water[center] = wcc * (1-Kev) + sdw - sediment[center] = scc + sds - sc = where(wcc>0, scc/wcc, 2*Kc) - fKc = Kc*svdw - #fKc = Kc*np.sin(Ka*angle)*svdw*wcc - #ds = where(sc > fKc, -Kd * scc, Ks * svdw) - ds = where(fKc>sc,(fKc-sc)*Ks,(fKc-sc)*Kdep)*wcc - self.flowrate[center] = svdw - self.scour[center] = ds - self.sedimentpct[center] = sc - self.capacity[center] = fKc - #rock[center] = rcc - ds - #rock[center] = where(rcc<0,0,rcc) # there isn't really a bottom to the rock but negative values look ugly - sediment[center] = scc + ds + sds - #print("sdw", sdw[10,15]) + scc = sediment[center] + rcc = rock[center] + water[center] = wcc * (1-Kev) + sdw + sediment[center] = scc + sds + sc = where(wcc > 0, scc / wcc, 2 * Kc) + fKc = Kc*svdw + ds = where(fKc > sc, (fKc - sc) * Ks, (fKc - sc) * Kdep) * wcc + self.flowrate[center] = svdw + self.scour[center] = ds + self.sedimentpct[center] = sc + self.capacity[center] = fKc + sediment[center] = scc + ds + sds + def flow(self, Kc, Ks, Kz, Ka, numexpr): + zeros = np.zeros + where = np.where + min = np.minimum + max = np.maximum + abs = np.absolute + arctan = np.arctan + sin = np.sin + + center = (slice( 1, -1,None),slice( 1, -1,None)) + rock = self.center + ds = self.scour[center] + rcc = rock[center] + rock[center] = rcc - ds * Kz + # there isn't really a bottom to the rock but negative values look ugly + rock[center] = where(rcc<0,0,rcc) - zeros = np.zeros - where = np.where - min = np.minimum - max = np.maximum - abs = np.absolute - arctan = np.arctan - sin = np.sin - - center = (slice( 1, -1,None),slice( 1, -1,None)) - #print("CentreSlice\n", np.array_str(center,precision=3), file=sys.stderr) - #up = (slice(None, -2,None),slice( 1, -1,None)) - #down = (slice( 2, None,None),slice( 1, -1,None)) - #left = (slice( 1, -1,None),slice(None, -2,None)) - #right = (slice( 1, -1,None),slice( 2,None,None)) - - #water = self.water - rock = self.center - #sediment = self.sediment - #height = rock + water - #sc = where(water>0, sediment/water, 0) ##!! this gives a runtime warning for division by zero - #sdw = zeros(water[center].shape) - #svdw = zeros(water[center].shape) - #sds = zeros(water[center].shape) - #angle = zeros(water[center].shape) - #print(height[center]) - #print(water[center]) - #for d in (up,down,left,right): - #if(numexpr and numexpr_available): - #hdd = height[d] - #hcc = height[center] - #dw = ne.evaluate('hdd-hcc') - #inflow = ne.evaluate('dw > 0') - #wdd = water[d] - #wcc = water[center] - #dw = ne.evaluate('where(inflow, where(wdd<dw, wdd, dw), where(-wcc>dw, -wcc, dw))/4.0') # nested where() represent min() and max() - #sdw = ne.evaluate('sdw + dw') - #scd = sc[d] - #scc = sc[center] - #rockd= rock[d] - #rockc= rock[center] - #sds = ne.evaluate('sds + dw * where(inflow, scd, scc)') - #svdw = ne.evaluate('svdw + abs(dw)') - #angle= ne.evaluate('angle + arctan(abs(rockd-rockc))') - #else: - #dw = (height[d]-height[center]) - #inflow = dw > 0 - #dw = where(inflow, min(water[d], dw), max(-water[center], dw))/4.0 - #sdw = sdw + dw - #sds = sds + dw * where(inflow, sc[d], sc[center]) - #svdw = svdw + abs(dw) - #angle= angle + np.arctan(abs(rock[d]-rock[center])) - - #if(numexpr and numexpr_available): - #wcc = water[center] - #scc = sediment[center] - #rcc = rock[center] - #water[center] = ne.evaluate('wcc + sdw') - #sediment[center] = ne.evaluate('scc + sds') - #sc = ne.evaluate('where(wcc>0, scc/wcc, 2000*Kc)') - #fKc = ne.evaluate('Kc*sin(Ka*angle)*svdw') - #ds = ne.evaluate('where(sc > fKc, -Kd * scc, Ks * svdw)') - #rock[center] = ne.evaluate('rcc - ds') - #rock[center] = ne.evaluate('where(rcc<0,0,rcc)') # there isn't really a bottom to the rock but negative values look ugly - #sediment[center] = ne.evaluate('scc + ds') - #else: - #wcc = water[center] - #scc = sediment[center] - ds = self.scour[center] - rcc = rock[center] - #water[center] = wcc + sdw - #sediment[center] = scc + sds - #sc = where(wcc>0, scc/wcc, 2*Kc) - #fKc = Kc*np.sin(Ka*angle)*svdw - #ds = where(sc > fKc, -Kd * scc, Ks * svdw) - rock[center] = rcc - ds * Kz - rock[center] = where(rcc<0,0,rcc) # there isn't really a bottom to the rock but negative values look ugly - #sediment[center] = scc + ds def rivergeneration(self, rainamount, rainvariance, userainmap, Kc, Ks, Kdep, Ka, Kev, Kspring, Kspringx, Kspringy, Kspringr, numexpr): self.init_water_and_sediment() self.rain(rainamount, rainvariance, userainmap) self.zeroedge(self.water) self.zeroedge(self.sediment) - #self.spring(Kspring, Kspringx, Kspringy, Kspringr) self.river(Kc, Ks, Kdep, Ka, Kev, numexpr) self.watermax = np.max(self.water) + def fluvial_erosion(self, rainamount, rainvariance, userainmap, Kc, Ks, Kdep, Ka, Kspring, Kspringx, Kspringy, Kspringr, numexpr): - #self.init_water_and_sediment() - #self.rain(rainamount, rainvariance, userainmap) - #self.zeroedge(self.water) - #self.zeroedge(self.sediment) - #self.spring(Kspring, Kspringx, Kspringy, Kspringr) self.flow(Kc, Ks, Kdep, Ka, numexpr) self.flowratemax = np.max(self.flowrate) self.scourmax = np.max(self.scour) self.scourmin = np.min(self.scour) self.sedmax = np.max(self.sediment) - print("DSMinMax", np.min(self.scour), np.max(self.scour)) + def analyze(self): - self.neighborgrid() - # just looking at up and left to avoid needless doubel calculations - slopes=np.concatenate((np.abs(self.left - self.center),np.abs(self.up - self.center))) - return '\n'.join(["%-15s: %.3f"%t for t in [ - ('height average', np.average(self.center)), - ('height median', np.median(self.center)), - ('height max', np.max(self.center)), - ('height min', np.min(self.center)), - ('height std', np.std(self.center)), - ('slope average', np.average(slopes)), - ('slope median', np.median(slopes)), - ('slope max', np.max(slopes)), - ('slope min', np.min(slopes)), - ('slope std', np.std(slopes)) - ]] - ) + self.neighborgrid() + # just looking at up and left to avoid needless doubel calculations + slopes=np.concatenate((np.abs(self.left - self.center),np.abs(self.up - self.center))) + return '\n'.join(["%-15s: %.3f"%t for t in [ + ('height average', np.average(self.center)), + ('height median', np.median(self.center)), + ('height max', np.max(self.center)), + ('height min', np.min(self.center)), + ('height std', np.std(self.center)), + ('slope average', np.average(slopes)), + ('slope median', np.median(slopes)), + ('slope max', np.max(slopes)), + ('slope min', np.min(slopes)), + ('slope std', np.std(slopes)) + ]] + ) + class TestGrid(unittest.TestCase): - def test_diffuse(self): - g=Grid(5) - g.peak(1) - self.assertEqual(g.center[2,2],1.0) - g.diffuse(0.1, numexpr=False) - for n in [(2,1),(2,3),(1,2),(3,2)]: - self.assertAlmostEqual(g.center[n],0.1) - self.assertAlmostEqual(g.center[2,2],0.6) - - def test_diffuse_numexpr(self): - g=Grid(5) - g.peak(1) - g.diffuse(0.1, numexpr=False) - h=Grid(5) - h.peak(1) - h.diffuse(0.1, numexpr=True) - self.assertEqual(list(g.center.flat),list(h.center.flat)) - - def test_avalanche_numexpr(self): - g=Grid(5) - g.peak(1) - g.avalanche(0.1, numexpr=False) - h=Grid(5) - h.peak(1) - h.avalanche(0.1, numexpr=True) - print(g) - print(h) - np.testing.assert_almost_equal(g.center,h.center) + def test_diffuse(self): + g = Grid(5) + g.peak(1) + self.assertEqual(g.center[2,2],1.0) + g.diffuse(0.1, numexpr=False) + for n in [(2,1),(2,3),(1,2),(3,2)]: + self.assertAlmostEqual(g.center[n],0.1) + self.assertAlmostEqual(g.center[2,2],0.6) + + + def test_diffuse_numexpr(self): + g = Grid(5) + g.peak(1) + g.diffuse(0.1, numexpr=False) + h = Grid(5) + h.peak(1) + h.diffuse(0.1, numexpr=True) + self.assertEqual(list(g.center.flat),list(h.center.flat)) + + + def test_avalanche_numexpr(self): + g = Grid(5) + g.peak(1) + g.avalanche(0.1, numexpr=False) + h = Grid(5) + h.peak(1) + h.avalanche(0.1, numexpr=True) + print(g) + print(h) + np.testing.assert_almost_equal(g.center,h.center) + if __name__ == "__main__": - import argparse - - parser = argparse.ArgumentParser(description='Erode a terrain while assuming zero boundary conditions.') - parser.add_argument('-I', dest='iterations', type=int, default=1, help='the number of iterations') - parser.add_argument('-Kd', dest='Kd', type=float, default=0.01, help='Diffusion constant') - parser.add_argument('-Kh', dest='Kh', type=float, default=6, help='Maximum stable cliff height') - parser.add_argument('-Kp', dest='Kp', type=float, default=0.1, help='Avalanche probability for unstable cliffs') - parser.add_argument('-Kr', dest='Kr', type=float, default=0.1, help='Average amount of rain per iteration') - parser.add_argument('-Kspring', dest='Kspring', type=float, default=0.0, help='Average amount of wellwater per iteration') - parser.add_argument('-Kspringx', dest='Kspringx', type=float, default=0.5, help='relative x position of spring') - parser.add_argument('-Kspringy', dest='Kspringy', type=float, default=0.5, help='relative y position of spring') - parser.add_argument('-Kspringr', dest='Kspringr', type=float, default=0.02, help='radius of spring') - parser.add_argument('-Kdep', dest='Kdep', type=float, default=0.1, help='Sediment deposition constant') - parser.add_argument('-Ks', dest='Ks', type=float, default=0.1, help='Soil softness constant') - parser.add_argument('-Kc', dest='Kc', type=float, default=1.0, help='Sediment capacity') - parser.add_argument('-Ka', dest='Ka', type=float, default=2.0, help='Slope dependency of erosion') - parser.add_argument('-ri', action='store_true', dest='rawin', default=False, help='use Blender raw format for input') - parser.add_argument('-ro', action='store_true', dest='rawout', default=False, help='use Blender raw format for output') - parser.add_argument('-i', action='store_true', dest='useinputfile', default=False, help='use an inputfile (instead of just a synthesized grid)') - parser.add_argument('-t', action='store_true', dest='timingonly', default=False, help='do not write anything to an output file') - parser.add_argument('-infile', type=str, default="-", help='input filename') - parser.add_argument('-outfile', type=str, default="-", help='output filename') - parser.add_argument('-Gn', dest='gridsize', type=int, default=20, help='Gridsize (always square)') - parser.add_argument('-Gp', dest='gridpeak', type=float, default=0, help='Add peak with given height') - parser.add_argument('-Gs', dest='gridshelf', type=float, default=0, help='Add shelve with given height') - parser.add_argument('-Gm', dest='gridmesa', type=float, default=0, help='Add mesa with given height') - parser.add_argument('-Gr', dest='gridrandom', type=float, default=0, help='Add random values between 0 and given value') - parser.add_argument('-m', dest='threads', type=int, default=1, help='number of threads to use') - parser.add_argument('-u', action='store_true', dest='unittest', default=False, help='perfom unittests') - parser.add_argument('-a', action='store_true', dest='analyze', default=False, help='show some statistics of input and output meshes') - parser.add_argument('-d', action='store_true', dest='dump', default=False, help='show sediment and water meshes at end of run') - parser.add_argument('-n', action='store_true', dest='usenumexpr', default=False, help='use numexpr optimizations') - - args = parser.parse_args() - print("\nInput arguments:") - print("\n".join("%-15s: %s"%t for t in sorted(vars(args).items())), file=sys.stderr) - - if args.unittest: - unittest.main(argv=[sys.argv[0]]) - sys.exit(0) - - if args.useinputfile: - if args.rawin: - grid = Grid.fromRaw(args.infile) - else: - grid = Grid.fromFile(args.infile) - else: - grid = Grid(args.gridsize) - - if args.gridpeak > 0 : grid.peak(args.gridpeak) - if args.gridmesa > 0 : grid.mesa(args.gridmesa) - if args.gridshelf > 0 : grid.shelf(args.gridshelf) - if args.gridrandom > 0 : grid.random(args.gridrandom) - - if args.analyze: - print('\nstatistics of the input grid:\n\n', grid.analyze(), file=sys.stderr, sep='' ) - t = getptime() - for g in range(args.iterations): - if args.Kd > 0: - grid.diffuse(args.Kd, args.usenumexpr) - if args.Kh > 0 and args.Kp > rand(): - grid.avalanche(args.Kh, args.usenumexpr) - if args.Kr > 0 or args.Kspring > 0: - grid.fluvial_erosion(args.Kr, args.Kc, args.Ks, args.Kdep, args.Ka, args.Kspring, args.Kspringx, args.Kspringy, args.Kspringr, args.usenumexpr) - t = getptime() - t - print("\nElapsed time: %.1f seconds, max memory %.1f Mb.\n"%(t,grid.maxrss), file=sys.stderr) - if args.analyze: - print('\nstatistics of the output grid:\n\n', grid.analyze(), file=sys.stderr, sep='') - - if not args.timingonly: - if args.rawout: - grid.toRaw(args.outfile, vars(args)) + import argparse + + parser = argparse.ArgumentParser(description='Erode a terrain while assuming zero boundary conditions.') + parser.add_argument('-I', dest='iterations', type=int, default=1, help='the number of iterations') + parser.add_argument('-Kd', dest='Kd', type=float, default=0.01, help='Diffusion constant') + parser.add_argument('-Kh', dest='Kh', type=float, default=6, help='Maximum stable cliff height') + parser.add_argument('-Kp', dest='Kp', type=float, default=0.1, help='Avalanche probability for unstable cliffs') + parser.add_argument('-Kr', dest='Kr', type=float, default=0.1, help='Average amount of rain per iteration') + parser.add_argument('-Kspring', dest='Kspring', type=float, default=0.0, help='Average amount of wellwater per iteration') + parser.add_argument('-Kspringx', dest='Kspringx', type=float, default=0.5, help='relative x position of spring') + parser.add_argument('-Kspringy', dest='Kspringy', type=float, default=0.5, help='relative y position of spring') + parser.add_argument('-Kspringr', dest='Kspringr', type=float, default=0.02, help='radius of spring') + parser.add_argument('-Kdep', dest='Kdep', type=float, default=0.1, help='Sediment deposition constant') + parser.add_argument('-Ks', dest='Ks', type=float, default=0.1, help='Soil softness constant') + parser.add_argument('-Kc', dest='Kc', type=float, default=1.0, help='Sediment capacity') + parser.add_argument('-Ka', dest='Ka', type=float, default=2.0, help='Slope dependency of erosion') + parser.add_argument('-ri', action='store_true', dest='rawin', default=False, help='use Blender raw format for input') + parser.add_argument('-ro', action='store_true', dest='rawout', default=False, help='use Blender raw format for output') + parser.add_argument('-i', action='store_true', dest='useinputfile', default=False, help='use an inputfile (instead of just a synthesized grid)') + parser.add_argument('-t', action='store_true', dest='timingonly', default=False, help='do not write anything to an output file') + parser.add_argument('-infile', type=str, default="-", help='input filename') + parser.add_argument('-outfile', type=str, default="-", help='output filename') + parser.add_argument('-Gn', dest='gridsize', type=int, default=20, help='Gridsize (always square)') + parser.add_argument('-Gp', dest='gridpeak', type=float, default=0, help='Add peak with given height') + parser.add_argument('-Gs', dest='gridshelf', type=float, default=0, help='Add shelve with given height') + parser.add_argument('-Gm', dest='gridmesa', type=float, default=0, help='Add mesa with given height') + parser.add_argument('-Gr', dest='gridrandom', type=float, default=0, help='Add random values between 0 and given value') + parser.add_argument('-m', dest='threads', type=int, default=1, help='number of threads to use') + parser.add_argument('-u', action='store_true', dest='unittest', default=False, help='perfom unittests') + parser.add_argument('-a', action='store_true', dest='analyze', default=False, help='show some statistics of input and output meshes') + parser.add_argument('-d', action='store_true', dest='dump', default=False, help='show sediment and water meshes at end of run') + parser.add_argument('-n', action='store_true', dest='usenumexpr', default=False, help='use numexpr optimizations') + + args = parser.parse_args() + print("\nInput arguments:") + print("\n".join("%-15s: %s"%t for t in sorted(vars(args).items())), file=sys.stderr) + + if args.unittest: + unittest.main(argv=[sys.argv[0]]) + sys.exit(0) + + if args.useinputfile: + if args.rawin: + grid = Grid.fromRaw(args.infile) + else: + grid = Grid.fromFile(args.infile) else: - grid.toFile(args.outfile) + grid = Grid(args.gridsize) + + if args.gridpeak > 0 : grid.peak(args.gridpeak) + if args.gridmesa > 0 : grid.mesa(args.gridmesa) + if args.gridshelf > 0 : grid.shelf(args.gridshelf) + if args.gridrandom > 0 : grid.random(args.gridrandom) + + if args.analyze: + print('\nstatistics of the input grid:\n\n', grid.analyze(), file=sys.stderr, sep='' ) + t = getptime() + for g in range(args.iterations): + if args.Kd > 0: + grid.diffuse(args.Kd, args.usenumexpr) + if args.Kh > 0 and args.Kp > rand(): + grid.avalanche(args.Kh, args.usenumexpr) + if args.Kr > 0 or args.Kspring > 0: + grid.fluvial_erosion(args.Kr, args.Kc, args.Ks, args.Kdep, args.Ka, args.Kspring, args.Kspringx, args.Kspringy, args.Kspringr, args.usenumexpr) + t = getptime() - t + print("\nElapsed time: %.1f seconds, max memory %.1f Mb.\n"%(t,grid.maxrss), file=sys.stderr) + if args.analyze: + print('\nstatistics of the output grid:\n\n', grid.analyze(), file=sys.stderr, sep='') + + if not args.timingonly: + if args.rawout: + grid.toRaw(args.outfile, vars(args)) + else: + grid.toFile(args.outfile) - if args.dump: - print("sediment\n", np.array_str(grid.sediment,precision=3), file=sys.stderr) - print("water\n", np.array_str(grid.water,precision=3), file=sys.stderr) - print("sediment concentration\n", np.array_str(grid.sediment/grid.water,precision=3), file=sys.stderr) + if args.dump: + print("sediment\n", np.array_str(grid.sediment,precision=3), file=sys.stderr) + print("water\n", np.array_str(grid.water,precision=3), file=sys.stderr) + print("sediment concentration\n", np.array_str(grid.sediment/grid.water,precision=3), file=sys.stderr) |