¿Extracción de elevación del archivo .HGT?

20

Quiero asignar una posición específica larga / lat en un mapa a la elevación de los archivos de datos SRTM3, pero no tengo idea de cómo encontrar el valor específico. Así que quiero un ejemplo de cómo puedo encontrar en N50E14.hgt elevación a 50 ° 24'58.888 "N, 14 ° 55'11.377" E.

MartinS
fuente
1
¿Que software estas usando? Hay algunas notas sobre el .hgtformato de archivo en la documentación de SRTM , pero una respuesta específica paso a paso depende del software que tenga disponible.
ungido el
1
No tengo software, soy programador de C # y estoy escribiendo mi propia aplicación. Puedo asignar long / lat a cada píxel y ahora quiero buscar la elevación a cada punto. El mejor formato de datos debería ser algo como el archivo CSV. Entonces en una fila puedo encontrar longitud, latitud, altitud. He buscado en la documentación de SRTM, pero aún no puedo imaginar cómo puedo proporcionar minería de datos en el archivo.
MartinS

Respuestas:

30

Formato de datos

Lo tomaré como un pequeño ejercicio sobre cómo programar un lector de datos. Echa un vistazo a la documentación :

Los datos SRTM se distribuyen en dos niveles: SRTM1 (para los EE. UU. Y sus territorios y posesiones) con datos muestreados a intervalos de un segundo de arco en latitud y longitud, y SRTM3 (para el mundo) muestreados a tres segundos de arco.

Los datos se dividen en mosaicos de latitud y longitud de uno por un grado en la proyección "geográfica", es decir, una presentación ráster con intervalos iguales de latitud y longitud en ninguna proyección, pero fácil de manipular y mosaico.

Los nombres de los archivos se refieren a la latitud y longitud de la esquina inferior izquierda del mosaico; por ejemplo, N37W105 tiene su esquina inferior izquierda a 37 grados de latitud norte y 105 grados de longitud oeste. Para ser más exactos, estas coordenadas se refieren al centro geométrico del píxel inferior izquierdo, que en el caso de los datos SRTM3 tendrá una extensión de aproximadamente 90 metros.

Los archivos de altura tienen la extensión .HGT y están firmados con enteros de dos bytes. Los bytes están en el orden "big-endian" de Motorola con el byte más significativo primero, directamente legible por sistemas como Sun SPARC, Silicon Graphics y computadoras Macintosh que utilizan procesadores Power PC. DEC Alpha, la mayoría de las PC y Macintosh construidas después de 2006 usan el orden Intel ("little-endian"), por lo que puede ser necesario un intercambio de bytes. Las alturas están en metros referenciados al geoide WGS84 / EGM96. A los vacíos de datos se les asigna el valor -32768.

Cómo proceder

Para su posición, 50 ° 24'58.888 "N 14 ° 55'11.377" E, ya encontró el mosaico correcto, N50E14.hgt. Veamos qué píxel le interesa. Primera latitud, 50 ° 24'58.888 "N:

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

segundos de arco Dividido por tres y redondeado al entero más cercano, se obtiene una fila de cuadrícula de 500. El mismo cálculo para la longitud da como resultado la columna de cuadrícula 1104.

La documentación de inicio rápido carece de información sobre cómo se organizan las filas y columnas en el archivo, pero en la documentación completa se afirma que

Los datos se almacenan en el orden principal de la fila (todos los datos de la fila 1, seguidos de todos los datos de la fila 2, etc.)

Es muy probable que la primera fila del archivo sea la más septentrional, es decir, si estamos interesados ​​en la fila 500 desde el borde inferior , en realidad tenemos que mirar la fila

1201 - 500 = 701

desde el principio si el archivo . Nuestra celda de cuadrícula es número

(1201 * 700) + 1104 = 841804

desde el inicio del archivo (es decir, omita 700 filas y, en la 701, tome la muestra 1104). Dos bytes por muestra significa que tenemos que omitir los primeros 1683606 bytes en el archivo y luego leer dos bytes para obtener nuestra celda de cuadrícula. Los datos son big-endian, lo que significa que debe intercambiar los dos bytes, por ejemplo, en plataformas Intel.

Programa de muestra

Un programa simplista de Python para recuperar los datos correctos se vería así (ver los documentos para el uso del módulo de estructura):

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)

Tenga en cuenta que la recuperación de datos eficiente debería verse un poco más sofisticada (por ejemplo, no abrir el archivo para todas y cada una de las muestras).

Alternativas

También puede usar un programa que pueda leer los archivos .hgt de fábrica. Pero eso es aburrido.

bhell
fuente
¡Golpéame y con más detalles para arrancar!
ungido el
Agradable explicar, te amo. Gracias por tu ayuda. Todos ustedes chicos.
MartinS
1
+1 Sí, las filas están ordenadas de norte a sur. Esto queda claro tan pronto como asigne uno de los archivos. Además, considere la obtención de alturas mediante interpolación bilineal entre los cuatro centros celulares que rodean la ubicación.
whuber
Gracias por esta información completa! Tengo una pregunta: cuando buscamos datos de elevación para 50 ° 24'58.888, ¿por qué restas la fila número 500 del borde inferior cuando las filas están ordenadas de norte a sur? ¡Gracias!
Georg
Si no me equivoco, creo que (j-1) será solo j. El valor j varía de 0 a 1200, por lo que no es necesario restar 1.
Malipivo
6

GDAL puede leer / escribir estos formatos ráster con el controlador SRTMHGT . Esto significa que puede ver el ráster con QGIS, ArcGIS o utilizar utilidades GDAL como gdallocationinfo para obtener valores de un punto, por ejemplo:

Convierta DMS a DD:

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

Luego, desde un shell, use gdallocationinfo file.hgt -wgs84 long lat:

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

La elevación es de 216 m.

Mike T
fuente
1
¿Qué hay de las ubicaciones en los lados sur u oeste?
Muhammet Ali Asan
2

Si usa QGIS, verifique si está instalado el complemento de Python "Point Sampling Tool". Lo encontrará en -> Mejoras (Python) -> Analizar.

Seleccione su capa de puntos de las posiciones requeridas, luego inicie el PST, elija el hgt (o cualquier archivo raster / polygone) y elija una nueva forma de punto para la salida.

Eso es todo :-)

  Chris
Chris Pallasch
fuente
0

La respuesta de Chris indica que es sencillo muestrear puntos de una capa en QGIS.

Sin embargo, dado que su respuesta a mi comentario aclara que está escribiendo su propio programa para leer valores de elevación de los .hgtarchivos, eche otro vistazo al PDF de inicio rápido en los documentos SRTM. Explica cómo se almacenan los datos de elevación. Para resumir:

  • Los archivos SRTM3 contienen una secuencia de valores enteros big-endian.
  • Cada valor entero representa una elevación "en metros referenciados al geoide WGS84 / EGM96", excepto los valores de -32768, que indican píxeles sin datos.
  • Hay 1201 líneas de 1201 muestras, por lo que debería haber 1442401 valores enteros en total.

Dice que puede convertir entre coordenadas lon / lat y píxeles, por lo que obtener la elevación es una cuestión de leer el valor entero del desplazamiento apropiado en el archivo. Dadas las coordenadas de píxeles xy en yrelación con la esquina superior izquierda de la escena, eso es básicamente offset = (y * 1201) + x. Pixel 0,0es el primer entero del archivo y pixel 1200,1200es el último entero del archivo.

ungido
fuente
1
Esto es correcto, pero le faltan algunos detalles cruciales proporcionados por la respuesta de bhell, incluido que las coordenadas están asociadas con los centros celulares . Así, por ejemplo, la esquina superior izquierda de N50E014.hgt se encuentra realmente en la longitud 13.99958 E, latitud 51.00042 N.
whuber