Skip to content

Commit

Permalink
fixed bug in net difference calculation, added copy elevation to curr…
Browse files Browse the repository at this point in the history
…ent mapset to prevent temporal framework error, and refactored code
  • Loading branch information
baharmon committed Dec 30, 2019
1 parent 7172ab5 commit 0825a1b
Showing 1 changed file with 50 additions and 116 deletions.
166 changes: 50 additions & 116 deletions grass7/raster/r.sim.terrain/r.sim.terrain.py
Original file line number Diff line number Diff line change
Expand Up @@ -544,6 +544,21 @@ def main():
expression="k_factor = {k_factor_value}".format(**locals()),
overwrite=True)

# copy the elevation raster if it is not in the current mapset
name = elevation.split("@")[0]
filename = gscript.read_command(
'g.list',
type='raster',
pattern=name,
mapset='.',
flags='m')
if not filename:
gscript.run_command(
'g.copy',
raster='{fullname},{name}'.format(fullname=elevation, name=name),
overwrite=True)
elevation = name

# create dynamic_evolution object
dynamics = DynamicEvolution(elevation=elevation,
mode=mode,
Expand Down Expand Up @@ -997,8 +1012,8 @@ def compute_difference(self, evolved_elevation, difference):
'r.mapcalc',
expression="{difference} = {evolved_elevation}-{elevation}".format(
difference=difference,
elevation=self.elevation,
evolved_elevation=evolved_elevation),
evolved_elevation=evolved_elevation,
elevation=self.elevation),
overwrite=True)
# gscript.write_command(
# 'r.colors',
Expand Down Expand Up @@ -1691,122 +1706,34 @@ def rainfall_event(self):
threads=self.threads,
fill_depressions=self.fill_depressions)

# determine mode and run model
if self.mode == 'simwe_mode':
(evolved_elevation, time, depth, erosion_deposition,
difference) = evol.erosion_deposition()
# remove relative timestamps from r.sim.water and r.sim.sediment
gscript.run_command(
'r.timestamp',
map=depth,
date='none')
gscript.run_command(
'r.timestamp',
map=erosion_deposition,
date='none')

elif self.mode == "usped_mode":
(evolved_elevation, time, depth, erosion_deposition,
difference) = evol.usped()

elif self.mode == "rusle_mode":
(evolved_elevation, time, depth, sediment_flux,
difference) = evol.rusle()

else:
raise RuntimeError(
'{mode} mode does not exist').format(mode=self.mode)

# register the evolved maps
gscript.run_command(
't.register',
type=raster,
input=self.elevation_timeseries,
maps=evolved_elevation,
start=evol.start,
increment=increment,
flags='i',
overwrite=True)
gscript.run_command(
't.register',
type=raster,
input=self.depth_timeseries,
maps=depth,
start=evol.start,
increment=increment,
flags='i',
overwrite=True)
try:
gscript.run_command(
't.register',
type=raster,
input=self.erdep_timeseries,
maps=erosion_deposition,
start=evol.start,
increment=increment,
flags='i',
overwrite=True)
except (NameError, CalledModuleError):
pass
try:
gscript.run_command(
't.register',
type=raster,
input=self.flux_timeseries,
maps=sediment_flux,
start=evol.start,
increment=increment,
flags='i', overwrite=True)
except (NameError, CalledModuleError):
pass
gscript.run_command(
't.register',
type=raster,
input=self.difference_timeseries,
maps=difference,
start=evol.start,
increment=increment,
flags='i',
overwrite=True)

# run the landscape evolution model
# as a series of rainfall intervals in a rainfall event
i = 1
i = 0
while i < iterations:

# update the elevation
evol.elevation = evolved_elevation
print evol.elevation

# update time
evol.start = time
print evol.start

# derive excess water (mm/hr) from rainfall rate (mm/hr)
# plus the depth (m) per rainfall interval (min)
gscript.run_command(
'r.mapcalc',
expression="{rain_excess}"
"={rain_intensity}"
"+{depth}"
"/1000."
"/{rain_interval}"
"*60.".format(
rain_excess=rain_excess,
rain_intensity=self.rain_intensity,
depth=depth,
rain_interval=self.rain_interval),
overwrite=True)

# update excess rainfall
rain_intensity = 'rain_intensity'
gscript.run_command(
'r.mapcalc',
expression="{rain_intensity} = {rain_excess}".format(
rain_intensity='rain_intensity',
rain_excess=rain_excess),
overwrite=True)
evol.rain_intensity = rain_intensity
if i > 0:
# derive excess water (mm/hr) from rainfall rate (mm/hr)
# plus the depth (m) per rainfall interval (min)
gscript.run_command(
'r.mapcalc',
expression="{rain_excess}"
"={rain_intensity}"
"+{depth}"
"/1000."
"/{rain_interval}"
"*60.".format(
rain_excess=rain_excess,
rain_intensity=self.rain_intensity,
depth=depth,
rain_interval=self.rain_interval),
overwrite=True)
# update excess rainfall
rain_intensity = 'rain_intensity'
gscript.run_command(
'r.mapcalc',
expression="{rain_intensity} = {rain_excess}".format(
rain_intensity='rain_intensity',
rain_excess=rain_excess),
overwrite=True)
evol.rain_intensity = rain_intensity

# determine mode and run model
if self.mode == "simwe_mode":
Expand Down Expand Up @@ -1894,6 +1821,13 @@ def rainfall_event(self):
name=['rain_excess'],
flags='f')

# update elevation
evol.elevation = evolved_elevation

# advance time
evol.start = time

# advance iterator
i = i+1

# compute net elevation change
Expand Down

0 comments on commit 0825a1b

Please sign in to comment.