Archive

QGIS

On my quest to create test data for spatial statistics, I’ve discovered income data for Austria per municipality on a news paper website:

Screenshot 2014-11-29 23.06.46

For further analysis, I decided to limit the area to Vienna and Lower Austria. Since the income data included GKZ “Gemeindekennzahl” IDs, it was possible to join them to municipalities extracted from OpenStreetMap using QuickOSM for QGIS. GRASS v.clean was used to clean the vector topology to the point where PySAL was able to compute spatial weights.

Using PySAL, I then computed income clusters: blue regions represent low clusters while red regions represent high clusters …

Municipality border data (c) OpenStreetMap and contributors Income data source: Statistik Austria via derStandard

Municipality border data (c) OpenStreetMap and contributors
Income data source: Statistik Austria via derStandard

The results show a statistically significant cluster of low income in the north west, in the area called Waldviertel, as well as a cluster of high income containing many of the municipalities surrounding Vienna, an area often referred to as the “Speckgürtel” (“bacon belt”).

As promised in my recent post “Experiments with Conway’s Game of Life”, I have been been looking into how to improve my first implementation. The new version which you can now find on Github is fully contained in one Python script which runs in the QGIS console. Additionally, the repository contains a CSV with the grid definition for a Gosper glider gun and the layer style QML.

Rather than creating a new Shapefile for each iteration like in the first implementation, this script uses memory layers to save the game status.

You can see it all in action in the following video:

(video available in HD)

Thanks a lot to Nathan Woodrow for the support in getting the animation running!

Sometimes there are still hick-ups causing steps to be skipped but overall it is running nicely now. Another approach would be to change the layer attributes rather than creating more and more layers but I like to be able to go through all the resulting layers after they have been computed.

This experiment is motivated by a discussion I had with Dr. Claus Rinner about introducing students to GIS concepts using Conway’s Game of Life. Conway’s Game of Life is a popular example to demonstrate cellular automata. Based on an input grid of “alive” and “dead” cells, new cell values are computed on each iteration based on four simple rules for the cell and its 8 neighbors:

  1. Any live cell with fewer than two live neighbours dies, as if caused by under-population.
  2. Any live cell with two or three live neighbours lives on to the next generation.
  3. Any live cell with more than three live neighbours dies, as if by overcrowding.
  4. Any dead cell with exactly three live neighbours becomes a live cell, as if by reproduction.

(Source: Wikipedia – Conway’s Game of Life)

Based on these simple rules, effects like the following “glider gun” can be achieved:

Gospers glider gun.gif
Gospers glider gun” by KieffOwn work. Licensed under CC BY-SA 3.0 via Wikimedia Commons.

There are some Game of Life implementations for GIS out there, e.g. scripts for ArcGIS or a module for SAGA. Both of these examples are raster-based. Since I couldn’t find any examples of raster manipulation like this in pyQGIS, I decided to instead implement a vector version: a Processing script which receives an input grid of cells and outputs the next iteration based on the rules of Game of Life. In the following screencast, you can see the Processing script being called repeatedly by a script from the Python console:

So far, it’s a quick and dirty first implementation. To make it more smooth, I’m considering adding spatial indexing and using memory layers instead of having Processing create a bunch of Shapefiles.

It would also be interesting to see a raster version done in PyQGIS. Please leave a comment if you have any ideas how this could be achieved.

Today’s post is inspired by a recent thread on the QGIS user mailing list titled “exporting text to Illustrator?”. The issue was that with the introduction of the new labeling system, all labels were exported as paths when creating an SVG. Unnoticed by almost everyone (and huge thanks to Alex Mandel for pointing out!) an option has been added to 2.4 by Larry Shaffer which allows exporting labels as texts again.

To export labels as text, open the Automatic Placement Settings (button in the upper right corner of the label dialog) and uncheck the Draw text as outlines option.

Screenshot 2014-09-20 21.03.26

Note that we are also cautioned that

For now the developers recommend you only toggle this option right
before exporting
and that you recheck it after.

Alex even recorded a video showcasing the functionality:

When mapping flows or other values which relate to a certain direction, styling these layers gets interesting. I faced the same challenge when mapping direction-dependent error values. Neighboring cell pairs were connected by two lines, one in each direction, with an associated error value. This is what I came up with:

srtm_errors_1200px

Each line is drawn with an offset to the right. The size of the offset depends on the width of the line which in turn depends on the size of the error. You can see the data-defined style properties here:

directed_error_style

To indicate the direction, I added a marker line with one > marker at the center. This marker line also got assigned the same offset to match the colored line bellow. I’m quite happy with how these turned out and would love to hear about your approaches to this issue.

srtm_errors_detail

These figures are part of a recent publication with my AIT colleagues: A. Graser, J. Asamer, M. Dragaschnig: “How to Reduce Range Anxiety? The Impact of Digital Elevation Model Quality on Energy Estimates for Electric Vehicles” (2014).

In my opinion, Stamen’s Toner-lite map is one of the best background maps to use together with colorful overlays. The only downsides of using it in QGIS are that the OpenLayers plugin can not provide the tiles at print resolution and that the projection is limited to Web Mercator. That’s why I’ve started to recreate the style for OSM Spatialite databases:

toner-lite

So far, there are styles for lines and polygons and they work quite well for the scale range between 1:1 and 1:250000. As always, you can download the styles from QGIS-resources on Github.

Using OSM data in QGIS is a hot topic but so far, no best practices for downloading, preprocessing and styling the data have been established. There are many potential solutions with all their advantages and disadvantages. To give you a place to start, I thought I’d share a workflow which works for me to create maps like the following one from nothing but OSM:

osm_google_100k

Getting the data

Raw OSM files can be quite huge. That’s why it’s definitely preferable to download the compressed binary .pbf format instead of the XML .osm format.

As a download source, I’d recommend Geofabrik. The area in the example used in this post is part of the region Pays de la Loire, France.

Preparing the data for QGIS

In the preprocessing step, we will extract our area of interest and convert the .pbf into a spatialite database which can be used directly in QGIS.

This can be done in one step using ogr2ogr:

C:\Users\anita_000\Geodata\OSM_Noirmoutier>ogr2ogr -f "SQLite" -dsco SPATIALITE=YES -spat 2.59 46.58 -1.44 47.07 noirmoutier.db noirmoutier.pbf

where the -spat option controls the area of interest to be extracted.

When I first published this post, I suggested a two step approach. You can find it here for future reference:

For the first step: extracting the area of interest, we need Osmosis. (For Windows, you can get osmosis from openstreetmap.org. Unpack to use. Requires Java.)

When you have Osmosis ready, we can extract the area of interest to the .osm format:

C:\Users\anita_000\Geodata\OSM_Noirmoutier>..\bin\osmosis.bat --read-pbf pays-de-la-loire-latest.osm.pbf --bounding-box left=-2.59 bottom=46.58 right=-1.44 top=47.07 --write-xml noirmoutier.osm

While QGIS can also load .osm files, I found that performance and access to attributes is much improved if the .osm file is converted to spatialite. Luckily, that’s easy using ogr2ogr:

C:\Users\anita_000\Geodata\OSM_Noirmoutier>ogr2ogr -f "SQLite" -dsco SPATIALITE=YES noirmoutier.db noirmoutier.osm

Finishing preprocessing in QGIS

In QGIS, we’ll want to load the points, lines, and multipolygons using Add SpatiaLite Layer:

Screenshot 2014-05-31 11.39.40

When we load the spatialite tables, there are a lot of features and some issues:

  • There is no land polygon. Instead, there are “coastline” line features.
  • Most river polygons are missing. Instead there are “riverbank” line features.

Screenshot 2014-05-31 11.59.58

Luckily, creating the missing river polygons is not a big deal:

  1. First, we need to select all the lines where waterway=riverbank.
    Screenshot 2014-05-31 13.14.00
  2. Then, we can use the Polygonize tool from the processing toolbox to automatically create polygons from the areas enclosed by the selected riverbank lines. (Note that Processing by default operates only on the selected features but this setting can be changed in the Processing settings.)
    Screenshot 2014-05-31 13.40.16

Creating the land polygon (or sea polygon if you prefer that for some reason) is a little more involved since most of the time the coastline will not be closed for the simple reason that we are often cutting a piece of land out of the main continent. Therefore, before we can use the Polygonize tools, we have to close the area. To do that, I suggest to first select the coastline using "other_tags" LIKE '%"natural"=>"coastline"%' and create a new layer from this selection (save selection as …) and edit it (don’t forget to enable snapping!) to add lines to close the area. Then polygonize.

Screenshot 2014-05-31 14.38.48

Styling the data

Now that all preprocessing is done, we can focus on the styling.

You can get the styles used in the map from my Github QGIS-resources repository:

  • osm_spatialite_googlemaps_multipolygon.qml … rule-based renderer incl. rules for: water, natural, residential areas and airports
  • osm_spatialite_googlemaps_lines.qml … rule-based renderer incl. rules for roads, rails, and rivers, as well as rules for labels
  • osm_spatialite_googlemaps_roadshields.qml … special label style for road shields
  • osm_spatialite_googlemaps_places.qml … label style for populated places such as cities and towns

qgis_osm_google_100k

