Я хочу назначить конкретную позицию long / lat на карте для отметки из файлов данных SRTM3, но не знаю, как найти конкретное значение. Поэтому я хочу получить пример того, как я могу найти в N50E14.hgt возвышение до 50 ° 24'58.888 "северной широты, 14 ° 55'11.377" восточной долготы.

Формат данных

Я возьму это как небольшое упражнение в том, как запрограммировать считыватель данных. Посмотрите на документацию :

Данные SRTM распределяются по двум уровням: SRTM1 (для США, их территорий и владений), данные которых отбираются с интервалами в одну угловую секунду по широте и долготе, а SRTM3 (для всего мира) отбирается по трем угловым секундам.

Данные делятся на фрагменты широты и долготы по одному градусу в «географической» проекции, то есть в растровом представлении с равными интервалами широты и долготы без какой-либо проекции, но с простой манипуляцией и мозаикой.

Имена файлов относятся к широте и долготе нижнего левого угла элемента мозаичного изображения - например, у нижнего левого угла N37W105 37 градусов северной широты и 105 градусов западной долготы. Точнее, эти координаты относятся к геометрическому центру нижнего левого пикселя, который в случае данных SRTM3 будет иметь протяженность около 90 метров.

Файлы высоты имеют расширение .HGT и подписаны двухбайтовыми целыми числами. Байты имеют порядок старшего байта в Motorola, причем сначала идет старший байт, который может быть непосредственно прочитан такими системами, как Sun SPARC, Silicon Graphics и Macintosh, использующими процессоры Power PC. DEC Alpha, большинство компьютеров и компьютеров Macintosh, выпущенных после 2006 года, используют порядок Intel («little-endian»), поэтому может потребоваться некоторая замена байтов. Высота указана в метрах относительно геоида WGS84 / EGM96. Пустоты данных присваиваются значения -32768.

Как продолжить

Для вашей позиции, 50 ° 24'58.888 "N 14 ° 55'11.377" E, вы уже нашли правильную плитку, N50E14.hgt. Давайте выясним, какой пиксель вас интересует. Первая широта, 50 ° 24'58.888 "N:

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

секунды дуги. Разделенный на три и округленный до ближайшего целого числа дает строку сетки 500. То же самое вычисление для долготы приводит к столбцу сетки 1104.

В документации по быстрому старту отсутствует информация о том, как строки и столбцы организованы в файле, но в полной документации указано, что

Данные хранятся в главном порядке строк (все данные для строки 1, затем все данные для строки 2 и т. Д.)

Первая строка в файле, скорее всего, самая северная, т. Е. Если мы заинтересованы в строке 500 от нижнего края , нам нужно взглянуть на строку

1201 - 500 = 701

с самого начала, если файл . Наша ячейка сетки является числом

(1201 * 700) + 1104 = 841804

от начала файла (то есть пропустить 700 строк, а в 701-м взять образец 1104). Два байта на образец означают, что мы должны пропустить первые 1683606 байтов в файле и затем прочитать два байта, чтобы получить нашу ячейку сетки. Данные с прямым порядком байтов, что означает, что вам нужно поменять местами два байта, например, на платформах Intel.

Пример программы

Упрощенная программа на Python для получения нужных данных будет выглядеть следующим образом (см. Документацию по использованию модуля 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)

Обратите внимание, что эффективный поиск данных должен выглядеть немного сложнее (например, не открывая файл для каждого образца).

альтернативы

Вы также можете использовать программу, которая может читать файлы .hgt из коробки. Но это скучно.

GDAL может читать / записывать эти растровые форматы с помощью драйвера SRTMHGT . Это означает, что вы можете просматривать растр с помощью QGIS, ArcGIS или использовать утилиты GDAL, такие как gdallocationinfo, для получения значений из точки, например:

Конвертировать DMS в DD:

• Широта: 50 ° 24'58,888 "N = 50 + (24/60) + (58,888 / 3600) = 50,4163577778
• Длинный: 14 ° 55'11,377 "E = 14 + (55/60) + (11,377 / 3600) = 14,9198269444

Затем из оболочки используйте gdallocationinfo file.hgt -wgs84 long lat:

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

Высота 216 м.

Если вы используете QGIS, проверьте, установлен ли плагин Python «Инструмент выборки точек». Вы найдете его в -> Улучшения (Python) -> Анализ.

Выберите свой точечный слой из требуемых позиций, затем запустите PST, выберите hgt (или любой другой растровый / полигонный файл) и выберите новую форму точки для вывода.

Это все :-)

  Chris

Ответ Криса указывает, что это просто для выборки точек из слоя в QGIS.

Однако, поскольку ваш ответ на мой комментарий уточняет, что вы пишете свою собственную программу для чтения значений высот из .hgtфайлов, по-другому взгляните на Quickstart PDF в документах SRTM. Это объясняет, как хранятся данные высот. Подвести итоги:

• Файлы SRTM3 содержат последовательность целочисленных значений с прямым порядком байтов.
• Каждое целочисленное значение представляет высоту «в метрах по отношению к геоиду WGS84 / EGM96», за исключением значений -32768, которые указывают пиксели без данных.
• Существует 1201 строка из 1201 выборок, поэтому в целом должно быть 1442401 целочисленных значений.

Вы говорите, что можете конвертировать между координатами lon / lat и пикселями, поэтому для получения высоты необходимо прочитать целочисленное значение из соответствующего смещения в файле. Учитывая координаты пикселей xи yотносительно верхнего левого угла сцены, это в основном offset = (y * 1201) + x. Пиксель 0,0- это первое целое число в файле, а пиксель 1200,1200- последнее целое число в файле.