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

Mathematische Grundlagen: Wegpunktprojektion

J

jahwe2000

Geocacher
Hallo,

ich engagiere mich (nicht uneigennützig) zur Zeit im GeoBeagle-Projekt (kurz: geocaching fürs Android-Handy-Betriebssystem).

Dort fehlt bisher eine Funktion, mittels der man vom aktuellen Punkt x meter in y Grad gehen kann. (Allgemein Wegpunktprojektion genannt).

Da GeoBeagle OpenSource ist, ist jeder eingeladen einzuspringen. Die Informatische Seite ist kein Problem, aber mir fehlt da das mathematische Fundament. Im Prinzip geht es ja darum, mit ein bisschen Sphärischer Geometrie die Koordinaten zu errechnen. Es gibt aber unzählige verschiedene Projektionsverfahren zur Projektion des Ellipsoiden auf eine 2D-Darstellung etc. Gibt es da nun ein Standardverfahren? Und wie sähe die Projektion (Heißt das wirklich noch so? Ich würde es intuitiv transposition nennen) mathematisch aus?

Hat jemand vielleicht einen Link mit weiterführenden Informationen? Oder noch besser: Gibt es Java-Bibliotheken, die das bereits können?

Ich freue mich über jede Art von (positiver ;-) ) Antwort.

Viele Grüße,
Philip
 
K

kiozen

Geomaster
Code: 
XY GPS_Math_Wpt_Projection(XY& pt1, double distance, double bearing)
{
    XY pt2;

    double d    = distance / 6378130.0;
    double lon1 = pt1.u;
    double lat1 = pt1.v;

    double lat2 = asin(sin(lat1) * cos(d) + cos(lat1) * sin(d) * cos(-bearing));
    double lon2 = cos(lat1) == 0 ? lon1 : fmod(lon1 - asin(sin(-bearing) * sin(d) / cos(lat1)) + PI, TWOPI) - PI;

    pt2.u = lon2;
    pt2.v = lat2;
    return pt2;
}


HTH

Oliver
 
pfeffer

pfeffer

Geowizard
diese Methode nimmt allerdings an, die Welt sei eine Kugel. Wenn man es auf einer Ellipse machen will, wird es nen ganzen Stück komplizierter:
Code: 
public class GeodeticCalculator
{
   static private final double TwoPi = 2.0 * Math.PI;
   
   /** Degrees/Radians conversion constant. */
   static private final double PiOver180 = Math.PI / 180.0;

