Latest Posts

(Not) Solving #quiztime with a spatial relationship query in Overpass-Turbo

2018-04-21 13:04

I saw a #quiztime on Twitter and it looked like a perfect geospatial puzzle so I had to give it a try: https://twitter.com/Sector035/status/987631780147679233

So we have something called “Lidl” which I know is a supermarket. And we have something else, with a potential name ending in “nden”. And that’s it.

If you have a global OSM database, this would be a simple query away. If you don’t have a global OSM database, the Overpass API let’s you query one. And overpass turbo can make it fun (not in this case though, I had a lot of trial and error until I found these nice examples in the Wiki).

I ended up looking for things within 50 meters to a Lidl named “-nden”

[out:json][timeout:180][maxsize:2000000000];

{{radius=50}}

node[name="Lidl"]->.lidls;

( 
   way(around.lidls:{{radius}})["name"~".*nden$"];
  node(around.lidls:{{radius}})["name"~".*nden$"];
);

(._;>;);  // whatever? 

Too far away from the Lidl.

But driving around the place in Google StreetView I could not find the spot and it does not look very much like the photo.

So I guess my query is fine but either the "thing" or the Lidl are not in OSM yet.

Oh well, I did learn something new (about Overpass) and it was fun. :)

Leave your thoughts

Specifying the read/open driver/format with GDAL/OGR

2018-01-20 12:01

For the commandline utilities you can’t. One possible workaround is using https://trac.osgeo.org/gdal/wiki/ConfigOptions#GDAL_SKIP and https://trac.osgeo.org/gdal/wiki/ConfigOptions#OGR_SKIP to blacklist drivers until the one you want is its first choice. A feature request exists at https://trac.osgeo.org/gdal/ticket/5917 but it seems that the option is not exposed as commandline option (yet?).

PS: If what you are trying to do is reading a .txt (or whatever) file as .csv and you get angry at the CSV driver only being selected if the file has a .csv extension, use CSV:yourfilename.txt

PPS: This post was motivated by not finding above information when searching the web. Hopefully this will rank high enough for *me* to find it next time. ;)

Leave your thoughts

z/x/y.ext tiles to mbtiles/gpkg

2017-11-27 20:11

I just had some trouble converting a set of tiles I scraped into a nice package. I tried both mb-util and tiles2gpkg_parallel.py. Both seemed to work but QGIS did not display anything. The “solution” was to convert the tiles from 512×512 pixels to 256×256. WTF?

tiles2gpkg_parallel.py also crashed at joining its intermediate files if used in Python 3.6. It needs “-tileorigin ul” for OSM-like tile names.

1 Comment

Turn a Raspberry into a no-fancyshmancy mpd server with Archlinux ARM

2017-08-05 13:08

Because my SD card keeps getting corrupted and I always start from scratch, here are my notes on how to turn a Raspberry Pi 2 with an always connected USB disk full of audio media into a nice little mpd server outputting via the 3.5mm jack. Probably works somewhat the same on a Raspberry Pi 3.

Install Archlinux ARM

Install https://archlinuxarm.org/platforms/armv7/broadcom/raspberry-pi-2#installation
And obviously update it right after, add your ssh key, whatever. Then:

Add packages to build a special mpd, sound stuff and common utils:
# pacman -Syu screen wget zip base-devel libmikmod unzip zziplib git doxygen boost alsa-utils hdparm ffmpeg htop libao audiofile libshout libmad faad2 libupnp libmms wavpack avahi libid3tag yajl libmpdclient

Get sound working

Add to /boot/config.cfg:

dtparam=audio=on  # https://archlinuxarm.org/platforms/armv7/broadcom/raspberry-pi-2#wiki
disable_audio_dither=1  # if you get white noise on low volume, https://www.raspberrypi.org/documentation/configuration/config-txt/audio.md
audio_pwm_mode=2  # better audio driver, https://www.raspberrypi.org/forums/viewtopic.php?f=29&t=136445

