Category Archives: GIS

Digitales Oberflächenmodell von Hamburg

Ach nee, der LGV hat ein bildbasiertes DOM von Hamburg im Transparenzportal veröffentlicht… Und sogar gleich zwei, eins von 2018, eins von 2020. Die Rasterweite ist 1 Meter, ob das wohl so vorlag oder für die Veröffentlichung gefiltert wurde?

Egal. Das ist ja großartig! Da werden eine Menge von Anwendungen ermöglicht (Sichtachsen! Verschattungen! Vermaschung! VR! AR!) und verschiedenste Akteure werden die Daten absolut feiern. Auch wenn es mit 1 Meter Auflösung wirklich mies grob ist, auf ein 1 Meter Gitter gerastert ist (nicht ausgedünnt, d. h. es ist teilweise stärker verfälscht und “daneben”) und “nur” bildbasiert (nicht gescannt) ist, geht da schon einiges mit.

Ausprobieren! Im Browser!

Achtung, frickelige Bedienung! Am besten den WASD-Möwen-Modus nutzen, mit Speed 1000. Oder mit einem Doppelklick irgendwo hinzoomen.

https://hamburg.datenatlas.de/DOM1_XYZ_HH_2018_04_30-colored.html
https://hamburg.datenatlas.de/DOM1_XYZ_HH_2020_04_30-colored.html

Datenaufbereitung als LAZ

Für 2018 liegen die Daten als 12768 einzelne XYZ-Kacheln vor, also als super ineffiziente Textdateien. Insgesamt sind es rund 22 Gigabyte. Für 2020 sind es stattdessen 827 größere Kacheln, aber ebenfalls in XYZ mit einem ähnlichem Platzbedarf.

Schönerweise gibt es freie Tools wie txt2las, was sie schnell und einfach ins super effiziente LAZ-Format umwandeln kann:

txt2las -i DOM1_XYZ_HH_2018_04_30/*.xyz -epsg 25832 -odir /tmp/laz/ -olaz

Hat bei mir ungefähr 5 Minuten gebraucht und da waren es nur noch 700 Megabyte. Das ZIP war übrigens mehr als 3 Gigabyte groß.

Zusammengefasst werden können die einzelnen Dateien mit lasmerge:

lasmerge -i /tmp/laz/*.laz -olaz -o DOM1_XYZ_HH_2018_04_30.laz

Interaktive 3D-Webanwendung

Dann noch schnell in den großartigen PotreeConverter von Markus Schütz geschmissen mit

PotreeConverter DOM1_XYZ_HH_2018_04_30.laz -o web/ --encoding BROTLI \
  --generate-page DOM1_XYZ_HH_2018_04_30 --title DOM1_XYZ_HH_2018_04_30

und 4 Minuten später ist die interaktive 3D-Webanwendung fertig (siehe unten), wegen der zusätzlichen Octree-Struktur jetzt bei rund 3 Gigabyte.

Punktwolke mit Farben aus Orthophoto einfärben

Die bereitgestellten Oberflächenmodelle sind so schlicht wie es nur geht, es sind reine XYZ-Daten ohne weitere Dimensionen wie Farbe o. ä.

Glücklicherweise gibt es ja auch die Orthophotos, eventuell wurde sogar dasselbe Bildmaterial genutzt? Da müsste mal jemand durch den Datenwust wühlen, die bei den DOPs werden die relevanten Metadaten nicht mitgeliefert…

Theoretisch könnte man sie also einfärben. Leider ist lascolor proprietär und kommt mit gruseligen, bösartigen Optionen, wenn man es wagt es “unlizenziert” zu nutzen (“Please note that the unlicensed version will (…) slightly change the LAS point order, and randomly add a tiny bit of white noise to the points coordinates once you exceed a certain number of points in the input file.”) und kann JPEG in GeoTIFF nicht lesen (so hab ich mir die DOPs aufbereitet). Eine Alternative ist das geniale PDAL. Mit einer Pipeline wie

{
    "pipeline": [
        "DOM1_XYZ_HH_2020_04_30.laz",
        {
            "type": "filters.colorization",
            "raster": "DOP20_HH_fruehjahrsbefliegung_2020.tif"
        },
        {
            "type": "writers.las",
            "compression": "true",
            "minor_version": "2",
            "dataformat_id": "3",
            "filename":"DOM1_XYZ_HH_2020_04_30-colored.laz"
        }
    ]
}

und

pdal pipeline DOM1_XYZ_HH_2020_04_30.laz+DOP20_HH_fruehjahrsbefliegung_2020_90.cog.tif.json

ist die Punktwolke innerhalb von Minuten coloriert und kann dann wie gehabt mit PotreeConverter in einen interaktiven 3D-Viewer gesteckt werden.

Das Ergebnis ist besser als erwartet, da es scheinbar tatsächlich die selben Bilddaten sind (für beide Jahre). Andererseits ist es auch nicht wirklich schick, da die DOPs nicht als True Orthophoto vorliegen und damit höhere Gebäude gekippt in den Bilder abgebildet sind. Sieht man hier schön am Planetarium.

DOM als GeoTIFF

Wer es lieber als GeoTIFF haben möchte, hat es etwas schwerer, denn GDAL kommt mit dieser Art von Kacheln (mit Lücken und in der bereitgestellten Sortierung) nicht gut klar. Mein Goto-Tool dafür ist GMT.

gmt xyz2grd $(gmt gmtinfo -I- *.xyz) -Vl -I1 -G/tmp/gmt.tif=gd:GTiff *.xyz

Das so erstellte GeoTIFF kann anschließend mit GDAL optimiert werden (ab GDAL 3.2.3/3.3):

gdal_translate -of COG -co COMPRESS=DEFLATE -co PREDICTOR=2 \
  --config GDAL_NUM_THREADS ALL_CPUS --config GDAL_CACHEMAX 50%  \
  -a_srs EPSG:25832 /tmp/gmt.tif DOM1_XYZ_HH_2020_04_30.tif

Hier mal im Vergleich mit dem DGM1 als Schummerungen:

Vermaschung als 3D-Modell

Leider habe ich keine gute Lösung für die 3D-Vermaschung gefunden. tin-terrain und dem2mesh kommen nicht mit so großen Datenmengen auf einmal klar und weiter hab ich nicht geschaut. Wer da was gutes weiß kann sich bei mir bei nächster Gelegenheit Kekse oder Bier abholen. ;)

Daten hinter den Bildern und den Viewern
Datenlizenz Deutschland Namensnennung 2.0, Freie und Hansestadt Hamburg, Landesbetrieb Geoinformation und Vermessung (LGV)

Von wo wurde der Hamburger Fernsehturm fotografiert

Vom großartigen @Fernsehturm_HH-Account inspiriert hab ich mal ein 4 Jahre altes Projektchen aufgewärmt: Verortete Fotos auf Flickr, die mit “Hamburg” und “Fernsehturm” getaggt sind. Alternativ als Linien vom (angeblichen) Aufnahmeort zum Turm.

Image

Ungefiltert und voller Mist, aber man kann z. B. toll sehen, wie einladend die Brückem vorm Kuhmühlenteich für Fotos sind.

Kostenlose Internetpunkteabsahnidee: Dasselbe mit dem Eiffelturm, dem Londoner Dildo u. ä. machen, Konkave Hülle außen drum, schick aufbereiten, als “Nur wo man $Wahrzeichen sieht, ist $Stadt” vermarkten, €€€€.

Interaktive Karte der Baugenehmigungen in Hamburg

Ich habe endlich mal ein über mehrere Jahre gereiftes (eher gealtertes und verfrickeltes) Projekt in einen einigermaßen vorzeigbaren Zustand gebracht: Eine interaktive Karte der Baugenehmigungen in Hamburg.

Je weniger transparent eine Fläche dargestellt ist, desto mehr Dokumente sind mit ihr verknüpft (ja, es ist ein Feature je Dokument D;). Eigentlich war die Seite anders aufgebaut, mit einem PDF-Viewer auf der rechten Seite. Aber da daten.transparenz.hamburg.de kein HTTPS kann (seid ihr auch so gespannt auf die UMTS-Auktion nächste Woche?), geht das aus Sicherheitsgründen nicht ohne ein Spiegeln der Daten oder einen Proxy.

Die Daten kommen größtenteils aus dem Transparenzportal. Für das Matching der angebenen Flurstücks-“IDs” zu den tatsächlichen Flurstücken war aber ein erheblicher Aufwand nötig. Das Drama ging bis hin zum Parsen aus PDFs, die mal so, mal so formatiert waren und natürlich auch voller Eingabefehler auf Behördenseite. Vielleicht schreibe ich da noch beizeiten mal einen Rant. TL;DR: Ohne die zugehörige Gemarkung ist mit einer Flurstücks-“ID” wie in den Daten angegeben keine räumliche Zuordnung möglich. In den veröffentlichten Daten stecken nur die Nenner der Flurstücksnummern, nicht aber die Gemarkungsnummern. Ziemlich absurd.

Das ganze ist nur ein Prototyp, vermutlich voller Fehler und fehlender Daten. Aber interessant und spaßig ist es, viel Freude also!

Es wäre noch eine MENGE zu tun, um das ganze rund zu machen. Falls du Lust hast, melde dich gerne. Es geht vom wilden Parsen, über Sonderregeln für kaputte Dokumente, zu Kartenstyling bis zur UI. Schön wäre es auch alles in einer anständigen Datenbank zu halten und nicht nur nach der räumlichen Dimension durchsuchen zu können.

Highlight current timeslice in a QGIS Atlas layout

Did this for an ex-colleague some months ago and forgot to share the how-to publically. We needed a visual representation of the current time in a layout that showed both a raster map (different layer per timeslice) and a timeseries plot of an aspect of the data (this was created outside QGIS).

Have lots of raster layers you want to iterate through. I have:

./ECMWF_ERA_40_subset/2019-01-01.tif
./ECMWF_ERA_40_subset/2019-01-02.tif
./ECMWF_ERA_40_subset/2019-01-03.tif
./ECMWF_ERA_40_subset/2019-01-04.tif
./ECMWF_ERA_40_subset/2019-01-05.tif
./ECMWF_ERA_40_subset/2019-01-06.tif
...

Create a new layer for your map extent. Draw your extent as geometry. Duplicate that geometry as many times as you have days. Alternatively you could of course have different geometries per day. Whatever you do, you need a layer with one feature per timeslice for the Atlas to iterate though. I have 30 days to visualise so I duplicated my extent 30 times.

Open the Field Calculator. Add a new field called date as string type (not as date type until some bug is fixed (sorry, did not make a note here, maybe sorting is/was broken?)) with an expression that represents time and orders chronologically if sorted by QGIS. For example: '2019-01-' || lpad(@row_number,2,0) (assuming your records are in the correct order if you have different geometries…)

Have your raster layers named the same way as the date attribute values.

Make a new layout.

For your Layout map check “Lock layers” and use date as expression for the “Lock layers” override. This will now select the appropriate raster layer, based on the attribute value <-> layer name, to display for each Atlas page.

Cool, if you preview the Atlas now you got a nice animation through your raster layers. Let’s do part 2:

In your layout add your timeseries graph. Give it a unique ID, e. g. “plot box”. Set its width and height via new variables (until you can get those via an expression this is needed for calculations below).

Create a box to visualise the timeslice. Set its width to map_get(item_variables('plot box'), 'plot_width') / @atlas_totalfeatures. For the height and y use/adjust this expression: map_get(item_variables('plot box'), 'plot_height'). For x comes the magic:

with_variable(
'days_total',
day(to_date(maximum("date"))-to_date(minimum("date")))+1,
-- number of days in timespan
-- +1 because we need the number of days in total
-- not the inbetween, day() to just get the number of days
with_variable(
'mm_per_day',
map_get(item_variables('plot box'), 'plot_width') / @days_total,
with_variable(
'days',
day(to_date(attribute(@atlas_feature, 'date'))-to_date(minimum("date"))),
-- number of days the current feature is from the first day
-- to_date because BUG attribute() returns datetime for date field
@mm_per_day * @days + map_get(item_variables('plot box'), 'plot_x')
)
)
)

This will move the box along the x axis accordingly.

Have fun!

How to get a dataset out of OpenStreetMap (OSM) and into QGIS easily

For small datasets just use the QuickOSM plugin, enter your key=value and it will load the data directly as QGIS layer(s).

For big datasets don’t run massive queries on the Overpass API but use prepared thematic or regional OpenStreetMap extracts.

For example, if you want all things tagged place=town globally, grab place_EPSG4326.zip from http://download.osmdata.xyz/ and run a filter on that locally.

Don’t be scared of processing a global OSM planet file either, Osmium makes filtering OSM data extremely easy:

osmium tags-filter in.osm.pbf n/place=town -o place=town.osm.pbf

Yes, a planet file is pretty big, but extracting specific features from that is not a big data problem and you must not be scared of it. Downloading the file will probably take you magnitudes longer than extracting something from it. For me it was 45 minutes for the download, then about 8 minutes for extracting on a seriously slow (~60MB/s) spinning metal hard disk drive (no SSD).

And yes, you can load OSM PBF directly into QGIS thanks to GDAL’s support for the format.

While we are at it: Don’t use GeoJSON for anything but data transfer and maybe storage. It is not an efficient format to power your layers and GIS analyses. (OSM PDF isn’t either.)

Make sure you actually do use spatial indexes


Ever ran some GIS analysis in QGIS and it took longer than a second? Chances are that your data did not have spatial indexes for QGIS to utilise and that it could have been magnitudes faster.

I realised just today, after years of using QGIS, that it did not automatically create a spatial index when saving a Shapefile. And because of that, lots of GIS stuff I did in the past, involving saving subsets of data to quick’n’dirty Shapefiles, was slower than necessary.

Sadly QGIS does not communicate lack of spatial indexing to the user in any way. I added a feature request to make Processing warn if no indexing is available.

An example: Running ‘Count points in polygon’ on 104 polygons with 223210 points:

  • Points in original GML file: 449 seconds
    • GML is not a format for processing but meant for data transfer, never ever think of using it for anything else
  • Points in ESRI Shapefile: 30 seconds
  • Points in GeoPackage: 3 seconds
  • Points in ESRI Shapefile with spatial index: 3 seconds
    • Same Shapefile as before but this time I had created a .qix index

So yeah, make sure you don’t only use a reasonable format for your data. And also make sure you do actually have an spatial index.

For Shapefiles, look for a .qix or .sbn side-car file somewhere in the heap of files. In QGIS you can create a spatial index for a vector layer either using the “Create spatial index” algorithm in Processing or the button of the same title in the layer properties’ source tab.

PS: GeoPackages always have a spatial index. That’s reason #143 why they are awesome.

Calculating the metric distance between latitudes with PostGIS

Postgres’ WITH clauses and window functions are so awesome.

-- generate values from -90 to 90, increment by 1
WITH values AS (
  SELECT generate_series AS latitude
  FROM generate_series(-90, 90)
),
-- create a geographic point each
points AS (
  SELECT ST_MakePoint(0, latitude)::geography AS point FROM values
)
SELECT
  -- latitude values of subsequent points
  format(
    '%s° to %s°',
    ST_Y(point::geometry),
    ST_Y(lag(point::geometry) OVER ())
  ) AS latitudes,
  -- geographic distance between subsequent points, formatted to kilometers
  format(
    '%s km',
    to_char(ST_Distance(point, lag(point) OVER ())/1000, '999D99')
  ) AS distance
FROM points
OFFSET 1  -- skip the first row, no lag there
;

    latitudes  |  distance  
 --------------+------------
  -89° to -90° |  111.69 km
  -88° to -89° |  111.69 km
  -87° to -88° |  111.69 km
  -86° to -87° |  111.69 km
  -85° to -86° |  111.69 km
  -84° to -85° |  111.68 km
  -83° to -84° |  111.68 km
  -82° to -83° |  111.67 km
  -81° to -82° |  111.67 km
  -80° to -81° |  111.66 km
  -79° to -80° |  111.66 km
  -78° to -79° |  111.65 km
  -77° to -78° |  111.64 km
  -76° to -77° |  111.63 km
  -75° to -76° |  111.62 km
  -74° to -75° |  111.61 km
  -73° to -74° |  111.60 km
  -72° to -73° |  111.59 km
  -71° to -72° |  111.58 km
  -70° to -71° |  111.57 km
  -69° to -70° |  111.56 km
  -68° to -69° |  111.54 km
  -67° to -68° |  111.53 km
  -66° to -67° |  111.51 km
  -65° to -66° |  111.50 km
  -64° to -65° |  111.49 km
  -63° to -64° |  111.47 km
  -62° to -63° |  111.45 km
  -61° to -62° |  111.44 km
  -60° to -61° |  111.42 km
  -59° to -60° |  111.40 km
  -58° to -59° |  111.39 km
  -57° to -58° |  111.37 km
  -56° to -57° |  111.35 km
  -55° to -56° |  111.33 km
  -54° to -55° |  111.31 km
  -53° to -54° |  111.30 km
  -52° to -53° |  111.28 km
  -51° to -52° |  111.26 km
  -50° to -51° |  111.24 km
  -49° to -50° |  111.22 km
  -48° to -49° |  111.20 km
  -47° to -48° |  111.18 km
  -46° to -47° |  111.16 km
  -45° to -46° |  111.14 km
  -44° to -45° |  111.12 km
  -43° to -44° |  111.10 km
  -42° to -43° |  111.08 km
  -41° to -42° |  111.06 km
  -40° to -41° |  111.04 km
  -39° to -40° |  111.03 km
  -38° to -39° |  111.01 km
  -37° to -38° |  110.99 km
  -36° to -37° |  110.97 km
  -35° to -36° |  110.95 km
  -34° to -35° |  110.93 km
  -33° to -34° |  110.91 km
  -32° to -33° |  110.90 km
  -31° to -32° |  110.88 km
  -30° to -31° |  110.86 km
  -29° to -30° |  110.84 km
  -28° to -29° |  110.83 km
  -27° to -28° |  110.81 km
  -26° to -27° |  110.80 km
  -25° to -26° |  110.78 km
  -24° to -25° |  110.77 km
  -23° to -24° |  110.75 km
  -22° to -23° |  110.74 km
  -21° to -22° |  110.72 km
  -20° to -21° |  110.71 km
  -19° to -20° |  110.70 km
  -18° to -19° |  110.69 km
  -17° to -18° |  110.67 km
  -16° to -17° |  110.66 km
  -15° to -16° |  110.65 km
  -14° to -15° |  110.64 km
  -13° to -14° |  110.63 km
  -12° to -13° |  110.63 km
  -11° to -12° |  110.62 km
  -10° to -11° |  110.61 km
  -9° to -10°  |  110.60 km
  -8° to -9°   |  110.60 km
  -7° to -8°   |  110.59 km
  -6° to -7°   |  110.59 km
  -5° to -6°   |  110.58 km
  -4° to -5°   |  110.58 km
  -3° to -4°   |  110.58 km
  -2° to -3°   |  110.58 km
  -1° to -2°   |  110.58 km
  0° to -1°    |  110.57 km
  1° to 0°     |  110.57 km
  2° to 1°     |  110.58 km
  3° to 2°     |  110.58 km
  4° to 3°     |  110.58 km
  5° to 4°     |  110.58 km
  6° to 5°     |  110.58 km
  7° to 6°     |  110.59 km
  8° to 7°     |  110.59 km
  9° to 8°     |  110.60 km
  10° to 9°    |  110.60 km
  11° to 10°   |  110.61 km
  12° to 11°   |  110.62 km
  13° to 12°   |  110.63 km
  14° to 13°   |  110.63 km
  15° to 14°   |  110.64 km
  16° to 15°   |  110.65 km
  17° to 16°   |  110.66 km
  18° to 17°   |  110.67 km
  19° to 18°   |  110.69 km
  20° to 19°   |  110.70 km
  21° to 20°   |  110.71 km
  22° to 21°   |  110.72 km
  23° to 22°   |  110.74 km
  24° to 23°   |  110.75 km
  25° to 24°   |  110.77 km
  26° to 25°   |  110.78 km
  27° to 26°   |  110.80 km
  28° to 27°   |  110.81 km
  29° to 28°   |  110.83 km
  30° to 29°   |  110.84 km
  31° to 30°   |  110.86 km
  32° to 31°   |  110.88 km
  33° to 32°   |  110.90 km
  34° to 33°   |  110.91 km
  35° to 34°   |  110.93 km
  36° to 35°   |  110.95 km
  37° to 36°   |  110.97 km
  38° to 37°   |  110.99 km
  39° to 38°   |  111.01 km
  40° to 39°   |  111.03 km
  41° to 40°   |  111.04 km
  42° to 41°   |  111.06 km
  43° to 42°   |  111.08 km
  44° to 43°   |  111.10 km
  45° to 44°   |  111.12 km
  46° to 45°   |  111.14 km
  47° to 46°   |  111.16 km
  48° to 47°   |  111.18 km
  49° to 48°   |  111.20 km
  50° to 49°   |  111.22 km
  51° to 50°   |  111.24 km
  52° to 51°   |  111.26 km
  53° to 52°   |  111.28 km
  54° to 53°   |  111.30 km
  55° to 54°   |  111.31 km
  56° to 55°   |  111.33 km
  57° to 56°   |  111.35 km
  58° to 57°   |  111.37 km
  59° to 58°   |  111.39 km
  60° to 59°   |  111.40 km
  61° to 60°   |  111.42 km
  62° to 61°   |  111.44 km
  63° to 62°   |  111.45 km
  64° to 63°   |  111.47 km
  65° to 64°   |  111.49 km
  66° to 65°   |  111.50 km
  67° to 66°   |  111.51 km
  68° to 67°   |  111.53 km
  69° to 68°   |  111.54 km
  70° to 69°   |  111.56 km
  71° to 70°   |  111.57 km
  72° to 71°   |  111.58 km
  73° to 72°   |  111.59 km
  74° to 73°   |  111.60 km
  75° to 74°   |  111.61 km
  76° to 75°   |  111.62 km
  77° to 76°   |  111.63 km
  78° to 77°   |  111.64 km
  79° to 78°   |  111.65 km
  80° to 79°   |  111.66 km
  81° to 80°   |  111.66 km
  82° to 81°   |  111.67 km
  83° to 82°   |  111.67 km
  84° to 83°   |  111.68 km
  85° to 84°   |  111.68 km
  86° to 85°   |  111.69 km
  87° to 86°   |  111.69 km
  88° to 87°   |  111.69 km
  89° to 88°   |  111.69 km
  90° to 89°   |  111.69 km

					

Open Layers View Tracker

I built this last year for some research and then swiftly forgot about releasing it to the public. Here it is now:

https://gitlab.com/Hannes42/OpenLayersViewTracker

Try it online: https://hannes42.gitlab.io/OpenLayersViewTracker/

Some awful but working JavaScript code to track a website user’s interaction with a Open Layers map. You can use this to do awesome user studies and experiments.

  • Runs client-side
  • You will get a polygon of each “view”!
  • You can look at them in the browser!
  • There are also timestamps! Hyperaccurate in milliseconds since unix epoch!
  • And a GeoJSON export!
  • This works with rotated views!
  • Written for Open Layers 4 using some version of JSTS, see the libs/ directory. No idea if it works with the latest versions or if Open Layers changed their API again.

Please do some funky research with it and tell me about your experiences! Apart from that, you are on your own.

There is a QGIS project with example data included. Check out the Atlas setup in the Print Layout!

Screenshot from a browser session

Resulting GeoJSON in QGIS

So if Time Manager supports the timestamp format you could interactively scroll around. I did not try, that plugin is so finicky.

Replaying what I looked at in a Open Layers web map, using QGIS Atlas

Average Earth from Space 2018

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 https://worldview.earthdata.nasa.gov/, grabbed Soumi VIIRS images via https://wiki.earthdata.nasa.gov/display/GIBS/GIBS+API+for+Developers 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. ;-)