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.

A hurricane map in QGIS, from geodata and hacks

2 Replies

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…