• Willkommen im Geoclub - dem größten deutschsprachigen Geocaching-Forum. Registriere dich kostenlos, um alle Inhalte zu sehen und neue Beiträge zu erstellen.

Schnittpunkt zweier Orthodrome?

M

muza

Geocacher
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

Wallraff

Geocacher
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

KoenigDickBauch

Geoguru
muza schrieb:
Weiss vielleicht jemand, wie ich das lösen könnte?
Zum Vergrößern anklicken....
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
 
OP
M

muza

Geocacher
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

KoenigDickBauch

Geoguru
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.
Zum Vergrößern anklicken....
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.
 
OP
M

muza

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

KoenigDickBauch

Geoguru
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

t31

Geowizard
@KoenigDickBauch

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

muza

Geocacher
KoenigDickBauch schrieb:
muza schrieb:
Wow, das ging ja schnell - dankeschön! :gott:
Werde es dann morgen probieren die Berechnung zu implementieren.
Zum Vergrößern anklicken....
Worin?
Zum Vergrößern anklicken....
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.
 
OP
M

muza

Geocacher
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

Wallraff

Geocacher
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
 
OP
M

muza

Geocacher
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

KoenigDickBauch

Geoguru
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).
 
OP
M

muza

Geocacher
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).
Zum Vergrößern anklicken....

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

muza

Geocacher
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. ;)
Zum Vergrößern anklicken....

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

KoenigDickBauch

Geoguru
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

KoenigDickBauch

Geoguru
muza schrieb:
In der Beschreibung von atan2(x,y) steht, dass es im Prinzip auch atan( x / y ) berechnet, aber diesmal vorzeichenbehaftet, oder?
Zum Vergrößern anklicken....
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!
Zum Vergrößern anklicken....
Nicht ganz! Schaue Beitrag zuvor.
 
OP
M

muza

Geocacher
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 :^^:
 
Du musst dich einloggen oder registrieren, um hier zu antworten.
Oben