Talk:Using GRASS GIS through Python and tangible interfaces (workshop at FOSS4G NA 2016)

TODO

• link TL presentation

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()
    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)