Zukünftig weltweit 30m SRTM-Daten verfügbar

Ich arbeite zur Zeit an einem Projekt mit Geheimhaltung, aber ich darf euch kurz mitteilen, wieso genaue Höhendaten interessant sind:

Ein moderner Hybrid benötigen Höhendaten, um zu kalkulieren, wann sie auf der Route wo aufladen können. Wenn man z.B. von einer Landstraße bergab fährt, und am Ende des Berges ein Dorf/eine Stadt ist, springt der Motor so früh an, dass er durch die Stadt mit dem Elektromotor fahren kann. So werden Emissionen (Lärme und Abgase) in der Stadt verhindert.

Das Projekt ist sehr interessant und noch vor wenigen Jahren wäre keiner auf die Idee gekommen, Höhendaten für sowas zu nutzen.

Das Beispiel habe ich mir in QGIS mal genauer angeschaut und einen DEM-Vergleich gebastelt:
DEM Vergleich

Das Höhenprofil in BRouter Web (CGIAR-CSI SRTM 90m) zeigt bis zu 9m Höhendifferenz bei eigentlich flacher Straße. Leider scheint sich auch beim neuen 30m SRTM das Gebäude auf die Umgebung auszuwirken, da hätte ich eine genauere Trennung erwartet.

Links: Gebäude in OpenStreetMap und StreetView. DEM: NASA SRTM 1 Arc-Second Global, CGIAR-CSI SRTM, EU-DEM

.

Sehr interssant…

zwei Dinge fallen mir natürlich auf:

Einmal das Alignment: wenn SRTM-90 so wie dokumentiert durch Mittelung von jeweils 9 Pixeln in SRTM-30 entstanden ist, dann müssten die Erhebungen in beiden Bildern den gleichen Schwerpunkt haben, es scheint aber eindeutig verschoben. Ich nehme an, dass Dir QGIS das Alignment aus den Meta-Daten der Dateien selbst gemacht hat und Du unschuldig bist? Also muss der Fehler entweder in QGIS sein oder in einem der beiden Datensätze.

Und EU-DEM scheint hier subtanziell anders zu sein. so, als hätten sie das Gebäude weg-gefiltert. Nach so einem Filter suche ich noch, der nicht einfach nur Gebäude in die Breite zerlaufen lässt, sondern solche spitzen Strukturen tatsächlich als “ungeografisch” erkennt und wegfiltert.

Die erwähnen nichts vom Filtern, nur vom gewichteten Mitteln: " based on SRTM and ASTER GDEM data fused by a weighted averaging approach" (steht da). Vielleicht hat an dieser Stelle einfach das flachere ASTER beim Gewichten gewonnen.

Die heruntergeladenen GeoTiffs hab ich einfach direkt als Rasterlayer in QGIS hinzugefügt. Zum Hinterlegen mit Bing per OpenLayers Plugin hab ich das Projekt-KBS auf Pseudo Mercator / EPSG:3857 mit Auto-Transformation gestellt, der Versatz ist aber auch ohne da.

Ich bin da schon etwas arglos rangegangen und hab mich darauf verlassen, dass QGIS das schon richtig macht, wobei ich mich nicht wirklich gut damit auskenne und eine Fehlbedienung nicht ausschließen will. Etwas wundern tu ich mich über das Ergebnis schon auch, kann ja mal schauen wie ich das verifizieren kann und evtl. auch mein Vorgehen nachvollziehbar dokumentieren.

Leicht süd-östlich zur Straße hin versetzt gibt es ja auch bei EU-DEM noch eine Erhebung, mit einem Maximalwert von 7.739999771118164 zwar deutlich abgeflacht, aber trotzdem noch erkennbar.

Ich hab’ das nochmal geprüft: Mit dem CGIAR-Datensatz der Version 4.1 (da wo die Kacheln 6001*6001 Pixel haben) passt das vom relativen Alignment gegenüber SRTM 30m.

