In the category “last night on Twitter”, a challenge I couldn’t resist: creating illuminated contours (aka Tanaka contours) in QGIS. Daniel P. Huffman started the thread by posting this great example:


This was quickly picked up by Hannes Kröger who blogged about his first attempt at reproducing the effect using QGIS and GIMP. Obviously, that left the challenge of finding a QGIS-only solution.

Everything that’s needed to create this effect is a DEM. As Hannes describes in his post, the DEM can then be used to compute the contour lines, e.g. with Raster | Extraction | Contour:

gdal_contour -a ELEV -i 100.0 C:\Users\anita\Geodata\misc\mt-st-helens\ C:/Users/anita/Geodata/misc/mt-st-helens/countours

Screenshot 2015-05-24 11.17.49


In order to be able to compute the brightness of the illuminated contours, we need to compute the orientation of every subsection of the contours. Therefore, we need to split the contour lines at each node. One way to do this is using v.split from the Processing toolbox:

Screenshot 2015-05-24 11.23.11

When we split the contours and visualize the result using arrows, we can see that they all wrap around the mountain in clockwise direction (light DEM cells equal higher elevation):


After the split, we can compute the orientation of the contour subsections using, for example, a user-defined function:

Screenshot 2015-05-24 19.09.12

This function can then be used in a Field calculator expression:

Screenshot 2015-05-24 19.11.53

Based on the orientation, we can then write an expression to control the contour line color. For example, if we want the sun to appear in the north west (-45°) we can use:

color_hsl( 0,0, 
  scale_linear( abs(
    ( CASE WHEN "azimuth"-45 < 0
      THEN "azimuth"-45+360 
      ELSE "azimuth"-45
    END )
  -180), 0, 180, 0, 100)

This will color the lines which are directly exposed to the sun white hsl(0,0,100) while the ones in the shadows will be black hsl(0,0,0).

Screenshot 2015-05-24 11.55.50

Use the Overlay layer blending mode to blend contours and DEM color:


The final step, to get as close to the original design as possible, is to create the effect of discrete elevation classes instead of a smooth color gradient. This can easily be achieved by changing the color interpolation mode of the DEM from Linear to Discrete:

Screenshot 2015-05-24 12.11.01

This leaves us with the following gorgeous effect:


As Hannes pointed out, another important aspect of Tanaka’s method is to also alter the contour line width. Lines in the sun or shadow should be wider (1 in this example) than those in orthogonal direction (0.2 in this example):

abs( abs(
  ( CASE WHEN "azimuth"-45 < 0
    THEN  "azimuth"-45+360
    ELSE  "azimuth"-45
  END )
-180) -90),
0, 90, 0.2, 1)



Today was the final day of #QGIS2015 the first joint QGIS conference and developer meeting. I had the pleasure to meet Time Manager co-developer Karolina Alexiou aka carolinux in person and give a talk including a hands-on workshop on Time Manager together. Time Manager makes it possible to explore spatio-temporal data by creating animations directly in QGIS.

The talk presents QGIS visualization tools with a focus on efficient use of layer styling to both explore and present spatial data. Examples include the recently added heatmap style as well as sophisticated rule-based and data-defined styles. The focus of this presentation is exploring and presenting spatio-temporal data using the Time Manager plugin. A special treat are time-dependent styles using expression-based styling which access the current Time Manager timestamp.

To download the example data and QGIS preojcts download

Today’s post is a short tutorial for creating trajectory animations with a fadeout effect using QGIS Time Manager. This is the result we are aiming for:

The animation shows the current movement in pink which fades out and leaves behind green traces of the trajectories.

About the data

GeoLife GPS Trajectories were collected within the (Microsoft Research Asia) Geolife project by 182 users in a period of over three years (from April 2007 to August 2012). [1,2,3] The GeoLife GPS Trajectories download contains many text files organized in multiple directories. The data files are basically CSVs with 6 lines of header information. They contain the following fields:

Field 1: Latitude in decimal degrees.
Field 2: Longitude in decimal degrees.
Field 3: All set to 0 for this dataset.
Field 4: Altitude in feet (-777 if not valid).
Field 5: Date – number of days (with fractional part) that have passed since 12/30/1899.
Field 6: Date as a string.
Field 7: Time as a string.

Data prep: PostGIS

