0

I have downloaded a large set of GridFloat (.flt, .hdr) DEM files from USGS NED (1") in order to implement my own elevation service on my website. I would like to be able to look up an elevation from this fileset, given latitude and longitude as inputs. I use Perl for my website development. The files have a conventional naming scheme, and I am able to get the appropriate tile filename using the lat/lng. Howevever, accessing the internals of the file is where I'm having an issue.

I know the file is in a fairly straightforward format (.flt, apparently called "Gridfloat"), but I could use some help figuring out the magic numbers for calculating where in the file I need to seek to for a given lat/lng, and how to handle byte order and so on so that I end up with an elevation. From what I understand, apparently row ordering can be an issue, as well as byte ordering. I am looking for a recipe that does not involve use of any third party libraries such as GDAL, which I think are overly complicated and slow for what I want to do. I think it should be possible to just open the file, seek to a position based on some calculation, read some bytes and then unpack them into the correct byte order. Here is an example .hdr file that accompanies floatn48w097_1.flt, I think it has the necessary info. There are a bunch of other files that come with the .zip, including .prj, but I believe those are for a commercial program like ArcInfo. I think everything I need should be in the following .hdr file.

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

What I'm really hoping for is a formula for calculating the row and column from the lat/lng, then another formula for translating the row/column into a position for seek, how many bytes to read, and how to convert those raw bytes into an integer (or whatever it is these files contain). I feel that this could be a very fast operation, without all the overhead involved with the larger libraries which seem to be focused on doing a lot of stuff that I don't need.

I don't need Perl code, just pseudocode showing the calculations for row/col offsets etc would be more than enough. I believe the files are binary format, a straightforward grid of 4-byte numbers. The file example that goes with the .hdr file above has a size of 52186176, and when you multiply the ncols by nrows (from the .hdr), you get 13046544. which divides nicely into the file size by 4. So I assume it's just a matter of getting the right formula for row/col based on lat/lng, and then getting the bytes swizzled into the right order. I've just not done this much.

I found some reference to the Gridfloat format here: coolutils.com/formats/flt so apparently the file consists of a grid of 64-bit floating point values.

Thanks!

Neil Gunton
  • 51
  • 1
  • 7
  • I'm just guessing that most of the Perl experts on this site aren't familiar with the internal format of a Gridfloat file. Is it plain text? Binary? Can you edit your question to show a portion of a file? – ThisSuitIsBlackNot Aug 21 '14 at 22:33
  • Hi, thanks for responding. I don't need Perl code, just pseudocode showing the calculations for row/col offsets etc would be more than enough. I believe the files are binary format, a straightforward grid of 4-byte numbers. The file example that goes with the .hdr file above has a size of 52186176, and when you multiply the ncols by nrows (from the .hdr), you get 13046544. which divides nicely into the file size by 4. So I assume it's just a matter of getting the right formula for row/col based on lat/lng, and then getting the bytes swizzled into the right order. I've just not done this much. – Neil Gunton Aug 21 '14 at 22:39
  • I found some reference to the Gridfloat format here: http://www.coolutils.com/formats/flt so apparently the file consists of a grid of 64-bit floating point values. – Neil Gunton Aug 21 '14 at 22:41
  • I would edit those details into your question. Take a look at [`pack`](http://perldoc.perl.org/functions/pack.html) and [`unpack`](http://perldoc.perl.org/functions/unpack.html). Perhaps something like `my @vals = do { local $/; open my $fh, '<', $file or die $!; unpack 'V*', <$fh> };`, which would give you an array of values to index into or iterate through. – ThisSuitIsBlackNot Aug 21 '14 at 22:53
  • Although frankly, I would almost always prefer to reuse existing code than reinvent the wheel, even if it comes with a lot of features I don't need. You may decide you want those features in the future, and code that's used by many people tends to be better tested. – ThisSuitIsBlackNot Aug 21 '14 at 22:57
  • Thanks, I have added info to the main question as you suggested. I do not want to read the entire file in each time, that would be very inefficient for what I'm doing (looking up thousands of points in real time to construct an elevation profile for a map route). It should be possible with just a simple file lookup, it's just that the row/col ordering and byte ordering is tying my brain in knots. I'll probably be able to work it out, was just wondering if anyone already had done so. Thanks again. – Neil Gunton Aug 21 '14 at 23:01
  • I did find an example of what I'm after here, someone posted some PHP code for reading an elevation from an SRTM format file (.hgt): https://groups.google.com/forum/#!topic/geonames/935qIjOAq_c That is slightly different, not sure if the byte ordering etc are similar or not, but it might be a starting point... – Neil Gunton Aug 21 '14 at 23:02
  • I think you misunderstood; you would only `unpack` the file into an array once. The array would be large (*at least* 400 MB on a 32-bit system for your sample file, based on a back-of-the-envelope calculation), but you could easily index into it, e.g. `my $elev = $array[42];` Once you have the actual elevations in an array like this, you don't have to worry about endianness (you already did that with `unpack`); you just have to calculate the offset into the array based on the cell size and coordinates of one of the corners. – ThisSuitIsBlackNot Aug 22 '14 at 00:55
  • Of course, since you know the size of each field, you could also seek directly to that point in the file, read 4 bytes, and use `unpack` on that, which would require significantly less memory. – ThisSuitIsBlackNot Aug 22 '14 at 01:06
  • This file is just one tile out of over 1,700. Unzipped, they take on the order of tens or maybe over 100 GB (not sure yet, still downloading the rest). It wouldn't be efficient to have to load the entire file into memory. I think I might have a solution, I'll post it as an answer below - but thanks for your help, much appreciated! – Neil Gunton Aug 22 '14 at 01:20

2 Answers2

1

Ok, I think I have an answer. The following is Perl routine, which seems to give back reasonable looking elevation values when tested with the USGS NED1 .flt files. The script takes latitude and longitude as command line arguments, looks up the file and indexes into the grid.

#!/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";

Hope this might be useful to anyone else trying to do the same kind of thing... I've tested this method now and it seems to produce good looking elevation profiles which are smoother than those from the 3" SRTM data.

Neil Gunton
  • 51
  • 1
  • 7
0

Neil put me on the right track but I think there's a few problems with his original answer. I've added some fixes and improvements including on-the-fly download of the needed tile from the 1/3 arc second (10 meter) dataset, proper parsing of the header file, and what I believe is corrected indexing.

This is still mostly illustrative and should be improved before production use, particularly, hanging on to the header information and the file handle for repeated queries.

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

halfer
  • 19,824
  • 17
  • 99
  • 186
biomiker
  • 3,108
  • 1
  • 29
  • 26