top logo
  • Page 1 of 3 ( 54 posts )
  • >>
  • Underdarktalks about »
  • gis

Last update:
Tue Apr 22 06:41:17 2014

A Django site.

This weekend, I had the pleasure to join Tim Sutton for the second edition of the QGIS Podcast. Every episode, the podcast aims to summarize the latest mailing list discussions and greatest new features.
This episode’s topics include: new CAD tools, usability and the new UX mailing list, new QGIS user groups (QUGs), point cloud support plans, and QGIS design.

If you would like to ask a question or suggest a topic, you can write to

FOSS4G 2014 is taking off

If you want to become an active part of this year’s FOSS4G, it’s now time to start thinking about your contributions!

FOSS4G 2014 will be taking place in Portland, Oregon, USA from Sept 8th-12th. Like last year in Nottingham, there will be a regular track for presentations as well as an academic track and a series of workshops.


If you are looking for inspiration, you might want the check out last year’s programme or read about the interesting story behind this years conference logo.

A QGIS 2.2 preview

With the major release of version 2.0, QGIS is once more returning to a fast release cycle. You can find the project road map on The QGIS 2.2 release is scheduled for Feb, 21st and we are already in feature freeze. This means that now is the time to get the nightly version, do some testing and report possible bugs before the new version is being shipped.

Like for version 2.0, the QGIS team has prepared a great visual change log listing many new features. For me, one of the feature highlights is the possibility to export maps with world files from Print Composer because it means that we can finally create high-resolution, georeferenced images comfortably from within the application.

Another feature which will help save a lot of time is the ability to invert color ramps. So far, we had to recreate the color ramp or use work-arounds involving expression-based color settings to achieve the same effect.


These are just my personal favorites. If you haven’t checked out the change log yet, I certainly encourage you to have a look and decide for yourself. Also, if you find the time, please help by testing and reporting any issues you encounter. This way, we can all help to make 2.2 another successful release.

Happy new year!

and thank you for a great 2013!

It has been a very busy year between writing my first book, going to FOSS4G, writing my first journal article and continuing to write this blog. The blog view counter shows a staggering 310,000 views for 2013.

The most popular posts of 2013 were:

  1. pgRouting 2.0 for Windows quick guide
  2. Vintage map design using QGIS
  3. Group Stats tutorial
  4. the Print Composer 2.0 series
  5. and Public transport isochrones with pgRouting

All the best for 2014!

OSM quality assessment with QGIS: network length

In my previous post, I presented a Processing model to determine positional accuracy of street networks. Today, I’ll cover another very popular tool to assess OSM quality in a region: network length comparison. Here’s the corresponding slide from my FOSS4G presentation which shows an example of this approach applied to OSM and OS data in the UK:


One building block of this tool is the Total graph length model which calculates the length of a network within specified regions. Like the model for positional accuracy, this model includes reprojection steps to ensure all layers are in the same CRS before the actual geoprocessing starts:


The final Compare total graph length model combines two instances of “Total graph length” whose results are then joined to eventually calculate the length difference (lenDIFF).


As usual, you can find the models on Github. If you have any questions, don’t hesitate to ask in the comments and if you find any issues please report them on Github.

OSM quality assessment with QGIS: positional accuracy

Over the last years, research on OpenStreetMap data quality has become increasingly popular. At this year’s FOSS4G, I had the honor to present some work we did at the AIT to assess OSM quality in Vienna, Austria. In the meantime, our paper “Towards an Open Source Analysis Toolbox for Street Network Comparison” has been published for early access. Thanks to the conference organizers who made this possible! I’ve implemented comparison tools found in related OSM literature as well as new tools for oneway street and turn restriction comparison using Sextante scripts and models for QGIS 1.8. All code is available on Github to enable collaboration. If you are interested in OSM data quality research, I’d like to invite you to give the tools a try.

Since most users probably don’t have access to QGIS 1.8 anymore, I’ll be updating the tools to QGIS 2.0 Processing. I’m starting today with the positional accuracy comparison tool. It is based on a method described by Goodchild & Hunter (1997). Here’s the corresponding slide from my FOSS4G presentation:


The basic idea is to evaluate the positional accuracy of a street graph by comparing it with a reference graph. To do that, we check how much of the graph lies within a certain tolerance (buffer) of the reference graph.

The processing model uses the following input: the two street graphs which should be compared, the size of the buffer (tolerance for positional accuracy), a polygon layer with analysis regions, and the field containing the region id. This is how the model looks in Processing modeler:


First, all layers are reprojected into a common CRS. This will have to be adjusted if the tool is used in other geographic regions. Then the reference graph is buffered and – since I found that dissolving buffers directly in the buffer tool can become very slow with big datasets – the faster difference tool is used to dissolve the buffers before we calculate the graph length inside the buffer (inbufLEN) as well as the total graph length in the analysis region (totalLEN). Finally, the two results are joined based on the region id field and the percentage of graph length within the buffered reference graph (inbufPERC) is calculated. A high percentage shows that both graphs agree very well geometrically.

The following image shows the tool applied to a sample of OpenStreetMap (red) and official data published by the city of Vienna (purple) at Wien Handelskai. OSM was used as a reference graph and the buffer size was set to 10 meters.


In general, both graphs agree quite well. The percentage of the official graph within 10 meters of the OSM graph is 93% in the 20th district. In the above image, we can see that links available in OSM are not contained in the official graph (mostly pedestrian/bike links) and there seem to be some connectivity issues as well in the upper right corner of the image.

In my opinion, Processing models are a great solution to document geoprocessing work flows and share them with others. If you want to collaborate on building more models for OSM-related analysis, just leave a comment bellow.

Interview on GIS Lounge

It has been a real pleasure to chat with Caitlin Dempsey at GIS Lounge about open source GIS and how I got hooked on QGIS.

In related news: It’s great to see the many great and creative contributions to the QGIS Map Showcase on Flickr! If you have some maps you are proud of, please share them with the community. If you would like to see your image reused on the QGIS website or in other QGIS marketing material, please choose an appropriate license for your image.

I’ve also started to work on a new Processing script that identifies similar line features. It currently uses a length comparison and the Hausdorff distance between two line features to calculate the similarity value, but more on that in a future post!

Getting ready for FOSS4G

It’s almost here, the biggest open source GIS event of the year: the FOSS4G 2013 in Nottingham. It’s going to be my first visit to FOSS4G and I’m looking forward to present a project I did together with two colleagues at AIT where we compared OpenStreetMap to the official Austrian street network using tools developed in QGIS 1.8 Sextante. The presentation is scheduled for the first day and it would be great to meet you there:

I also have the honor to present Victor Olaya’s Sextante/Processing in a workshop together with Paolo Cavallini on the 20th:

I guess I owe Victor a geobeer or two ;-)

See you in Nottingham! And for those who can’t make it to the UK: I’ll try to keep you posted if the conference wifi allows it.

Print Composer 2.0 – Take #7

Today’s post: More print composer overview magic!

Inverted Map Overviews

Thanks to the “Invert overview” option, we can now chose between highlighting the detail area (left example in the image) or blocking out the surrounding area (right example).


The “Lock layers for map item” option can come in very handy if you want to reduce the number of layers in the overview map while still keeping all layers of interest in the main map.

Mapping OGDWien Population Density

The city of Vienna provides both subdistrict geometries and population statistics. Mapping the city’s population density should be straightforward, right? Let’s see …

We should be able to join on ZBEZ and SUB_DISTRICT_CODE, check! But what about the actual population counts? Unfortunately, there is no file which simply lists population per subdistrict. The file I found contains four lines for each subdistrict: females 2011, males 2011, females 2012 and males 2012. That’s not the perfect format for mapping general population density.

A quick way to prepare our input data is applying pivot tables, eg. in Open Office: The goal is to have one row per subdistrict and columns for population in 2011 and 2012:

Export as CSV, add CSVT and load into QGIS. Finally, we can join geometries and CSV table:

A quick look at the joined data confirms that each subdistrict now has a population value. But visualizing absolute values results in misleading maps. Big subdistricts with only average density will overrule smaller but much denser subdistricts:

That’s why we need to calculate population density. This is easy to do using Field Calculator. The subdistrict file already contains area values but even if they were missing, we could calculate it using the $area operator: "pop2012" / ($area / 10000). The resulting population density in population per ha finally shows which subdistricts are the most densely populated:

One could argue that this is still no accurate representation of population density: Big parts of some subdistricts are actually covered by water – especially along the Danube – and therefore uninhabited. There are also big parks which could be excluded from the subdistrict area. But that’s going to be the topic of another post.

If you want to use my results so far, you can download the GeoJSON file from Github.

WFS to PostGIS in 1 Step

This is an update to my previous post “WFS to PostGIS in 3 Steps”. Thanks to Even Rouault’s comments and improvements to GDAL, it is now possible to load Latin1-encoded WFS (like the one by into PostGIS in just one simple step.

To use the following instructions, you’ll have to get the latest GDAL (

You only need to run SDKShell.bat to set up the environment and ogr2ogr is ready for action:

C:\Users\Anita>cd C:\release-1600-gdal-mapserver
C:\release-1600-gdal-mapserver>ogr2ogr -overwrite -f PostgreSQL PG:"user=myuser password=mypassword dbname=wien_ogd" "WFS:"

Thanks everyone for your comments and help!

QGIS Server on Windows7 Step-by-step

After my successful experiment with QGIS Server on Ubuntu, I took a shot at Windows7. These are my notes on installing QGIS Server following the instructions on the wiki.


Using OSGeo4W installer it is easy to install QGIS Server: Just mark qgis-server for installation from “Web” category (in the “Advanced” installation).

All other necessary packages will be selected automatically.

As mentioned on the wiki, the next step is to tell Apache which port number to use. Apache (2.2.14-4 from OSGeo4W) does not have any default IP/port set and it fails to start. To fix this, we need to edit the file


and change

Listen @apache_port_number@

to our needs, e.g. to listen on port 80

Listen 80

The last thing I had to do to get QGIS Server working was to copy two files
libeay32.dll and ssleay32.dll
from C:\OSGeo4W\apache\bin
to C:\OSGeo4W\apps\qgis\bin.

The GetCapabilities request should work now


Adding a QGIS project file

To add a project file to the server, we stay in C:\OSGeo4W\apps\qgis\bin. If we put a project file in this directory, it will be served by default (without having to pass the optional map parameter).

For this test, I added my vienna.qgs project file. This is how my QGIS bin folder looks like (notice the .dll files we copied from Apache/bin and the project file):

Next, we have to restart Apache to force QGIS Server to load the project file. The OSGeo4W installation provides a handy “Apache-Monitor” GUI to restart Apache. If it fails, try to reboot ;)

Let’s test the setup using “Add WMS Layer” in QGIS by adding the service URL and ticking “Ignore GetMap URI …” and “Ignore GetFeature URI …”.

The project layers are now available through WMS and can be loaded into your client.

QGIS "Add WMS Layer" dialog with my newly created WMS


That wasn’t bad. The wiki page was very helpful and I didn’t encounter any real problems. Editing a config file and copying a few .dlls is easy enough.

Since linking files is not one of Windows’ strengths, I’d expect a server with multiple projects to get quite messy. But it certainly works for home use and experiments.

Beautiful Global Projections – Adding Custom Projections to QGIS

This year we are celebrating Gerardus Mercator’s 500th birthday. We have all grown very accustomed to his Mercator projection but I want to take the chance to explore some alternatives:

Radical Cartography features an extensive projection reference compiled by Bill Rankin. One of the more exotic projections is “Van der Grinten I” by Alphons J. van der Grinten, 1898. It has a “pleasant balance of shape and scale distortion”. The “boundary is a circle” and “all parallels and meridians are circular arcs (spacing of parallels is arbitrary)”.

Using the name, we can try to find the projection definition on One of the definitions that works well in QGIS is “ESRI:53029 Sphere Van der Grinten I” with the following proj4 string:

+proj=vandg +lon_0=0 +x_0=0 +y_0=0 +R_A +a=6371000 +b=6371000 +units=m +no_defs

In QGIS Settings – Custom CRS, we can add this projection to the list of available CRS:

  1. Press the “Star” button to add a new empty entry”
  2. Add the Name and proj4 string
  3. Press the “Save” button to make the changes permanent

Custom projection dialog

Once this is done, “Van der Grinten I” can be selected for on-the-fly reprojection. I’ve been using NaturalEarth’s land and ocean dataset. The result might not be a perfect circle (due to the coarseness of the dataset) but I find it very appealing:

Van der Grinten projection

A Guide to Beautiful Reliefs in QGIS

This week Sourcepole released a new addition to the Raster Terrain Analysis plugin: a sophisticated Relief tool. (More info in their announcement) This plugin is shipped with QGIS (developer version, not in 1.7.3 release) by default but you might have to activate it in Plugin Manager:

The plugin dialog is quite self-explanatory. You can chose the elevation file, output path and any of the numerous raster formats. The z factor is a bit more mysterious. We will have a look at that in a second. The rest of the dialog is the relief color editor. Pressing Create automatically will give you a color gradient to start with.

Relief tool dialog

But what’s the z factor good for?

I’ve tried a few different settings using free NASA SRTM data and it seems that higher values lead to a smoother relief (Please ignore the water areas):

z factor = 100 z factor = 1000 z factor = 10000 z factor = 50000 z factor = 100000


As Marco noted in the comments: The z factor is used if the x/y units are different from the z unit.

  • If everything is in meters, use z factor 1.0 (default).
  • If x/y is in degree and z in meters, use z factor 111120.
  • If x/y degree and z is feet, use z factor 370400.

In the example above SRTM rasters are in WGS84 with heights in meters. That’s why the result using a z factor of 100000 looks so good.

In my opinion the results look great even with the coarse SRTM dataset I used. Looking forward to all the great QGIS maps we will see in the future.

Interstate Shield Styles in QGIS

This is a follow-up to my post “Google Maps”-Style Road Maps in QGIS and an answer to @mattwigway’s comment about how to create styles for e.g. US interstates.

interstate example map

To create a style like this, you can use the “Marker line” feature. The marker can be built using a blank SVG shield and additional text markers to add the road number on top.

creation of the interstate shield marker

This marker line can be put “on central point” to only show up once along the road.

creation of the interstate style

You can find interstate and other shield style on e.g. Wikimedia Commons. With Inkscape, you can remove the number and thus create a blank marker template.

Tweets to QGIS

The idea behind this post was to create a video of twitter activity using Time Manager. You can watch the results of my first test run here:

And this is how it’s done:

First, you have to collect some tweets with location information. The following command will collect tweets within a certain geographic region from the Twitter Stream API using curl. You need a Twitter user account to use the API. (Curl comes readily available with OSGeo4W install.)

curl -k -d @locations.txt -uuser:password > tweets.json

The contents of locations.txt is the geographic extent you are interested in, e.g. for Austria:


After collecting some data, you can load the tweets into QGIS. Executing the following lines in Python Console will add an in-memory point layer to the map. (I am only extracting coordinates and time stamp from the tweets, but you can access more information through the JSON object.)

import simplejson
from PyQt4.QtCore import *
from datetime import *


# create layer
vl = QgsVectorLayer("Point", "tweets", "memory")
pr = vl.dataProvider()

# add fields
pr.addAttributes( [ QgsField("t", QVariant.String) ] )

# create features
for line in f:
      fet.setAttributeMap({0:QVariant(str(datetime.strptime(j['created_at'],'%a %b %d %H:%M:%S +0000 %Y')))})



To use the result in Time Manager, you have to export the layer to e.g. Shapefile because it’s not possible to add query strings to in-memory layers.

If you are interested in learning more about PyQGIS, you can find a lot of useful material in the PyQGIS Cookbook.

Visualizing Global Connections

Today, I’ve been experimenting with data from They offer airport, airline and route data for download. The first idea that came to mind was to connect airports on a shared route by lines. This kind of visualization just looks much nicer if the connections are curved instead of simple straight lines.

Luckily, that’s pretty easy to do using PostGIS. After loading airport positions and route data, we can create the connection lines like this (based on [postgis-users] Great circle as a linestring):

UPDATE experimental.airroutes
SET the_geom =
(SELECT ST_Transform(ST_Segmentize(ST_MakeLine(
       ST_Transform(a.the_geom, 953027),
       ST_Transform(b.the_geom, 953027)
     ), 100000 ), 4326 )
 FROM experimental.airports a, experimental.airports b
 WHERE = airroutes.source_id
   AND = airroutes.dest_id

The CRS used in the query is not available in PostGIS by default. You can add it like this (source:

INSERT into spatial_ref_sys (srid, auth_name, auth_srid, proj4text, srtext) values ( 953027, 'esri', 53027, '+proj=eqdc +lat_0=0 +lon_0=0 +lat_1=60 +lat_2=60 +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m +no_defs ', 'PROJCS["Sphere_Equidistant_Conic",GEOGCS["GCS_Sphere",DATUM["Not_specified_based_on_Authalic_Sphere",SPHEROID["Sphere",6371000,0]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Equidistant_Conic"],PARAMETER["False_Easting",0],PARAMETER["False_Northing",0],PARAMETER["Central_Meridian",0],PARAMETER["Standard_Parallel_1",60],PARAMETER["Standard_Parallel_2",60],PARAMETER["Latitude_Of_Origin",0],UNIT["Meter",1],AUTHORITY["EPSG","53027"]]');

This is an example visualization (done in QGIS) showing only flight routes starting from Vienna International Airport:

Flight routes from Vienna International

Connections crossing the date line are currently more problematic. Lines would have to be split, otherwise this is what you’ll get:

Date line trouble

QGIS Server on Windows

Want to finally try QGIS Server on your Windows system? Thankfully, Till Adams has written a new HowTo for QGIS Server on Windows. Till shows how to install OSGeo4W version and how to get it running.

If you first want to see what QGIS Server and QGIS Webclient can do, check Webgis-Uster homepage by Andreas Neumann.

Open Government Data Wien in QGIS

Vienna’s Open Government Data initiative is publishing an increasing amount of Geodata and the best thing is: They’re publishing it via open standardized services! Both WMS and WFS are available and ready to be used in QGIS.

Let’s see how we can use the data available through their WFS using the district information layer “Bezirksgrenzen” as an example. The page lists a GML and a GeoJSON version of WFS. For now, we’ll use GeoJSON.

In QGIS, the layer can be loaded using “Add vector layer” – “Protocol” and inserting the GeoJSON url there. The encoding should be changed to ISO8859-15 to account for “Umlaute”.

The loaded GeoJSON layer "Bezirksgrenzen"

Now, we have geodata. Let’s add some attribute data too! Attribute data is available in CSV format. After downloading e.g. some information on the district population, check the file content and remove excessive header lines so that there is only one header line containing attribute names left. Then, you can load the CSV file into QGIS too (“Add vector layer”).

Last step: Joining geodata and attribute data! Go to the vector layer’s properties – Join tab and add the following join relation:

Joining GeoJSON and attribute layer

Now, the attribute table of the vector layer contains the additional CSV attributes – ready for further analysis. If you want to classify based on numerical CSV attributes, you’ll have to create a .csvt file first otherwise all fields are interpreted as texts.

Works great. Thumbs up for this great initiative!

QGIS 1.7 Released!

QGIS 1.7 has landed. After some delays due to a major infrastructure overhaul, a new version of QGIS is available for download. For a list of what’s new in 1.7 check the release announcement.


Powered by Django!