Archive

Tag Archives: GRASS

A few weeks ago, the city of Vienna released a great dataset: the so-called “Flächen-Mehrzweckkarte” (FMZK) is a polygon vector layer with an amazing level of detail which contains roads, buildings, sidewalk, parking lots and much more detail:

preview of the Flächen-Mehrzweckkarte

preview of the Flächen-Mehrzweckkarte

Now, of course we can use this dataset to create gorgeous maps but wouldn’t it be great to use it for analysis? One thing that has been bugging me for a while is routing for pedestrians and how it’s still pretty bad in many situations. For example, if I’d be looking for a route from the northern to the southern side of the square in the previous screenshot, the suggestions would look something like this:

Pedestrian routing in Google Maps

Pedestrian routing in Google Maps

… Great! Google wants me to walk around it …

Pedestrian routing on openstreetmap.org

Pedestrian routing on openstreetmap.org

… Openstreetmap too – but on the other side :P

Wouldn’t it be nice if we could just cross the square? There’s no reason not to. The routing graphs of OSM and Google just don’t contain a connection. Polygon datasets like the FMZK could be a solution to the issue of routing pedestrians over squares. Here’s my first attempt using GRASS r.walk:

Routing with GRASS r.walk

Routing with GRASS r.walk (Green areas are walk-friendly, yellow/orange areas are harder to cross, and red buildings are basically impassable.)

… The route crosses the square – like any sane pedestrian would.

The key steps are:

  1. Assigning pedestrian costs to different polygon classes
  2. Rasterizing the polygons
  3. Computing a cost raster for moving using r.walk
  4. Computing the route using r.drain

I’ve been using GRASS 7 for this example. GRASS 7 is not yet compatible with QGIS but it would certainly be great to have access to this functionality from within QGIS. You can help make this happen by supporting the crowdfunding initiative for the GRASS plugin update.

How do you objectively define and compute which parts of a network are in the center? One approach is to use the concept of centrality.

Centrality refers to indicators which identify the most important vertices within a graph. Applications include identifying the most influential person(s) in a social network, key infrastructure nodes in the Internet or urban networks, and super spreaders of disease. (Source: http://en.wikipedia.org/wiki/Centrality)

Researching this topic, it turns out that some centrality measures have already been implemented in GRASS GIS. thumbs up!

v.net.centrality computes degree, betweeness, closeness and eigenvector centrality.

As a test, I’ve loaded the OSM street network of Vienna and run

v.net.centrality -a input=streets@anita_000 output=centrality degree=degree closeness=closeness betweenness=betweenness eigenvector=eigenvector

grass_centrality

The computations take a while.

In my opinion, the most interesting centrality measures for this street network are closeness and betweenness:

Closeness “measures to which extent a node i is near to all the other nodes along the shortest paths”. Closeness values are lowest in the center of the network and higher in the outskirts.

Betweenness “is based on the idea that a node is central if it lies between many other nodes, in the sense that it is traversed by many of the shortest paths connecting couples of nodes.” Betweenness values are highest on bridges and other important arterials while they are lowest for dead-end streets.

(Definitions as described in more detail in Crucitti, Paolo, Vito Latora, and Sergio Porta. “Centrality measures in spatial networks of urban streets.” Physical Review E 73.3 (2006): 036125.)

Centrality: low values in pink, high values in green

Centrality: low values in pink, high values in green

Works great! Unfortunately, v.net.centrality is not yet part of the QGIS Processing GRASS toolbox. It would certainly be a great addition.

The following video by soerengebbert shows the latest developments in OS WPS applications: PyWPS (using wps-grass-bridge to integrate GRASS modules as WPS processes), GRASS 7, and QGIS using the qwps plugin by Horst Düster (all latest svn versions).

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

You can call almost any function in Python, using “grass.run_command(“function”, ….)”.

For your inspiration, you find GRASS Python scripts in their SVN.

%d bloggers like this: