Entfernung zwischen zwei GPS Koordinaten
Dominik Schmidt
Meine kleine GPS Simulation kommt langsam vorran. Ich kann schon die Daten
einlesen und in 3D darstellen (und drehen, etc). Das ganze schaut in etwa so
aus: http://www.logview.info/phpBB2/download.php?id=25 (hoffe der Link funzt
so ...).
Nun würde ich gerne die Gesamtstrecke und die Strecken zwischen jeweils zwei
Punkten berechnen. Mein Mathe in der Hinsicht ist aber etwas angestaubt und
deswegen habe ich gestern mal ein bisserl im Web gesucht. Gefunden habe ich
auch 2 Funktionen, aber jeweils mit dem Kommentar, dass das bei langen bzw.
kurzen Entfernungen nicht sauber funktioniert.
Deswegen mal die Frage hier ... Kann mir jemand eine Formel geben, die
universell einsatzbar ist und immer das richtige Ergebnis (=Entfernung)
liefert? Der Input wäre also z.B.:
N46.38.3242
E14.17.6446
N46.38.3250
E14.17.6470
Und die Funktion liefert dann 299m als Beispiel.
Und als weitere Frage ... Ich möchte die Daten eine Modellflugzeugs
visualisieren. Da ist es ja nun nicht ganz unerheblich, wie die aktuelle
Höhe ist. Also bräuchte ich als Erweiterung zur ersten Formel noch eine
Ergänzung und zwar mit Höhe. Der Input wäre dann z.B.:
N46.38.3242
E14.17.6446
123m
N46.38.3250
E14.17.6470
150m
Und das Ergebnis z.B. 400m.
Und last but not least ... Wie ist eurer Meinung nach die Höhenangabe des
GPS zu werten? Ist die Angabe halbwegs verlässlich, oder nicht? Da wir die
Höhe auch noch zusätzlich über einen Drucksensor erfassen, könnte ich da
evtl. Fehler ausgleichen.
So denn schon mal Danke!
Greetz Dominik
Fredie Kern
zur Streckenberechnung schlage ich Dir folgenden Weg vor.
1. Umrechnung der geographische (ellipsoidischen) Koordinaten
(N46.38.3242, E14.17.6470, h=150m) in kartesische (geozentrische) X,Y,Z
für beide Endpunkte der Strecke. Falls keine Höhenangabe verfügbar
einfach h=0 setzen; dann erhälts Du die Entfernung auf Meereshöhe.
2. Strecke aus X,Y,Z nach Pythagoras berechnen. s=\sqrt{X^2 + Y^2 + Z^2}
Gruß,
Fredie
Dominik Schmidt
news:3mdjvdF...@individual.net...
> 1. Umrechnung der geographische (ellipsoidischen) Koordinaten
> (N46.38.3242, E14.17.6470, h=150m) in kartesische (geozentrische) X,Y,Z
Gibt es da im Web irgendwo Infos wie die Umrechnung von statten geht? Habe
leider nix finden können.
Greetz Dominik
Gabriel Ebner
> Nun würde ich gerne die Gesamtstrecke und die Strecken zwischen jeweils
> zwei Punkten berechnen. Mein Mathe in der Hinsicht ist aber etwas
> angestaubt und deswegen habe ich gestern mal ein bisserl im Web gesucht.
> Gefunden habe ich auch 2 Funktionen, aber jeweils mit dem Kommentar, dass
> das bei langen bzw. kurzen Entfernungen nicht sauber funktioniert.
Besser die mit dem Arcussinus verwenden, auf große Entfernungen fallen
kleine Ungenauigkeiten nicht so auf.
http://en.wikipedia.org/wiki/Great_circle_distance
Gabriel.
Gabriel Ebner
v = a / sqrt(1 - e^2 * sin(lat)^2)
x = (v + h) * cos(lat) * cos(lon)
y = (v + h) * cos(lat) * sin(lon)
z = ((1 - e^2) * v + h) * sin(lat)
a und e^2 dem entsprechenden Ellipsoid anpassen (e^2 = 2*f - f^2)
Für WGS84 z.B:
a = 6378137
f = 1/298.257223563
e^2 = 2*f - f^2 = 0.00669437999014
Aber Achtung: die kartesische Entfernung ist mit Vorsicht zu genießen. Nur
für kurze Strecken mit Sichtverbindung brauchbar.
CK - Vertreibs & Handels GmbH
Dein Proggy sieht ja super aus! Ein konkretes Tool hab' ich leider nicht zur
Stelle, aber sieh Dir mal die site:
http://www.luftpiraten.de/navigation.html an, super Infos! Wenn Du die
Entfernung bezogen auf die Erdoberfläche haben willst, berechnest Du die
Bogenlänge b über dem Winkel alpha Dabei ist der Winkel alpha die Differenz
der Koordinaten (auf das Winkelformat achten! Minuten, Sekunden evtl.von
hddd°mm,mmmm' auf hddd,ddddd° umrechnen!) Radius r ist dabei 6371 km
(Erdradius).
b = PI * r * alpha / 180
erst einmal umrechnen (hddd°mm,mmmm' auf hddd,ddddd°):
Latitude1: N46°38,3242' = N46,63874°
Longitude 1: E14°17,6446' = E014,29408°
Latitude2: N46°38,3250' = N46,63875°
Longitude 2: E14°17,6470' = E014,29412°
diff. Latitude: 46,63874° - 46,63875° = -0,00001° (Differenz
geographische Breite, je mehr Nachkommastellen, umso besser :-) das
Vorzeichen ist nur wichtig, wenn Du auch die Richtung brauchst, falls Du
weiter als bis zum Äquator fliegst, geht auch Nord und Süd ein)
diff. Longitude: 014,29408° - 014,29412° = -0,00005° (Differenz
geographische Länge, je mehr Nachkommastellen, umso besser :-) das
Vorzeichen ist nur wichtig, wenn Du auch die Richtung brauchst, fliegst Du
westlicher als Greenwich oder bis hinter die Datumsgrenze geht sinngemäß Ost
und West mit ein)
alpha = sqr ( 0,00001° ^2 + 0,00005° ^2 ) = 0,00005099°
b = PI * 6371 km * 0,00005099° / 180 ~ 5 m
Mit der Höhenangabe des gps wäre ich vorsichtig (einerseits ist die
Genauigkeit nicht so berauschend, und Du musst nachsehen ob in deiner
Höhenangabe die geoidale Separation (Abplattung des Geoids) berücksichtigt
ist), bringt auch bei den Flughöhen nicht so viel.
Ich hab' einen gps-Datenlogger und einen gps-gestützten Autopiloten
programmiert, allerdings in PowerBasic (damits auch unter DOS oder WIN98SE
läuft). Das ganze bootet auf einem VIA nano-itx board von einem USB-Stick,
die Servosteuerung soll (leider noch nicht fertig) über den seriellen
Anschluss geschehen.
Viele Grüße
Mathias
Dominik Schmidt
> Dein Proggy sieht ja super aus!
Danke. Aber wie gesagt, da muss ich noch einiges dran pfeilen :-)
Vielen Dank für dein Beispiel. Ich denke damit werde ich es umsetzen können.
Die Frage ist nun ...
Wie genau ist die Berechnung? Geht die Formel sowohl mit kleinen wie auch
langen Distanzen?
>geoidale Separation
Hmm, was ist das nun wieder? Wo ersehe ich das im GPS Datenformat?
Greetz Dominik
PS @Mathias: Da du scheinbar tiefer in dem Thema drin bist ... Können wir
uns ein bisschen per Mail austauschen? Evtl. kannste ja auch mal mein Proggy
testen!?
Gabriel Ebner
> b = PI * r * alpha / 180
> alpha = sqr ( 0,00001° ^2 + 0,00005° ^2 ) = 0,00005099°
> b = PI * 6371 km * 0,00005099° / 180 ~ 5 m
Nordpol (N90, E0) - Nordpol (N90, E180):
alpha = sqrt(0°^2 + 180°^2) = 180°
b = pi * 6317km * 180° / 180 = pi * 6317km ~ 19845km
Ähm?
Gabriel.
Gabriel Ebner
> Wie genau ist die Berechnung? Geht die Formel sowohl mit kleinen wie auch
> langen Distanzen?
Je weiter die Punkte im Norden sind, und je weiter sie von einander entfernt
sind, desto ungenauer ist diese Formel.
>>geoidale Separation
> Hmm, was ist das nun wieder? Wo ersehe ich das im GPS Datenformat?
Im GPGGA-Datensatz, siehe
http://www.kowoma.de/gps/zusatzerklaerungen/NMEA.htm
Das ist Abweichung des Geoids (mittlerer Meeresspiegel) vom idealisierten
Rotationsellipsoid.
Gabriel.
Mathias Schröder
die geoidale Separation ist die Differenz zwischen dem Erdellipsoid (wegen
der Abplattung an den Polen) und dem MainSeaLevel (NormalNull). Wird im
NMEA183-Protokoll im $GPGGA-Satz mit ausgegeben. Üblicherweise geben die
gps-Geräte aber die Höhe über NN aus. An einem Austausch per mail hab' ich
auch Interesse, dann aber unter der mail: buerosc...@gmx.de , ist
nämlich meine richtige addy.
Viele Grüße
Mathias
Mathias Schröder
in BASIC sieht die Umrechnung so aus:
PI##=4*ATN(1)
ERDRADIUS=6371 'Erdradius in km.
LAT1##=14.29408
LON1##=46.63874
LAT2##=14.29412
LON2##=46.63875
DIFFLAT##=LAT2##-LAT1##
DIFFLON##=LON2##-LON1##
DIFFWIN##=SQR(DIFFLAT##^2+DIFFLON##^2)
DISTANZ#=PI##*ERDRADIUS/180*DIFFWIN##
ALFA##=ATN(-DIFFLON##/DIFFLAT##)
PRINT LAT1##,LON1##
PRINT LAT2##,LON2##
PRINT DIFFLAT##,DIFFLON##,DIFFWIN##
PRINT DISTANZ#;"km"
PRINT DISTANZ#/1.852;"nm"
PRINT 90+ALFA##*180/PI##;"grad (wahr)"
nicht schön (weil nur Textausgabe) aber funzt :-)
Viele Grüße
Mathias
Dominik Schmidt
> auch Interesse, dann aber unter der mail: buerosc...@gmx.de , ist
Ich melde mich gleich von @home :-)
Greetz
Mathias Schröder
Nordpol (N90, E0) - Nordpol (N90, E180):
alpha = sqrt(0°^2 + 180°^2) = 180°
b = pi * 6317km * 180° / 180 = pi * 6317km ~ 19845km
ist doch korrekt :-) , der Erdumfang ist überschläglich 40.000 km, also
braucht man von Pol zu Pol ca. 20.000 km. Berechnet wird mit der Formel der
Kreisbogen, also der Weg auf der Erdoberfläche (NN).
Viele Grüße
Mathias
Mathias Schröder
ich hoffe, ich nerve Euch jetzt nicht mit meinem BASIC-Geschreibsel, aber
wen's interessiert. Datenmitschnitt vom COM-Port 1 mit 4800 baud, ausgelesen
werden: $GPRMC, $GPGSA und $GPGGA aus NMEA183.
---------------------------------------------------------------------------------------------------------
DIM KOMMARMC(20)
DIM KOMMAGSA(20)
DIM KOMMAGGA(20)
DIM GPRMC$(20)
DIM GPGSA$(20)
DIM GPGGA$(20)
CLS
$COM 16384
OPEN "COM1:4800,n,8,1,DS,RS,CS,CD,PE,ME,FE" AS #1
LOGDATEINAME$=LEFT$(DATE$,2)+MID$(DATE$,4,2)+LEFT$(TIME$,2)+MID$(TIME$,4,2)+".LOG"
OPEN LOGDATEINAME$ FOR OUTPUT AS #2
WHILE NOT INSTAT
IF LOC(1) > 0 THEN
ComPortInput$ = INPUT$(LOC(1),#1)
FOR i=1 TO LEN(ComPortInput$)
PortInput$=MID$(ComPortInput$,i,1)
Eingang$=Eingang$+PortInput$
IF PortInput$="$" THEN
IF LEN(Eingang$)>3 THEN
IF LEFT$(Eingang$,5)="GPRMC" THEN
GPRMC$=LEFT$(Eingang$,LEN(Eingang$)-3)
IF LEFT$(Eingang$,5)="GPGSA" THEN
GPGSA$=LEFT$(Eingang$,LEN(Eingang$)-3)
IF LEFT$(Eingang$,5)="GPGGA" THEN
GPGGA$=LEFT$(Eingang$,LEN(Eingang$)-3)
END IF
Eingang$=""
END IF
NEXT i
END IF
j=1
LENGPRMC=LEN(GPRMC$)
FOR i=1 to LENGPRMC
IF MID$(GPRMC$,i,1)="," THEN KOMMARMC(j)=i : j=j+1
NEXT i
FOR i=1 TO j-1
IF KOMMARMC(i+1)>0 THEN
GPRMC$(i)=MID$(GPRMC$,KOMMARMC(i)+1,KOMMARMC(i+1)-KOMMARMC(i)-1)
NEXT i
j=1
LENGPGSA=LEN(GPGSA$)
FOR i=1 to LENGPGSA
IF MID$(GPGSA$,i,1)="," THEN KOMMAGSA(j)=i : j=j+1
NEXT i
FOR i=1 TO j-1
IF KOMMAGSA(i+1)>0 THEN
GPGSA$(i)=MID$(GPGSA$,KOMMAGSA(i)+1,KOMMAGSA(i+1)-KOMMAGSA(i)-1)
NEXT i
j=1
LENGPGGA=LEN(GPGGA$)
FOR i=1 to LENGPGGA
IF MID$(GPGGA$,i,1)="," THEN KOMMAGGA(j)=i : j=j+1
NEXT i
FOR i=1 TO j-1
IF KOMMAGGA(i+1)>0 THEN
GPGGA$(i)=MID$(GPGGA$,KOMMAGGA(i)+1,KOMMAGGA(i+1)-KOMMAGGA(i)-1)
NEXT i
UHRZEIT$=GPRMC$(1)
UHRZEITSTUNDE$=LEFT$(UHRZEIT$,2)
UHRZEITSTUNDE%=VAL(UHRZEITSTUNDE$)
UHRZEITMINUTE$=MID$(UHRZEIT$,3,2)
UHRZEITMINUTE%=VAL(UHRZEITMINUTE$)
UHRZEITSEKUNDE$=MID$(UHRZEIT$,5,2)
UHRZEITSEKUNDE%=VAL(UHRZEITSEKUNDE$)
UHRZEIT$=UHRZEITSTUNDE$+":"+UHRZEITMINUTE$+":"+UHRZEITSEKUNDE$
DATUM$=GPRMC$(9)
DATUMJAHR$="20"+MID$(DATUM$,5,2)
DATUMJAHR%=VAL(DATUMJAHR$)
DATUMMONAT$=MID$(DATUM$,3,2)
DATUMMONAT%=VAL(DATUMMONAT$)
DATUMTAG$=LEFT$(DATUM$,2)
DATUMTAG%=VAL(DATUMTAG$)
DATUM$=DATUMMONAT$+"-"+DATUMTAG$+"-"+DATUMJAHR$
IF GPRMC$(2)="A" THEN DATASTATUS$="ok"
IF GPRMC$(2)="V" THEN DATASTATUS$="invalid"
SATZAHL$=GPGGA$(7)
SATZAHL%=VAL(SATZAHL$)
IF GPGSA$(1)="A" THEN MODUS1$="automatisch 2D/3D"
IF GPGSA$(1)="M" THEN MODUS1$="manuell 2D/3D"
IF GPGSA$(2)="1" THEN MODUS2$="Fix not available"
IF GPGSA$(2)="2" THEN MODUS2$="2D"
IF GPGSA$(2)="3" THEN MODUS2$="3D"
LATIT1$=GPRMC$(3)
LATITUDE#=ROUND(VAL(LEFT$(LATIT1$,2))+VAL(MID$(LATIT1$,3,7))/60,7)
LATIT2$=GPRMC$(4)
IF LATIT2$="S" THEN LATITUDE#=-LATITUDE#
LONGI1$=GPRMC$(5)
LONGITUDE#=ROUND(VAL(LEFT$(LONGI1$,3))+VAL(MID$(LONGI1$,4,7))/60,7)
LONGI2$=GPRMC$(6)
IF LONGI2$="W" THEN LONGITUDE#=-LONGITUDE#
HOEHE1$=GPGGA$(9)
HOEHE2$=GPGGA$(10)
HOEHE#=VAL(HOEHE1$)
HOEHE3$=GPGGA$(11)
HOEHE4$=GPGGA$(12)
SEPARATION#=VAL(HOEHE3$)
SPEED$=GPRMC$(7)
SPEED#=VAL(SPEED$)*1.852
COURSE$=GPRMC$(8)
KURS#=VAL(COURSE$)
LOCATE 1,40 : PRINT "Uhrzeit : ";UHRZEIT$
IF UHRZEITALT$<>UHRZEIT$ AND SATZAHL%>3 THEN
SOUND 4444,.1
WRITE #2
,DATUMJAHR%,DATUMMONAT%,DATUMTAG%,UHRZEITSTUNDE%,UHRZEITMINUTE%,UHRZEITSEKUNDE%,SATZAHL%,LATITUDE#,LONGITUDE#,SPEED#,KURS#,HOEHE#
END IF
UHRZEITALT$=UHRZEIT$
LOCATE 1,1 : PRINT "Datum : ";DATUM$
LOCATE 2,1 : PRINT "Data status : ";DATASTATUS$;" "
LOCATE 2,40 : PRINT "Satelitenzahl :";SATZAHL%
LOCATE 3,1 : PRINT "Umschaltmodus : ";MODUS1$;" "
LOCATE 3,40 : PRINT "Navigation : ";MODUS2$;" "
LOCATE 4,1 : PRINT "Latitude :";LATITUDE#;" "
LOCATE 4,40 : PRINT "Longitude :";LONGITUDE#;" "
LOCATE 5,1 : PRINT "Hoehe ue. NN :";HOEHE#;HOEHE2$+" "
LOCATE 5,40 : PRINT "Geoidal sep. :";SEPARATION#;HOEHE4$+" "
LOCATE 6,1 : PRINT "Speed :";SPEED#;" Km/h "
LOCATE 6,40 : PRINT "Kurs :";KURS#;" Grad "
LOCATE 8,1
FOR i=1 TO 16
PRINT i,GPRMC$(i),GPGSA$(i),GPGGA$(i)
NEXT i
WEND
CLOSE 1
CLOSE 2
IF DATUM$<>"" THEN DATE$=DATUM$
IF UHRZEIT$<>"" THEN TIME$=UHRZEIT$
CLS
END
---------------------------------------------------------------------------------------------------------
Ist für PowerBASIC geschrieben, müsste aber auch mit kleinen Änderungen
unter QBASIC oder Ähnlichem laufen. Vorsicht, am Ende setzt das Programm die
Rechneruhr auf UTC!
Viele Grüße
Mathias
Gabriel Ebner
Das Problem ist nur, dass dies hier die Entfernung vom NORDpol zum NORDpol
ist, die eigentlich 0 sein sollte.
Gabriel.
Mathias Schröder
dann rechnest Du einfach mit:
Nordpol (N90, E0) - Nordpol (N90, E180):
alpha = sqrt(0°^2 + 360°^2) = 360°
b = pi * 6317km * 360° / 180 = pi * 6317km ~ 39690 km
und es stimmt auch wieder ;-)
Viele Grüße
Mathias
Mathias Schröder
sorry, das mit NORDpol-NORDpol hatte ich überlesen, bei den Breitengraden
kann man Nord und Süd vorzeichenbehaftet einsetzen, sinngemäß muss man
östliche und westliche Länge in der Berechnung mit Vorzeichen versehen.
Viele Grüße
Mathias