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


line stmt bran cond sub pod time code
1             package Geo::Elevation::HGT;
2              
3 1     1   69409 use 5.010;
  1         4  
4 1     1   5 use strict;
  1         2  
  1         20  
5 1     1   4 use warnings;
  1         2  
  1         34  
6 1     1   612 use POSIX ();
  1         6436  
  1         28  
7 1     1   7 use Carp;
  1         2  
  1         1041  
8             # set the version for version checking
9             our $VERSION = '0.04';
10             # file-private lexicals
11             my $grid_size; # .hgt grid size = 3601x3601 for 1-minute DEMs or 1201x1201 for 3-minute DEMs
12             my @DEMnames;
13             my @default_DEMs;
14             my @want_DEMnames;
15             my $status_descr;
16             my $url;
17             my $folder;
18             my $cache_folder;
19             my $debug;
20             my $switch;
21             my $cache;
22             my $fail;
23              
24             sub new {
25 0     0 1   my ($class, %params) = @_;
26 0           %params = (
27             url => "https://elevation-tiles-prod.s3.amazonaws.com/skadi", # info at https://registry.opendata.aws/terrain-tiles/
28             folder => "https://elevation-tiles-prod.s3.amazonaws.com/skadi",
29             cache_folder => "",
30             debug => 0,
31             %params
32             );
33 0           my $self = {};
34 0           while ( my($key,$value) = each %params ) {
35 0           $self->{$key} = $value;
36             }
37 0           bless $self, $class;
38 0           return $self;
39             }
40              
41             sub get_elevation_batch_hgt {
42 0     0 1   my ($self, $batch_latlon) = @_;
43 0           my @elegeh;
44 0           for my $latlon ( @$batch_latlon ) {
45 0           my ($lat, $lon) = @$latlon;
46 0           push (@elegeh, $self->get_elevation_hgt ($lat, $lon));
47             }
48 0           return \@elegeh;
49             }
50              
51             sub get_elevation_hgt {
52 0     0 1   my ($self, $lat, $lon) = @_;
53 0           $folder = $self->{folder};
54 0           $url = $self->{url};
55 0           $cache_folder = $self->{cache_folder};
56 0           $debug = $self->{debug};
57 0           $status_descr = "Memory";
58 0 0         say STDERR "get_elevation_hgt" if $debug;
59 0 0         say STDERR " Parameters: folder=>'$folder', url=>'$url', cache_folder=>'$cache_folder', debug=>$debug" if $debug;
60 0 0         say STDERR " Input lat lon: $lat $lon" if $debug;
61 0           my $flat = POSIX::floor $lat;
62 0           my $flon = POSIX::floor $lon;
63 0 0         my $ns = $flat < 0 ? "S" : "N";
64 0 0         my $ew = $flon < 0 ? "W" : "E";
65 0           my $lt = sprintf ("%02s", abs($flat));
66 0           my $ln = sprintf ("%03s", abs($flon));
67 0           my $DEMname = "$ns$lt$ew$ln";
68 0 0         say STDERR " Tile lat lon: $flat $flon" if $debug;
69             # read DEM unless already defined
70 0 0 0       say STDERR " Using DEM in memory: '$flat $flon'" if ($debug and defined $self->{DEMs}{$DEMname}{DEM});
71 0   0       $self->{DEMs}{$DEMname}{DEM} //= $self->_readDEM($DEMname); # //= Logical Defined-Or Assignment Operator
72 0           my $dem = $self->{DEMs}{$DEMname}{DEM};
73 0 0         say STDERR " Status: $status_descr" if $debug;
74 0           $self->{lat} = $lat;
75 0           $self->{lon} = $lon;
76 0           $self->{status_descr} = $status_descr;
77 0           $self->{switch} = $switch;
78 0           $self->{cache} = $cache;
79 0           $self->{fail} = $fail;
80 0           $self->{DEMname} = $DEMname;
81 0 0         if (ref($dem) eq "") {
82 0 0         say STDERR " No data in DEM: '$flat $flon' returning elevation 0" if $debug;
83 0           $self->{grid_size} = 0;
84 0           $self->{elevation} = 0;
85 0           return $self->{elevation};
86             }
87 0           $grid_size = sqrt (length ($$dem)/2); # grid size of DEM
88 0 0 0       unless ($grid_size == 3601 or $grid_size == 1201) {
89 0           croak "Unknown tile format for '$self->{DEMs}{$DEMname}{DEMpath}': grid size is '$grid_size', should be 3601 or 1201";
90             }
91             # the DEMs start in the NW corner with $grid_size - 1 intervals
92 0           my $ilat = (1 - ($lat - $flat)) * ($grid_size - 1);
93 0           my $ilon = ($lon - $flon) * ($grid_size - 1);
94 0 0         say STDERR " Grid size lat lon: $grid_size $ilat $ilon" if $debug;
95 0           $self->{grid_size} = $grid_size;
96 0           $self->{elevation} = _interpolate ($dem, $ilat, $ilon);
97 0           return $self->{elevation};
98             }
99              
100             sub _interpolate {
101 0     0     my ($f, $x, $y) = @_;
102 0           my $x1 = POSIX::floor $x;
103 0           my $x2 = POSIX::ceil $x;
104 0           my $y1 = POSIX::floor $y;
105 0           my $y2 = POSIX::ceil $y;
106 0           my $f11 = unpack ("s>*", substr ($$f, 2*($x1*$grid_size+$y1), 2)); # unpack signed big-endian 16-bit integer to elevation value
107 0           my $f21 = unpack ("s>*", substr ($$f, 2*($x2*$grid_size+$y1), 2)); # unpack signed big-endian 16-bit integer to elevation value
108 0           my $f12 = unpack ("s>*", substr ($$f, 2*($x1*$grid_size+$y2), 2)); # unpack signed big-endian 16-bit integer to elevation value
109 0           my $f22 = unpack ("s>*", substr ($$f, 2*($x2*$grid_size+$y2), 2)); # unpack signed big-endian 16-bit integer to elevation value
110 0 0         say STDERR " Grid corners: ($x1,$y1) ($x2,$y1) ($x1,$y2) ($x2,$y2)" if $debug;
111 0 0         say STDERR " Elevation at corners: $f11 $f21 $f12 $f22" if $debug;
112             # bilinear interpolation as per https://github.com/racemap/elevation-service/blob/master/hgt.js
113             # using the simplifying fact that ($x2-$x1)==1 and ($y2-$y1)==1
114 0           my $xx = $x - $x1;
115 0           my $yy = $y - $y1;
116 0           my $f1 = _avg ($f11, $f21, $xx);
117 0           my $f2 = _avg ($f12, $f22, $xx);
118 0 0         say STDERR " Interpolated elevation: "._avg ($f1, $f2, $yy) if $debug;
119 0           return _avg ($f1, $f2, $yy);
120             }
121              
122             sub _avg {
123 0     0     my ($f1, $f2, $x) = @_;
124 0           return $f1 + ($f2 - $f1) * $x;
125             }
126              
127             sub _readDEM {
128 1     1   658 use IO::Uncompress::AnyUncompress qw(anyuncompress $AnyUncompressError);
  1         73552  
  1         131  
129 1     1   781 use HTTP::Tiny;
  1         35921  
  1         640  
130 0     0     my ($self, $DEMname) = @_;
131 0           my $nslt = substr($DEMname,0,3);
132 0           my $path_to_hgt_gz = "$nslt/$DEMname.hgt.gz";
133 0           $switch = 0;
134 0           $cache = 0;
135 0           $fail = 0;
136 0           @default_DEMs = ("$path_to_hgt_gz", "$DEMname.zip");
137 0           @want_DEMnames = ("$DEMname.hgt.gz", "$DEMname.zip");
138 0           my $dem;
139 0 0         $status_descr = -d $folder ? "Folder" : "Url";
140 0 0         my $path = -d $folder ? _findDEM ($folder) : "$folder/$path_to_hgt_gz";
141 0 0         if (-e $path) {
142 0 0         say STDERR " Reading DEM '$path'" if $debug;
143 0           $self->{DEMs}{$DEMname}{DEMpath}=$path;
144 0 0         anyuncompress $path => \$dem or croak "anyuncompress failed on '$path': $AnyUncompressError";
145 0           return \$dem;
146             }
147 0 0         unless ($path =~ m/^https?:\/\//) {
148 0 0         say STDERR " DEM '$path' not found -> switch to '$url'" if $debug;
149 0           $status_descr .= "->Switch";
150 0           $switch = 1;
151             }
152 0 0         my $cache_path = $cache_folder ne "" ? _findDEM ($cache_folder) : undef;
153 0 0 0       if (defined $cache_path and -e $cache_path) {
154 0 0         say STDERR " Reading DEM from cache '$cache_path'" if $debug;
155 0           $self->{DEMs}{$DEMname}{DEMpath}=$cache_path;
156 0           $status_descr .= "->Cached";
157 0           $cache = 1;
158 0 0         anyuncompress $cache_path => \$dem or croak "anyuncompress failed on '$cache_path': $AnyUncompressError";
159 0           return \$dem;
160             }
161 0           $path = "$url/$path_to_hgt_gz";
162 0 0         say STDERR " Getting DEM '$path'" if $debug;
163 0 0         $status_descr .= "->Url" unless ($status_descr eq "Url");
164 0           my $response = HTTP::Tiny->new->get($path); # get gzip archive file .hgt.gz
165 0 0         unless ($response->{success}) {
166             # no success
167 0           carp " DEM '$path'- $response->{status} $response->{reason}. All of its elevations will read as 0";
168 0           $status_descr .= "->Failed";
169 0           $fail = 1;
170 0           return 0;
171             }
172 0 0 0       if ($cache_folder ne "" and -d $cache_folder) {
173 0 0         unless (-d "$cache_folder/$nslt") {mkdir "$cache_folder/$nslt", 0755}
  0            
174 0 0         open FILE, '>', "$cache_path" or croak "'$cache_path' error opening: $!";
175 0           binmode FILE;
176 0           print FILE $response->{content};
177 0           close FILE;
178 0           $status_descr .= "->Saved";
179 0 0         say STDERR " Saved DEM cache file '$cache_path'" if $debug;
180             }
181 0           $self->{DEMs}{$DEMname}{DEMpath}=$path;
182 0 0         anyuncompress \$response->{content} => \$dem or croak "anyuncompress failed on '$path': $AnyUncompressError";
183 0           return \$dem;
184             }
185              
186             sub _findDEM {
187 1     1   10 use File::Find;
  1         2  
  1         271  
188 0     0     my ($folder) = @_;
189 0           for my $default (@default_DEMs) {
190 0 0 0       say STDERR " Found default DEM '$folder/$default'" if ($debug and -e "$folder/$default");
191 0 0         return "$folder/$default" if (-e "$folder/$default");
192             }
193 0           splice (@DEMnames,0);
194 0           find(\&_wanted, $folder);
195 0 0 0       say STDERR " Found wanted DEM '$DEMnames[0]'" if ($debug and defined $DEMnames[0]);
196 0   0       return $DEMnames[0] // "$folder/$default_DEMs[0]";
197             }
198              
199             sub _wanted {
200 1     1   10 use List::Util 'any';
  1         3  
  1         177  
201 0     0     my $filename = $_;
202 0 0   0     push (@DEMnames, $File::Find::name) if any {$_ =~ m/^$filename$/i} @want_DEMnames; # add file to list if matched
  0            
203             }
204              
205             1; # End of Geo::Elevation::HGT
206              
207             __END__