Desejo atribuir uma posição longa / lat específica em um mapa à elevação dos arquivos de dados SRTM3, mas não tenho idéia de como encontrar o valor específico. Então, eu quero um exemplo de como posso encontrar na elevação N50E14.hgt a 50 ° 24'58.888 "N, 14 ° 55'11.377" E.

Formato de dados

Vou usar isso como um pequeno exercício de como programar um leitor de dados. Dê uma olhada na documentação :

Os dados SRTM são distribuídos em dois níveis: SRTM1 (para os EUA e seus territórios e posses) com dados amostrados em intervalos de um segundo de arco em latitude e longitude, e SRTM3 (para o mundo) amostrados em três segundos de arco.

Os dados são divididos em blocos de latitude e longitude de um por um grau na projeção "geográfica", ou seja, uma apresentação raster com intervalos iguais de latitude e longitude em nenhuma projeção, mas fáceis de manipular e mosaico.

Os nomes de arquivo se referem à latitude e longitude do canto inferior esquerdo do bloco - por exemplo, o N37W105 tem seu canto inferior esquerdo a 37 graus de latitude norte e 105 graus de longitude oeste. Para ser mais exato, essas coordenadas se referem ao centro geométrico do pixel inferior esquerdo, que no caso dos dados do SRTM3 terá cerca de 90 metros de extensão.

Os arquivos de altura têm a extensão .HGT e são assinados com números inteiros de dois bytes. Os bytes estão na ordem "big endian" da Motorola, com o byte mais significativo primeiro, diretamente legível por sistemas como os computadores Sun SPARC, Silicon Graphics e Macintosh usando processadores Power PC. DEC Alpha, a maioria dos computadores PC e Macintosh criados após 2006 usa a ordem Intel ("little-endian"), portanto, pode ser necessária alguma troca de bytes. As alturas estão em metros referenciados ao geóide WGS84 / EGM96. Vazios de dados recebem o valor -32768.

Como proceder

Para sua posição, 50 ° 24'58.888 "N 14 ° 55'11.377" E, você já encontrou o ladrilho correto, N50E14.hgt. Vamos descobrir em qual pixel você está interessado. Primeira latitude, 50 ° 24'58.888 "N:

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

segundos de arco. Dividido por três e arredondado para o número inteiro mais próximo, fornece uma linha de grade de 500. O mesmo cálculo para longitude resulta na coluna de grade 1104.

A documentação de início rápido carece de informações sobre como as linhas e colunas são organizadas no arquivo, mas na documentação completa , afirma-se que

Os dados são armazenados na ordem principal da linha (todos os dados da linha 1, seguidos por todos os dados da linha 2, etc.)

A primeira linha do arquivo é provavelmente a mais setentrional, ou seja, se estivermos interessados ​​na linha 500 a partir da borda inferior , na verdade temos que olhar para a linha

1201 - 500 = 701

desde o início, se o arquivo . Nossa célula da grade é o número

(1201 * 700) + 1104 = 841804

desde o início do arquivo (pule 700 linhas e, na 701ª, pegue a amostra 1104). Dois bytes por amostra significa que precisamos pular os primeiros 1683606 bytes no arquivo e depois ler dois bytes para obter nossa célula da grade. Os dados são big-endian, o que significa que você precisa trocar os dois bytes em, por exemplo, plataformas Intel.

Programa de exemplo

Um programa Python simplista para recuperar os dados corretos seria assim (consulte os documentos para uso do módulo struct):

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)

Observe que a recuperação eficiente de dados teria que parecer um pouco mais sofisticada (por exemplo, não abrir o arquivo para cada amostra).

Alternativas

Você também pode usar um programa que pode ler os arquivos .hgt imediatamente. Mas isso é chato.

O GDAL pode ler / gravar esses formatos de varredura com o driver SRTMHGT . Isso significa que você pode visualizar a varredura com QGIS, ArcGIS ou usar utilitários GDAL como gdallocationinfo para obter valores a partir de um ponto, por exemplo:

Converter DMS em DD:

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

Em seguida, a partir de um shell, use gdallocationinfo file.hgt -wgs84 long lat:

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

A altitude é de 216 m.

Se você usa o QGIS, verifique se o plug-in python "Point Sampling Tool" está instalado. Você o encontrará em -> Aprimoramentos (Python) -> Analisar.

Selecione a camada de pontos das posições necessárias, inicie o PST, escolha o hgt (ou qualquer arquivo de varredura / polígono) e escolha uma nova forma de ponto para a saída.

Isso é tudo :-)

  Chris

A resposta de Chris indica que é fácil amostrar pontos de uma camada no QGIS.

No entanto, como sua resposta ao meu comentário esclarece que você está escrevendo seu próprio programa para ler valores de elevação dos .hgtarquivos, verifique novamente o PDF de início rápido nos documentos do SRTM. Explica como os dados de elevação são armazenados. Para resumir:

• Os arquivos SRTM3 contêm uma sequência de valores inteiros big endian.
• Cada valor inteiro representa uma elevação "em metros referenciados ao geóide WGS84 / EGM96", exceto os valores de -32768, que indicam pixels sem dados.
• Existem 1201 linhas de 1201 amostras, portanto, deve haver 1442401 valores inteiros por completo.

Você diz que pode converter entre coordenadas lon / lat e pixels, portanto, obter a elevação é uma questão de ler o valor inteiro do deslocamento apropriado no arquivo. Dadas as coordenadas de pixel xe em yrelação ao canto superior esquerdo da cena, isso é basicamente offset = (y * 1201) + x. Pixel 0,0é o primeiro número inteiro no arquivo e pixel 1200,1200é o último número inteiro no arquivo.