Schnittpunkt zweier Orthodrome?

muza · 15 Dez. 2008

Servus!

Bin neu hier im Forum, habe aber direkt ein Problem - vielleicht kann mir ein kluger Kopf weiterhelfen. :gott:

Gegeben habe ich 4 Punkte A,B,C,D auf der Erdkugel (die Genauigkeit auf der Kugel reicht hier allemal) in Geokoordinaten, von denen immer 2 durch eine Strecke miteinander verbunden sind (A mit B und C mit D).
Die beiden Strecken AB und CD sind jeweils Orthodrome, liegen also beide auf einem Großkreis der Erdkugel.

Nun muss ich wissen, ob es einen Schnittpunkt der beiden Orthodrome gibt.

Habe mir überlegt, dass man die beiden Großkreise erhält, auf denen die Orthodrome liegen, indem man eine Ebene, die durch die beiden Endpunkte der Orthodrome und den Erdmittelpunkt geht, mit der Kugeloberfläche schneidet.
Problem ist allerdings, dass sich die Großkreise immer in zwei Punkten schneiden werden und ich nicht weiss, ob sie auf den Strecken AB und CD liegen.

Weiss vielleicht jemand, wie ich das lösen könnte?
Wenn weitere Infos benötigt werden, reiche ich diese natürlich nach.

Viele Dank schonmal und viele Grüße!

Wallraff · 15 Dez. 2008

Hallo,

Schnittpunkt

(Unter uns, ich hab's gerade ergoogelt und nicht selbst getestet.)

Ahso. Jetzt wissen wir noch nicht, ob der Schnittpunkt innerhalb oder außerhalb liegt.
Verflixt.

Wallraff

KoenigDickBauch · 15 Dez. 2008

muza schrieb:
Weiss vielleicht jemand, wie ich das lösen könnte?

Kochrezept:

Code:

1) Kordinaten nach X,Y,Z abbilden.
2) Jeweils mit den dem 0,0,0 ein ebene Aufspannen ABO und CDO 
   Sollten nicht alle auf einer geraden liegen.
3) Diese beiden Ebenen verschneiden und Geradengleichung ermitteln mit O als Ursprung.
4) Zwei Hilfspunkte erzeugen mit +-gemittelten Radius. Pab und Pcd
5) A,B,Pab auf Polarkoordinaten in der Ebene ABO und
   C,D,Pcd auf Polarkoordinaten in der Ebene CDO abbilden mit O jeweils als Zentrum.
6) Über die Winkel entscheiden ist Pab innerhalb der kürzeren Strecke von AB und 
   Pcd innerhalb der kürzeren Strecke von CD ist.
7) Pab und Pcd auf Lat,Lon,H abbilden.

KDB

muza · 15 Dez. 2008

Super, danke euch beiden - das hilft schonmal sehr weiter!

@ KoenigDickBauch:
Thx fürs Kochrezept

Bei Punkt 6 und 7 hapert es bei mir noch etwas am Verständniss:

6) Wie genau kann ich über einen Winkel prüfen, ob eine Punkt - etwa Pab - auf dem Othodrom ab liegt? Welcher Winkel wäre das in diesem Fall?

7) Lat und Lon verstehe ich, meinst Du mit H die Höhe über dem Meeresspiegel? In dem Fall wäre die wohl 0 bzw. der Radius, da von einer idealen Kugel ausgegangen wird ohne Berge und Täler.

KoenigDickBauch · 16 Dez. 2008

muza schrieb:
Super, danke euch beiden - das hilft schonmal sehr weiter!

@ KoenigDickBauch:
Thx fürs Kochrezept Bei Punkt 6 und 7 hapert es bei mir noch etwas am Verständniss:

6) Wie genau kann ich über einen Winkel prüfen, ob eine Punkt - etwa Pab - auf dem Othodrom ab liegt? Welcher Winkel wäre das in diesem Fall?

7) Lat und Lon verstehe ich, meinst Du mit H die Höhe über dem Meeresspiegel? In dem Fall wäre die wohl 0 bzw. der Radius, da von einer idealen Kugel ausgegangen wird ohne Berge und Täler.

Wenn wir die Punkte in Polarkoordinaten haben schauen wir uns die Winkel an. A hat 49° und B = 47° so haben wir zwei Bögen einmal mit 2° und der andere mit 358° hat nun Pab = 408° so liegt er auf dem kurzen Bogen. (Alle Winkel so berechnen das A<B und A<Pab und B<2*360 und Pab<2*360 ist. Ist AB <180° dann muss A<Pab<B anderfalls B<Pab)

Da die Erde eher ein Ellipsoid ist wirst, du mit den mittleren Radius in der Regel nicht die Oberfläche treffen, daher die Höhe h.

muza · 16 Dez. 2008

Wow, das ging ja schnell - dankeschön! :gott:

Werde es dann morgen probieren die Berechnung zu implementieren.

KoenigDickBauch · 16 Dez. 2008

muza schrieb:
Wow, das ging ja schnell - dankeschön!
Werde es dann morgen probieren die Berechnung zu implementieren.

Worin?

KoenigDickBauch · 16 Dez. 2008

Ich habe mal ein Script für Mopsos gemacht:

Code:

program Orthodrome;
var
  LatA,LonA: real;
  LatB,LonB: real;
  LatC,LonC: real;
  LatD,LonD: real;

  A1,A2,A3: real;
  B1,B2,B3: real;
  C1,C2,C3: real;
  D1,D2,D3: real;

  e1,e2,e3: real;
  f1,f2: real;
  g1,g2,g3: real;

  n: real;
  r: real;

  P1,P2,P3: real;
  Q1,Q2,Q3: real;

  LatP,LonP,HP: real;
  LatQ,LonQ,HQ: real;

begin
   r := 6371000.8;

   LatA:=49;
   LonA:=8;
   LatB:=49;
   LonB:=10;
   LatC:=50;
   LonC:=9;
   LatD:=48;
   LonD:=9;

   Writeln('A: ',WGS84(LatA,LonA));
   Writeln('B: ',WGS84(LatB,LonB));
   Writeln('C: ',WGS84(LatC,LonC));
   Writeln('D: ',WGS84(LatD,LonD));

   WGS84toXYZ(LatA,LonA,0,A1,A2,A3);
   WGS84toXYZ(LatB,LonB,0,B1,B2,B3);
   WGS84toXYZ(LatC,LonC,0,C1,C2,C3);
   WGS84toXYZ(LatD,LonD,0,D1,D2,D3);

   e1:=a2*b3-a3*b1;
   e2:=a3*b1-a1*b3;
   e3:=a1*b2-a2*b1;

   f1:=e1*c1+e2*c2+e3*c3;
   f2:=e1*d1+e2*d2+e3*d3;

   n:=(f2/f1);

   g1:=n*c1+d1;
   g2:=n*c2+d2;
   g3:=n*c3+d3;

   n := sqrt(g1*g1+g2*g2+g3*g3);

   g1:=g1/n;
   g2:=g2/n;
   g3:=g3/n;
   
   P1:=r*g1;
   P2:=r*g2;
   P3:=r*g3;

   Q1:=-P1;
   Q2:=-P2;
   Q3:=-P3;

   XYZtoWGS84(P1,P2,P3,LatP,LonP,HP);
   XYZtoWGS84(Q1,Q2,Q3,LatQ,LonQ,HQ);

   writeln;
   Writeln('P: ',WGS84(LatP,LonP));
   Writeln('Q: ',WGS84(LatQ,LonQ));
   writeln;
   writeln('HP: ',HP:7:0,' m');
   writeln('HQ: ',HQ:7:0,' m');
end.

Das Ergebnis ist:

Code:

