如何在 Perl 中从 USGS NED DEM GridFloat 文件中读取高程答案

【问题标题】：How to read elevation from USGS NED DEM GridFloat file in Perl如何在 Perl 中从 USGS NED DEM GridFloat 文件中读取高程
【发布时间】：2014-10-15 16:55:29
【问题描述】：

我从 USGS NED (1") 下载了大量 GridFloat (.flt, .hdr) DEM 文件，以便在我的网站上实现我自己的高程服务。我希望能够查找高程从这个文件集中，给定纬度和经度作为输入。我使用 Perl 进行网站开发。这些文件有一个传统的命名方案，我可以使用 lat/lng 获取适当的切片文件名。但是，访问文件是我遇到问题的地方。

我知道该文件是一种相当简单的格式（.flt，显然称为“Gridfloat”），但我可以使用一些帮助找出魔术数字来计算我需要在文件中寻找给定纬度的位置/lng，以及如何处理字节顺序等，以便我最终得到一个提升。据我了解，显然行排序以及字节排序可能是一个问题。我正在寻找一个不涉及使用任何第三方库（如 GDAL）的配方，我认为这对于我想做的事情来说过于复杂和缓慢。我认为应该可以只打开文件，根据一些计算寻找一个位置，读取一些字节，然后将它们解压缩为正确的字节顺序。这是一个伴随floatn48w097_1.flt的示例.hdr文件，我认为它有必要的信息。 .zip 附带了许多其他文件，包括 .prj，但我相信这些文件适用于 ArcInfo 等商业程序。我认为我需要的一切应该在下面的 .hdr 文件中。

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

我真正希望的是一个从 lat/lng 计算行和列的公式，然后是另一个将行/列转换为查找位置的公式，要读取多少字节，以及如何转换将这些原始字节转换为整数（或这些文件包含的任何内容）。我觉得这可能是一个非常快速的操作，而不会涉及到大型库所涉及的所有开销，这些库似乎专注于做很多我不需要的事情。

我不需要 Perl 代码，只需要显示行/列偏移等计算的伪代码就足够了。我相信这些文件是二进制格式，一个简单的 4 字节数字网格。上面 .hdr 文件附带的文件示例的大小为 52186176，当您将 ncols 乘以 nrows（来自 .hdr）时，您得到 13046544。它很好地将文件大小除以 4。所以我假设它是只需根据 lat/lng 获得正确的 row/col 公式，然后将字节调整为正确的顺序。我只是没做过这么多。

我在这里找到了对 Gridfloat 格式的一些参考：coolutils.com/formats/flt 显然该文件由 64 位浮点值网格组成。

谢谢！

【问题讨论】：

我只是猜测这个站点上的大多数 Perl 专家并不熟悉 Gridfloat 文件的内部格式。是纯文本吗？二进制？您可以编辑您的问题以显示文件的一部分吗？
您好，感谢您的回复。我不需要 Perl 代码，只需要显示行/列偏移量等计算的伪代码就足够了。我相信这些文件是二进制格式，一个简单的 4 字节数字网格。上面 .hdr 文件附带的文件示例的大小为 52186176，当您将 ncols 乘以 nrows（来自 .hdr）时，您得到 13046544。它很好地将文件大小除以 4。所以我假设它是只需根据 lat/lng 获得正确的 row/col 公式，然后将字节调整为正确的顺序。我只是没有做这么多。
我在这里找到了一些对 Gridfloat 格式的参考：coolutils.com/formats/flt 所以显然该文件由 64 位浮点值的网格组成。
我会将这些详细信息编辑到您的问题中。看看pack 和unpack。也许像my @vals = do { local $/; open my $fh, '<', $file or die $!; unpack 'V*', <$fh> }; 这样的东西会给你一个值数组来索引或迭代。
虽然坦率地说，我几乎总是更喜欢重用现有代码而不是重新发明轮子，即使它带有很多我不需要的功能。您可能会决定在未来需要这些功能，并且许多人使用的代码往往会经过更好的测试。

标签： perl elevation

【解决方案1】：

好的，我想我有答案了。以下是 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 数据更平滑的高程剖面。

【讨论】：

【解决方案2】：

尼尔让我走上了正轨，但我认为他的原始答案存在一些问题。我添加了一些修复和改进，包括从 1/3 角秒（10 米）数据集中即时下载所需的图块、正确解析头文件以及我认为更正的索引。

这仍然主要是说明性的，应该在生产使用之前进行改进，特别是保留标题信息和重复查询的文件句柄。

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

【讨论】：

谢谢halfer。不会再发生了。 :)