Записки программиста, обо всем и ни о чем. Но, наверное, больше профессионального.

2006-12-15

Расстояние между двумя точками, эллипсоид. Не новость конечно, но народ интересуется часто

Очередной раз столкнулся с задачей поточнее посчитать расстояние между двумя точками на местности. Некоторые скажут "боян"... да и хрен с ними. Короче, кто до сих пор не знает как посчитать расстояние на эллипсоиде - прошу:

Координаты в десятичных градусах, типа долгота lon: 37.53 широта lat: 55.40. Результат в метрах а заодно и азимут направления движения в градусах от 0 до 360.
Тем кто с GPS работает - нужно часто. Правда еще нужно уметь перевести координаты из NMEA записи в десятичные градусы, это для пятого класса задача... ладно, сначала переведем NMEA запись координат в градусы:


function nmea2DecDegree($nmeaCoord) { // decimal(9,7), convert nmea ddmm.mmmm (dddmm.mmmm) to decimal degree
$res = 0.0;
$pp = strpos($nmeaCoord, '.');
if ($pp === false || $pp <>

$d = substr($nmeaCoord, 0, $pp-2);
$m = substr($nmeaCoord, $pp-2);
$res = (double)$d + ((double)$m / 60.0);
return $res;
}


ну а теперь посчитаем расстояние в метрах и азимут в градусах:


function calcDistance ($StartLat, $StartLon, $EndLat, $EndLon) { // array, 'dist' => 'x' meters, 'bearing' => 'y' degree; from decimal degree coords
//~ 10 сантиметров точность, когда
//~ abs(lon1-lon2) <= 0.0000014
//~ abs(lat1-lat2) <= 0.0000008
$res = array('dist'=>0.0, 'bearing'=>0.0);
$StartLat = (double)str_replace(',', '.', $StartLat); $StartLon = (double)str_replace(',', '.', $StartLon);
$EndLat = (double)str_replace(',', '.', $EndLat); $EndLon = (double)str_replace(',', '.', $EndLon);

if ($StartLat == 0.0 || $StartLon == 0.0 || $EndLat == 0.0 || $EndLon == 0.0) return $res;
if ( abs($StartLon - $EndLon) <= 0.0000014 && abs($StartLat - $EndLat) <= 0.0000008 ) return $res;

$D2R = 0.01745329251994330; // Pi/180
$R2D = 57.29577951308230000; // 180/Pi
/*
// http://www.pcigeomatics.com/cgi-bin/pcihlp/PROJ%7CEARTH+MODELS%7CELLIPSOIDS
// ELLIPS Descriptor Semi-Major Axis (A) (metres) Semi-Minor Axis (B) (metres)
// 12 WGS 1984 6378137.000000 6356752.314245
// 15 Krassovsky 1940 6378245.000000 6356863.018800
*/
$a = 6378137.0; // Semi-major axis of ellipsoid in meters
// $b = 6356752.314245; // Semi-minor axis of ellipsoid
$e2 = 0.00673949674233346; // 2nd eccentricity squared Geocent_ep2 = (Geocent_a2 - Geocent_b2) / Geocent_b2

$fdLambda = ($StartLon - $EndLon) * $D2R;
$fdPhi = ($StartLat - $EndLat) * $D2R;
$fPhimean = ($StartLat + $EndLat) / 2.0 * $D2R;
$fTemp = 1 - $e2 * pow( sin($fPhimean), 2);
$fRho = $a * (1 - $e2) / pow($fTemp, 1.5);
$fNu = $a / sqrt(1 - $e2 * sin($fPhimean) * sin($fPhimean) );
$fz = 2 * asin(sqrt (
pow(sin( $fdPhi / 2.0), 2 ) + cos( $EndLat * $D2R )
* cos( $StartLat * $D2R ) * pow(sin($fdLambda / 2.0), 2)
) );

$fAlpha = asin(cos( $EndLat * $D2R ) * sin($fdLambda) / sin($fz) );
$fR = $fRho * $fNu / ($fRho * pow(sin($fAlpha), 2) + $fNu * pow(cos($fAlpha), 2));
$res['dist'] = $fz * $fR;

if ($res['dist'] <= 0.999) return $res;

$Bearing = abs($fAlpha * $R2D);
if (($StartLat <= $EndLat) and ($StartLon > $EndLon)) $Bearing = 360 - $Bearing;
elseif (($StartLat > $EndLat) and ($StartLon >= $EndLon)) $Bearing = 180 + $Bearing;
elseif (($StartLat > $EndLat) and ($StartLon < $EndLon)) $Bearing = 180 - $Bearing;

$res['bearing'] = $Bearing;
return $res;
}

Комментариев нет:

Отправить комментарий

Архив блога

Ярлыки

linux (241) python (191) citation (186) web-develop (170) gov.ru (159) video (124) бытовуха (115) sysadm (100) GIS (97) Zope(Plone) (88) бурчалки (84) Book (83) programming (82) грабли (77) Fun (76) development (73) windsurfing (72) Microsoft (64) hiload (62) internet provider (57) opensource (57) security (57) опыт (55) movie (52) Wisdom (51) ML (47) driving (45) hardware (45) language (45) money (42) JS (41) curse (40) bigdata (39) DBMS (38) ArcGIS (34) history (31) PDA (30) howto (30) holyday (29) Google (27) Oracle (27) tourism (27) virtbox (27) health (26) vacation (24) AI (23) Autodesk (23) SQL (23) Java (22) humor (22) knowledge (22) translate (20) CSS (19) cheatsheet (19) hack (19) Apache (16) Manager (15) web-browser (15) Никонов (15) functional programming (14) happiness (14) music (14) todo (14) Klaipeda (13) PHP (13) course (13) scala (13) weapon (13) HTTP. Apache (12) SSH (12) frameworks (12) hero (12) im (12) settings (12) HTML (11) SciTE (11) USA (11) crypto (11) game (11) map (11) HTTPD (9) ODF (9) Photo (9) купи/продай (9) benchmark (8) documentation (8) 3D (7) CS (7) DNS (7) NoSQL (7) cloud (7) django (7) gun (7) matroska (7) telephony (7) Microsoft Office (6) VCS (6) bluetooth (6) pidgin (6) proxy (6) Donald Knuth (5) ETL (5) NVIDIA (5) Palanga (5) REST (5) bash (5) flash (5) keyboard (5) price (5) samba (5) CGI (4) LISP (4) RoR (4) cache (4) car (4) display (4) holywar (4) nginx (4) pistol (4) spark (4) xml (4) Лебедев (4) IDE (3) IE8 (3) J2EE (3) NTFS (3) RDP (3) holiday (3) mount (3) Гоблин (3) кухня (3) урюк (3) AMQP (2) ERP (2) IE7 (2) NAS (2) Naudoc (2) PDF (2) address (2) air (2) british (2) coffee (2) fitness (2) font (2) ftp (2) fuckup (2) messaging (2) notify (2) sharepoint (2) ssl/tls (2) stardict (2) tests (2) tunnel (2) udev (2) APT (1) CRUD (1) Canyonlands (1) Cyprus (1) DVDShrink (1) Jabber (1) K9Copy (1) Matlab (1) Portugal (1) VBA (1) WD My Book (1) autoit (1) bike (1) cannabis (1) chat (1) concurrent (1) dbf (1) ext4 (1) idioten (1) join (1) krusader (1) license (1) life (1) migration (1) mindmap (1) navitel (1) pneumatic weapon (1) quiz (1) regexp (1) robot (1) science (1) serialization (1) spatial (1) tie (1) vim (1) Науру (1) крысы (1) налоги (1) пианино (1)