A: N 49° 0.0000 E 8° 0.0000
B: N 49° 0.0000 E 10° 0.0000
C: N 50° 0.0000 E 9° 0.0000
D: N 48° 0.0000 E 9° 0.0000

P: N 49° 1.2881 E 9° 0.0000
Q: S 49° 1.2881 W 171° 0.0000

HP:    5008 m
HQ:    5008 m

Im ersten Moment war ich überrascht, das es nicht genau 49° war. Aber klar, dieser Kreis ist ja geneigt und der Schnittpunkt muss was nördlicher liegen.

KDB

t31 · 16 Dez. 2008

@KoenigDickBauch

ich würde mir das gerne in Excel nachbauen - was macht denn WGS84toXYZ?

muza · 16 Dez. 2008

KoenigDickBauch schrieb:
muza schrieb:

Wow, das ging ja schnell - dankeschön!
Werde es dann morgen probieren die Berechnung zu implementieren.

Zum Vergrößern anklicken....

Worin?

Mache ein Projekt in Java, wo ich so einen Schnittpunkt brauche.
Es wäre auch sehr interessant was genau WGS84toXYZ und XYZtoWGS84 genau machen - Kann man das mit der Höhe auch weglassen?

Hab mich heute erstmal vergeblich abgemüht mit der Variante, die Wallraff gepostet hat, da sie mir einfacher erschien.
Um die Punkte nach dieser Methode zu berechnen, braucht man aber erstmal den nordpolnächsten Punkt der beiden Großkreise.
Um diesen nordpolnächsten Punkt zu erhalten braucht man widerrum erstmal den Kurswinkel, siehe hier, den man zumindest theoretisch aus den beiden Punkten berechnen kann, etwa wie hier.
Problem ist, dass bei cos(omega) omega positiv oder negativ sein kann, was man mit dem acos nicht rausbekommt.
Gibt es für den nordpolnächsten Punkt nicht direkt eine Formel, in die man die zwei Punkte der Orthodromen einsetzen kann?

@KDB:
Hab dein Skript mal in Java nachprogrammiert, aber irgendwo ist noch der Wurm drin.
Kannst du in deinem Skript vielleicht noch die Variablen A1-3, B1-3, C1-3, D1-3, P1-3 und Q1-3 ausgeben lassen?
Dann könnte ich vergleichen, an welcher stelle mein Programm einen Fehler hat.

muza · 16 Dez. 2008

Code:

public class Orthodrome 
{
	public static double r = 6371000.8;
	
	//Hilfsmethode
	public static double[] WGS84toXYZ(   double Lon, double Lat  )
	{
		double[] erg = new double[3];
		Lon += 180.0;
		Lat += 90.0;
		erg[0] = r * Math.sin( Math.toRadians( Lat ))  * Math.cos( Math.toRadians( Lon ));
		erg[1] = r * Math.sin( Math.toRadians( Lat ))  * Math.sin( Math.toRadians( Lon ));
		erg[2] = r * Math.cos( Math.toRadians( Lat ));
		return erg;
	}
	
	//Hilfsmethode
	public static double[] XYZtoWGS84 ( double x, double y, double z )
	{
		double[] erg = new double[2];
		erg[0] = Math.toDegrees(  Math.atan( y / x  )  ) - 180.0 ; 	//Lon
		erg[1] = Math.toDegrees( Math.acos( z / r  ) ) - 90.0  ; //Lat
		return erg; 
	}
	