   /**
    * Calculate the destination and final bearing after traveling a specified
    * distance, and a specified starting bearing, for an initial location. This
    * is the solution to the direct geodetic problem.
    * 
    * @param ellipsoid reference ellipsoid to use
    * @param start starting location
    * @param startBearing starting bearing (degrees)
    * @param distance distance to travel (meters)
    * @param endBearing bearing at destination (degrees) element at index 0 will
    *            be populated with the result
    * @return
    */
   static public TrackPoint calculateEndingGlobalCoordinates(Ellipsoid ellipsoid, TrackPoint start, double startBearing, double distance,
         double[] endBearing)
   {
      double a = ellipsoid.getSemiMajorAxis();
      double b = ellipsoid.getSemiMinorAxis();
      double aSquared = a * a;
      double bSquared = b * b;
      double f = ellipsoid.getFlattening();
      double phi1 = start.latDec  * PiOver180;
      double alpha1 = startBearing * PiOver180;
      double cosAlpha1 = Math.cos(alpha1);
      double sinAlpha1 = Math.sin(alpha1);
      double s = distance;
      double tanU1 = (1.0 - f) * Math.tan(phi1);
      double cosU1 = 1.0 / Math.sqrt(1.0 + tanU1 * tanU1);
      double sinU1 = tanU1 * cosU1;

      // eq. 1
      double sigma1 = Math.atan2(tanU1, cosAlpha1);

      // eq. 2
      double sinAlpha = cosU1 * sinAlpha1;

      double sin2Alpha = sinAlpha * sinAlpha;
      double cos2Alpha = 1 - sin2Alpha;
      double uSquared = cos2Alpha * (aSquared - bSquared) / bSquared;

      // eq. 3
      double A = 1 + (uSquared / 16384) * (4096 + uSquared * (-768 + uSquared * (320 - 175 * uSquared)));

      // eq. 4
      double B = (uSquared / 1024) * (256 + uSquared * (-128 + uSquared * (74 - 47 * uSquared)));

      // iterate until there is a negligible change in sigma
      double deltaSigma;
      double sOverbA = s / (b * A);
      double sigma = sOverbA;
      double sinSigma;
      double prevSigma = sOverbA;
      double sigmaM2;
      double cosSigmaM2;
      double cos2SigmaM2;

      for (;;)
      {
         // eq. 5
         sigmaM2 = 2.0 * sigma1 + sigma;
         cosSigmaM2 = Math.cos(sigmaM2);
         cos2SigmaM2 = cosSigmaM2 * cosSigmaM2;
         sinSigma = Math.sin(sigma);
         double cosSignma = Math.cos(sigma);

         // eq. 6
         deltaSigma = B
               * sinSigma
               * (cosSigmaM2 + (B / 4.0)
                     * (cosSignma * (-1 + 2 * cos2SigmaM2) - (B / 6.0) * cosSigmaM2 * (-3 + 4 * sinSigma * sinSigma) * (-3 + 4 * cos2SigmaM2)));

         // eq. 7
         sigma = sOverbA + deltaSigma;

         // break after converging to tolerance
         if (Math.abs(sigma - prevSigma) < 0.0000000000001) break;

         prevSigma = sigma;
      }

      sigmaM2 = 2.0 * sigma1 + sigma;
      cosSigmaM2 = Math.cos(sigmaM2);
      cos2SigmaM2 = cosSigmaM2 * cosSigmaM2;

      double cosSigma = Math.cos(sigma);
      sinSigma = Math.sin(sigma);

      // eq. 8
      double phi2 = Math.atan2(sinU1 * cosSigma + cosU1 * sinSigma * cosAlpha1, (1.0 - f)
            * Math.sqrt(sin2Alpha + Math.pow(sinU1 * sinSigma - cosU1 * cosSigma * cosAlpha1, 2.0)));

      // eq. 9
      double tanLambda = sinSigma * sinAlpha1 / (cosU1 * cosSigma - sinU1 * sinSigma * cosAlpha1);
      double lambda = Math.atan(tanLambda);

      // eq. 10
      double C = (f / 16) * cos2Alpha * (4 + f * (4 - 3 * cos2Alpha));

      // eq. 11
      double L = lambda - (1 - C) * f * sinAlpha * (sigma + C * sinSigma * (cosSigmaM2 + C * cosSigma * (-1 + 2 * cos2SigmaM2)));

      // eq. 12
      double alpha2 = Math.atan2(sinAlpha, -sinU1 * sinSigma + cosU1 * cosSigma * cosAlpha1);

      // build result
      double latitude = phi2 / PiOver180;
      double longitude = start.lonDec + (L / PiOver180);

      if ((endBearing != null) && (endBearing.length > 0))
      {
         endBearing[0] = alpha2 / PiOver180;
      }

      return new TrackPoint(latitude, longitude);
   }
Das hat MiK so in Cachewolf programmiert (GPL) nach diesen Formeln http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf

Gruß,
Pfeffer.
 
E

eremiljo

Geocacher
Wie groß ist denn die Abweichung der beiden Verfahren in den Geocaching-typischen Dimensionen? (Also Projektionsweite bis 1km?)
 
KoenigDickBauch

KoenigDickBauch

Geoguru
eremiljo schrieb:
Wie groß ist denn die Abweichung der beiden Verfahren in den Geocaching-typischen Dimensionen? (Also Projektionsweite bis 1km?)
Zum Vergrößern anklicken....

Code: 
49°
8°

Winkel 45
Entfernung 1000

WGS84:
N 49° 0.3815
E 8° 0.5799

KUGEL (Radius=6366707.02)
N 49° 0.3818
E 8° 0.5820

WGS84
77.75° 2.62m

Kugel
360° 2.61m

Bei reiner Projektion ist der Unterschied zu vernachlässigen. Wenn man einen passenden Radius für unsere Breiten nimmt.

Aber wenn damit dann gepeilt werden soll könnte sich der Fehler sehr stark bemerkbar machen.
 
satanklaus

satanklaus

Geomaster
Meine Frage passt ganz gut zu den letzten beiden Posts:

ich habe bei manchen Multis nicht immer alle erforderlichen Zahlen, um die nächste Koordinate zu berechnen. Manchmal hat man aber einen plausiblen Bereich, der sich aus anderen Infos ergibt (z.B. 2 <= x <= 5 ). Je nachdem, wo die unbekannt Zahl in der Formel auftaucht, hat sie einen mehr oder weniger großen Einfluss auf das Ergebnis bzw den Fehler. Wenn es z.B. die letzte Stelle betrifft kann man ja vielleicht per Suche in einem bestimmten Radius weiterkommen. Daher habe ich mich gefragt, wieviel eine Stelle Unterschied (also z.B. dd mm.mm2 statt dd mm.mm1) wohl in Metern ausmacht. Empirisch bin ich dem unterwegs mal nachgegangen, indem ich z.B. genau Richtung Norden gelaufen bin und die Koordinaten beobachtet habe. Das waren vielleicht 5-10m.

Wegpunkte mit allen Möglichkeiten anzulegen und damit das Areal einzugrenzen wollte ich nicht, daher die Frage:
Gibt es für sowas eine Faustformel?

Danke!
 
Du musst dich einloggen oder registrieren, um hier zu antworten.
Oben