¿Extrae la elevación del archivo .HGT?

18

Quiero asignar una posición lat / long específica en un mapa a la elevación desde 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.

    
pregunta Martin Ševic 11.12.2012 - 21:04

4 respuestas

27

Formato de datos

Lo tomaré como un pequeño ejercicio sobre cómo programar un lector de datos. Eche 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 en un grado en proyección "geográfica", es decir, una presentación de trama con intervalos iguales de latitud y longitud en ninguna proyección, pero fácil de manipular y en 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   los bytes están en el orden "big-endian" de Motorola con el byte más significativo primero,   legibles directamente por sistemas como Sun SPARC, Silicon Graphics y computadoras Macintosh que usan procesadores Power PC. DEC Alpha, la mayoría de las PC y las computadoras Macintosh construidas después de 2006 usan el pedido de 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ó la pieza correcta, N50E14.hgt. Averigüemos en qué píxel estás interesado. 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 número entero más cercano, se obtiene una fila de cuadrícula de 500. El mismo cálculo para los resultados de longitud en 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 mayor 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 al norte, 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 701a toma de muestra 1104). Dos bytes por muestra significa que debemos 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 necesita intercambiar los dos bytes, por ejemplo. Plataformas Intel.

Programa de ejemplo

Un programa de Python simplista para recuperar los datos correctos se vería así (vea los documentos para su 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 tendría que parecer un poco más sofisticada (por ejemplo, no abrir el archivo para cada muestra).

Alternativas

También puedes usar un programa que pueda leer los archivos .hgt de la caja. Pero eso es aburrido.

    
respondido por el bhell 11.12.2012 - 22:46
5

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

Convertir DMS a DD:

  • Lat: 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.

    
respondido por el Mike T 16.02.2015 - 05:08
1

Si utiliza QGIS, verifique si el complemento de python "Herramienta de muestreo de puntos" está instalado. Lo encontrarás en - > Mejoras (Python) - > Analizar.

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

Eso es todo :-)

  Chris
    
respondido por el Chris Pallasch 11.12.2012 - 22:06
0

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

Sin embargo, como su respuesta a mi comentario aclara que está escribiendo su propio programa para leer los valores de elevación de los archivos .hgt , vuelva a consultar Quickstart PDF 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 x y y relativas a la esquina superior izquierda de la escena, eso es básicamente offset = (y * 1201) + x . Pixel 0,0 es el primer entero en el archivo y pixel 1200,1200 es el último entero en el archivo.

    
respondido por el anoved 11.12.2012 - 22:52

Lea otras preguntas en las etiquetas