Make sure audio output is enabled, not muted, in # alsamixer.
# speaker-test -c6 -twav

Connect and mount external disk with all your media

Connect USB disk. Check UUID using ls -l /dev/disk/by-uuid/ or lsblk -f, add it to /etc/fstab at a mount point of your choice:
UUID=12341234-1234-1234-1234-123412341234 /media/egon ext4 defaults,nofail,x-systemd.device-timeout=1 0 2
and made sure that it can spin down by adding /etc/udev/rules.d/50-hdparm.rules (you might want to verify the hdparm call works first):
ACTION=="add|change", KERNEL=="sd[a-z]", ATTR{queue/rotational}=="1", RUN+="/usr/bin/hdparm -B 1 -S 60 -M /dev/%k"
Then check if things work with # mount -a

Compile mpd with zip and curl support (optional, if not wanted, just install mpd from the repos)

I compiled my own mpd because I wanted zip support and Archlinux’s does not ship with that.
$ wget https://aur.archlinux.org/cgit/aur.git/snapshot/mpd-git.tar.gz
$ tar xfvz mpd-git.tar.gz
$ cd mpd-git
$ nano PKGBUILD

Add ‘armv7h’ to the archs

and add some fancy configure options:

      --enable-zzip \
      --enable-mikmod \
      --enable-modplug \
      --enable-curl 

$ makepkg
# pacman -U the resulting package

Takes about 30 minutes.

Configure mpd

/etc/mpd.conf:

user "mpd"
pid_file "/run/mpd/mpd.pid"
db_file "/var/lib/mpd/mpd.db"
state_file "/var/lib/mpd/mpdstate"
playlist_directory "/var/lib/mpd/playlists"
music_directory    "/media/egon"

audio_output {
  type            "alsa"
  name            "default"
  mixer_type      "software"      # optional
}

# systemctl enable mpd
# systemctl start mpd

WLAN

TODO ;)

Once it all works, make an image of it so that next time installation is just dd.

Leave your thoughts

A hurricane map in QGIS, from geodata and hacks

2017-08-01 22:08

Cartography Inspirationator John Nelson made an awesome map of hurricanes and later posted detailed how-tos for ArcMap and recently ArcGIS Pro. Right after his first post I started rebuilding it in QGIS using draw effects for adding the colored outer glow instead of using image icons, adding a vignette on-the-fly and adjusting the background raster’s saturation on the fly. All in all, less manual work, more dynamic processing in QGIS. I quickly got frustrated though and gave up.

More than a year later (triggered by the new ArcGIS Pro how-to storymap) I revisited my draft and finished it. So here is the QGIS version. As no-one is paying me to write this, I currently cannot be arsed to make it as fancy as John’s posts. Sorry! :)

This project shows QGIS’ strength in features and struggle with performance. Also some bugs. I stopped working on this once I liked the looks. It is not optimized in any way. So this is just how I ended up doing it. You could do much better. I like proof of concepts. And short sentences.

