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; } } }