Advertisements

Archive

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.

Advertisements

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.

This is something I have been wanting to do for a long time: map which areas of Vienna have fast access to a certain kind of infrastructure. Now, I finally found time and data to perform this analysis. Data used is OSM road data (Cloudmade shapefile) for Austria and metro station coordinates for Vienna by Max Kossatz and Robert Harm.

Before importing the OSM roads into PostGIS, I cut out my area of interest and created a clean topology using GRASS v.clean.break. Once loaded into the database, assign_vertex_id() function does the rest and the network is ready for routing and distance calculations.
For the metro stations, I calculated the nearest network node using George MacKerron’s Nearest Neighbor function.

Catchments were calculated using driving_distance() function. It returns distance to a given metro station for all network nodes (up to a maximum distance). The result can be interpolated to show e.g. which areas are at most 1 km away from any metro station.

1 km catchments around metro stations in Vienna

Close-up look at the 1 km catchment zone border

Once set up, performing this analysis is reasonably fast. Instead of metro stations, any other infrastructure coverage can be analyzed easily. I could imagine this being really useful when looking for a new flat: “Find me an area close to work, a metro station and a highschool.”

The next great thing would be to have all data for calculation of transit travel times too. Yes, I’m looking at you Wiener Linien!

This year, three QGIS projects made it into Google Summer of Code. Congratulations to all successful students!

These are the accepted OSGEO projects:

QGIS

GRASS

  • GRASS WXGUI WMS service rendering
  • Completion of wxGUI Nviz extension for 3D data visualization in GRASS GIS
  • r.in.modis for GRASS GIS
  • Graphical User Interface for the hydrological tools r.stream* in GRASS GIS

pgRouting

  • Multi-modal public transit routing for pgRouting
  • Time Dependent \ Dynamic Shortest Path Algorithm Implementation for pgRouting

gvSIG

  • Add support to vector data formats for gvSIG Mini
  • Design and implement an API for tiled vectorial support of geo-location data services for gvSIG Mini
    Integration of GGL2 into gvSIG

Opticks

  • Development of a ship detection and classification toolkit for SAR imagery in Opticks
  • Astronomical processing tools for Opticks
  • Photography processing tools for Opticks

Geoserver

  • Enhancing Geoserver Authentication

uDig

  • Catalog View of uDig
  • OSM data mining and editing capabilities in uDig and Geotools

MapServer

  • INSPIRE View Service for MapServer

Others

  • Geoprocessing with Neo4j Spatial
  • PyOSSIM: Python bindings for OSSIM libraries

For a full list including student names check http://www.google-melange.com/gsoc/org/google/gsoc2011/osgeo.

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

QGIS and GRASS in Local Government Bushfire Hazard Mapping – A Case Study.

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()
%d bloggers like this: