File Coverage

blib/lib/Geo/Elevation/HGT.pm

Criterion	Covered	Total	%
statement	29	209	13.8
branch	0	94	0.0
condition	0	27	0.0
subroutine	10	21	47.6
pod	3	3	100.0
total	42	354	11.8

package Geo::Elevation::HGT;

use 5.010;
use strict;
use warnings;
use POSIX ();
use Carp;
# set the version for version checking
our $VERSION     = '0.08';
# file-private lexicals
my $grid_size; # .hgt grid size = 3601x3601 for 1-minute DEMs or 7201x7201 for 0.5-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 $bicubic;
my $bicubic_current;
my $bicubic_mixed;
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,
    bicubic => 0,
    bicubic_current => 0,
    bicubic_mixed => 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};
  $bicubic = $self->{bicubic};
  $bicubic_current = $self->{bicubic_current};
  $bicubic_mixed = $self->{bicubic_mixed};
  $status_descr = "Memory";
  say STDERR "get_elevation_hgt" if $debug;
  say STDERR "  Parameters: folder=>'$folder', url=>'$url', cache_folder=>'$cache_folder', debug=>$debug, bicubic=>$bicubic" 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 == 7201 or $grid_size == 1201) {
    croak "Unknown tile format for '$self->{DEMs}{$DEMname}{DEMpath}': grid size is '$grid_size', should be 3601 or 7201 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} = $bicubic ? _interpolate_bicubic ($dem, $ilat, $ilon) : _interpolate_bilinear ($dem, $ilat, $ilon);
  return $self->{elevation};
}

sub _interpolate_bicubic {
  use List::Util qw( min max );
  my ($f, $x, $y) = @_;
  my $x1 = POSIX::floor $x;
  my $x2 = POSIX::ceil  $x;
  my $x0 = max ($x1-1, 0);
  my $x3 = min ($x2+1, $grid_size-1);
  my $y1 = POSIX::floor $y;
  my $y2 = POSIX::ceil  $y;
  my $y0 = max ($y1-1, 0);
  my $y3 = min ($y2+1, $grid_size-1);

  my $f00 = unpack ("s>*", substr ($$f, 2*($x0*$grid_size+$y0), 2));    # unpack signed big-endian 16-bit integer to elevation value
  my $f01 = unpack ("s>*", substr ($$f, 2*($x0*$grid_size+$y1), 2));    # unpack signed big-endian 16-bit integer to elevation value
  my $f02 = unpack ("s>*", substr ($$f, 2*($x0*$grid_size+$y2), 2));    # unpack signed big-endian 16-bit integer to elevation value
  my $f03 = unpack ("s>*", substr ($$f, 2*($x0*$grid_size+$y3), 2));    # unpack signed big-endian 16-bit integer to elevation value

  my $f10 = unpack ("s>*", substr ($$f, 2*($x1*$grid_size+$y0), 2));    # unpack signed big-endian 16-bit integer to elevation value
  my $f11 = unpack ("s>*", substr ($$f, 2*($x1*$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 $f13 = unpack ("s>*", substr ($$f, 2*($x1*$grid_size+$y3), 2));    # unpack signed big-endian 16-bit integer to elevation value

  my $f20 = unpack ("s>*", substr ($$f, 2*($x2*$grid_size+$y0), 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 $f22 = unpack ("s>*", substr ($$f, 2*($x2*$grid_size+$y2), 2));    # unpack signed big-endian 16-bit integer to elevation value
  my $f23 = unpack ("s>*", substr ($$f, 2*($x2*$grid_size+$y3), 2));    # unpack signed big-endian 16-bit integer to elevation value

  my $f30 = unpack ("s>*", substr ($$f, 2*($x3*$grid_size+$y0), 2));    # unpack signed big-endian 16-bit integer to elevation value
  my $f31 = unpack ("s>*", substr ($$f, 2*($x3*$grid_size+$y1), 2));    # unpack signed big-endian 16-bit integer to elevation value
  my $f32 = unpack ("s>*", substr ($$f, 2*($x3*$grid_size+$y2), 2));    # unpack signed big-endian 16-bit integer to elevation value
  my $f33 = unpack ("s>*", substr ($$f, 2*($x3*$grid_size+$y3), 2));    # unpack signed big-endian 16-bit integer to elevation value

  say STDERR "  Grid corners: ($x0,$y0)  ($x1,$y1)  ($x2,$y2)  ($x3,$y3)" if $debug;
  # say STDERR "  Elevation at corners: $f11  $f21  $f12  $f22" if $debug;
  say STDERR "  Elevation at corners: $f00  $f01  $f02  $f03" if $debug;
  say STDERR "  Elevation at corners: $f10  $f11  $f12  $f13" if $debug;
  say STDERR "  Elevation at corners: $f20  $f21  $f22  $f23" if $debug;
  say STDERR "  Elevation at corners: $f30  $f31  $f32  $f33" if $debug;
  # Bicubic interpolation as per https://www.paulinternet.nl/?page=bicubic
  # using the simplifying fact that ($x2-$x1)==1 and ($y2-$y1)==1
  my $xx = $x - $x1;
  my $yy = $y - $y1;
  my $f0 = _bicubic ($f00, $f01, $f02, $f03, $yy);
  my $f1 = _bicubic ($f10, $f11, $f12, $f13, $yy);
  my $f2 = _bicubic ($f20, $f21, $f22, $f23, $yy);
  my $f3 = _bicubic ($f30, $f31, $f32, $f33, $yy);
  say STDERR "  bicubic interpolated elevation: "._bicubic ($f0, $f1, $f2, $f3, $xx) if $debug;
  return _bicubic ($f0, $f1, $f2, $f3, $xx);
}

sub _bicubic {
  my ($f0, $f1, $f2, $f3, $x) = @_;
  if ($bicubic_current) {
    return $f1 + $x*($f1 - $f0 + $x*(2.0*$f0 - 5.0*$f1 + 4.0*$f2 - $f3 + $x*(3.0*($f1 - $f2) + $f3 - $f0)));    # use slope of a line between the previous and the «current» point as the derivative at a point, that is (p1-p0) and (p3-p2)
  }
  elsif ($bicubic_mixed) {
    return $f1 + 0.75 * $x*($f2 - $f0 + $x*(2.0*$f0 - 5.0*$f1 + 4.0*$f2 - $f3 + $x*(3.0*($f1 - $f2) + $f3 - $f0)));    # use «current» and «next» point «mixed» slope
  }
  else {
    return $f1 + 0.5 * $x*($f2 - $f0 + $x*(2.0*$f0 - 5.0*$f1 + 4.0*$f2 - $f3 + $x*(3.0*($f1 - $f2) + $f3 - $f0)));    # use slope of a line between the previous and the «next» point as the derivative at a point, that is (p2-p0)/2 and (p3-p1)/2
  }
}

sub _interpolate_bilinear {
  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 "  bilinear 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 my $file_handle, '>', "$cache_path" or croak "'$cache_path' error opening: $!";
    binmode $file_handle;
    print $file_handle $response->{content};
    close $file_handle;
    $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__