source: subversion/applications/utils/distance_maps/distance_png.pl @ 16638

Last change on this file since 16638 was 16638, checked in by mungewell, 10 years ago

New tools for walking all highways and computing distance to each node

  • Property svn:executable set to *
File size: 10.0 KB
Line 
1#!/usr/bin/perl
2
3# All rights to this package are hereby disclaimed and its contents released
4# into the public domain by the author. Where this is not possible, you may
5# use this file under the same terms as Perl itself.
6
7# Originally written by Frederik Ramm <frederik@remote.org>, 14 July, 2009
8# Modified by Simon Wood <simon@mungewell.org>, 14 July 2009
9
10# Usage:
11# distance.pl SOURCE.OSM <target_lat> <target_lon> [<interpolation>]
12#
13# interpolation  0 = all nodes in SOURCE.OSM
14# interpolation -1 = junction/terminus nodes in SOURCE.OSM
15# interpolation >0 = all nodes in SOURCE.OSM, plus interpolated nodes seperated
16#                    by no more than 'X'm
17
18# NODES:
19#   * $nodes = { "1234" => { "lat" => 1.2, "lon" => "2.3", "nb" => [
20#       { "node" => ..., "len" => 1.2345 }, ... }, ... }
21#     
22# WAYS:
23#   * not stored, only processed to build the node neighbour
24#     array.
25#
26# -----------------------------------------------------------------------------
27
28use strict;
29use warnings;
30use bytes;
31
32use GD;
33use IO::Seekable;
34use Math::Trig qw(great_circle_distance deg2rad pi tan sec);
35
36my ($maxlat, $minlat, $maxlon, $minlon);
37my ($width, $height, $scale);
38my $start_time;
39my $nodes = {};
40
41my $output_file = "output.png";
42my $input_file = $ARGV[0];
43my $start_lat = $ARGV[1];
44my $start_lon = $ARGV[2];
45
46my $max_between = $ARGV[3] + 0;
47my $newnode = -1;
48
49# read into memory
50parse_file($input_file);
51analyze_graph();
52
53# must have an OSM node to start at, but only lat/lon given above, so
54# find nearest (caution
55my $startnode = find_nearest_node($start_lat, $start_lon);
56# the 2000 is the distance in metres from starting point
57flood_fill($startnode, 2000);
58draw_graph(20000);
59
60# -------------------------------------------------------------------------------
61# draw_graph()
62#
63# makes a PNG file from the graph.
64#
65# parameters:
66# * startnode - node reference to highlight in graph
67# -------------------------------------------------------------------------------
68
69sub draw_graph
70{
71    timer_start("draw_graph");
72    $width = 1024;
73    $scale = $width / ($maxlon-$minlon) / 10000;
74    $height = (ProjectF($maxlat) - ProjectF($minlat)) * 180 / pi * 10000 * $scale;
75    my $maxdist = shift;
76
77    my $pic = new GD::Image($width, $height);
78    my $gray = $pic->colorAllocate(128, 128, 128);
79    my $white = $pic->colorAllocate(255, 255, 255);
80    my $black = $pic->colorAllocate(0, 0, 0);
81
82    my $colors = [
83        $pic->colorAllocate(255, 0, 0),
84        $pic->colorAllocate(128, 0, 0),
85    ];
86    my $colindex = 0;
87
88    $pic->filledRectangle(0, 0, $width, $height, $white);
89
90    # draw grey lines
91    foreach my $node(values %$nodes)
92    {
93        my $pos1 = project($node);
94        $pic->filledRectangle($pos1->[0]-1, $pos1->[1]-1, $pos1->[0]+1, $pos1->[1]+1, $gray);
95        foreach my $other(@{$node->{nb}})
96        {
97            my $pos2 = project($other->{nd});
98            $pic->line(@$pos1, @$pos2, $gray);
99        }
100    }
101
102    $colindex=0;
103    foreach my $node(values %$nodes)
104    {
105        next unless defined($node->{dist});
106        if ($max_between>-1 || @{$node->{nb}}!=2) 
107        {
108            my $color = $colors->[$colindex + 1];
109
110            if (@{$node->{nb}} < 2 )
111            {
112                # select brighter color for terminus nodes
113                $color = $colors->[$colindex + 0];
114            }
115
116
117            my $pos1 = project($node);
118            $pic->filledRectangle($pos1->[0]-2, $pos1->[1]-2, $pos1->[0]+2, $pos1->[1]+2, $color);
119        }
120    }
121
122    open(PNG, ">$output_file") or die;
123    print PNG $pic->png;
124    close(PNG);
125    timer_stop();
126}
127
128# -------------------------------------------------------------------------------
129# analyze_graph()
130#
131# computes max/min lon/lat.
132# -------------------------------------------------------------------------------
133sub analyze_graph()
134{
135    timer_start("analyze_graph");
136    $minlon = 999; 
137    $minlat = 999;
138    $maxlat = -999;
139    $maxlon = -999;
140
141    foreach my $node(values %$nodes)
142    {
143        $maxlat = $node->{lat} if ($node->{lat} > $maxlat);
144        $minlat = $node->{lat} if ($node->{lat} < $minlat);
145        $maxlon = $node->{lon} if ($node->{lon} > $maxlon);
146        $minlon = $node->{lon} if ($node->{lon} < $minlon);
147    }
148    timer_stop();
149
150}
151
152# -------------------------------------------------------------------------------
153# find_nearest_node()
154#
155# finds and returns the node (reference) that is nearest to given position.
156# (does not use proper great circle formula, just EPSG4326 distance)
157# -------------------------------------------------------------------------------
158sub find_nearest_node
159{
160    timer_start("find_nearest_node");
161    my ($lat, $lon) = @_;
162    my $best_dist = 999999;
163    my $best_node;
164
165    foreach my $node(values %$nodes)
166    {
167        my $d = ($node->{lat}-$lat)**2 + ($node->{lon}-$lon)**2;
168        if ($d < $best_dist)
169        {
170            $best_dist = $d;
171            $best_node = $node;
172        }
173    }
174    timer_stop();
175    return $best_node;
176}
177
178sub timer_start
179{
180    my $current_task = shift;
181    $start_time = time();
182    $|=1;
183    print "$current_task...";
184}
185
186sub timer_stop
187{
188    my $elaps = time()-$start_time;
189    printf " %ds\n", $elaps;
190}
191
192sub project
193{
194    my $node = shift;
195    return [
196        $width - ($maxlon-$node->{lon})*10000*$scale, 
197        $height + (ProjectF($minlat)-ProjectF($node->{lat})) * 180/pi * 10000 * $scale 
198    ];
199}
200
201sub ProjectF
202{
203    my $Lat = shift() / 180 * pi;
204    my $Y = log(tan($Lat) + sec($Lat));
205    return($Y);
206}
207
208# -------------------------------------------------------------------------------
209# flood_fill()
210#
211# recursive flood-filling of the graph. starts at given start node and assigns
212# remaining dist; then visits all neighbours and updates their respective dists
213# with the remaining dist minus dist to get there.
214# -------------------------------------------------------------------------------
215sub flood_fill
216{
217    my ($startnode, $remain_dist) = @_;
218    timer_start("flood_fill");
219    my $queue = [ $startnode ];
220    $startnode->{dist} = 0;
221    while (my $node = shift @$queue)
222    {
223        my $d = $node->{dist};
224        foreach my $neigh(@{$node->{nb}})
225        {
226            next if ($neigh->{len}<=0);
227            my $newd = $d + $neigh->{len};
228            # we still have distance left to travel to that node
229            if (!defined($neigh->{nd}{dist}) || $newd < $neigh->{nd}{dist})
230            {
231                # the node has not been visited, or has been visited
232                # on a longer journey.
233                $neigh->{nd}{dist} = $newd;
234                push(@$queue, $neigh->{nd});
235            }
236        }
237    }
238    timer_stop();
239}
240
241sub parse_file
242{
243    my $filename = shift;
244    my $in = new IO::File $filename, "r" or die "cannot open $filename";
245    my $buffer;
246    my $edge = {};
247    my $waynodes;
248    my $hwtype;
249    my $weight_forward;
250    my $weight_backward;
251
252    timer_start("parsing ways");
253    while(<$in>)
254    {
255        $buffer .= $_;
256        while($buffer =~ /^\s*([^<]*?)\s*(<[^>]*>)(.*)$/)
257        {
258            $buffer = $3;
259            my $xmltag = $2;
260            if ($xmltag =~ /^\s*<way /)
261            {
262                $waynodes = [];
263                undef $hwtype;
264                $weight_forward=1;
265                $weight_backward=1;
266            }
267            elsif ($xmltag =~ /^\s*<\/osm/)
268            {
269                last;
270            }
271            elsif ($xmltag =~ /^\s*<relation/)
272            {
273                last;
274            }
275            elsif ($xmltag =~ /^\s*<tag k=['"]highway["'].*v=["'](.*)['"]/)
276            {
277                $hwtype = $1;
278            }
279            elsif ($xmltag =~ /^\s*<tag k=['"]oneway["'].*v=["'](.*)['"]/)
280            {
281                $weight_backward=-1;
282            }
283            elsif ($xmltag =~ /^\s*<nd ref=['"](\d+)["']/)
284            {
285                push(@$waynodes, $1);
286            }
287            elsif ($xmltag =~ /^\s*<\/way/)
288            {
289                next unless defined $hwtype;
290                # this is where you would check whether $hwtype contains a way
291                # you want considered, and do "next" if not
292                for (my $i=1; $i<scalar(@$waynodes); $i++)
293                {
294                    # you could also assign different weights to different kinds of
295                    # ways but for now we treat them all the same
296                    $edge->{$waynodes->[$i-1]}{$waynodes->[$i]} = $weight_forward;
297                    $edge->{$waynodes->[$i]}{$waynodes->[$i-1]} = $weight_backward;
298                }
299            }
300        }
301    }
302    timer_stop();
303
304    # ways now parsed.
305
306    $in->seek(0, SEEK_SET);
307    undef $buffer;
308
309    timer_start("parsing nodes");
310    while(<$in>)
311    {
312        $buffer .= $_;
313        while($buffer =~ /^\s*([^<]*?)\s*(<[^>]*>)(.*)$/)
314        {
315            $buffer = $3;
316            my $xmltag = $2;
317            if ($xmltag =~ /^\s*<node.*\sid=["'](\d+)['"].*lat=["']([0-9.Ee-]+)["'].*lon=["']([0-9.Ee-]+)["']/)
318            {
319                if (defined($edge->{$1}))
320                {
321                    # we want this node
322                    $nodes->{$1} = { lat=>$2, lon=>$3, id=>$1 };
323                }
324            }
325            elsif ($xmltag =~ /^\s*<way /)
326            {
327                last;
328            }
329        }
330    }
331    timer_stop();
332
333    $in->close();
334
335    # nodes now parsed, but lengths of edges still need to be computed.
336
337    timer_start("computing edge lengths");
338    foreach my $nid(keys %$edge)
339    {
340        my $node = $nodes->{$nid};
341        die "node $nid referenced in way but not present" unless defined $node;
342        foreach my $othernid(keys(%{$edge->{$nid}}))
343        {
344            my $other = $nodes->{$othernid};
345            die "node $othernid referenced in way but not present" unless defined $other;
346            my $dist = great_circle_distance(deg2rad($node->{lon}), deg2rad(90-$node->{lat}),
347                deg2rad($other->{lon}), deg2rad(90-$other->{lat}), 6378135); # gives metres
348            $nodes->{$nid}{nb} = [] unless defined $nodes->{$nid}{nb};
349            push(@{$nodes->{$nid}{nb}}, { nd => $other, len => $dist * $edge->{$nid}{$othernid} });
350        }
351    }
352    timer_stop();
353}
354
Note: See TracBrowser for help on using the repository browser.