	//Hauptprogramm
	public static void main(String[] args  )
	{
	   double LatA=49;
	   double LonA=8;
	   double LatB=49;
	   double LonB=10;
	   double LatC=50;
	   double LonC=9;
	   double LatD=48;
	   double LonD=9;

	   double[] A = WGS84toXYZ(LonA, LatA);
	   double[] B = WGS84toXYZ(LonB, LatB);
	   double[] C = WGS84toXYZ(LonC, LatC);
	   double[] D = WGS84toXYZ(LonD, LatD);
	   	   
	   double e1 = A[1] * B[2] - A[2] * B[0];
	   double e2 = A[2] * B[0] - A[0] * B[2];
	   double e3 = A[0] * B[1] - A[1] * B[0];

	   double f1 = e1 * C[0] + e2 * C[1] + e3 * C[2];
	   double f2 = e1 * D[0] + e2 * D[1] + e3 * D[2];

	   double n=(f2/f1);

	   double g1 = n * C[0] + D[0];
	   double g2 = n * C[1] + D[1];
	   double g3 = n * C[2] + D[2];

	   n = Math.sqrt(g1*g1+g2*g2+g3*g3);

	   g1=g1/n;
	   g2=g2/n;
	   g3=g3/n;
	   
	   double P1=r*g1;
	   double P2=r*g2;
	   double P3=r*g3;

	   double Q1=-P1;
	   double Q2=-P2;
	   double Q3=-P3;
	   
	   double[] P =XYZtoWGS84(P1,P2,P3);
	   double[] Q = XYZtoWGS84(Q1,Q2,Q3);
	   
	   System.out.println( "P - Lon:" + P[0] + " Lat:" + P[1]  );
	   System.out.println( "Q - Lon:" + Q[0] + " Lat:" + Q[1]  );
	}
}

Ausgabe

Code:

P - Lon:-171.0 Lat:49.02161758691602
Q - Lon:-171.0 Lat:-49.02161758691602

Also sinngemäß

Code:

P: N 49.02161758691602  W 171.0 
Q: S 49.02161758691602  W 171.0

Anscheinend wird der Längengrad falsch berechnet, da P und Q nicht gegenüberliegen auf der Kugel.
Hab die Formeln dazu auf dieser Seite unten gefunden.
Nach der Formel muss aber wohl zweimal die gleiche Länge rauskommen, da x / y immer gleich -x / -y ist.
Sieht jemand den Fehler?

Hier noch meine Werte für P und Q im Raum

Code:

P1: -4126500.8040500185
P2: -653573.5203482106
P3: -4809832.009672982
Q1: 4126500.8040500185
Q2: 653573.5203482106
Q3: 4809832.009672982

Wallraff · 17 Dez. 2008

Hallo,

Mist, da muss ich doch zugeben, dass ich die gefundene Formel nicht ganz überrissen hatte.
Sah ganz einfach aus, so wie der Geradenschnitt in der Ebene.

KDB kann ich natürlich nicht einfach für das saubere Programm loben ...

Das Haar in der Suppe= der Radius.
Warum ausgerechnet 6371000.8
Alles so elegant, aber XYZ auf dem Ellipsoid und der Schnittpunkt mit dem Standardradius beissen sich so unschön...?

BTW: Hast Du Urlaub, oder hast Du den Kollegen gesagt, Du bist auf Dienstreise und hast Dich eingeschlossen ?
Muss ich noch programmieren :/

Grüße Wallraff

muza · 17 Dez. 2008

Der Radius müsste doch eigentlich keine Rolle spielen, da er sich durch die Umrechnung zurück in Längen- und Breitengrad wieder wegkürzt zum Schluss.
Hmm.. bildet die WGS84toXYZ Funktion auf einen Ellipsoiden oder eine Kugel ab?
Weil mir würde eine Kugel erstmal reichen um die Berechnung einfacher zu halten.

KoenigDickBauch · 17 Dez. 2008

r ist Wurst
atan2(y,x)
auf einer Kugel ist alles viel leichter.

PS: Korrigiere r darf nicht 0 sein und sollte einen Punkt erzeugen, der in der Nähe der Oberfläche ist, damit die Numerik mit gleichen Größen arbeitet (Rundungsfehler).

muza · 17 Dez. 2008

KoenigDickBauch schrieb:
r ist Wurst
atan2(y,x)
auf einer Kugel ist alles viel leichter.

