Как прочитать высоту из файла GridFloat USGS NED DEM в Perl

Я загрузил большой набор файлов DEM GridFloat (.flt, .hdr) из USGS NED (1"), чтобы реализовать собственный сервис высот на своем веб-сайте. Я хотел бы иметь возможность искать высотные отметки из этого набора файлов. , учитывая широту и долготу в качестве входных данных. Я использую Perl для разработки своего веб-сайта. Файлы имеют обычную схему именования, и я могу получить соответствующее имя файла плитки, используя широту / долготу. Однако доступ к внутренней части файла находится там, где У меня проблема.

Я знаю, что файл имеет довольно простой формат (.flt, по-видимому, называется «Gridfloat»), но мне не помешала бы помощь в выяснении магических чисел для расчета, где в файле мне нужно искать для заданной широты/долготы, и как обрабатывать порядок байтов и т. д., чтобы я получил высоту. Насколько я понимаю, очевидно, что порядок строк может быть проблемой, а также порядок байтов. Я ищу рецепт, который не предполагает использования каких-либо сторонних библиотек, таких как GDAL, которые, как мне кажется, слишком сложны и медленны для того, что я хочу сделать. Я думаю, что должна быть возможность просто открыть файл, найти позицию на основе некоторых вычислений, прочитать несколько байтов, а затем распаковать их в правильном порядке байтов. Вот пример файла .hdr, который сопровождает floatn48w097_1.flt, я думаю, в нем есть необходимая информация. Вместе с .zip поставляется множество других файлов, включая .prj, но я полагаю, что они предназначены для коммерческой программы, такой как ArcInfo. Я думаю, что все, что мне нужно, должно быть в следующем файле .hdr.

ncols         3612
nrows         3612
xllcorner     -97.00166666667
yllcorner     46.99833333333
cellsize      0.000277777777778
NODATA_value  -9999
byteorder     LSBFIRST

На что я действительно надеюсь, так это на формулу для вычисления строки и столбца из широты/длины, а затем еще одну формулу для перевода строки/столбца в позицию для поиска, сколько байтов нужно прочитать и как преобразовать эти необработанные байты. в целое число (или что-то еще, что эти файлы содержат). Я чувствую, что это может быть очень быстрая операция без всех накладных расходов, связанных с большими библиотеками, которые, похоже, сосредоточены на выполнении множества вещей, которые мне не нужны.

Мне не нужен код Perl, просто псевдокод, показывающий расчеты для смещений строк/столбцов и т.д., был бы более чем достаточным. Я считаю, что файлы имеют двоичный формат, простую сетку 4-байтовых чисел. Пример файла, который идет с файлом .hdr выше, имеет размер 52186176, и когда вы умножаете ncols на nrows (из .hdr), вы получаете 13046544., что хорошо делится на размер файла на 4. Поэтому я предполагаю, что это просто вопрос получения правильной формулы для строки / столбца на основе широты / долготы, а затем получения байтов в правильном порядке. Я просто не так много сделал.

Я нашел ссылку на формат Gridfloat здесь: coolutils.com/formats/flt, так что, очевидно, файл состоит из сетки 64-битных значений с плавающей запятой.

Спасибо!


perl elevation
person Neil Gunton    schedule 21.08.2014    source источник
comment
Я просто предполагаю, что большинство экспертов по Perl на этом сайте не знакомы с внутренним форматом файла Gridfloat. Это обычный текст? Бинарный? Можете ли вы отредактировать свой вопрос, чтобы показать часть файла?   -  person ThisSuitIsBlackNot    schedule 22.08.2014
comment
Привет, спасибо за ответ. Мне не нужен код Perl, просто псевдокод, показывающий расчеты для смещений строк/столбцов и т.д., был бы более чем достаточным. Я считаю, что файлы имеют двоичный формат, простую сетку 4-байтовых чисел. Пример файла, который идет с файлом .hdr выше, имеет размер 52186176, и когда вы умножаете ncols на nrows (из .hdr), вы получаете 13046544., что хорошо делится на размер файла на 4. Поэтому я предполагаю, что это просто вопрос получения правильной формулы для строки / столбца на основе широты / долготы, а затем получения байтов в правильном порядке. Я просто не так много сделал.   -  person Neil Gunton    schedule 22.08.2014
comment
Я нашел ссылку на формат Gridfloat здесь: coolutils.com/formats/flt файл состоит из сетки 64-битных значений с плавающей запятой.   -  person Neil Gunton    schedule 22.08.2014
comment
Я бы отредактировал эти детали в вашем вопросе. Взгляните на pack и unpack. Возможно, что-то вроде my @vals = do { local $/; open my $fh, '<', $file or die $!; unpack 'V*', <$fh> };, что даст вам массив значений для индексации или повторения.   -  person ThisSuitIsBlackNot    schedule 22.08.2014
comment
Хотя, честно говоря, я почти всегда предпочел бы повторно использовать существующий код, чем изобретать велосипед, даже если он содержит множество функций, которые мне не нужны. Вы можете решить, что вам нужны эти функции в будущем, и код, который используется многими людьми, как правило, лучше тестируется.   -  person ThisSuitIsBlackNot    schedule 22.08.2014
comment
Спасибо, я добавил информацию к основному вопросу, как вы предложили. Я не хочу каждый раз читать весь файл, это было бы очень неэффективно для того, что я делаю (просматривая тысячи точек в реальном времени, чтобы построить профиль высот для маршрута карты). Это должно быть возможно с помощью простого поиска файлов, просто порядок строк / столбцов и порядок байтов связывают мой мозг узлами. Я, вероятно, смогу решить это, просто интересно, если кто-нибудь уже сделал это. Спасибо еще раз.   -  person Neil Gunton    schedule 22.08.2014
comment
Я нашел здесь пример того, что мне нужно, кто-то опубликовал PHP-код для чтения высоты из файла формата SRTM (.hgt): groups.google.com/forum/#!topic/geonames/935qIjOAq_c Это немного отличается, не уверен, что порядок байтов и т. д. аналогичен или нет, но это может быть отправной точкой...   -  person Neil Gunton    schedule 22.08.2014
comment
Я думаю, вы неправильно поняли; вы бы только unpack файл в массив один раз. Массив будет большим (не менее 400 МБ в 32-разрядной системе для файла примера, исходя из предварительных расчетов), но вы можете легко его проиндексировать, например. my $elev = $array[42]; Когда у вас есть фактические высоты в таком массиве, вам не нужно беспокоиться о порядке следования байтов (вы уже сделали это с unpack); вам просто нужно вычислить смещение в массиве на основе размера ячейки и координат одного из углов.   -  person ThisSuitIsBlackNot    schedule 22.08.2014
comment
Конечно, поскольку вы знаете размер каждого поля, вы также можете перейти непосредственно к этой точке в файле, прочитать 4 байта и использовать для этого unpack, что потребует значительно меньше памяти.   -  person ThisSuitIsBlackNot    schedule 22.08.2014
comment
Этот файл — всего лишь одна плитка из более чем 1700. В разархивированном виде они занимают порядка десятков или, может быть, более 100 ГБ (пока не уверен, все еще скачиваю остальное). Было бы неэффективно загружать весь файл в память. Я думаю, что у меня может быть решение, я опубликую его как ответ ниже, но спасибо за вашу помощь, очень признателен!   -  person Neil Gunton    schedule 22.08.2014

Ответы (2)


arrow_upward
1
arrow_downward

Хорошо, думаю, у меня есть ответ. Ниже приведена подпрограмма Perl, которая, кажется, возвращает разумно выглядящие значения высоты при тестировании с файлами USGS NED1 .flt. Сценарий принимает широту и долготу в качестве аргументов командной строки, ищет файл и индексирует его в сетке.

#!/usr/bin/perl

use strict;
use POSIX;
use Math::Round;

sub get_elevation
{
    my ($lat, $lng) = @_;

    my $lat_degree = ceil ($lat);
    my $lng_degree = floor ($lng);

    my $lat_letter = ($lat >= 0) ? 'n' : 's';
    my $lng_letter = ($lng >= 0) ? 'e' : 'w';

    my $lng_tilenum = abs($lng_degree);
    my $lat_tilenum = abs($lat_degree);

    my $tilename =  $lat_letter . sprintf('%02d', $lat_tilenum) . $lng_letter . sprintf('%03d',$lng_tilenum);
    my $path = "/data/elevation/ned1/$tilename/float${tilename}_1.flt";

    print "path = $path\n";

    die "No such file" if (!-e($path));

    my ($lat_fraction, $lat_integral) = modf (abs($lat));
    my $row = floor ((1 - $lat_fraction) * 3600);

    my ($lng_fraction, $lng_integral) = modf (abs($lng));
    my $col = floor ((1 - $lng_fraction) * 3600);

    open(FILE, "<$path");

    my $pos = (3612 * 4 * 6) + (3612 * 4 * $row) + (4 * 6) + ($col * 4);

    seek (FILE, $pos, SEEK_SET);

    my $buffer;
    read (FILE, $buffer, 4);

    close (FILE);

    my ($elevation) = unpack('f', $buffer);

    if ($elevation == -9999)
    {
        return 'undefined';
    }

    return $elevation;
}

my $lat = $ARGV[0];
my $lng = $ARGV[1];
my $elevation = get_elevation ($lat, $lng);

print "Elevation for ($lat, $lng) = $elevation meters (", $elevation * 3.28084, " feet)\n";

Надеюсь, это может быть полезно для всех, кто пытается сделать то же самое... Сейчас я протестировал этот метод, и, похоже, он дает хорошо выглядящие профили высот, которые более гладкие, чем профили из 3-дюймовых данных SRTM.

person Neil Gunton    schedule 22.08.2014

arrow_upward
0
arrow_downward

Нил поставил меня на правильный путь, но я думаю, что с его первоначальным ответом есть несколько проблем. Я добавил несколько исправлений и улучшений, включая загрузку на лету нужного тайла из набора данных 1/3 угловой секунды (10 метров), правильный разбор файла заголовка и то, что я считаю исправленной индексацией.

Это по-прежнему в основном иллюстративно и должно быть улучшено перед использованием в производстве, в частности, в зависимости от информации заголовка и дескриптора файла для повторяющихся запросов.

https://gist.github.com/biomiker/32fe34e1fa1bb49ae1135ab6652f596d

person biomiker    schedule 10.05.2017
comment
Спасибо, полумер. Больше не повторится. :) - person biomiker; 17.05.2017