Categories for

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


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

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 and to blacklist drivers until the one you want is its first choice. A feature request exists at 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

One SRTMGL1 GeoTIFF to rule them all

2017-01-10 12:01

So about half a year ago Lukas Martinelli asked about a global SRTM GeoTIFF. The SRTM elevation data is usually shared in many small tiles which can be ideal for some cases but annoying for others. I like downloading and processing big files so I took that as a challenge. It’s probably some mental thing. I never finished this blog post back then. It’s probably another mental thing. 8) Read on if you are interested in some GNU coreutils hackery as well as GDAL magick.

Here is how I did it. Endless thanks as usual to Even Rouault who fixed GDAL bugs and gave great hints. I learned more about GeoTIFFs and GDAL than I should have needed to know.

Downloading the source files

The source files are available at (warning, visiting this URL will make your browser cry and potentially render your system unresponsive). To download them you must create an NASA Earthdata Login at :\


                         U.S. GOVERNMENT COMPUTER

This US Government computer is for authorized users only.  By accessing this
system you are consenting to complete monitoring with no expectation of privacy.
Unauthorized access or use may subject you to disciplinary action and criminal

OMGOMGOMG. I am sure they would have used <blink> if accessibility guidelines allowed it. As I am not sure what their Terms of Service are, I will not give you a copy’n’paste ready line to download them all. wget or aria2 work well. You should get 14297 files with a total size of 98G.

Inspecting the ZIP files and preparing for GDAL’s vsizip

Each of those ZIP files has just a single file “.hgt” inside. GDAL specifically has support for the “SRTM HGT Format“, so we know it can read those files. We just need to extract them. (If you clicked that link you see that with GDAL 2.2 it will support directly loading the data from the ZIPs, that’s just Even being awesome.)

We don’t want to uncompress all those files just because we want to build a GeoTIFF from them later, do we? Luckily GDAL has a vsizip thingie which allows it to read files inside zip archives. Wicked! For this we need the “deep” paths to the files though, for example

