Coder | DevOp | Admin

Überschneidungslängen GPX Tracks

In einem touristischen Wegenetz gibt es benannte Routen von irgendwelchen a nach irgendwelchen b. Diese Wege wurden einzeln als gpx Track mit einem online GIS Editor und seiner autotrace Funktion erstellt und dann in einer Datei mit rund 40 tracks bereitgestellt. Aufgabe war es herauszufinden:

• Gesamtlänge in m des Netzes und der Routen
• Länge der Abschnitte, die einfach, doppelt, dreifach oder gar vierfach von Routen belegt sind

Diese Aufgabe ist im allgemeinen Fall nicht ganz einfach, da jede Route grundsätzlich aus einer beliebigen Punktfolge bestehen kann. Die Formulierung der Bedingungen ist schon nicht einfach. Man könnte formulieren, das die Punkte auf der Polygonfläche der Wege der Routen liegen muss (ca. 3-7m breit). Alternativ in geringem Abstand zueinander.

Geht es einfacher?

Die eher interessierte Frage nach dem "Wie habt ihr die Tracks erstellt?" brachte den Hinweis darauf, das die Daten mit einem auf OpenSource Software und OpenStreetmap Daten basierenden Autotracer online erzeugt wurden.

Dies führte mich zu der Hypothese, das die Daten der überlappenden Abschnitte auch mit gleichen Koordinaten dargestellt sein könnten, da zugrundeliegende Daten und vermutlich auch der Traceralgorithmus gleich oder mindestens ähnlich sind. Ein paar Tests verifizierten diese Hypothese. Beide Systeme ergaben gleiche Punktkoordinaten für die gleichen Abschnitte Ein Unterschied gab es in der Genauigkeit der Daten, mit 6 bzw. 7 Nachkommastellen waren die Koordinaten beschrieben.

Wie kann man das algorithmisch verwenden?

Mit den Koordinaten können die gleichen Punkte und Segmente identifiziert werden. Weil wir uns in einem kleinen lokalen Rahmen bewegen, würden grundsätzlich auch die gerundeten 6 Nachkommastellen reichen.

<trkseg>
    <trkpt lat="51.2543327" lon="11.6689334" />
    ...
</trkseg>

=> idPunkt : "51.25433311.668933"

Die Längen der Segmente können dann berechnet und in Gruppen aufsummiert werden.

Identitäten

Die Identität eines Punktes ist damit durch die Koordinaten gegeben, die Identität eines Abschnittes durch zwei bezogene Punkte. Eine Richtung ist durch die Punktreihenfolge definiert.

Koordinaten als Hash Key

Mit dem Aufbau eines Hashes, der als key für die Elemente die normierten, auf 6 Stellen gerundeten Koordinaten verwendet, kann man einfach die Segmente finden und zählen. Über die Endpunktkoordinaten kann die Länge berechnet werden.

Code

#! /usr/bin/perl
# (c) 2018 Ervin Peters
use strict;
use warnings;
use diagnostics;

use utf8;
use Math::Trig;
use XML::LibXML::Simple qw(XMLin);

use Encode;
use Term::Encoding qw(term_encoding);
my $encoding = term_encoding;
binmode STDIN,":encoding($encoding)";
binmode STDOUT, ":encoding($encoding)";

my $lines = {};
my $fileName = $1 || './20180508-test.gpx';

sub calcDist {
    my ( $self, $other ) = @_;
    my $dlong = $self->{lon} - $other->{lon};
    my $dlat  = $self->{lat} - $other->{lat};
    my $dx = ( 111120 * cos( deg2rad( $self->{lat} 
        + $dlat / 2 ) ) ) * $dlong;
    my $dy = $dlat * 111120;
    return sqrt( $dx * $dx + $dy * $dy );
}

sub checkLine{
    my ($p1,$p2) = @_;
    my $lineKey = substr($p1->{lat},0,10) 
        . substr($p1->{lon},0,10)
        . substr($p2->{lat},0,10)
        . substr($p2->{lon},0,10);

    if (defined($lines->{$lineKey})){
        # segment gibt es schon, hochzählen
        $lines->{$lineKey}->{cnt}++;
    }else{
        # entgegengesetzte Richtung
        my $lineKeyR = $p2->{lat} 
            . $p2->{lon} 
            . $p1->{lat} 
            . $p1->{lon};

        if (defined($lines->{$lineKeyR})){
            # segment gibt es schon, hochzählen
            $lines->{$lineKeyR}->{cnt}++;
        }else{
            $lines->{$lineKey} = {
                cnt =>1,
                length => calcDist($p1,$p2)
            };
        }
    }
}

my $trkLen ={};

sub checkPoints{
    my ($points,$trkName) = @_;

    my $p2;
    foreach my $p1 (@{$points}){
        if (defined($p2)){
            checkLine($p1,$p2);
            $trkLen->{$trkName} += calcDist($p1,$p2);
        }
        $p2 = $p1;
    }
}

my $gpx = XMLin($fileName);

while (my ($key,$val) = each %{$gpx}) {
  if ($key eq 'trk'){
    while(my($trkName,$trkData) = each %{$val}){
      $trkLen->{$trkName} = 0;
      while (my ($trkDataKey, $trkDataValue) 
                    = each %{$trkData}) {
        if ($trkDataKey eq 'trkseg'){
          if (ref($trkDataValue) eq 'HASH'){
            while (my ($segItemKey, $segItemData) 
                      = each %{$trkDataValue}) {
              if ($segItemKey eq 'trkpt'){
                checkPoints($segItemData,$trkName);
              }
            }
          }elsif (ref($trkDataValue) eq 'ARRAY'){
            foreach my $trkSeg (@{$trkDataValue}){
              while (my ($segItemKey, $segItemData) 
                        = each %{$trkSeg}) {
                if ($segItemKey eq 'trkpt'){
                  checkPoints($segItemData,$trkName);
                }
              }
            }
          }
        }
      }
    }
  }
}

my @len = (0,0,0,0,0,0);

while (my ($seg,$data) = each %{$lines}){
    $len[$data->{cnt}] += $data->{length};
}

while (my ($trk,$length) = each %{$trkLen}){
    $length = int($length);
    print "$trk : $length m\n";
}
print "\nBelegung:\n";
for (my $i = 1; $i < 6; $i++){
    my $length = int($len[$i]);
    print "$i-fach belegt : $length m\n" if ($length > 0);
}
print "fertig\n";
                        

Download: Überschneidungslängen.pl
Die Ausgabe sieht dann für meine Daten so aus:

muli $ ./calcs.pl
01. Verbindungsweg1 : 363 m
...
45 Verindungsweg12 : 8534 m

Belegung:
1-fach belegt : 199575 m
2-fach belegt : 14204 m
3-fach belegt : 490 m
4-fach belegt : 166 m
fertig            
muli $
                    

Dies ist die Untergrenze der Werte der Überlappungen, da der Autotracer gelegentlich nicht den Weg findet und manuell eine Brücke eingebaut wurde.