通過經緯度坐標計算距離的方法（經緯度距離計算）
最近在網上搜索“通過經緯度坐標計算距離的方法”，發現網上大部分都是如下的代碼：
#define PI 3.14159265
static double Rc = 6378137; // 赤道半徑
static double Rj = 6356725; // 極半徑
class JWD
{
public:
double m_Longitude, m_Latitude;
double m_RadLo, m_RadLa;
double Ec;
double Ed;
public:
JWD(double longitude, double latitude)
{
m_Longitude = longitude;
m_Latitude = latitude;
m_RadLo = longitude * PI/180.;
m_RadLa = latitude * PI/180.;
Ec = Rj + (Rc - Rj) * (90.-m_Latitude) / 90.;
Ed = Ec * cos(m_RadLa);
}
~JWD() {};
};
static JWD GetJWDB(JWD A, double x,double y)
{
double dx=x;
double dy=y;
double BJD = (dx/A.Ed + A.m_RadLo) * 180./PI;
double BWD = (dy/A.Ec + A.m_RadLa) * 180./PI;
JWD B(BJD, BWD);
return B;
}
void main()
{
double referla=30.0;
double referlo=60.0;
double dx=500.0;
double dy=60.0;
JWD A(referla,referlo),B(0.0,0.0);
B=GetJWDB(A,dx,dy);
cout < < " LA = " < < B.m_Latitude < < " LO= " < < B.m_Longitude < < endl;
}
上面這段與之類似的代碼是最容易通過搜索引擎找到的。大部分都是相互的轉載，從來沒有說明和注釋。給初學者造成了很大的困擾。特別是：
Ec = Rj + (Rc - Rj) * (90.-m_Latitude) / 90.;
Ed = Ec * cos(m_RadLa);
Ec、Ed這2個參數，有人還在CSDN上發帖詢問“函數中Ec，Ed是什么意思”。但從來沒有見到有人回答。
提問的鏈接：
http://www.gisforum.net/bbs/TopicOther.asp?t=5&BoardID=33&id=155609
http://bbs.csdn.net/topics/320024634
我剛開始接觸這個問題，在中文世界中也只能搜到這些Ctrl+C 到Ctrl+V的東西。全球最大的中文互聯網上也沒有任何解答。已經明白的人知道很簡單，但初學者肯定很疑惑。更何況現在LBS應用已經非常多了。肯定有非常多的人已經找到了更好的解答方式。 但對於我來說，用最常用的關鍵詞，最容易的找到了這些復制的答案，但疑點重重。本着好奇的心，經過一番了解，我把結果給大家分析一下。下面是C#的代碼：
public const double Ea = 6378137; // 赤道半徑 public const double Eb = 6356725; // 極半徑 private static void GetJWDB(double GLAT, double GLON, double distance, double angle, out double BJD, out double BWD) { double dx = distance * 1000 * Math.Sin(angle * Math.PI / 180.0); double dy = distance * 1000 * Math.Cos(angle * Math.PI / 180.0); //double ec = 6356725 + 21412 * (90.0 - GLAT) / 90.0; // 21412 是赤道半徑與極半徑的差 double ec = Eb + (Ea-Eb) * (90.0 - GLAT) / 90.0; double ed = ec * Math.Cos(GLAT * Math.PI / 180); BJD = (dx / ed + GLON * Math.PI / 180.0) * 180.0 / Math.PI; BWD = (dy / ec + GLAT * Math.PI / 180.0) * 180.0 / Math.PI; }
上面這個函數一看就是懂中文的人搞的，名字叫GetJWDB，取得經緯度。他根據輸入的經度、緯度、距離和一個角度，得到另外一個經緯度坐標，輸出參數為BJD、BWD。
1）angle * Math.PI / 180.0 作用是將角度轉換為弧度，經緯度坐標是角度值，計算時需要換為弧度。這里所有的計算都是用弧度。
2）函數以正北方(due north) 也就是指南針的方向為0度，順時針方向增加。如下圖，Distance距離如果是d的話，dx就是x軸方向的長度，即longitude經度方向的長度；dy就是y軸方向的長度，即latitude緯度方向的長度。
dx、dy的計算方式也可以是以正東(due east)方向為0度。那：dx=distance* Cos(θ)，dy=distance*Sin(θ)。區別是Cos 與Sin互換。（圖上未標）
圖1
3）Ea 表示赤道半徑，Eb表示極半徑，地球是一個近似球體，Ea與Eb稍微有點差距。ec的作用就是修正因為緯度不斷變化的球半徑長度。如果在GLAT=0，即在赤道上的時候，ec=Eb+(Ea-Eb)*(90-0)/90=Ea，那ec就剛好是赤道半徑Ea；如果在極點GLAT=90，ec=Eb+(Ea-Eb)*(90-90)/90=Eb，那ec 就剛好是極半徑Eb。
4）Ed是GLAT所在緯度的緯度圈的半徑，如下圖：
圖2
截面過球心，此時截面的面積最大，此圓叫球的大圓(Great Cycle)，沿着經線進行截面，得到的都是大圓(Great Cycle)。球面被不經過球心. 的截面所截得的圓. 叫做小圓。緯度圈所在的圓是一個小圓。地球半徑R，平均值R=6371.0Km， Ed如果用地球半徑R表示，那就是Ed=R*Cos(θ)，可以參看
《根據2個經緯度點，計算這2個經緯度點之間的距離(通過經度緯度得到距離)》這里提到的“A、D所在緯度圓圈的半徑（AO`）”。放到上面函數里，因為它不斷修正地球半徑ec，那就是ed = ec * Math.Cos(GLAT * Math.PI / 180);
5）按照地球經緯度的划分，如下圖：
每等分的緯度之間，經度線段的長度是一定的。 A段，B段長度是一樣的。
每一等分的經度之間，緯度線段的長度是從赤道向2極點減小的。C端，D段的長度是不一樣的。
圖3
參看上圖，那dx / ed 就相當於是在GLAT這個緯度上dx長度與總長度的占比，算出來應該是個經度跨度。如果這個經度跨度加上起始給定的經度就是最終的經度。
同理 dy/R就是在GLON這個經度上的dy長度與地球平均半徑R的占比，算出來應該是一個緯度跨度。如果這個緯度跨度加上起始給定的緯度就是最終的緯度。這里使用了R，取地球平均半徑。
dy/ec 就是用不斷修正的ec半徑替代了平均半徑R。
(dy / ec + GLAT * Math.PI / 180.0) 就是起始緯度加上distance距離的最終緯度，同時需要將該結果轉換為角度。 轉換角度方法是：弧度* 180.0 / Math.PI。
BWD = (dy / ec + GLAT * Math.PI / 180.0) * 180.0 / Math.PI;
(dx / ed + GLON * Math.PI / 180.0)就是起始經度加上distance距離的最終經度，同時需要將該結果轉換為角度。
BJD = (dx / ed + GLON * Math.PI / 180.0) * 180.0 / Math.PI;
這個根據一個經緯度坐標、距離然后求另外一個經緯度坐標的作用，主要就是確定一個最小外包矩形(Minimum bounding rectangle，簡稱MBR)。例如，我要找一個坐標點(lat,lon)的5公里范圍內的所有商戶信息、景點信息等。這個MBR就是一個最大的范圍，這個矩形是包含5公里范圍內所有這些有效信息的一個最小矩形。利用公式，求出四個方向0度、90度、180度、270度方向上的四個坐標點就可以得到這個MBR。
public const double Ea = 6378137; // 赤道半徑 public const double Eb = 6356725; // 極半徑 private static void GetlatLon(double LAT, double LON, double distance, double angle, out double newLon, out double newLat) { double dx = distance * 1000 * Math.Sin(angle * Math.PI / 180.0); double dy = distance * 1000 * Math.Cos(angle * Math.PI / 180.0); double ec = Eb + (Ea-Eb) * (90.0 - LAT) / 90.0; double ed = ec * Math.Cos(LAT * Math.PI / 180); newLon = (dx / ed + LON * Math.PI / 180.0) * 180.0 / Math.PI; newLat = (dy / ec + LAT * Math.PI / 180.0) * 180.0 / Math.PI; } public static void GetRectRange(double centorlatitude, double centorLogitude, double distance, out double maxLatitude, out double minLatitude, out double maxLongitude, out double minLongitude) { double temp = 0.0; GetlatLon(centorlatitude, centorLogitude, distance, 0, out temp, out maxLatitude); GetlatLon(centorlatitude, centorLogitude, distance, 180, out temp, out minLatitude); GetlatLon(centorlatitude, centorLogitude, distance, 90, out maxLongitude, out temp); GetlatLon(centorlatitude, centorLogitude, distance, 270, out minLongitude , out temp); }
這里得到的maxLatitude, minLatitude, maxLongitude, minLongitude就組成矩形的四個頂點。
網上另外的方法，
http://www.movable-type.co.uk/scripts/latlong.html
這里的“Destination point given distance and bearing from start point”一節介紹了。我直接把代碼貼上來：
這里GetRectRange 這個函數 也是以正北方（due north）為起始角度，順時針方向，求得maxLatitude, minLatitude, maxLongitude, minLongitude 這樣一個MBR。2種方法的誤差很小。我覺得都是可以用的公式。
/// <summary> /// where φ is latitude, λ is longitude, θ is the bearing (clockwise from north), /// δ is the angular distance d/R; d being the distance travelled, R the earth’s radius /// bearing 方位 0，90，180，270 /// </summary> private static void GetNewLatLon(double lat, double lon, double d, double bearing, out double lat2, out double lon2) { lat2 = 0.0; lon2 = 0.0; double R = 6378.137; var φ1 = ConvertDegreesToRadians(lat); var λ1 = ConvertDegreesToRadians(lon); var θ = ConvertDegreesToRadians(bearing); var φ2 = Math.Asin(Math.Sin(φ1) * Math.Cos(d / R) + Math.Cos(φ1) * Math.Sin(d / R) * Math.Cos(θ)); var λ2 = λ1 + Math.Atan2(Math.Sin(θ) * Math.Sin(d / R) * Math.Cos(φ1), Math.Cos(d / R) - Math.Sin(φ1) * Math.Sin(φ2)); λ2 = (λ2 + 3 * Math.PI) % (2 * Math.PI) - Math.PI; // normalise to -180..+180° lat2 = ConvertRadiansToDegrees(φ2); lon2 = ConvertRadiansToDegrees(λ2); }
如果有一個應用，表里存有100萬的數據，數據包含一個lat、lon的經緯度信息。就可以先根據輸入的經緯度和距離得到一個MBR，然后通過類似
SELECT Id
FROM IdInfoTable
WHERE latitude >= minLat AND latitude < maxLat
AND longitude >= minLon AND longitude < maxLon
的方式過濾掉大部分的數據，然后再通過《根據2個經緯度點，計算這2個經緯度點之間的距離(通過經度緯度得到距離)》這里提到的方法進行精細過濾。
完整代碼：
using System; using System.Collections.Generic; using System.Linq; using System.Text; namespace GetJWD { public class GetMBR { public double MaxLatitude; public double MinLatitude; public double MaxLongitude; public double MinLongitude; public double MaxLatitude2; public double MinLatitude2; public double MaxLongitude2; public double MinLongitude2; public GetMBR(double centorlatitude, double centorLogitude, double distance) { GetRectRange(centorlatitude, centorLogitude, distance, out MaxLatitude, out MinLatitude, out MaxLongitude, out MinLongitude); GetRectRange2(centorlatitude, centorLogitude, distance, out MaxLatitude2, out MinLatitude2, out MaxLongitude2, out MinLongitude2); } public const double Ea = 6378137; // 赤道半徑 public const double Eb = 6356725; // 極半徑 private static void GetlatLon(double LAT, double LON, double distance, double angle, out double newLon, out double newLat) { double dx = distance * 1000 * Math.Sin(angle * Math.PI / 180.0); double dy = distance * 1000 * Math.Cos(angle * Math.PI / 180.0); double ec = Eb + (Ea - Eb) * (90.0 - LAT) / 90.0; double ed = ec * Math.Cos(LAT * Math.PI / 180); newLon = (dx / ed + LON * Math.PI / 180.0) * 180.0 / Math.PI; newLat = (dy / ec + LAT * Math.PI / 180.0) * 180.0 / Math.PI; } public static void GetRectRange(double centorlatitude, double centorLogitude, double distance, out double maxLatitude, out double minLatitude, out double maxLongitude, out double minLongitude) { double temp = 0.0; GetlatLon(centorlatitude, centorLogitude, distance, 0, out temp, out maxLatitude); GetlatLon(centorlatitude, centorLogitude, distance, 180, out temp, out minLatitude); GetlatLon(centorlatitude, centorLogitude, distance, 90, out maxLongitude, out temp); GetlatLon(centorlatitude, centorLogitude, distance, 270, out minLongitude, out temp); } public static void GetRectRange2(double centorlatitude, double centorLogitude, double distance, out double maxLatitude, out double minLatitude, out double maxLongitude, out double minLongitude) { double temp = 0.0; GetNewLatLon(centorlatitude, centorLogitude, distance, 0, out maxLatitude, out temp); GetNewLatLon(centorlatitude, centorLogitude, distance, 180, out minLatitude, out temp); GetNewLatLon(centorlatitude, centorLogitude, distance, 90, out temp, out maxLongitude); GetNewLatLon(centorlatitude, centorLogitude, distance, 270, out temp, out minLongitude); } /// <summary> /// where φ is latitude, λ is longitude, θ is the bearing (clockwise from north), /// δ is the angular distance d/R; d being the distance travelled, R the earth’s radius /// bearing 方位 0，90，180，270 /// </summary> private static void GetNewLatLon(double lat, double lon, double d, double bearing, out double lat2, out double lon2) { lat2 = 0.0; lon2 = 0.0; double R = 6378.137; var φ1 = ConvertDegreesToRadians(lat); var λ1 = ConvertDegreesToRadians(lon); var θ = ConvertDegreesToRadians(bearing); var φ2 = Math.Asin(Math.Sin(φ1) * Math.Cos(d / R) + Math.Cos(φ1) * Math.Sin(d / R) * Math.Cos(θ)); var λ2 = λ1 + Math.Atan2(Math.Sin(θ) * Math.Sin(d / R) * Math.Cos(φ1), Math.Cos(d / R) - Math.Sin(φ1) * Math.Sin(φ2)); λ2 = (λ2 + 3 * Math.PI) % (2 * Math.PI) - Math.PI; // normalise to -180..+180° lat2 = ConvertRadiansToDegrees(φ2); lon2 = ConvertRadiansToDegrees(λ2); } public static double ConvertDegreesToRadians(double degrees) { return degrees * Math.PI / 180; } public static double ConvertRadiansToDegrees(double radian) { return radian * 180.0 / Math.PI; } } class Test { static void Main(string[] args) { double latorg = 22.54587746, lonorg = 114.12873077; var gpsdis = new GetMBR(latorg, lonorg, 5); Console.WriteLine("maxlat:{0}, minlat:{1}, maxlon:{2}, minlon:{3}", gpsdis.MaxLatitude, gpsdis.MinLatitude, gpsdis.MaxLongitude, gpsdis.MinLongitude); Console.WriteLine("maxlat:{0}, minlat:{1}, maxlon:{2}, minlon:{3}", gpsdis.MaxLatitude2, gpsdis.MinLatitude2, gpsdis.MaxLongitude2, gpsdis.MinLongitude2); } } }
謹以此文紀念那篇CSDN上因為 “本帖子已過去太久遠了，不再提供回復功能。”而永遠至今晚為止都還沒有答案的帖子！
如今LBS應用泛濫，JavaScript到處都能看到源碼，gitHub上sourceCode隨處可見的時代，希望此文能引導后人，少走我的彎路。如果有更好的方案，也歡迎留言。值此慶祝五一佳節！