Tach,
Das ist in der Tat auch eine hochkomplexe Sache! Für die Lösung musste ich mein ganzes Team heranziehen.
1. Wir beschaffen die Daten:
Wir holen die Innenstadt von Köln ab:
wget -O data.osm http://www.openstreetmap.org/api/0.6/map?bbox=6.938,50.920,6.978,50.951
2. Wir erzeugen einen Routing-Graphen:
Dazu benutzen wir das überaus hochaufwendige Skript “analyze.pl”:
#!/usr/bin/perl
use strict;
use utf8;
binmode STDIN, ':utf8';
binmode STDOUT, ':utf8';
# type, id and attributes of current osm object
my $type;
my $id;
my %attr;
# lon lat of nodeIds
my %lon;
my %lat;
# nodeIds of ways
my @nodeIds;
while (<>) {
# start of osm object
if (m#^\s*<(changeset|node|way|relation)\s+id=['"](\d+)['"]#) {
$type = $1;
$id = substr($type,0,1).$2;
%attr = ();
@nodeIds= ();
}
# lon lat of nodeIds
if (m#^\s*<node\s+id=['"]([^'"]*)['"]\s+lat=['"]([^'"]*)['"]\s+lon=['"]([^'"]*)['"]#) {
$lat{"n$1"} = $2;
$lon{"n$1"} = $3;
next;
}
# nodes of ways
if (m#^\s*<nd\s+ref=['"]([^'"]*)['"]#) {
push (@nodeIds, "n$1");
next;
}
# tag value
if (m#^\s*<tag\s+k=['"]([^'"]*)['"]\s+v=['"]([^'"]*)['"]#) {
# xml entities are left as an exercise to the students™
$attr{$1} = $2;
next;
}
# end of way
next unless m#^\s*</way># && @nodeIds>=2;
#-------------------------------------------------------------------------------------
# preprocess way
# - test if usable
# - determine speed
#-------------------------------------------------------------------------------------
# skip if osm way is not a real world way
# (may accept bus or train or younameit)
next unless $attr{'highway'};
# too dangerous
next if $attr{'name'} =~ /katze/i;
# compute available forward and backward speed for way (m/s)
my $speedForward = 1;
my $speedBackward = 1;
$speedForward=$speedBackward = 5 if $attr{'highway'} eq "pedestrian";
# oneway
$speedForward = 0 if $attr{'oneway'} eq '-1';
$speedBackward= 0 if $attr{'oneway'} eq 'yes';
#-------------------------------------------------------------------------------------
# generate graph entries
#-------------------------------------------------------------------------------------
# coordinates of graph nodes
foreach my $nid (@nodeIds) {
print "node\t$nid\t$lon{$nid}\t$lat{$nid}\n";
}
# process segments of ways resulting in graph arcs
for (my $index=0; $index<@nodeIds-1; $index++) {
# start and end node of segment
my $snode = $nodeIds[$index ];
my $enode = $nodeIds[$index+1];
# distance in degrees
my $dLon = ($lon{$snode}-$lon{$enode}) * cos ($lat{$snode} * 3.1415925 / 180.0);
my $dLat = ($lat{$snode}-$lat{$enode});
# distance scaled to meters
my $distance = sqrt($dLon*$dLon + $dLat*$dLat) / 360.0 * 40000000;
# forward arc
if ($speedForward>0) {
my $time = $distance / $speedForward;
print "arc\t$snode\t$enode\t$time\t$id\tf\t$index\n";
}
# backward arc
if ($speedBackward>0) {
my $time = $distance / $speedBackward;
print "arc\t$enode\t$snode\t$time\t$id\tb\t$index\n";
}
}
}
Wir rufen das Skript auf:
./analyze.pl <data.osm >graph
Und erhalten eine Textdarstellung des Routing-Graphen “graph”.
Der Graph enthält Knoten mit Id und Koordinaten:
node n743325638 6.9722186 50.9463142
und gerichtete Kanten mit Angabe der verbundenen Knoten, der Entfernung, und der Quelle der Daten (OSM-Way-Id, Richtung, Segmentnummer):
arc n743325638 n743325592 10.0407049925005 w59865404 f 0
3. Wir führen die Erreichbarkeitsanalyse durch:
Dazu nutzen wir das unfassbar hochkomplexe “reachable.pl”:
#!/usr/bin/perl
# ./reachable.pl {lat} {lon} {distance} <data.graph
use strict;
use utf8;
# process args
my $startLat = 1.0 * shift (@ARGV);
my $startLon = 1.0 * shift (@ARGV);
my $maxDist = 1.0 * shift (@ARGV);
# graph
my %lon;
my %lat;
my %neighbours; # node -> list of nodes
my %arcLength; # snode,enode -> distance
#-----------------------------------------------------------------------
# read graph
#-----------------------------------------------------------------------
warn "reading graph ...\n";
while (<>) {
chomp;
my @f = split /\t/;
# node
if ($f[0] eq 'node') {
my $node = $f[1];
my $lon = $f[2];
my $lat = $f[3];
$lon{$node} = $lon;
$lat{$node} = $lat;
next;
}
# arc
if ($f[0] eq 'arc') {
my $snode = $f[1];
my $enode = $f[2];
my $dist = $f[3] || 0.001; # must NOT be zero
# ignode multiple arcs between two nodes
next if defined $arcLength{$snode, $enode};
# no neighbours yet -> create list
unless ($neighbours{$snode}) {
$neighbours{$snode} = [];
}
# store arc
$arcLength{$snode, $enode} = $dist;
push (@{$neighbours{$snode}}, $enode);
next;
}
}
#-----------------------------------------------------------------------
# lookup start node
#-----------------------------------------------------------------------
warn "lookup start node ...\n";
my $bestDistance = 1e30;
my $startNode;
my $cos = cos ($startLat * 3.1415925 / 180.0);
foreach my $node (keys %lon) {
# distance in degrees
my $dLon = ($lon{$node}-$startLon) * $cos;
my $dLat = ($lat{$node}-$startLat);
# skip square root and scaling
my $distance = $dLon*$dLon + $dLat*$dLat;
next unless $distance < $bestDistance;
# found better node
$startNode = $node;
$bestDistance = $distance;
}
die unless $startNode;
# start node found
warn "start node $startNode @ lat $lat{$startNode} lon $lon{$startNode}.\n";
print "start\t$lon{$startNode}\t$lat{$startNode}\n";
#-----------------------------------------------------------------------
# reachability analysis
#-----------------------------------------------------------------------
my %bestDistance; # distance from startnode to node
my %predecessor; # only needed if way from start to certain node to be displayed
$bestDistance{$startNode} = 0;
my @nodesToBeAnalyzed = ($startNode);
# as long as nodes available
while (@nodesToBeAnalyzed) {
# order by distance and get nearest, terminate if nearest outside limit
@nodesToBeAnalyzed = sort {$bestDistance{$a} <=> $bestDistance{$b}} @nodesToBeAnalyzed;
my $thisNode = $nodesToBeAnalyzed[0];
my $thisDist = $bestDistance{$thisNode};
last if $thisDist>$maxDist;
shift (@nodesToBeAnalyzed);
# check arcs
my $hasOuterNeighbour = 0;
foreach my $thatNode (@{$neighbours{$thisNode}}) {
# dist to thatNode via thisNode
my $viaDist = $thisDist + $arcLength{$thisNode, $thatNode};
# reachable?
if ($viaDist <= $maxDist) {
print "arc\t$lon{$thisNode}\t$lat{$thisNode}\t$lon{$thatNode}\t$lat{$thatNode}\n";
} else {
print "noarc\t$lon{$thisNode}\t$lat{$thisNode}\t$lon{$thatNode}\t$lat{$thatNode}\n";
$hasOuterNeighbour=1;
}
# already known distance of thatNode
my $thatDist = $bestDistance{$thatNode};
# if first touch of thatNode, add to todo list
push (@nodesToBeAnalyzed, $thatNode) unless defined $thatDist;
# remember best distance to thatNode
$bestDistance{$thatNode} = $viaDist if !defined $thatDist || $viaDist<$thatDist;
}
print "inner\t$lon{$thisNode}\t$lat{$thisNode}\n" if $hasOuterNeighbour;
}
# list nodes on todo list
foreach my $thisNode (@nodesToBeAnalyzed) {
print "outer\t$lon{$thisNode}\t$lat{$thisNode}\n";
}
# done
Wir rufen es mit der Startposition und maximalen Entfernung auf:
./reachable.pl 50.936 6.948 500 <graph >reachable.txt
und erhalten die Erreichbarkeitsanalyse “reachable.txt”.
Diese enthält den Startknoten “start”, gerade noch erreichbare Knoten “inner” und nicht mehr erreichbare Knoten “outer”, jeweils mit Angabe der Koordinaten:
start 6.9479894 50.9358304
inner 6.9510076 50.9330590
outer 6.9394997 50.9355044
Außerdem die erreichbaren Kanten “arc” und nur teilweise erreichbare Kanten “noarc” jeweils mit Anfangs- und Endkoordinaten:
arc 6.9479894 50.9358304 6.9470897 50.9358754
noarc 6.9510076 50.9330590 6.9525011 50.9332148
4. Wir machen das Ergebnis der Analyse sichtbar:
Dazu legen wir eine HTML-Datei “reachable.htm” (selbstverständlich auch überaus komplex!) an:
<html>
<head>
<title>Reachable</title>
<script type="text/javascript" src="http://www.openlayers.org/api/OpenLayers.js"></script>
<script type="text/javascript" src="http://www.openstreetmap.org/openlayers/OpenStreetMap.js"></script>
<script type="text/javascript">
var map;
window.onload = function () {
map = new OpenLayers.Map ('map');
map.displayProjection = new OpenLayers.Projection('EPSG:4326');
map.addLayer (new OpenLayers.Layer.OSM.Mapnik('Mapnik'));
var vectors;
map.addLayer (vectors = new OpenLayers.Layer.Vector ('Reachability'));
OpenLayers.Request.GET({
async: false,
url: 'reachable.txt',
success: function(request) {
var lines = (request.responseText||'').split('\n');
for (var i in lines) {
var f=lines[i].split('\t');
switch (f[0]) {
case 'arc':
p0 = new OpenLayers.Geometry.Point (f[1], f[2]).
transform(this.map.displayProjection, this.map.getProjectionObject());
p1 = new OpenLayers.Geometry.Point (f[3], f[4]).
transform(this.map.displayProjection, this.map.getProjectionObject());
geometry = new OpenLayers.Geometry.LineString ([p0, p1]);
style = {
strokeColor: 'blue',
strokeWidth: 4,
strokeOpacity: 0.5
}
feature = new OpenLayers.Feature.Vector(geometry, null, style);
vectors.addFeatures([feature]);
break;
case 'noarc':
p0 = new OpenLayers.Geometry.Point (f[1], f[2]).
transform(this.map.displayProjection, this.map.getProjectionObject());
p1 = new OpenLayers.Geometry.Point (f[3], f[4]).
transform(this.map.displayProjection, this.map.getProjectionObject());
geometry = new OpenLayers.Geometry.LineString ([p0, p1]);
style = {
strokeColor: 'red',
strokeWidth: 4,
strokeOpacity: 0.5
}
feature = new OpenLayers.Feature.Vector(geometry, null, style);
vectors.addFeatures([feature]);
break;
case 'start':
point = new OpenLayers.Geometry.Point (f[1], f[2]).
transform(this.map.displayProjection, this.map.getProjectionObject());
style = {
pointRadius: 5,
strokeColor: 'black',
strokeWidth: 1,
strokeOpacity: 1,
fillColor: 'yellow',
fillOpacity: 0.8
}
feature = new OpenLayers.Feature.Vector(point, null, style);
vectors.addFeatures([feature]);
break;
case 'inner':
point = new OpenLayers.Geometry.Point (f[1], f[2]).
transform(this.map.displayProjection, this.map.getProjectionObject());
style = {
pointRadius: 3,
strokeColor: 'blue',
strokeWidth: 1,
strokeOpacity: 1,
fillColor: 'blue',
fillOpacity: 0.5
}
feature = new OpenLayers.Feature.Vector(point, null, style);
vectors.addFeatures([feature]);
break;
case 'outer':
point = new OpenLayers.Geometry.Point (f[1], f[2]).
transform(this.map.displayProjection, this.map.getProjectionObject());
style = {
pointRadius: 3,
strokeColor: 'red',
strokeWidth: 1,
strokeOpacity: 1,
fillColor: 'red',
fillOpacity: 0.5
}
feature = new OpenLayers.Feature.Vector(point, null, style);
vectors.addFeatures([feature]);
break;
}
}
}
});
map.zoomToExtent (vectors.getDataExtent());
};
</script>
</head>
<body style="margin: 0">
<h1>Reachable</h1>
<div id="map" style="background: gray; width: 100%; height: 75%;"></div>
</body>
</html>
Wir legen HTML-Seite “reachable.htm” und Erreichbarkeits-Analyse “reachable.txt” in ein Verzeichnis und rufen die Webseite auf.
5. Wir erfreuen uns des Ergebnisses.
6. Nachbemerkung:
Diese Erreichbarkeitsanalyse wurde im Rahmen des Verzehrs eines Stückes Almkäse in Verbindung mit dem Genuss eines Glases Rotwein entwickelt. Der Source Code ist öffentlich zugänglich.
Plüschige Grüße
Der Assistent (orange)