On Fri, Sep 19, 2008 at 03:42:11AM +0200, Florian Heer wrote:
Bei einigen (ca. 20%) der Tracks, stimmen die Werte aber nicht (Fehler weit über 1%),Fehler von 1% würde ich bei Deiner Formel durchaus erwarten. Nimm Vincenty's Algorithmus, der ist besser. Gibts auf [1] in JavaScript, meine Übersetzungen in C und Python habe ich angehängt. Lizenz ist LGPL.
[1] http://www.movable-type.co.uk/scripts/latlong-vincenty.html CU Sascha -- http://sascha.silbe.org/ http://www.infra-silbe.de/
# arch-tag: 413a3fa1-5a89-41f2-a905-0ad78cdd0267 from math import * ############################################################################# ## ## Copyright (C) 2002-2008 Chris Veness -- http://www.movable-type.co.uk/scripts/latlong-vincenty.html ## Copyright (C) 2007-2008 Sascha Silbe -- http://sascha.silbe.org/ ## ## This file may be distributed and/or modified under the terms of the ## GNU Lesser General Public License version 2 as published by the Free ## Software Foundation and appearing in the file COPYING.LESSER included in ## the packaging of this file. ## ############################################################################# # source of algorithm: # Vincenty direct formula - T Vincenty, "Direct and Inverse Solutions of Geodesics on the # Ellipsoid with application of nested equations", Survey Review, vol XXII no 176, 1975 # http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf # # based on JavaScript implementation by Chris Veness, available at # http://www.movable-type.co.uk/scripts/latlong-vincenty-direct.html # Python implementation by Sascha Silbe # # License for original JavaScript implementation: # JavaScript functions are shown below. You are welcome to re-use these scripts [without any warranty express or implied] # provided you retain my copyright notice and when possible a link to my website (under a LGPL license). If you have any # queries or find any problems, please contact me. # (C) 2002-2008 Chris Veness # WGS-84 ellipsiod a = 6378137. b = 6356752.3142 f = 1/298.257223563 def distance(p1, p2) : p1Lat = p1[0]*pi/180 p1Lon = p1[1]*pi/180 p2Lat = p2[0]*pi/180 p2Lon = p2[1]*pi/180 L = p2Lon - p1Lon U1 = atan((1-f) * tan(p1Lat)) U2 = atan((1-f) * tan(p2Lat)) sinU1 = sin(U1) cosU1 = cos(U1) sinU2 = sin(U2) cosU2 = cos(U2) Lambda = L LambdaP = 2*pi iterLimit = 20 while ((fabs(Lambda-LambdaP) > 1e-12) and iterLimit) : iterLimit -= 1 sinLambda = sin(Lambda) cosLambda = cos(Lambda) sinSigma = sqrt((cosU2*sinLambda) * (cosU2*sinLambda) + (cosU1*sinU2-sinU1*cosU2*cosLambda) * (cosU1*sinU2-sinU1*cosU2*cosLambda)) if (sinSigma==0) : return 0, 0 # co-incident points cosSigma = sinU1*sinU2 + cosU1*cosU2*cosLambda sigma = atan2(sinSigma, cosSigma) sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma cosSqAlpha = 1 - sinAlpha*sinAlpha if (cosSqAlpha == 0) : cos2SigmaM = 0 # equatorial line: cosSqAlpha=0 (Rule 6) else : cos2SigmaM = cosSigma - 2*sinU1*sinU2/cosSqAlpha C = f/16*cosSqAlpha*(4+f*(4-3*cosSqAlpha)) LambdaP = Lambda Lambda = L + (1-C) * f * sinAlpha * (sigma + C*sinSigma*(cos2SigmaM+C*cosSigma*(-1+2*cos2SigmaM*cos2SigmaM))) if (iterLimit==0) : # formula failed to converge return None, None uSq = cosSqAlpha * (a*a - b*b) / (b*b) A = 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq))) B = uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq))) deltaSigma = B*sinSigma*(cos2SigmaM+B/4*(cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)-B/6*cos2SigmaM*(-3+4*sinSigma*sinSigma)*(-3+4*cos2SigmaM*cos2SigmaM))) s = b*A*(sigma-deltaSigma) # initial bearing bearing = atan2(cosU2*sinLambda, cosU1*sinU2 - sinU1*cosU2*cosLambda)*180/pi # final bearing # bearing = atan2(cosU1*sinLambda, -sinU1*cosU2 + cosU1*sinU2*cosLambda)*180/pi while (bearing < 0) : bearing += 360 return s, bearing def destination(origin, bearing, distance) : lat1, lon1 = origin s = float(distance) alpha1 = bearing*pi/180 sinAlpha1 = sin(alpha1) cosAlpha1 = cos(alpha1) tanU1 = (1-f) * tan(lat1*pi/180) cosU1 = 1. / sqrt((1 + tanU1*tanU1)) sinU1 = tanU1*cosU1 sigma1 = atan2(tanU1, cosAlpha1) sinAlpha = cosU1 * sinAlpha1 cosSqAlpha = 1 - sinAlpha*sinAlpha uSq = cosSqAlpha * (a*a - b*b) / (b*b) A = 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq))) B = uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq))) sigma = s / (b*A) sigmaP = 2*pi while (abs(sigma-sigmaP) > 1e-12) : cos2SigmaM = cos(2*sigma1 + sigma) sinSigma = sin(sigma) cosSigma = cos(sigma) deltaSigma = B*sinSigma*(cos2SigmaM+B/4*(cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)-B/6*cos2SigmaM*(-3+4*sinSigma*sinSigma)*(-3+4*cos2SigmaM*cos2SigmaM))) sigmaP = sigma sigma = s / (b*A) + deltaSigma tmp = sinU1*sinSigma - cosU1*cosSigma*cosAlpha1 lat2 = atan2(sinU1*cosSigma + cosU1*sinSigma*cosAlpha1, (1-f)*sqrt(sinAlpha*sinAlpha + tmp*tmp)) Lambda = atan2(sinSigma*sinAlpha1, cosU1*cosSigma - sinU1*sinSigma*cosAlpha1) C = f/16*cosSqAlpha*(4+f*(4-3*cosSqAlpha)) L = Lambda - (1-C) * f * sinAlpha * (sigma + C*sinSigma*(cos2SigmaM+C*cosSigma*(-1+2*cos2SigmaM*cos2SigmaM))) # revAz = atan2(sinAlpha, -tmp) # final bearing lat2D = lat2*180/pi lon2D = lon1+L*180/pi while (lat2D > 90) : lat2D -= 180 while (lat2D < -90) : lat2D += 180 while (lon2D > 180) : lon2D -= 360 while (lon2D < -180) : lon2D += 360 return (lat2D, lon2D)
// arch-tag: a7d2a3ab-1c56-4ef5-9edd-0c4f1a5de4dd // ############################################################################# // ## // ## Copyright (C) 2002-2008 Chris Veness -- http://www.movable-type.co.uk/scripts/latlong-vincenty.html // ## Copyright (C) 2007-2008 Sascha Silbe -- http://sascha.silbe.org/ // ## // ## This file may be distributed and/or modified under the terms of the // ## GNU Lesser General Public License version 2 as published by the Free // ## Software Foundation and appearing in the file COPYING.LESSER included in // ## the packaging of this file. // ## // ############################################################################# #include <math.h> #include <stdint.h> #include "vincenty.h" #ifndef M_PI // M_PI was removed in C99 :-/ #define M_PI 3.14159265358979323846264338327 #endif static double const a = 6378137., b = 6356752.3142, f = 1/298.257223563; // WGS-84 ellipsiod double distance(double _p1Lon, double _p1Lat, double _p2Lon, double _p2Lat, double *bearing) { const double p1Lat = _p1Lat*M_PI/180; const double p1Lon = _p1Lon*M_PI/180; const double p2Lat = _p2Lat*M_PI/180; const double p2Lon = _p2Lon*M_PI/180; const double L = p2Lon - p1Lon; const double U1 = atan((1-f) * tan(p1Lat)); const double U2 = atan((1-f) * tan(p2Lat)); const double sinU1 = sin(U1), cosU1 = cos(U1); const double sinU2 = sin(U2), cosU2 = cos(U2); double lambda = L, lambdaP = 2*M_PI; uint8_t iterLimit = 20; double sinLambda = 0., cosLambda = 0., sinSigma = 0., cosSigma = 0., sigma = 0., sinAlpha, cosSqAlpha = 0.; double cos2SigmaM = 0., C; double uSq, A, B, deltaSigma, s; while ((fabs(lambda-lambdaP) > 1e-12) && --iterLimit) { sinLambda = sin(lambda); cosLambda = cos(lambda); sinSigma = sqrt((cosU2*sinLambda) * (cosU2*sinLambda) + (cosU1*sinU2-sinU1*cosU2*cosLambda) * (cosU1*sinU2-sinU1*cosU2*cosLambda)); if (sinSigma==0) return 0; // co-incident points cosSigma = sinU1*sinU2 + cosU1*cosU2*cosLambda; sigma = atan2(sinSigma, cosSigma); sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma; cosSqAlpha = 1 - sinAlpha*sinAlpha; cos2SigmaM = cosSigma - 2*sinU1*sinU2/cosSqAlpha; if (isnan(cos2SigmaM)) cos2SigmaM = 0; // equatorial line: cosSqAlpha=0 (§6) C = f/16*cosSqAlpha*(4+f*(4-3*cosSqAlpha)); lambdaP = lambda; lambda = L + (1-C) * f * sinAlpha * (sigma + C*sinSigma*(cos2SigmaM+C*cosSigma*(-1+2*cos2SigmaM*cos2SigmaM))); } if (iterLimit==0) return NAN; // formula failed to converge uSq = cosSqAlpha * (a*a - b*b) / (b*b); A = 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq))); B = uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq))); deltaSigma = B*sinSigma*(cos2SigmaM+B/4*(cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)- B/6*cos2SigmaM*(-3+4*sinSigma*sinSigma)*(-3+4*cos2SigmaM*cos2SigmaM))); s = b*A*(sigma-deltaSigma); if (bearing) // initial bearing *bearing = atan2(cosU2*sinLambda, cosU1*sinU2 - sinU1*cosU2*cosLambda)*180/M_PI; // final bearing // *bearing = atan2(cosU1*sinLambda, -sinU1*cosU2 + cosU1*sinU2*cosLambda)*180/M_PI; return s; } void destination(double originLon, double originLat, double bearing, double distance, double *destLon, double *destLat) { double lat1, lon1, s, alpha1, sinAlpha1, cosAlpha1; double tanU1, cosU1, sinU1, sigma1, sinAlpha, cosSqAlpha = 0., uSq, A, B, sigma, sigmaP; double cos2SigmaM = 0., sinSigma = 0., cosSigma = 0., deltaSigma, tmp, lat2, Lambda, C, L, lat2D, lon2D; lat1 = originLat; lon1 = originLon; s = distance; alpha1 = bearing*M_PI/180; sinAlpha1 = sin(alpha1); cosAlpha1 = cos(alpha1); tanU1 = (1-f) * tan(lat1*M_PI/180); cosU1 = 1. / sqrt((1 + tanU1*tanU1)); sinU1 = tanU1*cosU1; sigma1 = atan2(tanU1, cosAlpha1); sinAlpha = cosU1 * sinAlpha1; cosSqAlpha = 1 - sinAlpha*sinAlpha; uSq = cosSqAlpha * (a*a - b*b) / (b*b); A = 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq))); B = uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq))); sigma = s / (b*A); sigmaP = 2*M_PI; while (fabs(sigma-sigmaP) > 1e-12) { cos2SigmaM = cos(2*sigma1 + sigma); sinSigma = sin(sigma); cosSigma = cos(sigma); deltaSigma = B*sinSigma*(cos2SigmaM+B/4*(cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)-B/6*cos2SigmaM*(-3+4*sinSigma*sinSigma)*(-3+4*cos2SigmaM*cos2SigmaM))); sigmaP = sigma; sigma = s / (b*A) + deltaSigma; } tmp = sinU1*sinSigma - cosU1*cosSigma*cosAlpha1; lat2 = atan2(sinU1*cosSigma + cosU1*sinSigma*cosAlpha1, (1-f)*sqrt(sinAlpha*sinAlpha + tmp*tmp)); Lambda = atan2(sinSigma*sinAlpha1, cosU1*cosSigma - sinU1*sinSigma*cosAlpha1); C = f/16*cosSqAlpha*(4+f*(4-3*cosSqAlpha)); L = Lambda - (1-C) * f * sinAlpha * (sigma + C*sinSigma*(cos2SigmaM+C*cosSigma*(-1+2*cos2SigmaM*cos2SigmaM))); // revAz = atan2(sinAlpha, -tmp) # final bearing lat2D = lat2*180/M_PI; lon2D = lon1+L*180/M_PI; while (lat2D > 90) lat2D = lat2D - 180; while (lat2D < -90) lat2D = lat2D + 180; while (lon2D > 180) lon2D = lon2D - 360; while (lon2D < -180) lon2D = lon2D + 360; *destLon = lon2D; *destLat = lat2D; }
// arch-tag: 0ad91b8e-ebf5-45ae-850f-541cb4fa3cc0 // ############################################################################# // ## // ## Copyright (C) 2008 Sascha Silbe -- http://sascha.silbe.org/ // ## // ## This file may be distributed and/or modified under the terms of the // ## GNU Lesser General Public License version 2 as published by the Free // ## Software Foundation and appearing in the file COPYING.LESSER included in // ## the packaging of this file. // ## // ############################################################################# #ifndef __VINCENTY_H__ #define __VINCENTY_H__ double distance(double p1Lon, double p1Lat, double p2Lon, double p2Lat, double *bearing); void destination(double originLon, double originLat, double bearing, double distance, double *destLon, double *destLat); #endif
signature.asc
Description: Digital signature
_______________________________________________ Talk-de mailing list [email protected] http://lists.openstreetmap.org/listinfo/talk-de
