# Solutions

1. Compute topographic index using r.topidx.

```def run_slope(scanned_elev, env, **kwargs):
gscript.run_command('r.topidx', input=scanned_elev, output='topidx', env=env)
```

2. Compute topographic aspect (slope orientation) using r.slope.aspect and reclassify it into 8 main directions.

```def run_aspect(scanned_elev, env, **kwargs):
gscript.run_command('r.slope.aspect', elevation=scanned_elev, aspect='aspect', env=env)
rules = ['45:135:1', '135:225:2', '225:315:3', '315:45:4']
gscript.write_command('r.recode', input='aspect', output='aspect_class', rules='-', stdin='\n'.join(rules), env=env)
# set new color table: green - yellow - red
gscript.run_command('r.colors', map='aspect_class', color='random', env=env)
```

3. Show areas with concave profile and tangential curvature (concave forms have negative curvature).

```def run_curvatures(scanned_elev, env, **kwargs):
gscript.run_command('r.param.scale', input=scanned_elev, output='profile_curv', method='profc', size=11, env=env)
gscript.run_command('r.param.scale', input=scanned_elev, output='tangential_curv', method='crosc', size=11, env=env)
gscript.mapcalc("concave = if (profile_curv < 0 && tangential_curv < 0, 1, null())", env=env)
```

4. Derive peaks using either r.geomorphon or r.param.scale and convert them to points (using r.to.vect and v.to.points). From each of those points compute visibility with observer height of your choice a derive a cumulative viewshed layer where the value of each cell represents the number of peaks the cell is visible from (use r.series).

```def run_viewshed_peaks(scanned_elev, env, **kwargs):
gscript.run_command('r.geomorphon', dem=scanned_elev, forms='landforms',
search=16, skip=6, env=env)
gscript.mapcalc('peaks = if(landforms == 2, 1, null())', env=env)
gscript.run_command('r.to.vect', input='peaks', output='peaks_area', type='area', env=env)
gscript.run_command('v.to.points', input='peaks_area', output='peaks', type='centroid', flags='t', env=env)
coordinates = gscript.read_command('v.out.ascii', input='peaks', format='point', separator=',', env=env).strip()
if not coordinates:
return
i = 0
for coords in coordinates.splitlines():
print coords.split(',')[:2]
gscript.run_command('r.viewshed', input=scanned_elev, output='viewshed' + str(i),
coordinates=coords.split(',')[:2], observer_elevation=3, flags='b', env=env)
i += 1
gscript.run_command('r.series', input=['viewshed' + str(j) for j in range(i)], method='sum',
output='cumulative_viewshed', env=env)
gscript.run_command('r.colors', map='cumulative_viewshed', color='bcyr', env=env)
```

5. Find a least cost path between 2 points (for example from x=638360, y=220030 to x=638888, y=220388) where cost is defined as topographic index (trying avoid areas). Use r.topidx.

```import grass.script as gscript
def LCP(elevation, start_coordinate, end_coordinate, env):
gscript.run_command('r.topidx', input=scanned_elev, outout='topidx', env=env)
gscript.run_command('r.cost', input='topidx', output='cost', start_coordinates=start_coordinate,
outdir='outdir', flags='k', env=env)
gscript.run_command('r.colors', map='cost', color='gyr', env=env)
gscript.run_command('r.drain', input='cost', output='drain', direction='outdir',
drain='drain', flags='d', start_coordinates=end_coordinate, env=env)

if __name__ == '__main__':
elevation = 'elev_lid792_1m'
env = None
start = [638360, 220030]
end = [638888, 220388]
LCP(elevation, start, end, env)
```

6. Compute erosion with spatially variable landcover and soil erodibility (use rasters cfactorbare_1m and soils_Kfactor from the provided dataset). Reclassify the result into 7 classes based on severity of erosion and deposition:

```def run_usped(scanned_elev, env, **kwargs):
gscript.run_command('r.slope.aspect', elevation=scanned_elev, slope='slope', aspect='aspect', env=env)
gscript.run_command('r.watershed', elevation=scanned_elev, accumulation='flow_accum', threshold=1000, flags='a', env=env)
# topographic sediment transport factor
resolution = gscript.region()['nsres']
gscript.mapcalc("sflowtopo = pow(flow_accum * {res}.,1.3) * pow(sin(slope),1.2)".format(res=resolution), env=env)
# compute sediment flow by combining the rainfall, soil and land cover factors with the topographic sediment transport factor. We use a constant value of 270 for rainfall intensity factor
gscript.mapcalc("sedflow = 270. * {k_factor} * {c_factor} * sflowtopo".format(c_factor=cfactorbare_1m, k_factor=soils_Kfactor), env=env)
# compute divergence of sediment flow
gscript.run_command('r.divergence', magnitude='sedflow', direction='aspect', output='erosion_deposition', env=env)
colors = ['0% 100:0:100', '-100 magenta', '-10 red', '-1 orange', '-0.1 yellow', '0 200:255:200',
'0.1 cyan', '1 aqua', '10 blue', '100 0:0:100', '100% black']
gscript.write_command('r.colors', map='erosion_deposition',  rules='-', stdin='\n'.join(colors), env=env)
```

## Code from the TitanPad created by the workshop participants

```### Calculate Least Cost Path
import grass.script as gscript
def LCP(elevation, start_coordinate, end_coordinate, env):
gscript.run_command('r.slope.aspect', elevation=elevation, slope='slope', env=env)
gscript.run_command('r.cost', input='slope', output='cost', start_coordinates=start_coordinate,
outdir='outdir', flags='k', env=env)
gscript.run_command('r.colors', map='cost', color='gyr', env=env)
gscript.run_command('r.drain', input='cost', output='drain', direction='outdir',
drain='drain', flags='d', start_coordinates=end_coordinate, env=env)

if __name__ == '__main__':
elevation = 'elev_lid792_1m'
#    gscript.run_command('g.region', raster=elevation, flags='p')
env = None
start = [638469, 220070]
end = [638928, 220472]
LCP(scanned_elev, start, end, env)

##### ASPECT with reclassification into 8 classes
import grass.script as gscript

def run_aspect(scanned_elev, env, **kwargs):
# first we need to compute x- and y-derivatives
gscript.run_command('r.slope.aspect', elevation=scanned_elev, aspect='elev_aspect', env=env)
rules=['0:45:1', '45:90:2', '90:135:3', '135:180:4', '180:225:5', '225:270:6', '270:315:7', '315:360:8']
gscript.write_command('r.recode', input='elev_aspect', output='aspect_class', rules='-', stdin='\n'.join(rules), env=env)

if __name__ == '__main__':
import os
os.environ['GRASS_OVERWRITE'] = '1'
elevation = 'elev_lid792_1m'
run_aspect(scanned_elev=elevation, env=None)

## CENTER-POINT VIEWSHED
def centerPointViewshed(scanned_elev, env, **kwargs):
gscript.run_command('g.region', raster=scanned_elev)
# GET COORDINATES FOR CENTER OF INPUT RASTER
center = gscript.parse_command('g.region', raster=scanned_elev, flags='c',)
for k,v in center.iteritems():
if "north" in k:
north = k.split(" ")[-1]
if "east" in k:
east = k.split(" ")[-1]
print east,north
#ASSIGN COORDIATES AS east,north STRING AS PER r.viewshed REQUIREMENTS
centerPoint = [east,north]
print centerPoint
gscript.run_command('r.viewshed', input=scanned_elev, output=scanned_elev+'_viewshed',coordinates=centerPoint,observer_elevation="2.0")

if __name__ == '__main__':
elevation = 'elev_lid792_1m'
env = None
centerPointViewshed(scanned_elev=elevation, env=env)
## ...CENTER-POINT VIEWSHED

def run_waterflow(scanned_elev, env, **kwargs):
# first we need to compute x- and y-derivatives
gscript.run_command('r.slope.aspect', elevation=scanned_elev, dx='scan_dx', dy='scan_dy', env=env)
gscript.run_command('r.sim.water', elevation=scanned_elev, dx='scan_dx', dy='scan_dy',rain_value=150, depth='flow', env=env)

if __name__ == '__main__':
import os
os.environ['GRASS_OVERWRITE'] = '1'
elevation = 'elev_lid792_1m'
run_waterflow(scanned_elev=elevation, env=None)

### TOPINDEX CALCULATION
import grass.script as gscript

def run_topidx(scanned_elev, env, **kwargs):
# first we need to compute x- and y-derivatives
gscript.run_command('r.topidx', input='scanned_elev', output='topidx', env=env)

if __name__ == '__main__':
import os
os.environ['GRASS_OVERWRITE'] = '1'
elevation = 'elev_lid792_1m'
run_topidx(scanned_elev=elevation, env=None)

## PAC
import grass.script as gscript

def run_slope(scanned_elev, env, **kwargs):
gscript.run_command('r.slope.aspect', elevation=scanned_elev, slope='slope', env=env)

if __name__ == '__main__':
import os
elevation = 'elev_lid792_1m'
env = None
os.environ['GRASS_OVERWRITE'] = '1'
run_slope(scanned_elev=elevation, env=env)

##PAC

##
import grass.script as gscript

def run_solar_radiation(scanned_elev, env, **kwargs):
# convert date to day of year
import datetime
doy = datetime.datetime(2016, 5, 2).timetuple().tm_yday
# precompute slope and aspect
gscript.run_command('r.slope.aspect', elevation=scanned_elev, slope='slope', aspect='aspect', env=env, overwrite=True)
gscript.run_command('r.sun', elevation=scanned_elev, slope='slope', aspect='aspect', beam_rad='beam', step=1, day=doy, env=env, overwrite=True)
gscript.run_command('r.colors', map='beam', color='grey', flags='e')

if __name__ == '__main__':
elevation = 'elev_lid792_1m'
env = None

#Robert Dzur
def run_curvatures(scanned_elev, env, **kwargs):
gscript.run_command('r.param.scale', input=scanned_elev, output='landforms1',method='feature', size=9, env=env)
gscript.run_command('r.geomorphon', dem=scanned_elev, forms='landforms2',
search=16, skip=6, env=env)

##Damon
import grass.script as gscript

def run_lake(scanned_elev, env, **kwargs):
coordinates = [638830, 220150]
gscript.run_command('r.lake', elevation=scanned_elev, lake='output_lake',
coordinates=coordinates, water_level=120, env=env)

if __name__ == '__main__':
elevation = 'elev_lid792_1m'
env = None
run_lake(scanned_elev=elevation, env=env)

## PAC Solar
import grass.script as gscript

def run_solar_radiation(scanned_elev, env, **kwargs):
# convert date to day of year
import datetime
doy = datetime.datetime(2016, 5, 2).timetuple().tm_yday
# precompute slope and aspect
gscript.run_command('r.slope.aspect', elevation=scanned_elev, slope='slope', aspect='aspect', env=env)
gscript.run_command('r.sun', elevation=scanned_elev, slope='slope', aspect='aspect', beam_rad='beam', step=1, day=doy, env=env)
gscript.run_command('r.colors', map='beam', color='grey', flags='e')

if __name__ == '__main__':
import os
elevation = 'elev_lid792_1m'
env = None
os.environ['GRASS_OVERWRITE'] = '1'

#PAC Solar end

## Show Differences from Wiki Example

import grass.script as gscript
import os

def run_difference(real_elev, scanned_elev, env, **kwargs):
regression_params = gscript.parse_command('r.regression.line', flags='g', mapx=scanned_elev, mapy=real_elev, env=env)
gscript.mapcalc('{regression} = {a} + {b} * {before}'.format(a=regression_params['a'], b=regression_params['b'],
before=scanned_elev, regression='regression'), env=env)
gscript.mapcalc('{difference} = {regression} - {after}'.format(regression='regression', after=real_elev, difference='diff'), env=env)
gscript.write_command('r.colors', map='diff', rules='-', stdin="-100 black\n-20 red\n0 white\n20 blue\n100 black", env=env)

if __name__ == '__main__':
elevation = 'elev_lid792_1m'
sand_pile = 'sand_pile'
env = None
os.environ['GRASS_OVERWRITE'] = '1'
gscript.run_command('g.region', raster=elevation, flags='p')
gscript.run_command('r.surf.fractal', output=sand_pile)
run_difference(real_elev=elevation, scanned_elev=sand_pile, env=env)

## End Show Differences

#Topographic index (for real?)
import grass.script as gscript

def run_topoi(scanned_elev, env, **kwargs):
gscript.run_command('r.topidx', input=scanned_elev, output='topoi', env=env, overwrite=True)

if __name__ == '__main__':
elevation = 'elev_lid792_1m'
env = None
run_topoi(scanned_elev=elevation, env=env)

#CAS
def run_watershed_slope(scanned_elev, env, **kwargs):
gscript.run_command('r.watershed', elevation=scanned_elev, accumulation='flow_accum',
basin='watersheds', threshold=1000, env=env)
gscript.run_command('r.slope.aspect', elevation=scanned_elev, slope='slope', env=env)
gscript.run_command('r.stats.zonal', base='watersheds', cover='slope', method='average',
output='watersheds_slope', env=env)
gscript.run_command('r.colors', map='watersheds_slope', color='bgyr', env=env, flags='w')

#Nick B

# TASSIA
def run_usped(scanned_elev, env, **kwargs):
gscript.run_command('r.slope.aspect', elevation=scanned_elev, slope='slope', aspect='aspect', env=env)
gscript.run_command('r.watershed', elevation=scanned_elev, accumulation='flow_accum', threshold=1000, flags='a', env=env)
# topographic sediment transport factor
resolution = gscript.region()['nsres']
gscript.mapcalc("sflowtopo = pow(flow_accum * {res}.,1.3) * pow(sin(slope),1.2)".format(res=resolution), env=env)
# compute sediment flow by combining the rainfall, soil and land cover factors with the topographic sediment transport factor. We use a constant value of 270 for rainfall intensity factor
gscript.mapcalc("sedflow = 270. * {k_factor} * {c_factor} * sflowtopo".format(c_factor=0.05, k_factor=0.1), env=env)
# compute divergence of sediment flow
gscript.run_command('r.divergence', magnitude='sedflow', direction='aspect', output='erosion_deposition', env=env)
colors = ['0% 100:0:100', '-100 magenta', '-10 red', '-1 orange', '-0.1 yellow', '0 200:255:200',
'0.1 cyan', '1 aqua', '10 blue', '100 0:0:100', '100% black']
gscript.write_command('r.colors', map='erosion_deposition',  rules='-', stdin='\n'.join(colors), env=env)

#Robert Dzur - ACC
def run_accum(scanned_elev, env, **kwargs):
gscript.run_command('r.watershed', elevation=scanned_elev, accumulation='elev_acc', threshold=1000)

#Sam Sifleet
def run_viewshed(scanned_elev, env, **kwargs):
coordinates = [638830, 220150]
gscript.run_command('r.viewshed', input=scanned_elev, output='viewshed', coordinates=coordinates, observer_elevation=1.75, flags='b', env=env)
gscript.run_command('r.colors', map='viewshed', color='grey')

#Jen Lishman
# from Wiki Least Cost Path Example
import grass.script as gscript
def LCP(elevation, start_coordinate, end_coordinate, env):
gscript.run_command('r.slope.aspect', elevation=elevation, slope='slope', env=env)
gscript.run_command('r.cost', input='slope', output='cost', start_coordinates=start_coordinate,
outdir='outdir', flags='k', env=env)
gscript.run_command('r.colors', map='cost', color='gyr', env=env)
gscript.run_command('r.drain', input='cost', output='drain', direction='outdir',
drain='drain', flags='d', start_coordinates=end_coordinate, env=env)
if __name__ == '__main__':
elevation = 'elev_lid792_1m'
env = None
start = [638469, 220070]
end = [638928, 220472]
LCP(elevation, start, end, env)

##PAC Cost`
import grass.script as gscript

def LCP(scanned_elevation, start_coordinate, end_coordinate, env):
gscript.run_command('r.slope.aspect', elevation=scanned_elevation, slope='slope', env=env)
gscript.run_command('r.cost', input='slope', output='cost_pac', start_coordinates=start_coordinate,
outdir='outdir', flags='k', env=env)
gscript.run_command('r.colors', map='cost_pac', color='ramp', env=env)
gscript.run_command('r.drain', input='cost_pac', output='drain', direction='outdir',
drain='drain', flags='d', start_coordinates=end_coordinate, env=env)

if __name__ == '__main__':
import os
os.environ['GRASS_OVERWRITE'] = '1'
elevation = 'elev_lid792_1m'
env = None
start = [638469, 220020]
end = [638928, 220472]
LCP(elevation, start, end, env)
## PAC COST
```