.hgt umwandeln für arcgis

Wenn du dir nicht sicher bist, in welchem der anderen Foren du die Frage stellen sollst, dann bist du hier im Forum für allgemeine Fragen sicher richtig.
Antworten
apogis
User
Beiträge: 5
Registriert: Dienstag 25. Mai 2010, 13:12

Hallo,
erstmal muss ich mich im Vorraus schonmal entschuldigen, da ich den Thread: 'An alle Schüler & Studenten mit Informatikproblemen' gerade gesehen habe; ich habe nämlich leider auch keine Ahnung von Python, aber auch niemanden, den ich sonst um Hilfe fragen könnte.
Bei mir geht es nicht um eine Hausaufgabe, sondern einen HiWi-Job, bei dem ich für den Geomorphologie-Lehrstuhl ein Höhenmodell erstellen & daraus die Lage verschiedener Terrassen & daraus die abgelaufenen Prozesse bestimmen soll, leider kann ich die Daten, die ich dafür erhalten habe nicht mal anschauen geschweige den irgendwie damit arbeiten. Die Daten sind in einer .hgt Datei, die 1.244.950.160 Bytes groß ist, Höhen sind dort wohl in 16bit signed integer enthalten & ich würde sie gerne in einem GIS darstellen.
ArcGIS hat zwar ein Tool dafür, aber nur für Daten, die eine Auflösung von 3' (2,884,802 bytes) bzw. 1' (25,934,402 bytes) haben. Meine Daten sind höher aufgelöst, also muss ich mir dafür wohl selber etwas zurechtschreiben. Hat aber bis jetzt nicht geklappt.
Zu dem Thema habe ich hier: http://stackoverflow.com/questions/3574 ... nary-files ein Python-Script gefunden & versucht es nachzuschreiben, sieht bei mir bis jetzt so aus:

Code: Alles auswählen

>>>from struct import unpack, calcsize 
>>>def read_row( row, row_lenght): 
...  format = "!h" 
...  for i in range(0, row_lenght): 
...   offset = i * calcsize(format) 
...   (height,) = unpack(format, row[offset : offset+calcsize(format)]) 
... 
da hört das Script in dem Beispiel auf & da ich noch nie was mit Python gemacht habe, weiß ich nicht was ich jetzt noch tun muss.
Ich hätte die Datei danach gerne als .grd vorliegen, aber bis jetzt hab ich ja noch gar keine Parameter eigegeben. Werde ich die erst eigeben müssen, wenn ich das Script in Arcgis ausführe?
& muss ich jetzt noch sagen, dass ich die Datei speichern will? & wie kann ich das Script speichern (das ist wahrscheinlich die peinlichste Frage von allen...)?

Wäre großartig wenn einer von euch irgendwelche Tipps hat, oder mir sagen kann wo ich soetwas nachlesen könnte.

Verzweifelte Grüße
Zuletzt geändert von Anonymous am Dienstag 25. Mai 2010, 14:43, insgesamt 1-mal geändert.
Grund: Quelltext in Code-Tags gesetzt.
frabron
User
Beiträge: 306
Registriert: Dienstag 31. März 2009, 14:36

Wo kommt die hgt-Datei denn her? Denn das ist doch das Format für SRTM Daten. Die sind halt nur spezifiziert für 1, bzw. 3 Grad Auflösung, d.h. die meisten Tools werden das nicht einlesen können. Bevor du dich hier abquälst: Bekommst du die Daten nicht in einem anderen Format?

