File Coverage

blib/lib/Geo/OSM/Imager.pm
Criterion Covered Total %
statement 15 17 88.2
branch n/a
condition n/a
subroutine 6 6 100.0
pod n/a
total 21 23 91.3


line stmt bran cond sub pod time code
1             package Geo::OSM::Imager;
2              
3 1     1   710 use 5.008001;
  1         3  
4 1     1   5 use strict;
  1         2  
  1         21  
5 1     1   14 use warnings;
  1         2  
  1         50  
6              
7             our $VERSION = "0.02";
8              
9 1     1   466 use Math::Trig;
  1         9192  
  1         136  
10 1     1   609 use Geo::Ellipsoid;
  1         3104  
  1         26  
11 1     1   196 use LWP::UserAgent;
  0            
  0            
12             use HTTP::Request::Common;
13             use List::Util qw(max);
14             use List::MoreUtils qw(minmax);
15             use Imager;
16              
17             =encoding utf-8
18              
19             =head1 NAME
20              
21             Geo::OSM::Imager - simplifies plotting onto OpenStreetMap tiles
22              
23             =head1 SYNOPSIS
24              
25             my $g=Geo::OSM::Imager->new(ua => 'MyApplication');
26             my $image=$g->init(\@points);
27             ...
28             my ($x,$y)=$g->latlon2xy($lat,$lon);
29             $image->circle(x=>$x,y=>$y,r=>50,color=>$blue);
30             ...
31             $image->circle($g->latlon2hash($lat,$lon),r=>50,color=>$blue);
32             ...
33             $image->write(file => 'test.png');
34              
35             =head1 DESCRIPTION
36              
37             This module sets up an Imager object made of OpenStreetMap tiles, for
38             drawing of geographic data.
39              
40             Beware of over-using OpenStreetMap tile servers, and see the usage
41             policy at https://operations.osmfoundation.org/policies/tiles/ .
42              
43             Be hesitant about drawing straight lines over long distances, as map
44             projections will cause distortion. Over more than a few hundred
45             metres, the author prefers to break the line into a series of points
46             and plot individual line segments.
47              
48             =cut
49              
50             =head1 USAGE
51              
52             =over
53              
54             =item new()
55              
56             Creates a new Geo::OSM::Imager object. Takes an optional hash of
57             parameters:
58              
59             maxx - maximum X size of the image, in pixels.
60             maxy - maximum Y size of the image, in pixels.
61              
62             The image will generally be between 50% and 100% of this size.
63              
64             margin - fractional margin around bounding points
65             marginlat - fractional latitude margin around bounding point
66             marginlon - fractional longitude margin around bounding points
67              
68             The fraction of the latitude/longitude span to leave as space around
69             the matter to be plotted. With a margin of zero, points will be
70             plotted right at the edges of the image. A margin of 1/7 works well,
71             and is the default. marginlat and marginlon allow you to define this
72             separately for latitude and longitude.
73              
74             tileage - minimum age to expire tiles
75              
76             The number of seconds after which a tile may be considered "old" and
77             re-downloaded. Tileserver usage policy forbids an expiry of less than
78             one week (604800s), which is the default.
79              
80             tiledir - directory for the tile cache
81              
82             The directory in which to store tiles; it must exist.
83              
84             tilesize - size of tiles
85              
86             The pixel size of each tile. Leave at its default of 256 unless you
87             know what you're doing.
88              
89             tileurl - base URL for downloading tiles
90              
91             The base URL for downloading files. If you are using a local
92             tileserver, or a public tileserver other than OpenStreetMap, set it
93             here.
94              
95             ua - user-agent
96              
97             Tileserver usage policy requires a "Valid HTTP User-Agent identifying
98             application".
99              
100             =cut
101              
102             sub new {
103             my ($pkg,%p)=@_;
104             my $self={margin => 1/7,
105             maxx => 2400,
106             maxy => 2400,
107             tileage => 604800,
108             tiledir => "$ENV{HOME}/Maps/OSM",
109             tilesize => 256,
110             tileurl => 'https://tile.openstreetmap.org/',
111             };
112             if (%p) {
113             foreach my $param (qw(margin marginlat marginlon maxx maxy tileage tiledir tilesize tileurl ua)) {
114             if (exists $p{$param}) {
115             $self->{$param}=$p{$param};
116             }
117             }
118             }
119             $self->{marginlat} ||= $self->{margin};
120             $self->{marginlon} ||= $self->{margin};
121             $self->{lwp}=LWP::UserAgent->new(agent => $self->{ua});
122             unless (defined $self->{ua}) {
123             die "Need a user-agent to access OpenStreetMap tile servers";
124             }
125             bless($self,$pkg);
126             return $self;
127             }
128              
129             =item init()
130              
131             Checks bounds and sets up the image. Pass an arrayref of points, each
132             of which can be either an arrayref [lat,lon] or a hashref including
133             lat and lon keys (or "latitude", "long", "longitude").
134              
135             These need not be the same points you're going to plot, though that's
136             obviously the easiest approach.
137              
138             Returns the Imager object.
139              
140             =cut
141              
142             sub init {
143             my ($self,$points)=@_;
144             my @series;
145             foreach my $p (@{$points}) {
146             if (ref $p eq 'ARRAY') {
147             map {push @{$series[$_]},$p->[$_]} (0,1);
148             } elsif (ref $p eq 'HASH') {
149             my $lat;
150             foreach my $latname (qw(lat latitude)) {
151             $lat ||= $p->{$latname};
152             }
153             my $lon;
154             foreach my $lonname (qw(lon long longitude)) {
155             $lon ||= $p->{$lonname};
156             }
157             if (defined $lat && defined $lon) {
158             push @{$series[0]},$lat;
159             push @{$series[1]},$lon;
160             }
161             }
162             }
163             my @minmax=map {[minmax(@{$_})]} @series;
164             $self->{geo}=Geo::Ellipsoid->new(units => 'degrees');
165              
166             my %bounds=(lon => [undef,undef],
167             lat => [undef,undef]);
168             $bounds{lat}[0]=$minmax[0][0]-($minmax[0][1]-$minmax[0][0])*$self->{marginlat};
169             $bounds{lat}[1]=$minmax[0][1]+($minmax[0][1]-$minmax[0][0])*$self->{marginlat};
170             $bounds{lon}[0]=$minmax[1][0]-($minmax[1][1]-$minmax[1][0])*$self->{marginlon};
171             $bounds{lon}[1]=$minmax[1][1]+($minmax[1][1]-$minmax[1][0])*$self->{marginlon};
172              
173             my $longdist=max(
174             $self->{geo}->to($bounds{lat}[0],$bounds{lon}[0],$bounds{lat}[0],$bounds{lon}[1]),
175             $self->{geo}->to($bounds{lat}[1],$bounds{lon}[0],$bounds{lat}[1],$bounds{lon}[1]),
176             ); # metres
177             my $longscale=$longdist/$self->{maxy}; # metres/pixel
178             my $latdist=max(
179             $self->{geo}->to($bounds{lat}[0],$bounds{lon}[0],$bounds{lat}[1],$bounds{lon}[0]),
180             $self->{geo}->to($bounds{lat}[0],$bounds{lon}[1],$bounds{lat}[1],$bounds{lon}[1]),
181             );
182             my $latscale=$latdist/$self->{maxx};
183             my $scale=max($longscale,$latscale); # make sure it fits, use wider scale
184             $self->{zoomlevel}=int(
185             log(
186             cos(
187             deg2rad(
188             ($bounds{lat}[0]+$bounds{lat}[1])/2
189             )
190             )*6378137.0*2*pi/$scale
191             )/log(2)-8
192             );
193              
194             while ($self->{zoomlevel} > 18) {
195             $self->{zoomlevel}--;
196             foreach my $mode (qw(lat lon)) {
197             my $mean=($bounds{$mode}[0]+$bounds{$mode}[1])/2;
198             foreach my $nn (0,1) {
199             $bounds{$mode}[$nn]+=($bounds{$mode}[$nn]-$mean);
200             }
201             }
202             }
203              
204             $self->{xmax}=int(($self->getTileNumber($bounds{lat}[1],$bounds{lon}[1]))[0]+.9999999);
205             $self->{xmin}=int(($self->getTileNumber($bounds{lat}[0],$bounds{lon}[0]))[0]);
206             $self->{ymax}=int(($self->getTileNumber($bounds{lat}[0],$bounds{lon}[0]))[1]+.9999999);
207             $self->{ymin}=int(($self->getTileNumber($bounds{lat}[1],$bounds{lon}[1]))[1]);
208              
209             my $img=Imager->new(xsize => $self->{tilesize}*($self->{xmax}-$self->{xmin}+1), ysize => $self->{tilesize}*($self->{ymax}-$self->{ymin}+1),
210             channels => 4);
211             mkdir "$self->{tiledir}/$self->{zoomlevel}";
212             foreach my $x ($self->{xmin}..$self->{xmax}) {
213             mkdir "$self->{tiledir}/$self->{zoomlevel}/$x";
214             foreach my $y ($self->{ymin}..$self->{ymax}) {
215             my $stub="$self->{zoomlevel}/$x/$y.png";
216             my $dl=1;
217             if (-e "$self->{tiledir}/$stub") {
218             my $fa=(stat("$self->{tiledir}/$stub"))[9];
219             if (time-$fa < $self->{tileage}) {
220             $dl=0;
221             }
222             }
223             if ($dl) {
224             my $rq=HTTP::Request->new(GET => "$self->{tileurl}/$stub");
225             my $rp=$self->{lwp}->request($rq);
226             if ($rp->is_success) {
227             open OUT,">$self->{tiledir}/$stub" or die "Can't open $self->{tiledir}/$stub for writing\n";
228             binmode OUT;
229             print OUT $rp->content;
230             close OUT;
231             } else {
232             die "Couldn't fetch $self->{tileurl}/$stub\n";
233             }
234             }
235             my $i=Imager->new;
236             $i->read(file => "$self->{tiledir}/$self->{zoomlevel}/$x/$y.png");
237             $img->rubthrough(left => $self->{tilesize}*($x-$self->{xmin}),
238             top => $self->{tilesize}*($y-$self->{ymin}),
239             src => $i);
240             }
241             }
242             $self->{offsetx}=$self->{offsety}=$self->{img}=0;
243             my $xclipmax=int(($self->latlon2xy($bounds{lat}[1],$bounds{lon}[1]))[0]);
244             my $xclipmin=int(($self->latlon2xy($bounds{lat}[0],$bounds{lon}[0]))[0]+.9999999);
245             my $yclipmax=int(($self->latlon2xy($bounds{lat}[0],$bounds{lon}[0]))[1]+.9999999);
246             my $yclipmin=int(($self->latlon2xy($bounds{lat}[1],$bounds{lon}[1]))[1]);
247             $self->{img}=$img->crop(left => $xclipmin,
248             top => $yclipmin,
249             width => $xclipmax-$xclipmin,
250             height => $yclipmax-$yclipmin);
251             $self->{offsetx}=-$xclipmin;
252             $self->{offsety}=-$yclipmin;
253             return $self->{img};
254             }
255              
256             =item image()
257              
258             Returns the Imager object.
259              
260             =cut
261              
262             # let the user plot onto it
263             sub image {
264             my ($self)=@_;
265             unless (exists $self->{img}) {
266             die "Not yet initialised.\n";
267             }
268             return $self->{img};
269             }
270              
271             =item zoom()
272              
273             Returns the zoom level of the initialised object. See
274             L and
275             L
276             for more.
277              
278             =cut
279              
280             sub zoom {
281             my ($self)=@_;
282             unless (exists $self->{img}) {
283             die "Not yet initialised.\n";
284             }
285             return $self->{zoomlevel};
286             }
287              
288             sub getTileNumber {
289             my ($self,$lat,$lon) = @_;
290             my $zoom = $self->{zoomlevel};
291             my $xtile = ($lon+180)/360 *2**$zoom;
292             my $ytile = (1 - log(tan(deg2rad($lat)) + sec(deg2rad($lat)))/pi)/2 *2**$zoom;
293             return ($xtile, $ytile);
294             }
295              
296             =item latlon2xy($lat,$lon)
297              
298             Given a (latitude, longitude) coordinate pair, returns the (x, y)
299             coordinate pair needed to plot onto the Imager object.
300              
301             =cut
302              
303             sub latlon2xy {
304             my ($self,$lat,$lon) = @_;
305             unless (exists $self->{img}) {
306             die "Not yet initialised.\n";
307             }
308             my ($x,$y)=$self->getTileNumber($lat,$lon);
309             $x=($x-$self->{xmin})*$self->{tilesize}+$self->{offsetx};
310             $y=($y-$self->{ymin})*$self->{tilesize}+$self->{offsety};
311             return ($x,$y);
312             }
313              
314             =item latlon2hash($lat,$lon)
315              
316             Given a (latitude, longitude) coordinate pair, returns a list of the
317             form ('x', $x, 'y', $y) for use with many Imager plotting functions.
318              
319             =cut
320              
321             sub latlon2hash {
322             my ($self,$lat,$lon) = @_;
323             my ($x,$y)=$self->latlon2xy($lat,$lon);
324             return ('x',$x,'y',$y);
325             }
326              
327             =item segment($lat1,$lon1,$lat2,$lon2,$step)
328              
329             Given two (latitude, longitude) coordinate pairs and a step value,
330             returns an arrayref of (latitude, longitude) coordinate pairs
331             interpolating the route on a great circle. This is generally worth
332             doing when distances exceed around 100 miles or high precision is
333             wanted.
334              
335             A positive step value is the length of each segment in metres. A
336             negative step value is the number of divisions into which the overall
337             line should be split.
338              
339             =cut
340              
341             sub segment {
342             my ($self,$lat1,$lon1,$lat2,$lon2,$step)=@_;
343             my @out=[$lat1,$lon1];
344             my ($r,$b)=$self->{geo}->to($lat1,$lon1,$lat2,$lon2);
345             if ($step<0) {
346             $step=-$r/$step;
347             }
348             my $ra=0;
349             while ($ra<$r) {
350             $ra+=$step;
351             push @out,[$self->{geo}->at($lat1,$lon1,$ra,$b)];
352             $out[-1][1]=$self->constrain($out[-1][1],180);
353             }
354             push @out,[$lat2,$lon2];
355             return @out;
356             }
357              
358             sub constrain {
359             my ($self,$angle,$range)=@_;
360             while ($angle>$range) {
361             $angle-=$range*2;
362             }
363             while ($angle<-$range) {
364             $angle+=$range*2;
365             }
366             return $angle;
367             }
368              
369             =back
370              
371             =head1 OTHER CONSIDERATIONS
372              
373             Note that you need not draw directly onto the supplied object: you can
374             create a new transparent image using the width and height of the one
375             provided by the module, draw onto that, and copy the results with a
376             rubthrough or compose command. See L for
377             more.
378              
379             =head1 BUGS
380              
381             Won't work to span +/- 180 degrees longitude.
382              
383             =head1 LICENSE
384              
385             Copyright (C) 2017 Roger Bell_West.
386              
387             This library is free software; you can redistribute it and/or modify
388             it under the same terms as Perl itself.
389              
390             =head1 AUTHOR
391              
392             Roger Bell_West Eroger@firedrake.orgE
393              
394             =head1 SEE ALSO
395              
396             L
397              
398             =cut
399              
400             1;