source: subversion/applications/utils/distance_maps/distance.pl @ 29705

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

changed format of output to OSM

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