1. 前言
??在我們平時(shí)使用美團(tuán)朽缴，餓了么等app進(jìn)行訂餐祟蚀，或者使用貓眼進(jìn)行訂電影票的時(shí)候，都有一個(gè)距離的排序科乎，表明該家店距離我們當(dāng)前的位置壁畸，這種基于地理位置的服務(wù)，統(tǒng)一被稱為L(zhǎng)BS（Location Based Service），而LBS的實(shí)現(xiàn)則是借助于GIS捏萍，WC（無(wú)線通信）等信息技術(shù)來實(shí)現(xiàn)太抓。而今天我們所要討論的就是這個(gè)距離的實(shí)現(xiàn)。
GIS令杈，Geographic information system腻异，地理信息系統(tǒng)。
2. 計(jì)算方式
??由于地球是一個(gè)橢圓形这揣，我們?cè)谟?jì)算的時(shí)候有點(diǎn)麻煩悔常，所以我們更常用的方式是將地球作為一個(gè)球形來計(jì)算，而計(jì)算球面上任意兩點(diǎn)之間的距離的公式通常有兩種：Great-circle distance和Haversine formula给赞，而目前大多數(shù)公司都是用的是Haversine公式机打，原因可以參考：
Great-circle distance公式用到了大量余弦函數(shù)， 而兩點(diǎn)間距離很短時(shí)（比如地球表面上相距幾百米的兩點(diǎn)）片迅，余弦函數(shù)會(huì)得出0.999…的結(jié)果残邀， 會(huì)導(dǎo)致較大的舍入誤差。而Haversine公式采用了正弦函數(shù)柑蛇，即使距離很小芥挣，也能保持足夠的有效數(shù)字。
而有關(guān)這兩者的介紹可以參考維基百科：Haversine formula 維基百科耻台，Great-circle distance 維基百科空免。而最終該公式的形式為：
至于為什么是這種形式，其實(shí)目前網(wǎng)上有許多推導(dǎo)公式盆耽，感興趣的可以看一下推導(dǎo)過程蹋砚，順便回憶一下自己當(dāng)年學(xué)過的數(shù)學(xué)知識(shí)：
1. 關(guān)于已知兩點(diǎn)經(jīng)緯度求球面最短距離的公式推導(dǎo)
2. 根據(jù)經(jīng)緯度計(jì)算兩點(diǎn)之間的距離的公式推導(dǎo)過程以及google.maps的測(cè)距函數(shù)
而如果要考慮到高度的影響的話，可以參考：https://stackoverflow.com/questions/3694380/calculating-distance-between-two-points-using-latitude-longitude-what-am-i-doi
??另外摄杂，還有一種方式是 Vincenty's formulae坝咐，該方式也是用于計(jì)算球體表面兩點(diǎn)之間距離的方式，而它所基于的就是地球是扁球體的形狀析恢，因此這種方式比假設(shè)地球是球體的方式應(yīng)該更加準(zhǔn)確墨坚，但實(shí)現(xiàn)起來比較麻煩。感興趣的可以查看下維基百科：Vincenty's formulae 維基百科
3. Java實(shí)現(xiàn)
接下來映挂，我們來看一下該公式的Java實(shí)現(xiàn)：
public final class DistanceUtils {
/**
* 地球半徑,單位 km
*/
private static final double EARTH_RADIUS = 6378.137;
/**
* 根據(jù)經(jīng)緯度泽篮，計(jì)算兩點(diǎn)間的距離
*
* @param longitude1 第一個(gè)點(diǎn)的經(jīng)度
* @param latitude1 第一個(gè)點(diǎn)的緯度
* @param longitude2 第二個(gè)點(diǎn)的經(jīng)度
* @param latitude2 第二個(gè)點(diǎn)的緯度
* @return 返回距離 單位千米
*/
public static double getDistance(double longitude1, double latitude1, double longitude2, double latitude2) {
// 緯度
double lat1 = Math.toRadians(latitude1);
double lat2 = Math.toRadians(latitude2);
// 經(jīng)度
double lng1 = Math.toRadians(longitude1);
double lng2 = Math.toRadians(longitude2);
// 緯度之差
double a = lat1 - lat2;
// 經(jīng)度之差
double b = lng1 - lng2;
// 計(jì)算兩點(diǎn)距離的公式
double s = 2 * Math.asin(Math.sqrt(Math.pow(Math.sin(a / 2), 2) +
Math.cos(lat1) * Math.cos(lat2) * Math.pow(Math.sin(b / 2), 2)));
// 弧長(zhǎng)乘地球半徑, 返回單位: 千米
s = s * EARTH_RADIUS;
return s;
}
public static void main(String[] args) {
double d = getDistance(116.308479, 39.983171, 116.353454, 39.996059);
System.out.println(d);
}
}
由于平時(shí)我們用到數(shù)學(xué)函數(shù)的地方不多，所以這里我們來簡(jiǎn)單介紹下用到的幾個(gè)函數(shù)：
Math.pow(x,y) //這個(gè)函數(shù)是求x的y次方
Math.toRadians //將一個(gè)角度測(cè)量的角度轉(zhuǎn)換成以弧度表示的近似角度
Math.sin //正弦函數(shù)
Math.cos //余弦函數(shù)
Math.sqrt //求平方根函數(shù)
Math.asin //反正弦函數(shù)
由于三角函數(shù)中特定的關(guān)聯(lián)關(guān)系袖肥，Haversine公式的最終實(shí)現(xiàn)方式可以有多種咪辱，比如借助轉(zhuǎn)角度的函數(shù)atan2：
public static double getDistance2(double longitude1, double latitude1,
double longitude2, double latitude2) {
double latDistance = Math.toRadians(longitude1 - longitude2);
double lngDistance = Math.toRadians(latitude1 - latitude2);
double a = Math.sin(latDistance / 2) * Math.sin(latDistance / 2)
+ Math.cos(Math.toRadians(longitude1)) * Math.cos(Math.toRadians(longitude2))
* Math.sin(lngDistance / 2) * Math.sin(lngDistance / 2);
double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));
return c * EARTH_RADIUS;
}
返回的單位是千米振劳，如果想返回米椎组，可以修改地球半徑的單位從千米到米，并且由于該結(jié)果是double類型的历恐，所以還可以借助Math.round方法進(jìn)行四舍五入為long類型寸癌，然后精確到米：
// ......
// 弧長(zhǎng)乘地球半徑(6378137)
s = s * EARTH_RADIUS;
// 返回類型: long专筷，單位: 米
return Math.round(s * 10000) / 10000;
接下來說幾點(diǎn)概念：
3.1 地球半徑
??由于地球不是一個(gè)完美的球體，所以并不能用一個(gè)特別準(zhǔn)確的值來表示地球的實(shí)際半徑蒸苇，不過由于地球的形狀很接近球體磷蛹，用[6357km] 到 [6378km]的范圍值可以涵蓋需要的所有半徑。并且通常情況下溪烤，地球半徑有幾個(gè)常用值：
- 極半徑味咳，從地球中心至南極或北極的距離， 相當(dāng)于6356.7523km檬嘀；
- 赤道半徑槽驶，從地球中心到赤道的距離，大約6378.137km鸳兽；
- 平均半徑掂铐，6371.393km，表示地球中心到地球表面所有各點(diǎn)距離的平均值揍异；
- RE全陨，地球半徑，有時(shí)被使用作為距離單位, 特別是在天文學(xué)和地質(zhì)學(xué)中常用衷掷，大概距離是6370.856km辱姨；
所以我們通過地球半徑進(jìn)行計(jì)算的時(shí)候，通常情況下戚嗅，我們可以使用上面的每一個(gè)值都可以進(jìn)行計(jì)算炮叶，不過或多或少都會(huì)有誤差的，但這樣的誤差是也是允許存在的渡处。這里參考自維基百科：維基百科-地球半徑
4. MySQL實(shí)現(xiàn)
同樣镜悉，在MySQL中實(shí)現(xiàn)該功能，計(jì)算公式還是通過Haversine公式医瘫。不過在Google Map中侣肄，已經(jīng)提供了相應(yīng)的實(shí)現(xiàn)方式，我們先來看一下醇份。
4.1 Google Map實(shí)現(xiàn)
首先稼锅，我們需要先創(chuàng)建表結(jié)構(gòu)：
CREATE TABLE `markers` (
`id` INT NOT NULL AUTO_INCREMENT PRIMARY KEY ,
`name` VARCHAR( 60 ) NOT NULL ,
`address` VARCHAR( 80 ) NOT NULL ,
`lat` FLOAT( 10, 6 ) NOT NULL ,
`lng` FLOAT( 10, 6 ) NOT NULL
) ENGINE = MYISAM ;
當(dāng)然存儲(chǔ)引擎可以是InnoDB。然后僚纷，進(jìn)行初始化數(shù)據(jù)：
INSERT INTO `markers` (`id`, `name`, `address`, `lat`, `lng`) VALUES ('1','Heir Apparel','Crowea Pl, Frenchs Forest NSW 2086','-33.737885','151.235260');
INSERT INTO `markers` (`id`, `name`, `address`, `lat`, `lng`) VALUES ('2','BeeYourself Clothing','Thalia St, Hassall Grove NSW 2761','-33.729752','150.836090');
INSERT INTO `markers` (`id`, `name`, `address`, `lat`, `lng`) VALUES ('3','Dress Code','Glenview Avenue, Revesby, NSW 2212','-33.949448','151.008591');
INSERT INTO `markers` (`id`, `name`, `address`, `lat`, `lng`) VALUES ('4','The Legacy','Charlotte Ln, Chatswood NSW 2067','-33.796669','151.183609');
INSERT INTO `markers` (`id`, `name`, `address`, `lat`, `lng`) VALUES ('5','Fashiontasia','Braidwood Dr, Prestons NSW 2170','-33.944489','150.854706');
INSERT INTO `markers` (`id`, `name`, `address`, `lat`, `lng`) VALUES ('6','Trish & Tash','Lincoln St, Lane Cove West NSW 2066','-33.812222','151.143707');
INSERT INTO `markers` (`id`, `name`, `address`, `lat`, `lng`) VALUES ('7','Perfect Fit','Darley Rd, Randwick NSW 2031','-33.903557','151.237732');
INSERT INTO `markers` (`id`, `name`, `address`, `lat`, `lng`) VALUES ('8','Buena Ropa!','Brodie St, Rydalmere NSW 2116','-33.815521','151.026642');
INSERT INTO `markers` (`id`, `name`, `address`, `lat`, `lng`) VALUES ('9','Coxcomb and Lily Boutique','Ferrers Rd, Horsley Park NSW 2175','-33.829525','150.873764');
INSERT INTO `markers` (`id`, `name`, `address`, `lat`, `lng`) VALUES ('10','Moda Couture','Northcote Rd, Glebe NSW 2037','-33.873882','151.177460');
然后就可以根據(jù)經(jīng)緯度值矩距，然后基于Haversine公式來查詢數(shù)據(jù)，假設(shè)我們要查詢latitude=37.38714,longitude=-122.083235怖竭，范圍在25英里內(nèi)的前20條數(shù)據(jù)：
SELECT id, ( 3959 * acos( cos( radians(37) ) * cos( radians( lat ) ) * cos( radians( lng ) - radians(-122) ) + sin( radians(37) ) * sin( radians( lat ) ) ) ) AS distance
FROM markers
HAVING distance < 25
ORDER BY distance
LIMIT 0 , 20;
而如果我們要查詢公里锥债，將3959英里也就是地球半徑，修改為6371即可。
Google Maps地址：Creating a Store Locator on Google Maps php&MySQL
4.2 st_distance函數(shù)
??MySQL其實(shí)在很早就提供了這種存儲(chǔ)經(jīng)緯度及相關(guān)運(yùn)算的功能哮肚，這種數(shù)據(jù)類型叫做空間數(shù)據(jù)類型登夫，而對(duì)應(yīng)的索引被稱為空間索引，但由于MySQL之前的版本對(duì)InnoDB支持的并不是太好允趟，所以使用的并不多恼策。不過MySQL5.6和MySQL5.7對(duì)此進(jìn)行了優(yōu)化，添加了st_distance等相關(guān)函數(shù)來支持經(jīng)緯度相關(guān)的計(jì)算潮剪。
??這里只來看一下st_distance函數(shù)的使用涣楷，其他相關(guān)的函數(shù)我會(huì)專門寫一篇文章來學(xué)習(xí)。我們還是拿上面Google Maps所建的表來測(cè)試抗碰，來按照距離進(jìn)行查詢：
SELECT
s.*,
(st_distance(point(lng, lat), point(-122.083235, 37.38714) ) * 111195) AS distance
FROM
markers s
ORDER BY
distance
其中总棵，point是MySQL的空間數(shù)據(jù)類型改含，先不多說這塊。就這樣捍壤，我們只需要通過st_distance函數(shù)就計(jì)算出了我們所需要查詢的結(jié)果，不過這里需要說一下：
- st_distance 函數(shù)返回的單位是degrees鹃觉，也就是空間單位的度數(shù)，我們?nèi)绻獙egrees轉(zhuǎn)換為米或者千米的話盗扇，需要乘以
EARTH_RADIUS * PI/180， EARTH_RADIUS 也就是地球半徑疗隶，至于是米還是千米，就看該變量的單位斑鼻。- 該運(yùn)算其實(shí)就相當(dāng)于對(duì)地球半徑進(jìn)行弧度與角度的轉(zhuǎn)換，也就是
Math.toRadians坚弱，而上面我們寫的111195其實(shí)是一個(gè)有誤差的值蜀备，該值就是通過該計(jì)算得出的結(jié)果；我們可以簡(jiǎn)單看一下toRadians實(shí)現(xiàn)：
public static double toRadians(double angdeg) {
return angdeg / 180.0 * PI;
}
這里的轉(zhuǎn)換參考自：Stackoverflow - Get Distance in Meters instead of degrees in Spatialite
其實(shí)杠氢，MySQL有提供直接查詢結(jié)果是米的函數(shù)：
st_distance_sphere，并且該函數(shù)的計(jì)算結(jié)果要比st_distance轉(zhuǎn)換為米的結(jié)果更精確脂凶。不過該函數(shù)是MySQL5.7之后才引入的宪睹，5.7之前還是需要通過計(jì)算轉(zhuǎn)換成米。更多可參考官方文檔地址：
MySQL 5.7 ST_Distance_Sphere(g1, g2 [, radius])
5 Geohash算法
??Geohash是目前比較主流的范圍搜索的算法艰猬，比如說搜索附近500米內(nèi)的地點(diǎn)這種問題。Geohash算法將二維的經(jīng)緯度編碼為一個(gè)字符串埋市，每個(gè)字符串代表了某一矩形區(qū)域冠桃，也就是說，這個(gè)矩形區(qū)域內(nèi)所有的點(diǎn)（經(jīng)緯度坐標(biāo)）都共享相同的GeoHash字符串道宅，這樣在查詢的時(shí)候就可以對(duì)該字符串做索引食听，然后根據(jù)該字符串進(jìn)行過濾。
Geohash算法的最大用途其實(shí)就是附近地址搜索了污茵，不過樱报，從geohash的編碼算法中可以看出它的一個(gè)缺點(diǎn)，也就是邊界問題：雖然兩個(gè)地點(diǎn)距離很近泞当，但恰好位于分界點(diǎn)的兩側(cè)迹蛤，這樣geohash字符串就會(huì)不相同，然后匹配的時(shí)候就會(huì)有問題襟士。不過要解決這個(gè)問題也很簡(jiǎn)單盗飒，就是計(jì)算的時(shí)候，計(jì)算出8個(gè)分別分布在周圍8個(gè)區(qū)域的地點(diǎn)陋桂。
??在實(shí)際應(yīng)用中逆趣，可以先根據(jù)Geohash篩選出附近的地點(diǎn)，然后再算出距離附近地點(diǎn)的距離嗜历。而如果要計(jì)算Geohash宣渗，可以通過 spatial4j 工具包來實(shí)現(xiàn)，GeohashUtils.encodeLatLon(lat, lon)梨州，默認(rèn)精度是12位痕囱，其中l(wèi)ucene就使用了spatial4j工具包來計(jì)算距離。
<dependency>
<groupId>org.locationtech.spatial4j</groupId>
<artifactId>spatial4j</artifactId>
<version>0.7</version>
</dependency>
有管GeoHash算法暴匠，可參考：
1. Geohash - 維基百科
2. GeoHash介紹-核心原理解析
3. Github-Java實(shí)現(xiàn)Geohash算法- github.com/GongDexing/Geohash
6. 其他
其實(shí)實(shí)現(xiàn)距離的方式有好多種咐蝇，比如說：
- mysql sql查詢
- mysql+geohash
- mysql 空間索引 (MySQL5.7版本以上)
- PostgreSQL/mongodb + geohash
- redis+geohash
- Lucene/Solr/ES + Spatial/geohash
??并且，這種基于搜索排序的功能其實(shí)正是Lucene這種搜索引擎和非關(guān)系型數(shù)據(jù)庫(kù)所擅長(zhǎng)的巷查。而對(duì)MySQL而言有序，一直以來MySQL在GIS上的功能支持都比較弱，并且僅有MyISAM引擎支持岛请，不過MySQL5.7之后發(fā)生了改變旭寿，提供了InnoDB引擎的GIS支持。所以崇败，針對(duì)MySQL的這塊功能肩祥，等接下來專門來學(xué)習(xí)一下混狠。
本文參考自：
1. 幾個(gè)地理位置信息處理方案的對(duì)比和分析
2. 空間索引 - 各數(shù)據(jù)庫(kù)空間索引使用報(bào)告
3. 美團(tuán)技術(shù)團(tuán)隊(duì)-地理空間距離計(jì)算優(yōu)化
