Makale Özeti

GPS alıcıların ucuzlaması, özellikle akıllı telefonların birçoğunun GPS alıcısı bulundurması konuma duyarlı (locatin based) uygulamaların ciddi bir şekilde artmasına yol açtı. Bu yazımızda koordinatları verilen iki nokta arasındaki uzaklığın nasıl hesaplanacağını öğreneceğiz. Aynı zamanda çok fazla teknik detaylara inmeden GPS sistemleri için kullanılan yeryüzü modelinden bahsedeceğim.

Makale

GPS alıcıların ucuzlaması, özellikle akıllı telefonların birçoğunun GPS alıcısı bulundurması konuma duyarlı (locatin based) uygulamaların ciddi bir şekilde artmasına yol açtı. Bu yazımızda koordinatları verilen iki nokta arasındaki uzaklığın nasıl hesaplanacağını öğreneceğiz. Aynı zamanda çok fazla teknik detaylara inmeden GPS sistemleri için kullanılan yeryüzü modelinden bahsedeceğim.

Malumunuz dünyamız küre şeklinde olmayıp kutuplardan basık, ekvatordan hafifçe şişkin bir şekle sahiptir. Basit düzey hesaplar/uygulamalar için dünyanın şekli ellipsoid olarak kabul edilip hesaplama yapılabilirse de dünyanın yüzey şeklindeki farklılıklar (dağlar, ovalar, vadiler, denizler vs) göz önünde bulundurulursa bu basit varsayımın birçok uygulama için kabul edilmeyecek seviyede hata oranına sahip olduğu görülecektir.

Dünyanın birebir modelini çıkarmak hem coğrafi şartlardan dolayı hem de ortaya çıkacak astronomik büyüklükte veriyi saklayıp makul bir sürede işleyebilecek sistemlerin olmayışından dolayı (şimdilik) mümkün görünmemektedir. Bu sebeple yeryüzündeki engebelerin yumuşatılıp hesaplanması kolay bir hale getirilmiş geoit [1] dediğimiz matematiksel yeryüzü modelleri türetilmiştir.

Yukarıdaki şekilde örnek bir yeryüzü kesitinde numaralandırılmış alanlar sırasıyla şunları temsil etmektedir: 1-Okyanuslar, 2- Referans elipsoit, 3- Yerel şakül/çekül, 4- Kıta, 5- Geoit.

Tüm dünya tarafından kullanılan matematiksel yeryüzü modelleri olduğu gibi yerel olarak kullanılan matematiksel modellerde vardır. Örneğin Kuzey Amerika için North American Datum of 1983 (NAD83), Afrika için Clarke (1880) sistemi kullanılmaktadır [2]. Uluslararası düzeyde World Geodetic System 84 (WGS84) modeli kabul görmüş ve şuan ki Küresel Konumlama Sistemi’nde (GPS) bu model kullanılmaktadır [3]. Bu model küresel çapta ± 1 metrelik bir doğruluk hassasiyetine sahiptir [4].

1974 yılında Birleşik Devletler Hava Kuvvetleri (U.S. Air Force) için çalışan Polonyalı jeofizikçi (jeodezist) Thaddeus Vincenty, yeryüzündeki iki nokta arasındaki uzaklığı çok yüksek derece doğrulukta hesaplayan bir formül geliştirdi. Bu formül kullanılan matematiksel dünya modeli üzerinde 0.5mm hassasiyetinde sonuç vermektedir [5]. Bu dönem boyunca (2010 bahar) Mobile and Wireless Systems Lab dersi kapsamında mobil sistemler için geliştirdiğimiz uygulamada bu formülün Chriss Veness tarafından yazılan JavaScript kodunu C# diline çevirip projemizde kullandım. Yoğun matematiksel hesaplamalar performans açısından soru işareti oluşturabilir. Bizim projemizde bu tür yoğun hesaplama gerektiren işlemleri Azure Cloud Computing tarafında çalıştırıp mobil istemciye sadece sonuçları gönderdik. Bu nedenle performans açısından bir sorunla karşılaşmadık.

        /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  */

        /* Vincenty Inverse Solution of Geodesics on the Ellipsoid (c) Ozcan ILIKHAN 2010                 */

        /* Based on javascript version developed by Chris Veness 2010.                                    */

        /*                                                                                                */

        /* from: Vincenty inverse 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                                             */

        /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  */

 

        ///<summary>

        /// Calculates geodetic distance between two points specified by latitude/longitude using

        /// Vincenty inverse formula for ellipsoids (based on WGS-84).

        ///</summary>

        ///<param name="lat1"> Latitude of first point in decimal degrees</param>

        ///<param name="lon1"> Longitude of first point in decimal degrees</param>

        ///<param name="lat2"> Latitude of second point in decimal degrees</param>

        ///<param name="lon2"> Longitude of second point in decimal degrees</param>

        ///<returns>Distance in meters between points</returns>

 

        public static double VincentyDistance(double lat1, double lon1, double lat2, double lon2)

        {

            double a = 6378137, b = 6356752.3142,  f = 1/298.257223563;  // WGS-84 ellipsoid params

            double L = DegToRad(lon2 - lon1);

            double U1 = Math.Atan((1-f) * Math.Tan(DegToRad(lat1) ));

            double U2 = Math.Atan((1 - f) * Math.Tan(DegToRad(lat2)));

            double sinU1 = Math.Sin(U1), cosU1 = Math.Cos(U1);

            double sinU2 = Math.Sin(U2), cosU2 = Math.Cos(U2);

            

            double cosSigma;

            double sigma;

            double sinAlpha;

            double cosSqAlpha;

            double cos2SigmaM;

            double sinLambda;

            double sinSigma;

            double cosLambda;

 

            double lambda = L, lambdaP, iterLimit = 100;

           

            do

            {

                sinLambda = Math.Sin(lambda);

                cosLambda = Math.Cos(lambda);

                sinSigma = Math.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 = Math.Atan2(sinSigma, cosSigma);

                sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma;

                cosSqAlpha = 1 - sinAlpha*sinAlpha;

                cos2SigmaM = cosSigma - 2*sinU1*sinU2/cosSqAlpha;

               

                if (Double.IsNaN(cos2SigmaM))

                    cos2SigmaM = 0;  // equatorial line: cosSqAlpha=0 (§6)

 

                double 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)));

 

            } while (Math.Abs(lambda-lambdaP) > 1e-12 && --iterLimit>0);

 

            if (iterLimit == 0)

                return Double.NaN;  // formula failed to converge

 

            double uSq = cosSqAlpha * (a*a - b*b) / (b*b);

            double A = 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)));

            double B = uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq)));

            double deltaSigma = B*sinSigma*(cos2SigmaM+B/4*(cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)-

                B/6*cos2SigmaM*(-3+4*sinSigma*sinSigma)*(-3+4*cos2SigmaM*cos2SigmaM)));

            double s = b*A*(sigma-deltaSigma);

          

            return s;

        }

 

        ///<summary>

        /// Converts decimal degrees to radians.

        ///</summary>

        private static double DegToRad(double deg)

        {

            return (deg * Math.PI / 180.0);

        }

 

        ///<summary>

        /// Converts radians to decimal degrees.

        ///</summary>

        private static double RadToDeg(double rad)

        {

            return (rad / Math.PI * 180.0);

        }

 

