Perl script for downloading and extracting USGS elevation data

gridFloatElevation.pl

#!/usr/bin/env perl

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

my $debug = 0;

my $dataDir = "/Users/miker/elevation/data";
my $usgsUrl = "https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/13/GridFloat";
# Note, FTP = "ftp://rockyftp.cr.usgs.gov/vdelivery/Datasets/Staged/NED/13/GridFloat"

sub download_tile {
    my $tilename = $_[0];

    # Note, using older dataset for simplicity
    my $filename = "$tilename.zip";
    
    my $cmd = "curl -f -o $dataDir/$tilename.zip $usgsUrl/$filename && unzip -d $dataDir/$tilename $dataDir/$tilename.zip";
    print "Downloading: $usgsUrl/$filename\n";
    print "To:          $dataDir/$tilename\n";
    system($cmd);
}

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

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

    my $lat_tilenum = abs(ceil ($lat));
    my $lng_tilenum = abs(floor($lng));

    my $tilename =  $lat_letter . sprintf('%02d', $lat_tilenum) . $lng_letter . sprintf('%03d',$lng_tilenum);

    print("$lat, $lng => $tilename\n") if $debug;
    
    my $path = "$dataDir/$tilename/float${tilename}_13.flt";
    
    if (! -e $path) {
			# Attempt download
			print ("$lat, $lng => $tilename");
			download_tile($tilename);
    }
    
    if (!-e $path) {
			print "Tile not found: $path\n";
			return undef;
    }

    # parse header file
    my $header_path = $path;
    $header_path =~ s/\.flt$/\.hdr/;

    warn("Parsing $header_path\n") if $debug;
    open(my $header_file, $header_path) || die "Can't read $header_path";
    my %header;
    while(<$header_file>) {
			my ($key, $val) = split;
			print "$key => $val\n" if $debug;
			$header{$key} = $val;
    }

    my $ncols =  $header{ncols} || die "invalid header (ncols)";
    my $nrows =  $header{nrows} || die "invalid header (nrows)";
    my $xllcorner  =   $header{xllcorner} || die "invalid header (xllcorner)";
    my $yllcorner   =  $header{yllcorner} || die "invalid header (yllcorner)";
    my $cellsize  =    $header{cellsize} || die "invalid header (cellsize)";
    my $NODATA_value =   $header{NODATA_value} || die "invalid header (NODATA_value)";
    my $byteorder =    $header{byteorder} || die "invalid header (byteorder)";

    # TODO: Handle $byteorder = MSBFIRST.  
    die "Unexpected byteorder $byteorder\n" if $byteorder ne 'LSBFIRST';
    
    my $bytesPerDatum = 4;
    
    my $height = $nrows * $cellsize;
    my $width  = $ncols * $cellsize;
    
    print ("Tile size = $width x $height\n") if $debug;
    
    my $x0 = $xllcorner;
    my $y0 = $yllcorner + ( $nrows * $cellsize );
    
    print ("Origin = $y0, $x0\n") if $debug;
    
    print ("Expected size = ", ($nrows * $ncols * $bytesPerDatum), "\n") if $debug;
	
    open(FILE, "<$path");

    # TODO:  Make sure this works in all hemispheres.
    my $row = round(( $y0 - $lat ) / $cellsize);
    my $col = round(( $lng - $x0 ) / $cellsize);

    # Is this sound edge case handling or hack that misses the point?
    $row = $nrows - 1 if $row >= $nrows;
    $col = $ncols - 1 if $col >= $ncols;

    if ($row < 0 || $row >= $nrows) {
			print("Error: invalid row $row\n");
			return undef;
    }

    if ($col < 0 || $col >= $ncols) {
			print("Error: invalid col $col\n");
			return undef;
    }

    my $pos = ($row * $bytesPerDatum * $ncols) + ($col * $bytesPerDatum);

    print("Looking for data at $row, $col => $pos\n") if $debug;

    seek (FILE, $pos, SEEK_SET);

    my $buffer;
    read (FILE, $buffer, $bytesPerDatum);  # 4 = 32 bit or 5 = 64 bit?

    close (FILE);

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

    if ($elevation == $NODATA_value)
    {
        return undef;
    }

    return $elevation;
}

if ($ARGV[0] && $ARGV[1]) {
    my $lat = $ARGV[0];
    my $lng = $ARGV[1];

    my $e = get_elevation($lat, $lng);
    print $e, "\n" if defined $e;
}
else {
    print "Usage:  $0 lat lng\n";
    exit(1);
}