PS: Korrigiere r darf nicht 0 sein und sollte einen Punkt erzeugen, der in der Nähe der Oberfläche ist, damit die Numerik mit gleichen Größen arbeitet (Rundungsfehler).

Hast du vielleicht eine Idee, warum bei mir nicht der richtige Punkt rauskommt?

KoenigDickBauch · 17 Dez. 2008

muza schrieb:
Hast du vielleicht eine Idee, warum bei mir nicht der richtige Punkt rauskommt?

JA Hint: steht im letzten Beitrag von mir.

muza · 17 Dez. 2008

KoenigDickBauch schrieb:
muza schrieb:

Hast du vielleicht eine Idee, warum bei mir nicht der richtige Punkt rauskommt?

Zum Vergrößern anklicken....

JA Hint: steht im letzten Beitrag von mir.

Okay, habs

In der Beschreibung von atan2(x,y) steht, dass es im Prinzip auch atan( x / y ) berechnet, aber diesmal vorzeichenbehaftet, oder?
Klappt jetzt jedenfalls - DANKE!

Code:

	//Hilfsmethode
	public static double[] XYZtoWGS84 ( double x, double y, double z )
	{
		double[] erg = new double[2];
		erg[0] = Math.toDegrees(  Math.atan2( y , x  )  ); 	//Lon
		erg[1] = Math.toDegrees( Math.acos( z / r  ) ) - 90.0  ; //Lat
		return erg; 
	}

Nachtrag:
An welcher Stelle könnte man jetzt am besten Prüfen, ob einer der gefundenen Schnittpunkte tatsächlich innerhalb der zwei Orthodromen liegt?
Hab mir überlegt, dass wenn Q der Schnittpunkt ist, der Abstand AB gleich AQ + BQ sein müsste.
Wie würde das über die Polarkoordinaten funktionieren?

KoenigDickBauch · 17 Dez. 2008

Und dann möchte ich auf zwei Tippfehler aus meinem Script hinweisen, die wie der Fehlerteufel es will, sich bei diesem Beispiel fast aufheben. Hier dein Code, suche nach einem Symmetriefehler (Vektorprodukt für Normvektor).

Code:

      double e1 = A[1] * B[2] - A[2] * B[0];
      double e2 = A[2] * B[0] - A[0] * B[2];
      double e3 = A[0] * B[1] - A[1] * B[0];

und hier muss es so heißen:

Code:

      double n=- (f2/f1);

KoenigDickBauch · 17 Dez. 2008

muza schrieb:
In der Beschreibung von atan2(x,y) steht, dass es im Prinzip auch atan( x / y ) berechnet, aber diesmal vorzeichenbehaftet, oder?

es heißt atan2(y,x) da es ja mit atan( y / x ) übereinstimmen soll. Ersetze / durch ,. Und noch besser es wird auf x=0 überprüft. Mach mal ein Beispiel, so daß die Pole als Ergebnis heraus kommen.

muza schrieb:
Klappt jetzt jedenfalls - DANKE!

Nicht ganz! Schaue Beitrag zuvor.

muza · 17 Dez. 2008

Stimmt, so siehts symmetrischer aus

Code:

	   double e1 = A[1] * B[2] - A[2] * B[1];
	   double e2 = A[2] * B[0] - A[0] * B[2];
	   double e3 = A[0] * B[1] - A[1] * B[0];

Gibt nun zusammen mit dem neuen (-)

Code:

P - Lon:-171.0 Lat:49.00432103313932
Q - Lon:8.999999999999993 Lat:-49.00432103313933

Sieht doch ganz gut aus.
Werde dann mal noch etwas mit anderen Koordinaten testen :^^:

Schnittpunkt zweier Orthodrome?

Geocacher

Geocacher

Geoguru

Geocacher

Geoguru

Geocacher

Geoguru

Geoguru

Geowizard

Geocacher

Geocacher

Geocacher

Geocacher

Geoguru

Geocacher

Geoguru

Geocacher

Geoguru

Geoguru

Geocacher