This post shows how to quickly and easily create a small QGIS plugin for counting the number of features within a vector layer.

To get started, you will need QGIS and Qt Designer (to design the user interface) installed. If you are on Windows, I suggest WinPython which provides Qt Designer and Spyder (a Python IDE).

The great thing about creating plugins for QGIS: There is a plugin for that! It’s called Plugin Builder. And while you are at it, also install Plugin Reloader. Reloader is very useful for plugin developers because it lets you quickly reload your plugin without having to restart QGIS every time you make changes to the code.

installPluginBuilder

Plugin Builder will create all the files we need for our plugin. Just start it and select a name for your plugin class (one word in CamelCase), as well as a name for the plugin itself and the plugin menu entry (can be multiple words). Once you press Ok, you’re asked to select a folder to store the plugin. You can save directly to the QGIS plugin folder ~\.qgis2\python\plugins.

pluginBuilder

Next, open the newly created folder (in my case ~\.qgis2\python\plugins\BuilderTest). Amongst other files, it contains the user interface file ui_buildertest.ui. Our plugin will count the number of features in a vector layer. Therefore, it needs a combobox which allows the user to select a layer. Open the .ui file in Qt Designer and add a combobox to the dialog. Change the object name of the combobox to layerCombo. We’ll later use this name in the plugin code to add items to the combobox. Save the dialog and close Qt Designer.

qtDesigner

Now, we need to compile the .ui and the resources.qrc file to turn the dialog and the icon into usable Python code. This is done on the command line. On Windows, I suggest using the OSGeo4W Shell. Navigate to the plugin folder and run:

pyuic4 -o ui_buildertest.py ui_buildertest.ui
pyrcc4 -o resources_rc.py resources.qrc

Note: Using the latest version of Plugin Builder, you only need to compile resources.qrc because the .ui is now loaded dynamically. Furthermore, you should use the following command to compile the resources:

pyrcc4 -o resources.py resources.qrc 

If you enable and run the plugin now, you will already see the dialog but the combobox will be empty. To populate the combobox, we need to write a few lines of code in buildertest.py. First, we’ll fetch all loaded layers and add all vector layers to the combobox. Then, we’ll add code to compute and display the number of features in the selected layer. To achieve this, we expand the run() method:

def run(self):        
    # show the dialog
    self.dlg.show()

    layers = QgsMapLayerRegistry.instance().mapLayers().values()
    for layer in layers:
        if layer.type() == QgsMapLayer.VectorLayer:
            self.dlg.layerCombo.addItem( layer.name(), layer ) 
         
    # Run the dialog event loop
    result = self.dlg.exec_()
    # See if OK was pressed
    if result == 1:
        # do something useful 
        index = self.dlg.layerCombo.currentIndex()
        layer = self.dlg.layerCombo.itemData(index)
        QMessageBox.information(self.iface.mainWindow(),"hello world","%s has %d features." %(layer.name(),layer.featureCount()))

When you are done with the code, you can use Plugin Reloader to load the new version. When you start the plugin now, the combobox will be populated with the names of the vector layers in your current project. And on pressing Ok, the plugin will compute and display the number of features.

builderTEst

builderTestResult

For more information on PyQGIS and more code samples I warmly recommend the PyQGIS Cookbook. Have fun!

Extracting POIs from OpenStreetMap is reasonably simple using Overpass API. A very convenient way to construct the query is to use a query builder which allows you to select the area of interest and builds queries for different servers.

xapi_query_builder

Of course you can fine-tune the query further. For example, you can add multiple key-value pairs to the query. I used the following query to select all Billa supermarkets:

www.overpass-api.de/api/xapi?*[shop=supermarket][name=Billa][bbox=15.96725,48.0432,16.79947,48.40915]

Note the * in the query? It means that I’m querying all kinds of features: nodes, ways, and relations.

Save the server response to a .osm file. This file can be loaded into QGIS using simple drag-and-drop or Add Vector Layer. A dialog will open where you can select the type of features you want to load from the file. You can simply use Select All and OK to load everything.

add_osm_file

My supermarket POIs came in two types: points and multipolygons. To style them both with nice supermarket SVG icons, I decided to use a Centroid fill with the SVG marker for the polygon layer:

osm_pois_billa

Open data and open source GIS … nice :-)

Today was the last day of the Vienna code sprint which brought together OSGeo developers from many projects. It’s been a great week thanks to organizers and sponsors!

The QGIS team was extremely busy working on the project’s web infrastructure (e.g. new plugins.qgis.org website) as well as hunting down and fixing bugs.

Check out some impressions on twitter.

qgis23_vienna

More pictures on the official blog: vienna2014.sprint.osgeo.org