Das “+1” in 6001 ist eine Erleichterung für die Interpolation, so dass man die Randzeilen immer doppelt hat. Das Pixel der Spalte 0 hat seinen Schwerpunkt (nicht seinen linken Rand) dann genau auf dem Längengrad, der die Kachel begrenzt. Bei dem CGIAR Daten der Vesion 4.0 gab es das glaubich nicht, da waren es 6000*6000 Pixel.

Jedenfalls liegt der gemeinsame Mittelpunkt der beiden hellsten Pixel in allen Datensätzen am identischen Punkt, und der trifft in das Gebäude auch in Google-Earth leicht (ca 30m) südlich, so wie im linken Bild von ikonor zu sehen.

Das gleiche gilt auch für den 90m Datensatz, den man bei http://earthexplorer.usgs.gov/ laden kann, aber hier sind seltsamerweise dann die Werte anders und auch das hab’ ich mir angeschaut, und es scheint so zu sein dass:

  • bei CGIAR 4.1 die 90m-Pixel tatsaechlich durch Mittelung der jeweils 9 Pixel entstanden sind

  • bei USGS aber der Wert des jeweils zentralen Pixels genommen wurde.

so dass

  • diese beiden 90m Datensaetze also eine unterschiedliche Charakteristik haben und für Routingzwecke nicht unbedingt austauschbar sind.
  • man mit dem 30m Datensatz aber je nach Wahl des eigenen Filters beides wieder nachbauen kann

Kommt man hier nicht mit einer Häufigkeitsanalyse über ein 5x5 oder 7x7 Raster am weitesten? Ein einzelnes Hochhaus erzeugt nur sehr wenige Abweichungen, dafür aber große. Ein Berg oder ähnlich erzeugt viele Abweichungen, im Schnitt dann eher kleine. Eine Klippe erzeugt auch ein sehr markantes Profil. Damit könnte man automatisch die Daten Klassifizieren.

Und dann stellt man den Hochpassfilter entsprechend starkt ein, je nachdem, was die vorherige Klassifizierung ergeben hat.

Ich hab jetzt ne ganze Weile rumgesucht und -probiert, kann mir das Ganze aber leider immer noch nicht wirklich erklären.
Bisherige Infos hab ich mal unter meiner User-Seite im Wiki gesammelt.

Ich hatte das CGIAR v4.1 GeoTIFF vom JRC (IT) Mirror verwendet, nur mit diesem tritt der Versatz gegenüber Original SRTM auf, mit den GeoTIFF und ArcASCII von CGIAR-CSI (USA) nicht. Alle sind v4.1 laut readme.txt, jedoch das von JRC mit 6000 Pixeln und geradem Ursprung, die anderen mit 6001 Pixeln und ungeradem Ursprung.

Bei gleichem Schwerpunkt der GeoTIFFs (AREA_OR_POINT=Area) und anderem Usprung erscheint mir der Versatz logisch. Ich verstehe nur nicht mehr, warum die anderen bei unterschiedlicher Schwerpunktdefinition übereinstimmen.

Beim GeoTIFF von cgiar.org mit AREA_OR_POINT=Area verstehe ich das auch so, beim ArcASCII würde ich aber “xllcorner” als Schwerpunkt auf dem linken Rand interpretieren (PixelIsArea vs. PixelIsPoint).

Da hat das GeoTIFF AREA_OR_POINT=Point in den Metadaten gesetzt, würde ich also als Schwerpunkt auf dem Rand interpretieren und verstehe nicht, warum das Alignment trotzdem mit AREA_OR_POINT=Area bei CGIAR übereinstimmt (beide +1 Pixel und ungerader Ursprung).

Da hat es wohl zwei Quellen mit unterschiedlichen Varianten bei Version 2 gegeben (Sampling vs. Averaging).

Bei der 90m Version 3.0 (SRTM Plus) vom November 2013 gibt es beide Varianten (averaged und sub-sampled). Löcher sind dort mit ASTER GDEM gefüllt, sodass ein Vergleich mit CGIAR-CSI evtl. schon interessant wäre. Beim laufenden 30m Release steht zwar auch was von “void filled”, aber keine Details und ich bin mir nicht sicher, ob das derselbe Datenbestand ist.