Daha öncede değindiğim gibi Vincenty formülü çok yüksek derecede doğruluk oranına sahiptir. Bu derece yüksek bir hassasiyet oranına ihtiyaç duymuyorsanız ve hesaplamanın hızlı olması sizin için önemli ise Haversine formülünü [6] kullanmanızı tavsiye ederim.

Haversine formülü küresel trigonometrideki genel bir formülün özel bir durumudur. Aşağıdaki şekilde görüldüğü gibi bir küre üzerindeki u, v, ve w noktalarını birleştiren büyük daireler (great circles) küresel bir üçgen oluşturur.

Bu küresel üçgenin açılarını ve kenar uzunluklarını çözmek için çeşitli teoremler kullanılır. Haversine formülü aşağıdaki verilen kosinüs teoreminden türetilmiştir.

Bu kuraldan yola çıkarak haversin kuralı, , elde edilir. İki nokta arasındaki uzaklığı ise şu formülle elde ederiz:

Formüldeki haversin fonksiyonu açılımı:  haversin(θ) = sin2(θ/2) = (1−cos(θ))/2

İşlemi biraz daha basitleştirmek için haversin foksiyonunun veya sinüs foksiyonun tersi alınarak denkemi çözebiliriz:

Aslında ben matematiği çok severim ama bu yazının asıl konusunun programlama olduğunu göz önünde bulundurarak bu formüllerin nasıl türetildiğinin ve formüldeki sembollerin detaylarına girmiyorum. Merak edenler küresel trigonometri (spherical trigonometry) ile alakalı kaynaklara müracaat edebilirler [7], [8].

Inverse haversine kullanarak iki nokta arası uzaklığı veren fonksiyonumuz:

 
 public double HaversineDistance( double lat1, double lon1, double lat2, double lon2)
{
	// Radius of earth in meters 
	double R = 6371 * 1000;
 
	double deltaLat = DegToRad(lat2 - lat1);
	double deltaLon = DegToRad(lon2 - lon1);
	double a = Math.Sin(deltaLat / 2) * Math.Sin(deltaLat / 2) +
			Math.Cos(DegToRad(lat1)) * Math.Cos(DegToRad(lat2)) *
                        Math.Sin(deltaLon / 2) * Math.Sin(deltaLon / 2);
	double c = 2 * Math.Asin( Math.Min(1, Math.Sqrt(a)) );
	double d = R * c;
	return d;
}

Bu makalede GPS koordinatları verilen iki nokta arasındaki mesafenin nasıl hesaplanacağını ayrıca konu ile alalı genel bazı coğrafi bilgileri öğrenmiş olduk. Bir sonraki yazıda görüşmek dileğiyle...

Özcan İLİKHAN

PhD Student
Department of Computer Sciences
University of Wisconsin-Madison
http://pages.cs.wisc.edu/~ilikhan 

Kaynaklar:

[1]- Geoid, http://en.wikipedia.org/wiki/Geoid

[2]- Geodetic Systems, http://www.ngs.noaa.gov/PUBS_LIB/Geodesy4Layman/TR80003B.HTM

[3]- World Geodetic System, http://en.wikipedia.org/wiki/WGS84

[4]- The Global WGS84 Coordinate System, http://vgn.dm.gov.ae/DMEGOV/DVRS/DVRS_docs/THE_GLOBAL_WGS_84_COORDINATE.pdf

[5]- Thaddeus Vincenty, http://en.wikipedia.org/wiki/Thaddeus_Vincenty

[6]- Haversine formula, http://en.wikipedia.org/wiki/Haversine_formula

[7]- Wolfram's mathworld: Spherical Trigonometry, http://mathworld.wolfram.com/SphericalTrigonometry.html

[8]- Green, R. M. Spherical Astronomy. New York: Cambridge University Press, 1985, http://www.amazon.com/exec/obidos/ASIN/0521317797/ref=nosim/weisstein-20