Difference between revisions of "Talk:Using GRASS GIS through Python and tangible interfaces (workshop at FOSS4G NA 2016)"

From GRASS-Wiki
Jump to navigation Jump to search
(→‎TODO: done, link added)
 
(6 intermediate revisions by 2 users not shown)
Line 1: Line 1:
=TODO=
* link TL presentation
= Solutions =
= Solutions =


Line 31: Line 28:


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).
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).
<source lang="python">
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)
</source>


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.
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.
<source lang="python">
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)
</source>


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:
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:
<source lang="python">
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)
</source>
== Code from the TitanPad created by the workshop participants ==
<source lang="python">
### 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
    run_solar_radiation(scanned_elev='elevation', env=env)
#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'
    run_solar_radiation(scanned_elev=elevation, env=env)
#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
</source>

Latest revision as of 02:30, 24 March 2017

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
    run_solar_radiation(scanned_elev='elevation', env=env)

#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'
    run_solar_radiation(scanned_elev=elevation, env=env)

#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