Data sources (I host a lightning-fast mirror of Natural Earth at https://www.datenatlas.de/geodata/public/sources/www.naturalearthdata.com/ if downloads are 404 again…):

That’s right, nothing but pure, unaltered geodata!

  1. Load the NE2_LR_LC_SR_W_DR.tif raster. Lower the saturation to about -50.
  2. Set your project’s projection to EPSG:3031 so that you get a nice polar viewpoint.
  3. Rotate the canvas by 150° because that’s what John did.
  4. Load the graticules. Use a rule-based style for ‘"direction" = 'W' or "direction" = 'E' or "direction" is NULL‘ so that you get all longitudes but only the equator from the latitudes. Set the layer transparency to 60 or something like that.
  5. Load the coastline.
  6. Add a vignette using this trick.
  7. Load the Allstorms.ibtracs_all_points.v03r09.shp Shapefile and then:
    1. Create a rule-based style with rules for the storm categories. I think I used this.
    2. Use white markers without outlines. Set their sizes using an expression on the wind speed or like me, manually to e.g. 0.7, 0.8, 1, 1, 1.6, 2 millimeters. This is something to play around with until it looks good.
    3. I used transparency for the markers of the lowest classes, 85% and 70%, the others are not transparent. This is something to play around with until it looks good. If I recall correctly I used transparency on the layer level here to keep the bubbly looks.
    4. By now you should be quite annoyed at how slow the rendering process is. >:) But wait, it gets much worse! :o)
    5. For each of the classes, enable draw effects on the markers.
    6. Set the source to use Addition blend mode. Set the source to be somewhat transparent, I used 50%, 80%, 80%, 50%, 40%, 0%. This is something to play around with until it looks good.
    7. Add outer glow and choose an appropriate color (Hint: Use a lot of saturation). Then play around with the spread, blur radius and transparency until it looks good. For some reason I ended up using 1mm/3/95%, 1mm/4/95%, 1mm/3/80%, 1.2mm/4/50%, 1.5mm/4/40%, 3mm/1/50%. Only now that I post this I realise how weirdly inconsistent this is and a quick test shows how irrelevant the blur radius changes are (except for the highest class). Oh well. It’s not fun to interate if you have to click so much and rendering that 40-50 seconds…
  8. That’s it! Done!

This was both fun and incredibly annoying. QGIS has the features but lacks in speed for this funky project (no wonder, blending 300,000 glowing points is not that nice). Here is a realtime video of how it rendered on my machine (take away some seconds from manually enabling the layers after another):

PS: Oh god this WordPress style sucks…

2 Comments

Flowers in QGIS?

2017-07-30 21:07

The other day I was working on visualisation of some intermediate research stuff and ended up with something looking like a bouquet of flowers.

Twitter liked it so here is a how-to.

  1. Have some lines that meet in a shared point. They should have a shared ID per group (in my case they were MultiLineStrings).
  2. Color them per ID.
  3. Turn the background black.
  4. Use an Arrow style for the lines. Set the Head thickness to half the Arrow width. Set the Head length to whatever you consider fancy.
  5. Remove the Outline (set No Pen).
  6. Set the Feature blending mode to Multiply (or to Screen or Dodge or Addition if you prefer fireworks to flowers).
  7. Set an appropriate color scheme. For flowers I think RdYlGn works great (that’s how I realised what my random tinkering had lead to) or PiYg or simply Spectral, for fireworks random colors.

That’s it! Now play around and have some fun!

Leave your thoughts

Initial setup failed. Cannot continue. Error: Couldn’t run mojosetup

2017-05-06 13:05

If you are trying to install a game downloaded from gog.com but you get something like

$ ./gog_shadow_tactics_blades_of_the_shogun_2.2.0.3.sh
Verifying archive integrity… All good.
Uncompressing Shadow Tactics: Blades of the Shogun (GOG.com) 100%
Collecting info for this system…
Operating system: linux
CPU Arch: x86_64
trying mojosetup in bin/linux/x86_64
USING en_US

PANIC
Initial setup failed. Cannot continue.

Error: Couldn’t run mojosetup

then your file is corrupted or incomplete. You will have to re-download it. Ignore the “bin/linux/x86_64” path, that references something inside the installer, not something on your system.

Petition GOG to add checksums to their download pages or better yet, have reasonable download options that support resuming and HTTPS and what not…

Leave your thoughts

FOSSGIS 2017

2017-03-28 18:03

Die FOSSGIS 2017 in Passau war grandios. Ich bin sooo froh, dass ich mich auf den weiten Weg gemacht hatte. Die Liste von Vortragsaufzeichnungen, die ich selbst noch anschauen will, ist lang… Ausprobieren muss ich unbedingt (mal wieder) ein aktuelles GRASS GIS, GVSIG CE (das Poster hat Lust gemacht), osmium, die ganzen Vector Tiles Tools uvm.

