Skip to content

Instantly share code, notes, and snippets.

@biomiker
Created May 10, 2017 23:42
Show Gist options
  • Save biomiker/32fe34e1fa1bb49ae1135ab6652f596d to your computer and use it in GitHub Desktop.
Save biomiker/32fe34e1fa1bb49ae1135ab6652f596d to your computer and use it in GitHub Desktop.
Perl script for downloading and extracting USGS elevation data
#!/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);
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment