Ciąg dlaszy wpisu Rejestrowanie wysokości przez odbiorniki GPS

Ciąg dalszy bo ,,ciąg technologiczny'' FIT →TCX →GPX →srtm.py →lepszy_GPX→Endomondo/Strava pozbawia niestety przesyłany plik GPX danych, które nie są wspierane przez format GPX (takie jak tętno dla przykładu).

Problem można rozwiązać dodając do pliku TCX poprawione wysokości z pliku lepszy_GPX:

#!/usr/bin/perl
use XML::LibXML;
use XML::LibXML::XPathContext;
use Getopt::Long;

binmode(STDOUT, ":utf8");

my $Gpx_file; my $Tcx_file;

## gpx -- plik GPX z poprawionymi wysokościami
## tcx -- oryginalny plik TCX (w którym mają być poprawione wysokości)
GetOptions( "gpx=s" => \$Gpx_file, "tcx=s" => \$Tcx_file, ) ;

my $parser = XML::LibXML->new();

for my $file2parse ("$Gpx_file") {

  my $doc = $parser->parse_file($file2parse);
  my $root = $doc->documentElement();

  my $xc = XML::LibXML::XPathContext->new( $root );
  $xc->registerNs('gpx', 'http://www.topografix.com/GPX/1/0');

  foreach my $t ($xc->findnodes('//gpx:trkpt')) {

    my $xpe = XML::LibXML::XPathContext->new( $t );
    $xpe->registerNs('gpx', 'http://www.topografix.com/GPX/1/0');
    $node++;
    $gmttime = ($xpe->findnodes('gpx:time'))[0]->textContent();
    $altitude = ($xpe->findnodes('gpx:ele'))[0]->textContent();
    $GpxPoints{"$gmttime"} = $altitude;
  }
}

## for $e (keys %GpxPoints) { print "$e => $GpxPoints{$e}\n"; }

for my $file2parse ("$Tcx_file") {
  # http://www.garmin.com/xmlschemas/TrainingCenterDatabase/v2
  my $doc = $parser->parse_file($file2parse);

  my $root = $doc->documentElement();

  my $xc = XML::LibXML::XPathContext->new( $root );
  $xc->registerNs('tcx',
      'http://www.garmin.com/xmlschemas/TrainingCenterDatabase/v2');

  foreach my $t ($xc->findnodes('//tcx:Trackpoint')) {

    my $xpe = XML::LibXML::XPathContext->new( $t );
    $xpe->registerNs('tcx',
       'http://www.garmin.com/xmlschemas/TrainingCenterDatabase/v2');

    $gmttime = ($xpe->findnodes('tcx:Time'))[0]->textContent();

    foreach my $xpa ( $xpe->findnodes('tcx:AltitudeMeters') ) {
       $oldAltitude = $xpa->textContent();
       ## jeżeli istnieje w pliku GPX punkt z tym samym stemplem 
       ## czasu zmień zawartość elementu tcx:AltitudeMeters
       if (exists $GpxPoints{$gmttime} ) {
          ## replace content of tcx:AltitudeMeters
          $xpa->removeChildNodes();
          $xpa->appendText("$GpxPoints{$gmttime}");
          $changedNodesNo++;
       } else {
          ## jeżeli nie istnieje usuń cały węzeł
          my $parent = $t->parentNode;
          $parent->removeChild( $t );
          $droppedNodes .= "$gmttime;";
          $droppedNodesNo++;
       }
    }
  }

### ###
print "<?xml version='1.0' encoding='UTF-8' standalone='no' ?>\n";
print $root->toString;
### ###
}

print STDERR "Zmieniono: $changedNodesNo / Usunięto: $droppedNodesNo\n";
print STDERR "($droppedNodes)\n";

Ponieważ skrypt srtm.py nie tworzy pliku GPX zawierającego wszystkie punkty z pliku TCX (z jakiś powodów niewielka część jest pomijana); warunek exists $GpxPoints{$gmttime} sprawdza czy istnieje w pliku GPX punkt z podanym stemplem czasowym. Jeżeli istnieje zmieniana jest wysokość, jeżeli nie to punkt jest usuwany.

Rejestrowanie wysokości przez odbiorniki GPS jak wiadomo odbywa się z błędem. Do tej pory mi to wisiało, ale problem stanął, gdy zachciało mi się dodać nachylenie do danych GPS (wyświetlanych jako napisy generowane autorskimi skryptami z pliku GPX/TCX), które to nachylenie nie zgadzało się wielokrotnie ze stanem faktycznym. Szczególnie tajemnicze były różnice na plus/minus 300% i więcej...

W trakcie prób rozwiązania tego problemu dowiedziałem się tego i owego na temat rejestracji wysokości. I o tym jest ten wpis. Natomiast problem z obliczaniem nachylenia został nierozwiązany. Wartości nachylenia pokazywane w trakcie jazdy przez urządzenia, takie jak Garmin Edge 500 są wiarygodnie, co by świadczyło, że zaimplementowano w nich jakiś całkiem sprytny algorytm wygładzający. Szukałem co na ten temat wie google, ale nic nie znalazłem. Próbowałem wygładzać wysokość/nachylenie w R za pomocą funkcji loess (a także średniej ruchomej) -- rezultaty były kiepskie.

Znalazłem natomiast informacje w jaki sposób można poprawić dokładność danych o wysokości, zarejestrowaną (niedokładnie) przez urządzenie GPS. Otóż można albo skorzystać z usługi Google Elevation (GE) albo użyć danych SRTM, przy czym minusem GE jest dzienny limit 2500 zapytań.