Selbst gewagt habe ich einen Lightning Talk über Interaktive Visualisierung von Geodaten in Jupyter Notebooks (Youtube) sowie einen Vortrag zu meinem Projekt GeoPackages der freien Hamburger Geodaten (Youtube) anzufertigen.

Der LT kam so extrem gut an, dass ich nächstes Mal wohl um einen richtigen Vortrag oder auch Workshop kaum herum komme. :o)

Für den Geopackage-Vortrag hatte ich leider die Daten und Skripte zuhause gelassen und musste daher etwas improvisieren… Trotzdem war kam auch er gut an und ich habe großartigen Input bekommen, z.B. dass es ein tolles neues QGIS-Plugin für GML Application Schema Gedöns gibt und einen GMLAS-Treiber in GDAL. Danke!

Bei spontanen QGIS-Anwender- und Vereinstreff habe ich nachgefragt, wie es eigentlich mit QGIS an den Hochschulen aussieht und wie ich meinen Arbeitgeber vielleicht mal auf den Weg von Esri/IDRISI zu QGIS bringen kann. Da war ausser Claas Leiners Lehre in Kassel wenig bekannt. Vielleicht starte ich mal eine kleine Recherche, um etwas Einblick in die Landschaft zu bekommen. Wäre doch klasse, wenn sich mehr Unis von proprietärer Software entsagen mögen!

Leave your thoughts

What does it look like if you move all countries onto the same location?

2017-01-22 12:01

Sunday pre-lunch Python fun: What does it look like if you move all countries onto the same location?

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
import fiona
from shapely.geometry import *
from fiona.crs import from_string
from fiona.transform import transform_geom
from shapely.affinity import translate
 
def get_biggest_polygon(multipolygon):
    assert isinstance(geometry, MultiPolygon)
 
    max_area = 0
    biggest_polygon = None
    for polygon in geometry:
        if polygon.area > max_area:
            max_area = polygon.area
            biggest_polygon = polygon
    return biggest_polygon
 
def project_locally(geometry, from_crs):
    """centered on the centroid of the geometry"""
    # ugly because i map/shape back and forth, maybe try shapely instead
    lat = geometry.centroid.y
    lon = geometry.centroid.x
 
    to_crs = from_string("+proj=aeqd  +R=6371000 +lat_0={lat} +lon_0={lon}".format(lat=lat, lon=lon))
    reprojected_geometry = transform_geom(
        from_crs, 
        to_crs, 
        mapping(geometry)
    )
 
    return shape(reprojected_geometry)
 
with fiona.open("ne_10m_admin_0_countries.shp") as countries:
    centered_polygons = []
 
    for country in countries:
        geometry = shape(country['geometry'])
 
        # only use the biggest part of each country, otherwise everything sucks
        if isinstance(geometry, MultiPolygon):
            polygon = get_biggest_polygon(geometry)
        else:
            polygon = geometry
 
        # project nicely
        polygon = project_locally(polygon, countries.crs)
 
        # centering on 0,0 is simply moving the geometry by MINUS its x/y
        dx = -polygon.centroid.x
        dy = -polygon.centroid.y
        translated_polygon = translate(polygon, dx, dy)
 
        centered_polygons.append(translated_polygon)
 
with open("/tmp/outfile.wkt", "w") as sink:
    for polygon in centered_polygons:
        sink.write(polygon.wkt+"\n")

Use this in any way you like but please share your creations and code as well. :)

Some rough explanation: For each country I check if it is a multipolygon and if so, use only its biggest “sub”-polygon in the next steps. I then project the WGS84 coordinates to an Azimuthal Equidistant projection centered on the centroid of the polygon. That new geometry gets shifted to sit on the origin of the system. I collect all those polygons and write them as plain WKT to a file. Styling was done in QGIS.

And for the smart folk, the same but without local projection:

2 Comments