naja, die Erklärung ist wohl einfach, dass sie das beim JRC mit dem Alignment nicht so genau genommen haben.

hier gibt’s noch eine Analyse zu den Fragen: http://imagico.de/map/srtm1_de.php

Ich habe in meiner Instanz von brouter-web jetzt mal die 30m Daten für vier 5*5 Kacheln eingesetzt (N35/W10 - N40/E10) das ist so ein Streifen von Portugal bis Korsika, jedenfalls ist Mallorca dabei.

Als Filter habe ich keinen diskreten Filter auf dem Raster benutzt, sondern die Interpolation auf diesem Raster ersetzt: Statt wie bisher 2-dimensional linear zwischen Pixeln mitteln (also die Flaechenanteile nehmen von einem Rechteck in Pixelgroesse, was man ueber das Raster schiebt), nehme ich jetzt eine Kreisscheibe von 50m radius und gewichte die Pixel nach den Flaechenanteilen der Ueberlappung).

Aus den 50m hab’ ich zwar keine Wissenschaft gemacht, scheint mir aber plausibel, weil bei dem bisherigen Verfahren mit 90m Daten ist der effektive Einflussradius ja je nach Alignment mit dem Gitter irgendwo zwischen 45m un 90m. Aber vielleicht kann ich wirklich noch weiter runter mit dem Radius, ohne Probleme mit dem Rauschen zu bekommen.

Man kann in BRouter-Web den Unterschied zwischen beiden Datensaetzen gut anschauen, indem man per Profil-Upload das Flag
“assign forceSecondaryData 1” im global-Kontext setzt, dann wird auf den “lastrun-backup” umgeschaltet.

Man findet durchaus Stellen, wo es einen signifikanten Unterschied macht, z.B. hier:

http://brouter.de/brouter-web/#zoom=14&lat=39.92777&lon=3.11565&layer=OpenStreetMap&lonlats=3.111706,39.929024|3.12398,39.931656&nogos=&profile=trekking&alternativeidx=0&format=geojson

Über den Earth Explorer lässt sich inzwischen auch Europa herunterladen: http://earthexplorer.usgs.gov/
Als Veröffentlichungsdatum steht der 23. September 2014, was aber nicht ganz stimmen kann.

Die 90m-Daten gibt es in einer “Void Filled”-Version und in einer “Non-Void Filled”-Version. Die 30m-Daten nur in der Version mit Löchern (“Some tiles may still contain voids”).

Steht irgendwo, ob es auch 30m-Daten mit geflickten Löchern geben wird? Ich finde nichts, aber das “may still contain voids” macht mir Hoffnung…

Grüße, Max

Kannst Du Dir selber basteln, indem Du überall da, wo im 30m Datensatz “void” steht, den Wert aus dem 90m-void-filled Datensatz nimmst.

Zu void-filled bei 30m hab ich auch keine weiteren Infos gesehen. Auf der NASA SRTM Homepage gibt es jedoch einen Eintrag “NASADEM: What’s Next with SRTM”, nach dem eine Überarbeitung der SRTM Daten unter dem Namen NASADEM stattfindet.

Ich habe ein Tool Leaflet-raw-DEM angefangen, mit dem Ziel Höhendaten in Rohform im Browser anzeigen und vergleichen zu können:

Das ist allerdings noch nicht ausgereift und mit Vorsicht zu genießen, im Moment mache ich aber mit anderen Dingen weiter.

Seit kurzem sind zumindest über den Earthexplorer die 30m SRTM-Daten auch für den Nahen Osten verfügbar.

Wurde die Lizenz von usgs schon überprüft. Mir scheint da das alte Problem zu sein, dass die Namensnennung wollen. “we request that the following statements be used when citing, copying, or reprinting data:…”
Reicht da eine generelle Nennung, oder eine im Änderungssatz?
Gibt es einen direkten Zugriff auf die Bilder für JOSM. In Afrika gibt es nämlich noch große weiße Flecken ohne nutzbare Luftbilder, da wäre das schon gut.

