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

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