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

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

changed format of output to OSM

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