File Coverage

blib/lib/Geo/Elevation/HGT.pm

Criterion	Covered	Total	%
statement	26	160	16.2
branch	0	76	0.0
condition	0	24	0.0
subroutine	9	18	50.0
pod	3	3	100.0
total	38	281	13.5

package Geo::Elevation::HGT;

use 5.010;
use strict;
use warnings;
use POSIX ();
use Carp;
# set the version for version checking
our $VERSION     = '0.04';
# file-private lexicals
my $grid_size; # .hgt grid size = 3601x3601 for 1-minute DEMs or 1201x1201 for 3-minute DEMs
my @DEMnames;
my @default_DEMs;
my @want_DEMnames;
my $status_descr;
my $url;
my $folder;
my $cache_folder;
my $debug;
my $switch;
my $cache;
my $fail;

sub new {
  my ($class, %params) = @_;
  %params = (
    url => "https://elevation-tiles-prod.s3.amazonaws.com/skadi",    # info at https://registry.opendata.aws/terrain-tiles/
    folder => "https://elevation-tiles-prod.s3.amazonaws.com/skadi",
    cache_folder => "",
    debug => 0,
    %params
  );
  my $self = {};
  while ( my($key,$value) = each %params ) {
    $self->{$key} = $value;
  }
  bless $self, $class;
  return $self;
}

sub get_elevation_batch_hgt {
  my ($self, $batch_latlon) = @_;
  my @elegeh;
  for my $latlon ( @$batch_latlon ) {
    my ($lat, $lon) = @$latlon;
    push (@elegeh, $self->get_elevation_hgt ($lat, $lon));
  }
  return \@elegeh;
}

sub get_elevation_hgt {
  my ($self, $lat, $lon) = @_;
  $folder = $self->{folder};
  $url = $self->{url};
  $cache_folder = $self->{cache_folder};
  $debug = $self->{debug};
  $status_descr = "Memory";
  say STDERR "get_elevation_hgt" if $debug;
  say STDERR "  Parameters: folder=>'$folder', url=>'$url', cache_folder=>'$cache_folder', debug=>$debug" if $debug;
  say STDERR "  Input lat lon: $lat  $lon" if $debug;
  my $flat = POSIX::floor $lat;
  my $flon = POSIX::floor $lon;
  my $ns = $flat < 0 ? "S" : "N";
  my $ew = $flon < 0 ? "W" : "E";
  my $lt = sprintf ("%02s", abs($flat));
  my $ln = sprintf ("%03s", abs($flon));
  my $DEMname = "$ns$lt$ew$ln";
  say STDERR "  Tile lat lon: $flat  $flon" if $debug;
  # read DEM unless already defined
  say STDERR "  Using DEM in memory: '$flat  $flon'" if ($debug and defined $self->{DEMs}{$DEMname}{DEM});
  $self->{DEMs}{$DEMname}{DEM} //= $self->_readDEM($DEMname);  # //= Logical Defined-Or Assignment Operator
  my $dem = $self->{DEMs}{$DEMname}{DEM};
  say STDERR "  Status: $status_descr" if $debug;
  $self->{lat} = $lat;
  $self->{lon} = $lon;
  $self->{status_descr} = $status_descr;
  $self->{switch} = $switch;
  $self->{cache} = $cache;
  $self->{fail} = $fail;
  $self->{DEMname} = $DEMname;
  if (ref($dem) eq "") {
    say STDERR "  No data in DEM: '$flat  $flon' returning elevation 0" if $debug;
    $self->{grid_size} = 0;
    $self->{elevation} = 0;
    return $self->{elevation};
  }
  $grid_size = sqrt (length ($$dem)/2);    # grid size of DEM
  unless ($grid_size == 3601 or $grid_size == 1201) {
    croak "Unknown tile format for '$self->{DEMs}{$DEMname}{DEMpath}': grid size is '$grid_size', should be 3601 or 1201";
  }
  # the DEMs start in the NW corner with $grid_size - 1 intervals
  my $ilat = (1 - ($lat - $flat)) * ($grid_size - 1);
  my $ilon =      ($lon - $flon)  * ($grid_size - 1);
  say STDERR "  Grid size lat lon: $grid_size  $ilat  $ilon" if $debug;
  $self->{grid_size} = $grid_size;
  $self->{elevation} = _interpolate ($dem, $ilat, $ilon);
  return $self->{elevation};
}

sub _interpolate {
  my ($f, $x, $y) = @_;
  my $x1 = POSIX::floor $x;
  my $x2 = POSIX::ceil  $x;
  my $y1 = POSIX::floor $y;
  my $y2 = POSIX::ceil  $y;
  my $f11 = unpack ("s>*", substr ($$f, 2*($x1*$grid_size+$y1), 2));    # unpack signed big-endian 16-bit integer to elevation value
  my $f21 = unpack ("s>*", substr ($$f, 2*($x2*$grid_size+$y1), 2));    # unpack signed big-endian 16-bit integer to elevation value
  my $f12 = unpack ("s>*", substr ($$f, 2*($x1*$grid_size+$y2), 2));    # unpack signed big-endian 16-bit integer to elevation value
  my $f22 = unpack ("s>*", substr ($$f, 2*($x2*$grid_size+$y2), 2));    # unpack signed big-endian 16-bit integer to elevation value
  say STDERR "  Grid corners: ($x1,$y1)  ($x2,$y1)  ($x1,$y2)  ($x2,$y2)" if $debug;
  say STDERR "  Elevation at corners: $f11  $f21  $f12  $f22" if $debug;
  # bilinear interpolation as per https://github.com/racemap/elevation-service/blob/master/hgt.js
  # using the simplifying fact that ($x2-$x1)==1 and ($y2-$y1)==1
  my $xx = $x - $x1;
  my $yy = $y - $y1;
  my $f1 = _avg ($f11, $f21, $xx);
  my $f2 = _avg ($f12, $f22, $xx);
  say STDERR "  Interpolated elevation: "._avg ($f1, $f2, $yy) if $debug;
  return _avg ($f1, $f2, $yy);
}

sub _avg {
  my ($f1, $f2, $x) = @_;
  return $f1 + ($f2 - $f1) * $x;
}

sub _readDEM {
  use IO::Uncompress::AnyUncompress qw(anyuncompress $AnyUncompressError);
  use HTTP::Tiny;
  my ($self, $DEMname) = @_;
  my $nslt = substr($DEMname,0,3);
  my $path_to_hgt_gz = "$nslt/$DEMname.hgt.gz";
  $switch = 0;
  $cache = 0;
  $fail = 0;
  @default_DEMs = ("$path_to_hgt_gz", "$DEMname.zip");
  @want_DEMnames = ("$DEMname.hgt.gz", "$DEMname.zip");
  my $dem;
  $status_descr = -d $folder ? "Folder" : "Url";
  my $path = -d $folder ? _findDEM ($folder) : "$folder/$path_to_hgt_gz";
  if (-e $path) {
    say STDERR "  Reading DEM '$path'" if $debug;
    $self->{DEMs}{$DEMname}{DEMpath}=$path;
    anyuncompress $path => \$dem or croak "anyuncompress failed on '$path': $AnyUncompressError";
    return \$dem;
  }
  unless ($path =~ m/^https?:\/\//) {
    say STDERR "  DEM '$path' not found -> switch to '$url'" if $debug;
    $status_descr .= "->Switch";
    $switch = 1;
  }
  my $cache_path = $cache_folder ne "" ? _findDEM ($cache_folder) : undef;
  if (defined $cache_path and -e $cache_path) {
    say STDERR "  Reading DEM from cache '$cache_path'" if $debug;
    $self->{DEMs}{$DEMname}{DEMpath}=$cache_path;
    $status_descr .= "->Cached";
    $cache = 1;
    anyuncompress $cache_path => \$dem or croak "anyuncompress failed on '$cache_path': $AnyUncompressError";
    return \$dem;
  }
  $path = "$url/$path_to_hgt_gz";
  say STDERR "  Getting DEM '$path'" if $debug;
  $status_descr .= "->Url" unless ($status_descr eq "Url");
  my $response = HTTP::Tiny->new->get($path);    # get gzip archive file .hgt.gz
  unless ($response->{success}) {
    # no success
    carp "  DEM '$path'- $response->{status} $response->{reason}. All of its elevations will read as 0";
    $status_descr .= "->Failed";
    $fail = 1;
    return 0;
  }
  if ($cache_folder ne "" and -d $cache_folder) {
    unless (-d "$cache_folder/$nslt") {mkdir "$cache_folder/$nslt", 0755}
    open FILE, '>', "$cache_path" or croak "'$cache_path' error opening: $!";
    binmode FILE;
    print FILE $response->{content};
    close FILE;
    $status_descr .= "->Saved";
    say STDERR "  Saved DEM cache file '$cache_path'" if $debug;
  }
  $self->{DEMs}{$DEMname}{DEMpath}=$path;
  anyuncompress \$response->{content} => \$dem or croak "anyuncompress failed on '$path': $AnyUncompressError";
  return \$dem;
}

sub _findDEM {
  use File::Find;
  my ($folder) = @_;
  for my $default (@default_DEMs) {
    say STDERR "  Found default DEM '$folder/$default'" if ($debug and -e "$folder/$default");
    return "$folder/$default" if (-e "$folder/$default");
  }
  splice (@DEMnames,0);
  find(\&_wanted, $folder);
  say STDERR "  Found wanted DEM '$DEMnames[0]'" if ($debug and defined $DEMnames[0]);
  return $DEMnames[0] // "$folder/$default_DEMs[0]";
}

sub _wanted {
  use List::Util 'any';
  my $filename = $_;
  push (@DEMnames, $File::Find::name) if any {$_ =~ m/^$filename$/i} @want_DEMnames;     # add file to list if matched
}

1; # End of Geo::Elevation::HGT

__END__