Wenn mich meine Rechenkünste nicht im Stich gelassen haben, ist die Ausdehnung deiner Datei ja 24.949,45 x 24.949,45 x 2, wobei mich da schon die Nachkommastelle stört( siehe http://www2.jpl.nasa.gov/srtm/faq.html zur Berechnung )

Zum Code selber kann ich leider nix beitragen, ausser das in dem Beispiel von einer bekannten row, bzw. height length ausgegangen wird (das ist das, was ich oben berechnet hab). Kennst du die?
BlackJack

Also ich würde das ja nicht in 2-Byte-Pärchen einlesen und mit `struct` umwandeln, sondern grössere "Streifen" mittels `array`-Modul aus der Standardbibliothek oder besser noch `numpy`.

Vielleicht gibt's ja auch eine Bibliothek mit der man effizient Ausschnitte von so einer grossen "raw"-Datei in `numpy`-Arrays einlesen kann!? Kennt da jemand etwas?
apogis
User
Beiträge: 5
Registriert: Dienstag 25. Mai 2010, 13:12

Die Daten, die man normalerweise bekommt sind auch in dieser Auflösung, aber das ist für die Fragestellung meines Chefs zu grob, deswegen haben wir diese besser aufgelösten Daten besorgt (wurden auch von der NASA erstellt & soweit ich weiß beim testen eines neuen Rechners (oder so) erzeugt & jetzt schon für ein oder zwei Paper weitergegeben). Bei einem der Autoren hab ich auch schon nach Hinweisen für die Umwandlung gefragt, die müssen das ja auch irgendwie gemacht haben, allerdings hab ich noch keine Antwort & versuche alle Möglichkeiten auszuschöpfen.
In der Headerdateisteht, dass die "Data file dimensions" 15670 und 19862 sind, was der Dateigröße durch 4 entspricht, ich weiß leider auch nicht warum es nicht die Hälfte ist bzw. was das dann zu bedeuten hat.
An welcher Stelle im Script müsste ich die row lenght den eintragen? Wenn ich jetzt mal davon ausgehe, dass die 15670 ist?
Ich hab auch immernoch nicht herausgefunden, wie man so ein script speichern kann... in python25, falls das einen Unterschied macht.

vielen Dank aber auf jeden Fall für deine Antwort! ich hab hier überhaupt niemanden mit dem ich mich irgendwie austauschen könnte...
apogis
User
Beiträge: 5
Registriert: Dienstag 25. Mai 2010, 13:12

Könnte es sein, dass die Daten dann als 4 Byte Integer vorliegen?
frabron
User
Beiträge: 306
Registriert: Dienstag 31. März 2009, 14:36

Naja, wenn die Dimension 15760*19862 ist, ist die Kachel halt nicht quadratisch, wie sonst, sondern rechteckig :D

Das Snippet von stackoverflow geht von quadratischen Kacheln aus, das wirst du so nicht verwenden können. Da bei dir der Faktor 4 anstatt 2 ist, gehe ich mal davon aus, dass die Daten nicht als Integer, 16bit vorliegen, sondern evt. 32bit, oder float. Aber dazu kenne ich mich wirklich nicht genug mit Binärformaten aus, um das sagen zu können.

Das Snippet bei Stackoverflow iteriert eigentlich nur über die Zeilen und Spalten der Datei und liest dabei immer zwei bytes ein (glaube ich), um die Höhe an der Stelle zu ermitteln. Soweit ich das Format verstanden habe, ist das eine Matrix(x,y) mit den z-Werten an einer Stelle x,y, die du auspacken musst. Aber wie Blackjack schon andeutet, ist das wohl nicht der beste Ansatz. Die Datei ist ja immerhin recht gross.
The NASA SRTM data files are in Big-Endian format, so depending on what platform you are reading the data on, you may have to do a conversion from Big-Endian to Little-Endian.
Hier auch aufpassen http://en.wikipedia.org/wiki/Endianness


Nachtrag:
Schau dir mal das Grass GIS Skript zum import an. Das ist ein Bash-Skript, das die Datei in eine Datei umwandelt, die sich mittels gdal einlesen lässt. Beachte besonders NROWS, NCOLS

Code: Alles auswählen

LL_LATITUDE=`echo $TILE  | cut -b2-3`
LL_LONGITUDE=`echo $TILE | cut -b5-7`

#are we on the southern hemisphere? If yes, make LATITUDE negative.
NORTH=`echo $TILE  | sed 's+.hgt++g' | cut -b1`
if [ "$NORTH" = "S" ] ; then
   LL_LATITUDE=`echo $LL_LATITUDE | awk '{printf "%.10f", $1 * -1 }'`
fi

#are we west of Greenwich? If yes, make LONGITUDE negative.
EAST=`echo $TILE  | sed 's+.hgt++g' | cut -b4`
if [ "$EAST" = "W" ] ; then
   LL_LONGITUDE=`echo $LL_LONGITUDE | awk '{printf "%.10f", $1 * -1 }'`
fi

# Calculate Upper Left from Lower Left
ULXMAP=`echo $LL_LONGITUDE | awk '{printf "%.1f", $1}'`
# SRTM90 tile size is 1 deg:
ULYMAP=`echo $LL_LATITUDE  | awk '{printf "%.1f", $1 + 1.0}'`


if [ $GIS_FLAG_1 -ne 1 ] ; then
  echo "BYTEORDER M
LAYOUT BIL
NROWS 1201
NCOLS 1201
NBANDS 1
NBITS 16
BANDROWBYTES 2402
TOTALROWBYTES 2402
BANDGAPBYTES 0
PIXELTYPE SIGNEDINT
NODATA -32768
ULXMAP $ULXMAP
ULYMAP $ULYMAP
XDIM 0.000833333333333
YDIM 0.000833333333333"> "$TILE.hdr"
apogis
User
Beiträge: 5
Registriert: Dienstag 25. Mai 2010, 13:12

Oh man, es ist ja noch schlimmer als ich dachte :P
Dass das Snippet von quadratischen Kacheln ausgeht wusste ich nicht, hab mich schon so gefreut überhaupt was brauchbares gefunden zu haben...
Dann werde ich mich wohl mal mit 'numpy' auseinandersetzen...
Ich hab aber scheinbar auch nichtmal alle notwendigen Informationen, ob es integer oder float sind weiß ich zB nicht, habs auch nur aus wegen dem Betrag bei stackoverflow angenommen...
Neben dem Scripten versuche ich auch es mit SAGA-GIS zu öffnen, da gibt es direkt die Möglichkeit Binary Raw Data einzulesen, aber mit den Parametern die ich bis jetzt hab kommt da auch nichts raus, was einem Höhenmodell ähnelt...

zum Nachtrag:
Oh, danke für das Script! Ich schaus mal durch.
frabron
User
Beiträge: 306
Registriert: Dienstag 31. März 2009, 14:36

Wenn du GRASS installiert hast: Eine Kopie vom Skript anlegen und mit den Werten für NROWS, NCOLS, NBITS(32?), BANDROWBYTES(4faches von Breite?), TOTALROWBYTES(4faches von Breite?) herumspielen.
apogis
User
Beiträge: 5
Registriert: Dienstag 25. Mai 2010, 13:12

Oh wow, ich habe ein Höhenmodell! :D
Ich habe in SAGA-GIS Zeilen- und Spaltenangaben vertauscht, weil ich dachte die Daten wären von einem anderen Ausschnitt... & sie waren 4byte floating, falls jemand auch mal in die Verlegenheit kommt mit diesen Daten zu Arbeiten.

Vielen dank auf jeden Fall frabron & BlackJack :)
Antworten