Latest Posts

Replicating a media-hyped color by numbers Etsy map in 10 minutes

2019-01-24 17:01

Thats beautiful… how long it took?

Well, that looks like QGIS’ random colors applied to

So I fired up QGIS, extracted the region from, realised those rivers came without a corresponding basin, extracted the region from, set the map background to black, set the rivers to render in white, set the rivers’ line width to correspond to their UP_CELLS attribute, put the basins on top, colored them randomly by BASIN_ID, set the layer rendering mode to Darken and that was it.

I should open an Etsy store.

Yes, I realise that replicating things is easier than creating them. But seriously, this is just a map of features colored by category and all the credit should go to


But Hannes, that original has some gradients!

Ok, then set the rivers not to white but a grey and the basin layer rendering mode to Overlay instead of Darken.

This product incorporates data from the HydroSHEDS database which is © World Wildlife Fund, Inc. (2006-2013) and has been used herein under license. WWF has not evaluated the data as altered and incorporated within, and therefore gives no warranty regarding its accuracy, completeness, currency or suitability for any particular purpose. Portions of the HydroSHEDS database incorporate data which are the intellectual property rights of © USGS (2006-2008), NASA (2000-2005), ESRI (1992-1998), CIAT (2004-2006), UNEP-WCMC (1993), WWF (2004), Commonwealth of Australia (2007), and Her Royal Majesty and the British Crown and are used under license. The HydroSHEDS database and more information are available at

Leave your thoughts

Average Earth from Space 2018

2019-01-21 13:01

Due to popular request, I now provide a high-resolution, digital download as an experiment. Any profit from this will go towards bug-fixes in GDAL, the free and open-source software I used for this project.

I collected all the daily images of the satellite-based Soumi VIIRS sensor and calculated the per-pixel average over the whole year of 2018.

You can explore it in full resolution in an interactive webmap.

WGS84 / Equirectangular
Web Mercator

You can see fascinating patterns, e.g. “downstream” of islands in the ocean or following big rivers like the Amazon. Be wary though as snow and clouds look the same here.

It’s also fun to look at the maximum to see cloudless regions (this image is highly exaggerated and does not make too much sense anyways):

Many thanks to Joshua Stevens for motivation and infos on data availability! The initial idea was of course inspired by Charlie ‘vruba’ Lyod‘s Cloudless atlas works. I explored imagery on, grabbed Soumi VIIRS images via and processed it in GDAL. Poke me repeatedly to blog about the process please.

Full resolution, georeferenced imagery (12288×6144 resp. 8192×8192) in PNG available on request.

I said average but of course it is the median. The average or mean does not make much sense to use nor does it look any good. ;-)


Data-defined images in QGIS Atlas

2019-01-09 18:01

Say you want to display a feature specific image on each page of a QGIS Atlas.

In my example I have a layer with two features:

"type": "FeatureCollection",
"features": [
"type": "Feature",
"properties": {
"image_path": "/tmp/1.jpeg"
"geometry": {
"type": "Polygon",
"coordinates": [[[9,53],[11,53],[11,54],[9,54],[9,53]]]
"type": "Feature",
"properties": {
"image_path": "/tmp/2.jpeg"
"geometry": {
"type": "Polygon",
"coordinates": [[[13,52],[14,52],[14,53],[13,53],[13,52]]]

And I also have two JPEG images, named “1.jpeg” and “2.jpeg” are in my /tmp/ directory, just as the “image_path” attribute values suggest.

The goal is to have a map for each feature and its image displayed on the same page.

Create a new print layout, enable Atlas, add a map (controlled by Atlas, using the layer) and also an image.

For the “image source” enable the data-defined override and use attribute(@atlas_feature, 'image_path') as expression.

That’s it, now QGIS will try to load the image referenced in the feature’s “image_path” value as source for the image on the Atlas page. Yay kittens!

Leave your thoughts


2018-12-21 15:12

Fergus Falls liegt im Westen Minnesotas, zwischen den Bundesstaaten Wisconsin und North Dakota, am nördlichen Rand der USA. Die Jahresdurchschnittstemperatur nahe der Grenze zu Kanada liegt bei plus 3 Grad Celsius, im Winter bei minus 20 Grad. Von den Hochhäusern New Yorks und den Stränden San Franciscos sind es mehr als 2200 Kilometer bis nach Fergus Falls. Von El Paso, der mexikanischen Grenze, braucht ein Auto 22 Stunden. Wie viele Mexikaner zieht es in diese Gegend? Wer stellt hier so ein Schild auf?

Der SPIEGEL hat in einem ersten Schritt seine Dokumentationsabteilung das Manuskript der Relotius-Geschichte noch einmal in Stichproben prüfen lassen – und musste feststellen, dass bei der Verifikation tatsächlich nicht sauber gearbeitet wurde. (Mehr zu den Sicherheitsmechanismen des Hauses lesen Sie hier.)

Ein paar Detail-Beispiele, die auffielen:

So sind es von Fergus Falls nicht 2200 Kilometer nach New York, wie es im Text steht, sondern nur 1888.

Ich *hasse* diesen Zahlenfetischismus (passt genau in diese gefühlsmanipulative Art von Reportagen), daher war mein Interesse geweckt (während im Hintergrund gerade “for king and country” aus den Lautsprechern schallte.

Warum sollten es exakt 1888 Kilometer sein? Was hatte die Spiegel-Dokumentation denn da gemessen? Rathaus zu Rathaus? Zentroid zu Zentroid der administrativen Ortsgrenzen? Die kürzeste Distanz?

Da im nächsten Satz des Artikels eine Reisedauer mit dem Auto angegeben wurde, scheint es naheliegend, dass der Autor auch für die Distanzen eine solche Reise angenommen hatte.

Und siehe da, Google Maps gibt für die Reise mit dem Auto von Fergus Falls (Minnesota 56537) nach New York 2218 Kilometer, alternativ 2335 Kilometer an.

Graphhopper bestätigt die 2218 Kilometer (exakt!!!11).

Bing bietet Routen von 2208 bis 2372 Kilometer an.

Yandex rundet sinnvoll auf 2300 bzw. 2400 Kilometer

Vier verschiedene Routenplaner bestätigen also eine Entfernung von grob 2200 Kilometern (oder mehr) und damit die Angabe im Artikel des Autors.

Was hat die Spiegel-Dokumentation also dann wohl gemacht?

Auf Wikipedia gibt es repräsentative Koordinaten für beide Orte, 46°16′59″N 96°04′39″W für Fergus Falls (irgendwo), 40°42′46″N 74°00′21″W für New York City (Rathaus).

Betrachtet man die Erde als Spheroid läge die direkte Distanz zwischen diesen bei etwa 1878 Kilometern.

Mit einem besser passenden Ellipsoid (WGS84) wären es 1882 Kilometer.

Die Spiegel-Dokumentation hat also wohl irgendwelche Koordinaten, relativ zentral in den beiden Orten, gewählt und die direkte Distanz zwischen ihnen berechnet. Immerhin geodätisch!

Fazit: Ach, watt weiß ich, dieses Detail im Faktencheck der Spiegel-Dokumentation ist definitiv verständnisloser Aktionismus. Distanzen sind ein kompliziertes Konzept. Ist eine Reisedistanz gemeint? Diese kann durch unterschiedliche Routen oder Umleitungen stark abweichen. Ist eine direkte Distanz gemeint? Falls ja, wie genau kann das Ergebnis jetzt sein? Was sind die exakten Punkte, zwischen denen gemessen wurde und welche Art von Genauigkeit bzw. Unsicherheit haben diese? Wann ist es sinnvoll eine Stadt mit einer Koordinate zu beschreiben?

1888 Kilometer ist in jedem Fall nichts als ein blindes Vortäuschen von Exaktheit.

PS: Zahlen vergöttern die, die sie nicht verstehen. Mir schwindelte.

Leave your thoughts

Fake chromatic aberration and dynamic label shadows in QGIS

2018-12-02 13:12

Welcome to Part 43 of “Fun with Weird Hacks in QGIS”!

Fake Chromatic Aberration

Get some lines or polygons! I used German administrative polygons and a “Outline: Simple line” style.

Add a “Geometry Generator” to the same layer and set it to “LineString / MultiLineString”. If you used polygons, you will need to use “boundary($geometry)” to get the border lines of your polygons.

Translate (shift) the geometries radially away from the map center, with an amount of translation based on their distance to the map center. This will introduce ugly artifacts cool glitches for big geometries. Note the magic constant of 150 I divided the distances by, I cannot be arsed to turn this into percentages so you will have to figure out what works for you. Also, if your map is in a different CRS than the layer you do this with, you will need to transform the coordinates (I do that for the labels below).


Color those lines in some fancy 80s color like #00FFFF.

We only want the color to appear in the edges of the map, so set the “Feature Blending” mode of the layer to “Lighten”. This will make sure the white lines do not get darker/colored.

I forgot to take an image for that and then it was lunch time. :x

Now do the same but for distances the other way around and color that in something else (like foofy #FF00FF).


Oh, can you see it already? Move the map around! Ohhh!

For the final touch, use the layer’s “Draw Effects” to replace the “Source” with a “Blur”. Be aware that the “Blur type” can quite strongly influence the look and find a setting for the “Blur strength” that works for you. I used “Guassian blur (quality)” with a strength of 2.

Dynamic Label Shadows

Get some data to label! I used ne_10m_populated_places_simple. Label it with labels placed “Offset from Point” without any actual offset. This is just to make sure calculations on the geometry’s location make sense to affect the labeling later.

Pick a wicked cool font like Lazer 84!

Add a “Buffer” to the labels and pick an appropriate color (BIG BOLD #FF00FF works well again).

Time for magic! Add a “Shadow” to the labels and use an appropriate color (I used the other color from earlier again, #00FFFF).

We want to make the label shadows be further away from the label if the feature is further away from the map center. So OVERRIDE the offset with code that does exactly that. Note another magic constant (relating to meters in EPSG:25832) and that I needed to transform coordinates here (my map is in EPSG:25832 while this layer is EPSG:4326).

		$geometry, 'EPSG:4326', 'EPSG:25832'

Cool. But that just looks weird. Time for another magical ingredient while OVERRIDING the angle at which the label is placed. You guessed it, radially away from the map center!

			$geometry, 'EPSG:4326', 'EPSG:25832'

By the way, it does look best if your text encoding is introducing bad character marks and missing umlauts!

That’s it, move the map around and be WOWed by the sweet effects and the time it takes to render :o)

And now it is your turn to apply this to some MURICA geodata and rake in meaningless internet points that make you feel good! Sad but true ;)

And also TODO: Make the constants percentages instead, that should make sure it works on any projection and any scale.

Leave your thoughts

Making a Star Wars hologram in QGIS

2018-11-11 21:11

I like doing silly things in QGIS.

So I wanted to make a Star Wars hologram (you know, that “You’re my only hope” Leia one) showing real geodata. What better excuse for abusing QGIS’ Inverted Polygons and Raster Fills. So, here is what I did:

  1. Find some star wars hologram leia image.
  2. Crudly remove the princess (GIMP’s Clone and Healing tools work nicely for this).
  3. In QGIS create an empty Polygon-geometry Scratchpad layer and set the renderer to Inverted Polygons to fill the whole canvas.
  4. Set the Fill to a Raster image Fill and load your image.
  5. Load some geodata, style it accordingly and rejoice.
    1. Get some geodata. I used Natural Earth’s countries, populated places and tiny countries (to have some stuff in the oceans), all in 110m.
    2. Select a nice projection, I used “+proj=ortho +lat_0=39 +lon_0=139 +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m +no_defs”.
    3. I used a three layer style for the countries:
      • A Simple Line outline with color #4490f3, stroke width 0.3mm and the Dash Dot Line stroke style.
      • A Line Pattern Fill with a spacing of 0.8mm, color #46a8f3 and a stroke width of 0.4mm.
      • And on top of those, for some noisiness a black Line Pattern Fill rotated 45° with a spacing of 1mm and stroke width 0.1mm.
      • Then the Feature Blending Mode Dodge to the Layer Rendering and aha!
      • More special effects come from Draw Effects, I disabled Source and instead used a Blur (Gaussian, strength 2) to lose the crispness and also an Outer Glow (color #5da6ff, spread 3mm, blur radius 3) to, well, make it glow.

    4. I used a three layer style for the populated places:
      • I used a Simple Marker using the “cross” symbol and a size of 1.8mm. The Dash Line stroke style gives a nice depth effect when in Draw Effects the Source is replaced with a Drop Shadow (1 Pixel, Blur radius 2, color #d6daff).
      • A Blending Mode of Addition for the layer itself makes it blend nicely into the globe.
    5. I used a three layer style for the tiny countries:
      • I used a white Simple Marker with a size of 1.2mm and a stroke color #79c7ff, stroke width 0.4mm.
      • Feature Blending Mode Lighten makes sure that touching symbols blob nicely into each other.

  6. You can now export the image at your screen resolution (I guess) using Project -> Import/Export or just make a screenshot.
  7. Or add some more magic with random offsets and stroke widths in combination with refreshing the layers automatically at different intervals:
Leave your thoughts

Building QGIS with debugging symbols

2018-09-29 14:09

As I keep searching the web for way too long again and again, I hope this post will be #1 next I forget how to build QGIS with debugging symbols.

Add CMAKE_BUILD_TYPE=Debug to the cmake invocation.


cmake -G "Unix Makefiles" ../ \

For a not as safe but more performant compilation, you can use RelWithDebInfo. I just found out today but will use that in the future rather than the full-blown Debug. See for some background.

On Archlinux, also add options=(debug !strip) in your PKGBUILD to have them not stripped away later.

Leave your thoughts

Writing one WKT file per feature in QGIS

2018-09-18 13:09

Someone in #qgis just asked about this so here is a minimal pyqgis (for QGIS 3) example how you can write a separate WKT file for each feature of the currently selected layer. Careful with too many features, filesystems do not like ten thousands of files in the same directory. I am writing them to /tmp/ with $fid.wkt as filename, adjust the path to your liking.

layer = iface.activeLayer()
features = layer.getFeatures()

for feature in features:
  geometry = feature.geometry()
  wkt = geometry.asWkt()
  fid = feature.attribute("fid")
  filename = "/tmp/{n}.wkt".format(n=fid)
  with open(filename, "w") as output:
1 Comment

zstandard compression in GeoTIFF

2018-08-23 18:08

I ran GDAL 2.4’s gdal_translate (GDAL 2.4.0dev-333b907 or GDAL 2.4.0dev-b19fd35e6f-dirty, I am not sure) on some GeoTIFFs to compare the new ZSTD compression support to DEFLATE in file sizes and time taken.

Hardware was a mostly idle Intel(R) Xeon(R) CPU E3-1245 V2 @ 3.40GHz with fairly old ST33000650NS (Seagate Constellation) harddisks and lots of RAM.

A small input file was DGM1_2x2KM_XYZ_HH_2016-01-04.tif with about 40,000 x 40,000 pixels at around 700 Megabytes.
A big input file was srtmgl1.003.tif with about 1,3000,000 x 400,000 pixels at 87 Gigabytes.
Both input files had been DEFLATE compressed at the default level 6 without using a predictor (that’s what the default DEFLATE level will make them smaller here).

gdal_translate -co NUM_THREADS=ALL_CPUS -co PREDICTOR=2 -co TILED=YES -co BIGTIFF=YES --config GDAL_CACHEMAX 6144 was used all the time.
For DEFLATE -co COMPRESS=DEFLATE -co ZLEVEL=${level} was used, for ZSTD -co COMPRESS=ZSTD -co ZSTD_LEVEL=${level}

Mind the axes, sometimes I used a logarithmic scale!

Small file



Big file




It has been some weeks since I really looked at the numbers, so I am making the following up spontaneously. Please correct me!
Those numbers in the findings below should be percentages (between the algorithms, to their default values, etc), but my curiosity was satisfied. At the bottom is the data, maybe you can take it to present a nicer evaluation? ;)

ZSTD is powerful and weird. Sometimes subsequent levels might lead to the same result, sometimes a higher level will be fast or bigger. On low levels it is just as fast as DEFLATE or faster with similar or smaller sizes.

A <700 Megabyte version of the small file was accomplished within a minute with DEFLATE (6) or half a minute with ZSTD (5). With ZSTD (17) it got down to <600 Megabyte in ~5 Minutes, while DEFLATE never got anywhere near that.
Similarly for the big file, ZSTD (17) takes it down to 60 Gigabytes but it took almost 14 hours. DEFLATE capped at 65 Gigabytes. The sweet spot for ZSTD was at 10 with 4 hours for 65 Gigabytes (DEFLATE took 11 hours for that).

In the end, it is hard to say what default level ZSTD should take. For the small file level 5 was amazing, being even smaller than and almost twice as fast as the default (9). But for the big file the gains are much more gradual, here level 3 or level 10 stand out. I/O might be to blame?

Yes, the machine was not stressed and I did reproduce those weird ones.

Raw numbers

Small file

Algorithm Level Time [s] Size [Bytes] Size [MB] Comment
ZSTD 1 18 825420196 787
ZSTD 2 19 783437560 747
ZSTD 3 21 769517199 734
ZSTD 4 25 768127094 733
ZSTD 5 31 714610868 682
ZSTD 6 34 720153450 687
ZSTD 7 40 729787784 696
ZSTD 8 42 729787784 696
ZSTD 9 51 719396825 686 default
ZSTD 10 63 719394955 686
ZSTD 11 80 719383624 686
ZSTD 12 84 712429763 679
ZSTD 13 133 708790567 676
ZSTD 14 158 707088444 674
ZSTD 15 265 706788234 674
ZSTD 16 199 632481860 603
ZSTD 17 287 621778612 593
ZSTD 18 362 614424373 586
ZSTD 19 549 617071281 588
ZSTD 20 834 617071281 588
ZSTD 21 1422 616979884 588
DEFLATE 1 25 852656871 813
DEFLATE 2 26 829210959 791
DEFLATE 3 32 784069125 748
DEFLATE 4 31 758474345 723
DEFLATE 5 39 752578464 718
DEFLATE 6 62 719159371 686 default
DEFLATE 7 87 710755144 678
DEFLATE 8 200 705440096 673
DEFLATE 9 262 703038321 670

Big file

Algorithm Level Time [m] Size [Bytes] Size [MB] Comment
ZSTD 1 70 76132312441 72605
ZSTD 2 58 75351154492 71860
ZSTD 3 63 73369706659 69971
ZSTD 4 75 73343346296 69946
ZSTD 5 73 72032185603 68695
ZSTD 6 91 72564406429 69203
ZSTD 7 100 71138034760 67843
ZSTD 8 142 71175109524 67878
ZSTD 9 175 71175109524 67878 default
ZSTD 10 235 69999288435 66757
ZSTD 11 406 69999282203 66757
ZSTD 12 410 69123601926 65921
ZSTD 13 484 69123601926 65921
ZSTD 14 502 68477183815 65305
ZSTD 15 557 67494752082 64368
ZSTD 16 700 67494752082 64368
ZSTD 17 820 64255634015 61279
ZSTD 18 869 63595433364 60649
ZSTD 19 1224 63210562485 60282
ZSTD 20 2996 63140602703 60216
ZSTD 21 lolno
DEFLATE 1 73 87035905568 83004
DEFLATE 2 76 85131650648 81188
DEFLATE 3 73 79499430225 75817
DEFLATE 4 77 75413492394 71920
DEFLATE 5 92 76248511117 72716
DEFLATE 6 129 73901542836 70478 default
DEFLATE 7 166 73120114047 69733
DEFLATE 8 407 70446588490 67183
DEFLATE 9 643 70012677124 66769
Leave your thoughts

Lots of CityGML files to a single obj to draco with localising coordinates to 0,0,0

2018-08-22 21:08

I used the CityGML LoD2 model of Hamburg here, from (If you do this for Hamburg, you might want to not include Neuwerk’s tiles ;))


python2 -i -o LoD2-DE_HH_2017-12-31.zip_obj/

converts each CityGML file to a obj (takes about half an hour on my system).


mesh_cat LoD2-DE_HH_2017-12-31.zip_obj/*.obj -o

turns them into one single obj.

mesh_filter -center -scale 0.001 -rot 90 -1 0 0 LoD2-DE_HH_2017-12-31.zip_localised.obj

translates “so center of mass is at (0,0,0)” and scales to something smaller and rotates so that z is up. trimesh2’s documentation on the rot parameter is insufficient, I think my solution means “rotate by 90 degrees -1 times around x, 0 times around y, 0 times around z” or something like that. Look into the source if that helps you.


draco_encoder -i LoD2-DE_HH_2017-12-31.zip_localised.obj -o LoD2-DE_HH_2017-12-31.zip_localised.obj.drc

makes it tiny (8 Megabyte).

I did not do any error checking, used a epsilon of 1 in CityGML2OBJs and just wanted to play around.

Leave your thoughts