Calculating Profile Information in GRASS and Python

This post describes how to calculate raster profile information (altitude changes to be exact) for a road vector layer from a DEM. DEM source is NASA’s free SRTM data.

The following script is based on scw’s answer to my related question on gis.stackexchange. This script is supposed to be run from inside a GRASS console (as I haven’t figured out yet how to import grass into python otherwise).

The idea behind this is to create points along the input road vectors and then sample the DEM values at those points. The resulting 3D points are then evaluated in a python function and the results inserted back into GRASS. There, the new values (altitude changes up and down hill) are joined back to the road elements.

import grass.script as grass
from pysqlite2 import dbapi2 as sqlite
from datetime import datetime

class AltitudeCalculator():
    def __init__(this, point_file):
        this.file = open(point_file,'r')
        this.prev_elevation = None
        this.prev_id = None
        this.increase = 0
        this.decrease = 0
        this.id = None
        this.output = ''
    
    def calculate(this):
        for line in this.file:
            line = line.split('|')
        
            #x = float(line[0])
            #y = float(line[1])
            elev = float(line[2]) 
            this.id = int(line[3])
        
            if this.prev_id:
                if this.prev_id == this.id:
                    # continue this road id
                    change = elev - this.prev_elevation 
                    if change > 0:
                        this.increase += change
                    else:
                        this.decrease += change   
                else:
                    # end of road id, save current data and reset
                    this.addOutputLine()
                    this.reset()
                    
            this.prev_id = this.id
            this.prev_elevation = elev
        
        this.addOutputLine() # last entry  

    def reset(this):
        this.increase = 0
        this.decrease = 0 
        
    def addOutputLine(this):
        this.output += "%i,%i,%i\n" % (this.prev_id, int(round(this.increase)), int(round(this.decrease)))
        
    def save(this, file_name):
        # create .csv file
        file = open(file_name,'w')
        file.write('road_id,sum_up,sum_down\n'+this.output) 
        file.close()
        # create .csvt file describing the columns
        file = open(file_name+'t','w')
        file.write('"Integer","Integer","Integer"')
        file.close()
        
def main():
    road_file = '/home/user/maps/roads.shp'
    elevation_file = '/home/user/maps/N48E016.hgt'
    point_file = '/home/user/maps/osm_road_pts_3d.txt'
    alt_diff_file = '/home/user/maps/alt_diff.csv'
    db_file = '/home/user/grass/newLocation/test/sqlite.db'
    
    print('db.connect: '+str(datetime.now()))
    grass.run_command("db.connect", driver='sqlite', database=db_file)
    
    # import files    
    print('import files: '+str(datetime.now()))
    grass.run_command("r.in.gdal", input=elevation_file, output='srtm_dem', overwrite=True)
    grass.run_command("v.in.ogr", dsn=road_file, output='osm_roads', min_area=0.0001, snap=-1, overwrite=True)
    
    # set resolution to match DEM
    print('set resolution to match DEM: '+str(datetime.now()))
    grass.run_command("g.region", rast='srtm_dem')

    # convert to points
    print('convert to points: '+str(datetime.now()))
    grass.run_command("v.to.points",  input='osm_roads', output='osm_road_pts', type='line', dmax=0.001, overwrite=True, flags='ivt')
    
    # extract elevation at each point
    print('extract elevation at each point: '+str(datetime.now()))
    grass.run_command("v.drape", input='osm_road_pts', output='osm_road_pts_3d', type='point', rast='srtm_dem', method='cubic', overwrite=True)
    
    # export to text files
    print('export to text files: '+str(datetime.now()))
    grass.run_command("v.out.ascii", input='osm_road_pts_3d', output=point_file, overwrite=True)
    
    # calculate height differences
    print('calculate height differences: '+str(datetime.now()))
    calculator = AltitudeCalculator(point_file)
    calculator.calculate()
    calculator.save(alt_diff_file)
    
    # import height differences into grass
    print('import height differeces into grass: '+str(datetime.now()))
    grass.run_command("db.droptable", table='alt_diff', flags='f')
    grass.run_command("db.in.ogr", dsn=alt_diff_file, output='alt_diff')
    
    # create an index
    print('create an index: '+str(datetime.now()))
    connection = sqlite.connect(db_file) 
    cursor = connection.cursor()
    sql = 'create unique index alt_diff_road_id on alt_diff (road_id)'
    cursor.execute(sql)
    connection.close()
    
    # join
    print('join: '+str(datetime.now()))
    grass.run_command("v.db.join", map='osm_roads', layer=1, column='cat', otable='alt_diff', ocolumn='road_id')
    
    print('Done: '+str(datetime.now()))

if __name__ == "__main__":
    main()