Zaczniejmy od prostszego przypadku, tj. dodania/zmiany wysokości z modelu SRTM. W tym celu pobieramy/instalujemy pakiet srtm.py. Pakiet zawiera m.in narzędzie pn. gpxelevations:

## Dodaje dane SRTM, zapisuje do pliku PLIK_SRTMS.gpx
gpxelevations -s -o PLIK.gpx -f PLIK_SRTMS.gpx
  
## Wygładza dane i zapisuje do pliku PLIK_SRTMS.gpx
gpxelevations -o 20160511.gpx -f PLIK_SRTMS.gpx

Dane SRTM można też dodać/zamienić posługując się okienkową aplikacją pn. GPSPrune, jak ktoś lubi klikać ale nie lubi wiersza poleceń.

Zakładając, że dane są w pliku GPX pobranym z urządzenia, poniższy skrypt wygeneruje plik CSV zawierający m.in. wysokość oryginalną, wysokość z modelu SRTM oraz wysokość wygładzoną:

#!/bin/bash

FILE=`basename $1 .gpx`

if [ -f "$FILE.gpx" ] ; then
  ## Dodanie lepszych wysokości
  gpxelevations -o $FILE.gpx -f ${FILE}_SRTM.gpx
  gpxelevations -s -o $FILE.gpx -f ${FILE}_SRTMS.gpx

  ## Zamiana na CSV (skrypt gpx2csv.pl zamieszczono dalej)
  gpx2csv.pl ${FILE}.gpx > ${FILE}.csv && \
  gpx2csv.pl ${FILE}_SRTM.gpx > ${FILE}_SRTM.csv && \
  gpx2csv.pl ${FILE}_SRTMS.gpx > ${FILE}_SRTMS.csv

  ## Scalenie w jeden plik CSV
  paste -d ';' ${FILE}.csv ${FILE}_SRTM.csv ${FILE}_SRTMS.csv | \
  awk -F';' ' BEGIN{ print "daytime;lat;long;ele;srtm;srtms"; };\
     {print $1 ";" $2 ";" $3 ";" $4 ";" $8 ";" $12}' > ${FILE}_ALLE.csv

else
  echo "*** USAGE: $0 PLIK.gpx ***"
fi

Teraz można porównać wysokość oryginalną, wysokość SRTM oraz wygładzoną na wykresie:

library(reshape)
require(ggplot2)

args = commandArgs(trailingOnly = TRUE);

if (length(args)==0) { stop("Podaj nazwę pliku CSV", call.=FALSE) }

fileBase <- gsub(".csv", "", args[1]);
outFile <- paste (fileBase, ".pdf", sep = "");

d <- read.csv(args[1], sep = ';',  header=T, na.string="NA");

ggplot(d, aes(x = as.POSIXct(daytime, format="%Y-%m-%dT%H:%M:%SZ"))) +
  geom_line(aes(y = ele, colour = 'zarejestrowana', group = 1), size=.5) +
  geom_line(aes(y = srtm, colour = 'srtm', group = 1), size=.5) +
  geom_line(aes(y = srtms, colour = "srtm (wygładzona)", group = 1), size=.5) +
  ylab(label="Wysokość [mnpm]") +
  xlab(label="czas") +
  labs(colour = paste("Plik:", fileBase )) +
  theme(legend.position="top") +
  theme(legend.text=element_text(size=12));

ggsave(file=outFile, width=12, height=5)

## Uruchomienie:
## Rscript gps_vs_srtm.R PLIK.csv

Oryginalne dane z urządzenia są systematycznie zaniżone.

Dodanie danych z Google Elevation jest równie proste. Poniższy skrypt -- jako przykład -- pobiera wysokość dla punktu o współrzędnych podanych jako argumenty wywołania (szerokość/długość):

use LWP::Simple; 
use JSON qw( decode_json );

my $geocodeapi =
 "https://maps.googleapis.com/maps/api/elevation/json?locations";

 ## szerokość = $ARGV[0] ; długość = $ARGV[1]
 my $url = $geocodeapi . "=$ARGV[0],$ARGV[1]";
 my $json = get($url);

 my $d_json = decode_json( $json );

 if ( $d_json->{status} eq 'OK' ) {
   $elevationG = $d_json->{results}->[0]->{elevation};
   $resolutionG = $d_json->{results}->[0]->{resolution};
   $latG = $d_json->{results}->[0]->{location}->{lat};
   $lngG = $d_json->{results}->[0]->{location}->{lng};
 }

 print "$elevationG\n";

Drugi z wykresów zawiera dane pobrane z Google Elevation Service.

Na zakończenie jeszcze skrypt zamieniający gpx na csv:

#!/usr/bin/perl
# Wykorzystanie gpx2csv.pl PLIK.gpx > PLIK.csv
#
use XML::DOM;
my $parser = new XML::DOM::Parser;

for my $file2parse (@ARGV) {
  my $doc = $parser->parsefile ($file2parse);
  for my $t ( $doc->getElementsByTagName ("trk") ) {

    for my $p ( $t->getElementsByTagName ("trkpt") ) {  $node++;
        my $latitude = $p->getAttributeNode ("lat")->getValue() ;
        my $longitude = $p->getAttributeNode ("lon")->getValue() ;

	$gmttime = $p->getElementsByTagName ("time")->
	   item(0)->getFirstChild->getNodeValue();
	
	 $altitude = $p->getElementsByTagName ("ele")->
	   item(0)->getFirstChild->getNodeValue();
	printf "%s;%s;%s;%s\n", $gmttime, $latitude, $longitude, $altitude;
       }
   }
}