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