The files inside the archives are always simply the same filename minus the “.SRTMGL1” and the .zip extension. Perfect, now we know all the files from inside the ZIPs! Right? Nope. Unfortunately some (17) of the archives do include the “.SRTMGL1” bit in their files, for example N38E051.SRTMGL1.hgt… >:(

That’s why we need to look into each zip file and determine the filename inside. (Alternatively *you* could extract them all after all, you with your fancy, huge SSD.)

We can simply use `unzip` with the `-t` switch to make it show what’s inside (and as a “free” benefit it will also check the file’s integrity for us). `tee` is used here to save the output to a file. This will take a while as it needs to read through all the files…

$ unzip -t "*" | tee unzip-t.log


testing: N29E000.hgt OK
No errors detected in compressed data of

testing: N27E033.hgt OK
No errors detected in compressed data of

testing: N45E066.hgt OK
No errors detected in compressed data of

14297 archives were successfully processed.
14297 archives were successfully processed.

Yay, no errors!

In the output we see the path to the archive and we see the name of the file inside. With some sed and grep we can easily construct “deep” paths.

First we remove the status lines, the blank lines and the very last full status line with `grep`:

$ cat unzip-t.log | grep -v -e "No errors" -e '^$' -e 'successfully processed'

testing: N12E044.hgt OK
testing: N29W002.hgt OK
testing: N45E116.hgt OK

Then we concatenate every consecutive two lines into one line with `paste` (I LOVE THIS 1 TRICK!):

$ cat unzip-t.log | grep -v -e "No errors" -e '^$' -e 'successfully processed' | paste - -

Archive: testing: N12E044.hgt OK
Archive: testing: N29W002.hgt OK
Archive: testing: N45E116.hgt OK

Finally we strip things away with sed (I don’t care that you would do it differently, this is how I quickly did it with my flair of hammering) and direct the output into a new file called `zips`:

$ cat unzip-t.log | grep -v -e "No errors" -e '^$' -e 'successfully processed' | paste - - | sed -e 's#Archive:[ ]*##' -e 's#\t[ ]*testing: #/#' -e 's#[ ]*OK##' > zips

Yay, we have all the actual inside-zip paths now!

If you are curious, you can try gdalinfo with those now. You need to prepend /vsizip/ to the path to make it read inside zip files.

$ gdalinfo /vsizip/

Driver: SRTMHGT/SRTMHGT File Format
Files: /vsizip/
Size is 3601, 3601
Coordinate System is:
        SPHEROID["WGS 84",6378137,298.257223563,
Origin = (43.999861111111109,13.000138888888889)
Pixel Size = (0.000277777777778,-0.000277777777778)
Corner Coordinates:
Upper Left  (  43.9998611,  13.0001389) ( 43d59'59.50"E, 13d 0' 0.50"N)
Lower Left  (  43.9998611,  11.9998611) ( 43d59'59.50"E, 11d59'59.50"N)
Upper Right (  45.0001389,  13.0001389) ( 45d 0' 0.50"E, 13d 0' 0.50"N)
Lower Right (  45.0001389,  11.9998611) ( 45d 0' 0.50"E, 11d59'59.50"N)
Center      (  44.5000000,  12.5000000) ( 44d30' 0.00"E, 12d30' 0.00"N)
Band 1 Block=3601x1 Type=Int16, ColorInterp=Undefined
  NoData Value=-32768
  Unit Type: m

Hooray! GDAL reads the hgt file from inside its zip!

We need to prepend all the paths with `/vsizip/` so let’s do that:

$ sed 's#^#/vsizip/##' zips > vsizips


Unfortunately those 17 misnamed files from earlier need special treatment… This is what my notes say, not sure what it was supposed to do. If you are recreating this all, please just ask me and I will look at it again. For now, let’s just pretend that this leads to a file called `hgtfiles` in which all the paths are perfect and all the files are perfect.

grep -v 'SRTMGL1.hgt$' vsizips > vsizipswithoutmisnamedfiles
mkdir misnamedfiles
cd misnamedfiles/
nano ../listofmisnamedfiles # insert the paths to those 17 zips here #TODO
while read filename; do unzip "../${filename}"; done < ../listofmisnamedfiles
rename 's/SRTMGL1.//' *.hgt
cd ..
find misnamedfiles/ -type f > listofmisnamedfileshgt
cat listofmisnamedfileshgt vsizipswithoutmisnamedfiles > hgtfiles

As I said above though, go with a recent GDAL and this is all much easier. Even even included a check for the different filenames inside, how can you not like that guy!

Building a Virtual Raster Table

Virtual Raster Tables (VRT) are some kind of files that pretend to be rasters. They are awesome. Here we use a VRT that simply turns all the small rasters we have into one big ass raster.

Ok, ready to create a Virtual Raster Table!

$ time gdalbuildvrt -input_file_list hgtfiles srtmgl1.003.vrt

0...10...20...30...40...50...60...70...80...90...100 - done.

real 0m22.679s
user 0m19.250s
sys 0m3.105s


You could go ahead and use this for your work/leisure now. But remember, it is tens of gigabytes of data so if you do not use it at a 1:1 scale things will not be fun and might fry your cat. You want overviews/pyramids.

Turning the VRT into a GeoTIFF (Optional)

Let’s make a HUGEGEOTIFF because that’s cool! You don’t have to, instead you could build overviews for the .vrt file.

We want it quick to read and small so I used DEFLATE, TILED and the horizontal predictor. I ran this on a weak i7 with 2G of RAM and can’t remember what the worst bottleneck was. Probably CPU.

$ time gdal_translate -co NUM_THREADS=ALL_CPUS --config PREDICTOR 2 -co COMPRESS=DEFLATE -co TILED=YES -co BIGTIFF=YES srtmgl1.003.vrt srtmgl1.003.vrt.tif

Input file size is 1296001, 417601
0...10...20...30...40...50...60...70...80...90...100 - done.

real 231m57.961s
user 503m23.044s
sys 3m14.600s

If you are courageous you can load that file in your GIS now. But again, unless you watch it at a 1:1 scale or something close to that it WILL not be much fun and potentially expose weaknesses in your GIS and fill up your system’s memory and crash and you had no unsaved work open, didn’t you?

Building Overviews

Overviews aka pyramids are usually about 1/3 the size of the full image. If they are not, you probably used different compression settings. This was the step that I expected to be just boring to wait for, but it turned out the most problematic. I tried all available free tools but none worked properly with an input this big. With GDAL we found a workaround after a while.

gdaladdo has problems building multiple overview levels with a file this big… The trick is to built the overviews sequentially, not in one invocation. The best solution at the moment was building overviews for overviews for overviews and so on until you reach a comfortable size. Something like:

gdaladdo -ro srtmgl1.003.tif 2
gdaladdo -ro srtmgl1.003.tif.ovr 2
gdaladdo -ro srtmgl1.003.tif.ovr.ovr 2
gdaladdo -ro srtmgl1.003.tif.ovr.ovr.ovr 2
gdaladdo -ro srtmgl1.003.tif.ovr.ovr.ovr.ovr 2 4 8 16 32 64 128 256 512

I used the following parameters:


Apparently gdaladdo automatically makes them tiled, which is good. I used to find out.

Creating overviews took just about 3 hours with this weird trick. So, in total the processing just takes half a day.


Enjoy! I impose no license/copyright/whatever on this derivative work.


Hillshading? Slope? You do it!

Found a better, faster way? You blog it!

We could add the overviews to the main image with `gdal_translate srtmgl1.003.tif srtmgl1.003.withovr.tif -co COPY_SRC_OVERVIEWS=YES [other options]` says Even.

Leave your thoughts

Fancy Free File Formats For Hamburg’s Open Geodata

2016-10-02 16:10

While this is about data of the city of Hamburg, Germany, I decided to post in English as GeoPackage Propaganda should be accessible. ;)

I am casually working on converting open geodata released via the Transparenzportal Hamburg to more usable GeoPackages, GeoTIFFs and similar formats with free and open-source tools like GDAL and GMT where possible. This includes things like orthophotos, ALKIS, addresses, DEM, districts etc. You can get a list of most available source data here but there are some datasets “hidden” in other categories as well. The data is usually released in GML or as gridded files (e.g. JPEG or XYZ tiles/files). While this is pretty much perfect as source formats, working with them is cumbersome. My goal is to make this data more accessible for anyone in tools like QGIS.

For now you can find the 20cm orthophotos for 2013-2015 and the 1m DEM in Mind the licenses, see the readme file for a bit of info. More to come, I want to plan the pipeline a bit better first though. There should be a full script & log from source to GeoPackage for each file.

I will also provide mirrors of the source files. If you want to collaborate, please contact me. Apart from Hamburg I will also add free/open datasets for the whole of Germany, things related to (nearby) bathymetry and some global ones. If you want Shapefiles, ECW, MrSID or similar, you can pay me for converting.

Leave your thoughts

Merge GML files into one SQLite or Spatialite file

2015-05-20 21:05

For example the buildings in Hamburg, Germany:

for file in *.xml
 if [ -f "${layer}".spatialite ]
  ogr2ogr -f "SQLite" -update -append "${layer}".spatialite "${file}" "${layer}" -dsco SPATIALITE=YES
  ogr2ogr -a_srs EPSG:25832 -f "SQLite" "${layer}".spatialite "${file}" "${layer}" -dsco SPATIALITE=YES

Remove the -dsco SPATIALITE=YES and change the output filename for SQLite. QGIS can work with both.

$ bash AX_Gebaeude

Be aware that Spatialite is much more sensitive to geometric problems. You might get things like

ERROR 1: sqlite3_step() failed:
ax_gebaeude.GEOMETRY violates Geometry constraint [geom-type or SRID not allowed] (19)
ERROR 1: ROLLBACK transaction failed: cannot rollback - no transaction is active
ERROR 1: Unable to write feature 1712 from layer AX_Gebaeude.

ERROR 1: Terminating translation prematurely after failed translation of layer AX_Gebaeude (use -skipfailures to skip errors)

but on the other hand, you get spatial indexing which makes queries or high zoom interaction much quicker.

Be aware that if you try to merge files into a Shapefile and fields are getting truncated, those fields will only be filled with data for the first file you merge. On the later files OGR will try to match the input field names to the merged file’s fieldnames, notice the difference and discard them. If you still want to convert to Shapefiles, check out the -fieldTypeToString IntegerList,StringList options.

Leave your thoughts