Ich möchte der Höhe aus SRTM3-Datendateien eine bestimmte Long / Lat-Position auf einer Karte zuweisen, habe jedoch keine Ahnung, wie der bestimmte Wert ermittelt werden soll. Ich möchte ein Beispiel dafür, wie ich in N50E14.hgt eine Höhe von 50 ° 24'58.888 "N, 14 ° 55'11.377" E finden kann.

Datei Format

Ich nehme es als kleine Übung, wie man einen Datenleser programmiert. Schauen Sie sich die Dokumentation an :

SRTM-Daten werden auf zwei Ebenen verteilt: SRTM1 (für die USA und ihre Gebiete und Besitztümer) mit Daten, die in Intervallen von einer Bogensekunde in Längen- und Breitengraden abgetastet werden, und SRTM3 (für die Welt), die mit drei Bogensekunden abgetastet werden.

Die Daten werden in einer "geografischen" Projektion in Breiten- und Längengrade unterteilt, dh in einer Rasterpräsentation mit gleichen Intervallen von Breiten- und Längengraden in keiner Projektion, die jedoch einfach zu manipulieren und zu mosaikieren ist.

Dateinamen beziehen sich auf den Breiten- und Längengrad der linken unteren Ecke der Kachel - z. B. liegt die linke untere Ecke von N37W105 bei 37 Grad nördlicher Breite und 105 Grad westlicher Länge. Genauer gesagt beziehen sich diese Koordinaten auf das geometrische Zentrum des unteren linken Pixels, das im Fall von SRTM3-Daten eine Ausdehnung von etwa 90 Metern hat.

Höhendateien haben die Erweiterung .HGT und sind mit zwei Byte Ganzzahlen signiert. Die Bytes sind in der "Big-Endian" -Reihenfolge von Motorola angegeben, wobei das höchstwertige Byte zuerst angezeigt wird und von Systemen wie Sun SPARC, Silicon Graphics und Macintosh-Computern mit Power PC-Prozessoren direkt gelesen werden kann. DEC Alpha, die meisten PCs und Macintosh-Computer, die nach 2006 gebaut wurden, verwenden die Intel-Reihenfolge ("Little-Endian"), so dass möglicherweise ein gewisser Byte-Austausch erforderlich ist. Die Höhenangaben in Metern beziehen sich auf das Geoid WGS84 / EGM96. Datenlücken wird der Wert -32768 zugewiesen.

Wie geht es weiter?

Für Ihre Position, 50 ° 24'58.888 "N 14 ° 55'11.377" E, haben Sie bereits die richtige Kachel gefunden, N50E14.hgt. Lassen Sie uns herausfinden, an welchem ​​Pixel Sie interessiert sind. Erster Breitengrad, 50 ° 24'58.888 "N:

24'58.888" = (24 * 60)" + 58.888" = 1498.888"

Bogensekunden. Geteilt durch drei und auf die nächste ganze Zahl gerundet ergibt sich eine Gitterreihe von 500. Dieselbe Berechnung für den Längengrad ergibt die Gitterspalte 1104.

In der Schnellstart-Dokumentation fehlen Informationen darüber, wie Zeilen und Spalten in der Datei organisiert sind. In der vollständigen Dokumentation heißt es jedoch, dass

Die Daten werden in der Hauptreihenfolge der Zeilen gespeichert (alle Daten für Zeile 1, gefolgt von allen Daten für Zeile 2 usw.).

Die erste Zeile in der Datei ist sehr wahrscheinlich die nördlichste, dh wenn wir an Zeile 500 vom unteren Rand interessiert sind , müssen wir uns die Zeile tatsächlich ansehen

1201 - 500 = 701

von anfang an ob die datei . Unsere Gitterzelle ist die Nummer

(1201 * 700) + 1104 = 841804

vom Dateianfang an (dh 700 Zeilen überspringen und im 701. Beispiel 1104 nehmen). Zwei Bytes pro Sample bedeuten, dass wir die ersten 1683606 Bytes in der Datei überspringen und dann zwei Bytes lesen müssen, um unsere Gitterzelle zu erhalten. Die Daten sind Big-Endian, was bedeutet, dass Sie die zwei Bytes auf z. B. Intel-Plattformen austauschen müssen.

Beispielprogramm

Ein einfaches Python-Programm zum Abrufen der richtigen Daten sieht folgendermaßen aus ( Informationen zur Verwendung des struct-Moduls finden Sie in den Dokumenten ):

import struct

def get_sample(filename, n, e):
    i = 1201 - int(round(n / 3, 0))
    j = int(round(e / 3, 0))
    with open(filename, "rb") as f:
        f.seek(((i - 1) * 1201 + (j - 1)) * 2)  # go to the right spot,
        buf = f.read(2)  # read two bytes and convert them:
        val = struct.unpack('>h', buf)  # ">h" is a signed two byte integer
        if not val == -32768:  # the not-a-valid-sample value
            return val
        else:
            return None

if __name__ == "__main__":
    n = 24 * 60 + 58.888
    e = 55 * 60 + 11.377
    tile = "N50E14.hgt"  # Or some magic to figure it out from position
    print get_sample(tile, n, e)

Beachten Sie, dass ein effizienter Datenabruf etwas ausgefeilter aussehen müsste (z. B. das Öffnen der Datei nicht für jede einzelne Probe).

Alternativen

Sie können auch ein Programm verwenden, mit dem die .hgt-Dateien sofort gelesen werden können. Aber das ist langweilig.

GDAL kann diese Rasterformate mit dem SRTMHGT-Treiber lesen / schreiben . Dies bedeutet, dass Sie das Raster mit QGIS, ArcGIS anzeigen oder GDAL-Dienstprogramme wie gdallocationinfo verwenden können , um Werte von einem Punkt abzurufen , z.

DMS in DD konvertieren:

• Lat: 50 ° 24'58.888 "N = 50 + (24/60) + (58.888 / 3600) = 50.4163577778
• Long: 14 ° 55'11.377 "E = 14 + (55/60) + (11.377/3600) = 14.9198269444

Verwenden Sie dann aus einer Shell gdallocationinfo file.hgt -wgs84 long lat:

$ gdallocationinfo N50E14.hgt -wgs84 14.9198269444 50.4163577778
Report:
  Location: (1104P,700L)
  Band 1:
    Value: 216

Die Höhe beträgt 216 m.

Wenn Sie QGIS verwenden, überprüfen Sie, ob das Python-Plugin "Point Sampling Tool" installiert ist. Sie finden es unter -> Verbesserungen (Python) -> Analysieren.

Wählen Sie Ihren Punktlayer für die gewünschten Positionen aus, starten Sie die PST, wählen Sie das HGT (oder eine beliebige Raster- / Polygondatei) und wählen Sie eine neue Punktform für die Ausgabe aus.

Das ist alles :-)

  Chris

Chris 'Antwort zeigt an, dass es einfach ist, Punkte aus einer Ebene in QGIS abzutasten.

Da Ihre Antwort auf meinen Kommentar jedoch verdeutlicht, dass Sie Ihr eigenes Programm zum Lesen von Höhenwerten aus den .hgtDateien schreiben, lesen Sie die Quickstart-PDF-Datei in den SRTM-Dokumenten. Es wird erklärt, wie die Höhendaten gespeichert werden. Zusammenfassen:

• SRTM3-Dateien enthalten eine Folge von Big-Endian-Integer-Werten.
• Jeder ganzzahlige Wert stellt eine Höhe "in Metern, bezogen auf das Geoid WGS84 / EGM96" dar, mit Ausnahme von Werten von -32768, die keine Datenpixel angeben.
• Es gibt 1201 Zeilen mit 1201 Abtastwerten, daher sollten insgesamt 1442401 Ganzzahlwerte vorhanden sein.

Sie sagen, Sie können zwischen Lon / Lat-Koordinaten und Pixeln konvertieren. Um die Höhe zu ermitteln, müssen Sie den ganzzahligen Wert aus dem entsprechenden Versatz in der Datei lesen. Bei gegebenen Pixelkoordinaten xund yrelativ zur oberen linken Ecke der Szene ist das im Grunde genommen offset = (y * 1201) + x. Pixel 0,0ist die erste Ganzzahl in der Datei und Pixel 1200,1200ist die letzte Ganzzahl in der Datei.