PS: Typisch USA, dass die Goggle Maps darüber gelegt haben…

edit:
Habe was über wms gefunden: https://lta.cr.usgs.gov/sites/default/files/WMShelp_0.pdf
Scheinbar muss man sich anmelden und jeweils eine Karte bestellen.

Das sind keine Karten, sondern nur die Höhendaten im 30m-Gitter. Aus denen werden dann z.B. Höhenlinien generiert.

Hier noch eine Meldung vom LP DAAC (auch in der Wochennotiz erwähnt):
https://lpdaac.usgs.gov/nasa_shuttle_radar_topography_mission_srtm_global_1_arc_second_data_released_over_middle_east

Beim LP DAAC gibt es auch einen direkten HTTP Download im HGT-Format (siehe auch SRTM im Wiki):
http://e4ftl01.cr.usgs.gov/SRTM/SRTMGL1.003/2000.02.11/

(via LP DAAC SRTMGL1 > Data Access > Data Pool > Direct HTTP Access > SRTM
oder https://lpdaac.usgs.gov/ Menü: Tools > Data Access > Data Pool > …)

Bin mir aber nicht sicher, ob das jetzt exakt die gleichen Daten sind wie vom USGS (im Earthexplorer separat unter NASA LPDAAC Collections).

Ich greif den Faden hier nochmal auf, weil ich es ja immer noch nicht geschafft habe, die SRTM-1 Daten für BRouter produktiv zu kriegen. Was mir noch gefehlt hatte war ein funktionierender Merge zwischen den SRTM-1 DAten und den “guten alten” CGIAR 4.1 Daten. Woran ich dabei fast verzweifelt war ist diese (völlig sinnlose…) 50-Degree Grenze in den Daten: oberhalb 50 Nord und unterhalb 50 Süd ist der longitudinale Pixel-Abstand verdoppelt. Sinnlos ist das deswegen, weil es vom Datenvolumen ja nur ca 10% ausmacht, und da hätte man ja auch fünfe gerade sein lassen können. Vielleicht war’s aber eher auch ein Rausch-Problem.

Fast verzweifelt war ich an der 50 Grad Grenze, weil sie ja in den CGIAR nicht sichtbar war. Sie war aber trotzdem da, denn wenn man sich die >50 Tiles dort genau anschaut sieht man, dass die “ungeraden” Spallten nur Mittelwerte der geraden Spalten sind. Die geraden Spalten wiederum sind die 3x3 Mittel der 1sec*2sec Pixel aus den Ursprungsdaten. Wenn man das mal verstanden hat, ist es auch für die >50 Grad Tiles ganz einfach, die CGIAR-4.1 und die SRTM1 Daten aufeinanderzulegen, und dann findet man in ca 99% eine Übereinstimmung (± 1/2m), wenn man für die SRTM1 Daten das 3x3 Mittel rechnet. Und das verbleibende 1% sind ja genau die Problemfälle, für die ich es bei den CGIAR-4.1 Daten belassen möchte.

Jetzt hab’ ich immerhin geschafft, den Merge zu machen und neue Tiles zu berechnen mit einer etwas breiteren Überlappung und ohne Artefakte an der 50-Grad Grenze.

Aber jetzt hab’ ich noch ein IT-Probem: Der Strato-Level-2 V-Server, der das alles rechnet hat ein extremes Storage-Problem, das den wöchtentlichen BRouter-Daten-Refresh schon jetzt grenzwertig lang werden lässt. Ich hab’ mir zwar schon Mühe gegeben, die neuen Höhendaten in ein kompaktes Format zu speichern (deutlich kompakter als die *bil.zip Daten vom earth-explorer), trotzdem is es auch dieser Infrastruktur zur Zeit nicht vertretbar, statt der 15 GB CGIAR-4.1 Daten jetzt ca 70 GB SRTM1-Daten im wöchentlichen Lauf zu lesen…