Writing one WKT file per feature in QGIS

1 Reply

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:
    output.write(wkt)

zstandard compression in GeoTIFF

2 Replies

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

DEFLATE

ZSTD

Big file

DEFLATE

ZSTD

Findings

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

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

2 Replies

I used the CityGML LoD2 model of Hamburg here, http://daten-hamburg.de/geographie_geologie_geobasisdaten/3d_stadtmodell_lod2/LoD2-DE_HH_2017-12-31.zip from http://suche.transparenz.hamburg.de/dataset/3d-stadtmodell-lod2-de-hamburg3. (If you do this for Hamburg, you might want to not include Neuwerk’s tiles ;))

1. https://github.com/tudelft3d/CityGML2OBJs/

python2 CityGML2OBJs.py -i LoD2-DE_HH_2017-12-31.zip/ -o LoD2-DE_HH_2017-12-31.zip_obj/

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

2. http://gfx.cs.princeton.edu/proj/trimesh2/

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

turns them into one single obj.

mesh_filter LoD2-DE_HH_2017-12-31.zip.obj -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.

3. https://google.github.io/draco/

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.

Mein Fahrrad!

Leave a reply

Vielleicht googelt ja jemand beim Gebrauchtfahrradkauf und meldet sich dann netterweise bei mir.
In der Nacht vom 27. auf den 28. Juli 2018 gestohlen:

Rabeneick
City-/Trekkingrad
Herrenrad (Diamantrahmen, nicht Trapez)
28 Zoll
52 cm Rahmengröße
Schwarz
8 Gang Shimano-Nabenschaltung
Felgenbremsen Vorne und Hinten
Rücktritt
Nabendynamo
Gepäckträger
Rahmennummer PBTDA 452596 266550 (oder schreibt man die anders? PBTDA452596266550 / PBTDA 452596266550 / 452596266550)
(Am linken Griff ein leicht hervorstehendes Endstück)

QGIS: “no result set” with PostGIS

Leave a reply

If QGIS tells you “no result set” while you are playing around with PostGIS query layers. If your query works fine in the DB Manager but “Load as Layer” fails (and not silently with “invalid PostgreSQL layer”).

Try turning it off and on again. Restart QGIS. It just might save you many minutes of unreasonable frustration.

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

Leave a reply

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. :)

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

Leave a reply

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. ;)

More structural detail in hillshading with this one weird trick

1 Reply

I wanted to have more structural detail in my QGIS live hillshading so I blended a duplicate layer with the sun angle set to directly overhead (90°) and the blending mode as darken.

It seems to make steep structures more prominent.

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

1 Reply

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.

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

Leave a reply

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.