Since any kind of GIS operation on text files will be quite inefficient, I decided to load the data into a PostGIS database. This table of millions of GPS points can then be sliced into appropriate chunks for exploration, for example, a day in Beijing:

    trajectories.t_datetime + interval '1 day' as t_to_datetime,
   FROM geolife.trajectories
   WHERE st_dwithin(trajectories.geom,
   AND trajectories.t_datetime >= '2008-11-11 00:00:00'
   AND trajectories.t_datetime < '2008-11-12 00:00:00'

Trajectory viz: a fadeout effect for point markers

The idea behind this visualization is to show both the current movement as well as the history of the trajectories. This can be achieved with a fadeout effect which leaves behind traces of past movement while the most recent positions are highlighted to stand out.

Map tiles by Stamen Design, under CC BY 3.0. Data by OpenStreetMap, under ODbL.

Map tiles by Stamen Design, under CC BY 3.0. Data by OpenStreetMap, under ODbL.

This effect can be created using a Single Symbol renderer with a marker symbol with two symbol layers: one layer serves as the highlights layer (pink) while the second layer represents the traces (green) which linger after the highlights disappear. Feature blending is used to achieve the desired effect for overlapping markers.

Screenshot 2015-05-06 23.52.40

The highlights layer has two expression-based properties: color and size. The color fades to white and the point size shrinks as the point ages. The age can be computed by comparing the point’s t_datetime timestamp to the Time Manager animation time $animation_datetime.

This expression creates the color fading effect:


and this expression makes the point size shrink:



I’m currently preparing this and a couple of other examples for my Time Manager workshop at the upcoming 1st QGIS conference in Nødebo. The workshop materials will be made available online afterwards.


[1] Yu Zheng, Lizhu Zhang, Xing Xie, Wei-Ying Ma. Mining interesting locations and travel sequences from GPS trajectories. In Proceedings of International conference on World Wild Web (WWW 2009), Madrid Spain. ACM Press: 791-800.
[2] Yu Zheng, Quannan Li, Yukun Chen, Xing Xie, Wei-Ying Ma. Understanding Mobility Based on GPS Data. In Proceedings of ACM conference on Ubiquitous Computing (UbiComp 2008), Seoul, Korea. ACM Press: 312-321.
[3] Yu Zheng, Xing Xie, Wei-Ying Ma, GeoLife: A Collaborative Social Networking Service among User, location and trajectory. Invited paper, in IEEE Data Engineering Bulletin. 33, 2, 2010, pp. 32-40.

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

Pedestrian routing on

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

Over the last couple of weeks, Karolina has been very busy improving and expanding Time Manager. This post is to announce the 1.6 release of Time Manager which brings you many fixes and exciting new features.

Screenshot 2015-03-25 17.58.38

What’s this feature interpolation you’re talking about?

Interpolation is really helpful if you have multiple observations of the same (moving) real-world object at different points in time and you want to visualize the movement between the observations. This can be used to visualize animal paths, vehicle tracks, or any other movement in space.

The following example shows a simple layer which contains 12 point features (3 for each id value).

Screenshot 2015-03-25 17.50.55

Using Time Manager interpolation, it is easy to create animations with interpolated positions between observations:


How is it done?

When you open the Time Manager 1.6 Settings | Add layer dialog, you will find a new option for interpolation settings. This first version supports linear interpolation of point features but more options might be added in the future. Note how the id attribute is specified to let Time Manager know which features belong to the same real-world object.

Screenshot 2015-03-25 17.43.08

For the interpolation, Time Manager creates a new layer which contains the interpolated features. You can see this layer in the layer list.

Screenshot 2015-03-25 17.46.13

I’m really looking forward to seeing all the great animations this feature will enable. Thanks Karolina for making this possible!

You probably remember my Game of Life posts from last year: Experiments with Conway’s Game of Life & More experiments with Game of Life where I developed a vector-based version of GoL.

Richard Wen and Claus Rinner at Ryerson University now published a raster-based version.

Here’s a screenshot of the script in action:

Screenshot 2015-03-08 20.04.07

The code is hosted on Github and I’m sure there will be many other ideas which can build on code snippets to read and write raster cell values.

For more info, please visit the GIS at Ryerson blog.

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:

Researching this topic, it turns out that some centrality measures have already been implemented in GRASS GIS. thumbs up! computes degree, betweeness, closeness and eigenvector centrality.

As a test, I’ve loaded the OSM street network of Vienna and run -a input=streets@anita_000 output=centrality degree=degree closeness=closeness betweenness=betweenness eigenvector=eigenvector


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, is not yet part of the QGIS Processing GRASS toolbox. It would certainly be a great addition.


Get every new post delivered to your Inbox.

Join 7,949 other followers

%d bloggers like this: