1.场景
随着智能手机和传感器技术的发展，LBS（Location based service）类的应用也逐渐多了起来。几乎每一个软件或多或少都要牵扯上点位置的事。这篇文章主要来讲解下怎么快速的搜寻某个位置点周边的数据。
将问题模型化就是，给一组数据，每个数据包括了这个数据的位置信息（经纬度）和其他信息，给一个搜寻点（经度和纬度）和搜寻范围（radius），在最少的计算和时间内将在搜寻点的搜寻范围内的数据找出来。
2.目前的解决方案
（1）圆形区域一一匹配搜索。
将list中的数据一个一个同搜索位置进行距离计算。几乎没有人这样用，每一个数据都需要进行开根号操作，数据量大的时候不可想象。
（2）先对数据库中的经度（longitude）和纬度（latitude）添加B+树索引，然后搜索以radius的两倍为边长的正方形，如下图：
（3）使用目前现成的拥有地理位置搜索的数据库，比如mongoDB，HBase等等。
（4）使用geohash数据结构，见下面介绍
3.geohash介绍
（1）感性认识GeoHash
首先来点感性认识，http://openlocation.org/geohash/geohash-js/ 提供了在地图上显示geohash编码的功能。
1）GeoHash将二维的经纬度转换成字符串，比如下图展示了北京9个区域的GeoHash字符串，分别是WX4ER，WX4G2、WX4G3等等，每一个字符串代表了某一矩形区域。也就是说，这个矩形区域内所有的点（经纬度坐标）都共享相同的GeoHash字符串，这样既可以保护隐私（只表示大概区域位置而不是具体的点），又比较容易做缓存，比如左上角这个区域内的用户不断发送位置信息请求餐馆数据，由于这些用户的GeoHash字符串都是WX4ER，所以可以把WX4ER当作key，把该区域的餐馆信息当作value来进行缓存，而如果不使用GeoHash的话，由于区域内的用户传来的经纬度是各不相同的，很难做缓存。
2）字符串越长，表示的范围越精确。如图所示，5位的编码能表示10平方千米范围的矩形区域，而6位编码能表示更精细的区域（约0.34平方千米）
4.png
3）字符串相似的表示距离相近（特殊情况后文阐述），这样可以利用字符串的前缀匹配来查询附近的POI信息。如下两个图所示，一个在城区，一个在郊区，城区的GeoHash字符串之间比较相似，郊区的字符串之间也比较相似，而城区和郊区的GeoHash字符串相似程度要低些。
城区
郊区
通过上面的介绍我们知道了GeoHash就是一种将经纬度转换成字符串的方法，并且使得在大部分情况下，字符串前缀匹配越多的距离越近，回到我们的案例，根据所在位置查询来查询附近餐馆时，只需要将所在位置经纬度转换成GeoHash字符串，并与各个餐馆的GeoHash字符串进行前缀匹配，匹配越多的距离越近。
（2）GeoHash算法的步骤
下面以北海公园为例介绍GeoHash算法的计算步骤
7.png
2.1. 根据经纬度计算GeoHash二进制编码地球纬度区间是[-90,90]， 北海公园的纬度是39.928167，可以通过下面算法对纬度39.928167进行逼近编码:1）区间[-90,90]进行二分为[-90,0),[0,90]，称为左右区间，可以确定39.928167属于右区间[0,90]，给标记为1；2）接着将区间[0,90]进行二分为 [0,45),[45,90]，可以确定39.928167属于左区间 [0,45)，给标记为0；3）递归上述过程39.928167总是属于某个区间[a,b]。随着每次迭代区间[a,b]总在缩小，并越来越逼近39.928167；4）如果给定的纬度x（39.928167）属于左区间，则记录0，如果属于右区间则记录1，这样随着算法的进行会产生一个序列1011100，序列的长度跟给定的区间划分次数有关。
首先将纬度范围(-90, 90)平分成两个区间(-90, 0)、(0, 90)， 如果目标纬度位于前一个区间，则编码为0，否则编码为1。
由于39.92324属于(0, 90)，所以取编码为1。
然后再将(0, 90)分成 (0, 45), (45, 90)两个区间，而39.92324位于(0, 45)，所以编码为0。
以此类推，直到精度符合要求为止，得到纬度编码为1011 1000 1100 0111 1001。
经度也用同样的算法，对(-180, 180)依次细分，得到116.3906的编码为1101 0010 1100 0100 0100
2.2. 组码通过上述计算，纬度产生的编码为10111 00011，经度产生的编码为11010 01011。偶数位放经度，奇数位放纬度，把2串编码组合生成新串：11100 11101 00100 01111。最后使用用0-9、b-z（去掉a, i, l, o）这32个字母进行base32编码，首先将11100 11101 00100 01111转成十进制，对应着28、29、4、15，十进制对应的编码就是wx4g。同理，将编码转换成经纬度的解码算法与之相反，具体不再赘述。
8.png
（3）GeoHash Base32编码长度与精度
下表摘自维基百科：http://en.wikipedia.org/wiki/Geohash可以看出，当geohash base32编码长度为8时，精度在19米左右，而当编码长度为9时，精度在2米左右，编码长度需要根据数据情况进行选择。
（4）GeoHash算法
上文讲了GeoHash的计算步骤，仅仅说明是什么而没有说明为什么？为什么分别给经度和维度编码？为什么需要将经纬度两串编码交叉组合成一串编码？本节试图回答这一问题。如图所示，我们将二进制编码的结果填写到空间中，当将空间划分为四块时候，编码的顺序分别是左下角00，左上角01，右下脚10，右上角11，也就是类似于Z的曲线，当我们递归的将各个块分解成更小的子块时，编码的顺序是自相似的（分形），每一个子快也形成Z曲线，这种类型的曲线被称为Peano空间填充曲线。这种类型的空间填充曲线的优点是将二维空间转换成一维曲线（事实上是分形维），对大部分而言，编码相似的距离也相近， 但Peano空间填充曲线最大的缺点就是突变性，有些编码相邻但距离却相差很远，比如0111与1000，编码是相邻的，但距离相差很大。
10.png
除Peano空间填充曲线外，还有很多空间填充曲线，如图所示，其中效果公认较好是Hilbert空间填充曲线，相较于Peano曲线而言，Hilbert曲线没有较大的突变。为什么GeoHash不选择Hilbert空间填充曲线呢？可能是Peano曲线思路以及计算上比较简单吧，事实上，Peano曲线就是一种四叉树线性编码方式。
（5）使用注意点
1）由于GeoHash是将区域划分为一个个规则矩形，并对每个矩形进行编码，这样在查询附近POI信息时会导致以下问题，比如红色的点是我们的位置，绿色的两个点分别是附近的两个餐馆，但是在查询的时候会发现距离较远餐馆的GeoHash编码与我们一样（因为在同一个GeoHash区域块上），而较近餐馆的GeoHash编码与我们不一致。这个问题往往产生在边界处。
12.png
解决的思路很简单，我们查询时，除了使用定位点的GeoHash编码进行匹配外，还使用周围8个区域的GeoHash编码，这样可以避免这个问题。2）我们已经知道现有的GeoHash算法使用的是Peano空间填充曲线，这种曲线会产生突变，造成了编码虽然相似但距离可能相差很大的问题，因此在查询附近餐馆时候，首先筛选GeoHash编码相似的POI点，然后进行实际距离计算。
4.实践
geohash
<?php
class Geohash{
private $coding="0123456789bcdefghjkmnpqrstuvwxyz";
private $codingMap=array();
public function Geohash()
{
for($i=0; $i<32; $i++)
{
$this->codingMap[substr($this->coding,$i,1)]=str_pad(decbin($i), 5, "0", STR_PAD_LEFT);
}
}
public function decode($hash)
{
$binary="";
$hl=strlen($hash);
for($i=0; $i<$hl; $i++)
{
$binary.=$this->codingMap[substr($hash,$i,1)];
}
$bl=strlen($binary);
$blat="";
$blong="";
for ($i=0; $i<$bl; $i++)
{
if ($i%2)
$blat=$blat.substr($binary,$i,1);
else
$blong=$blong.substr($binary,$i,1);
}
$lat=$this->binDecode($blat,-90,90);
$long=$this->binDecode($blong,-180,180);
// $lat=$this->binDecode($blat,2,54);
// $long=$this->binDecode($blong,72,136);
$latErr=$this->calcError(strlen($blat),-90,90);
$longErr=$this->calcError(strlen($blong),-180,180);
$latPlaces=max(1, -round(log10($latErr))) - 1;
$longPlaces=max(1, -round(log10($longErr))) - 1;
$lat=round($lat, $latPlaces);
$long=round($long, $longPlaces);
return array($lat,$long);
}
public function encode($lat,$long)
{
$plat=$this->precision($lat);
$latbits=1;
$err=45;
while($err>$plat)
{
$latbits++;
$err/=2;
}
$plong=$this->precision($long);
$longbits=1;
$err=90;
while($err>$plong)
{
$longbits++;
$err/=2;
}
$bits=max($latbits,$longbits);
$longbits=$bits;
$latbits=$bits;
$addlong=1;
while (($longbits+$latbits)%5 != 0)
{
$longbits+=$addlong;
$latbits+=!$addlong;
$addlong=!$addlong;
}
$blat=$this->binEncode($lat,-90,90, $latbits);
$blong=$this->binEncode($long,-180,180,$longbits);
$binary="";
$uselong=1;
while (strlen($blat)+strlen($blong))
{
if ($uselong)
{
$binary=$binary.substr($blong,0,1);
$blong=substr($blong,1);
}
else
{
$binary=$binary.substr($blat,0,1);
$blat=substr($blat,1);
}
$uselong=!$uselong;
}
$hash="";
for ($i=0; $i<strlen($binary); $i+=5)
{
$n=bindec(substr($binary,$i,5));
$hash=$hash.$this->coding[$n];
}
return $hash;
}
private function calcError($bits,$min,$max)
{
$err=($max-$min)/2;
while ($bits--)
$err/=2;
return $err;
}
private function precision($number)
{
$precision=0;
$pt=strpos($number,'.');
if ($pt!==false)
{
$precision=-(strlen($number)-$pt-1);
}
return pow(10,$precision)/2;
}
private function binEncode($number, $min, $max, $bitcount)
{
if ($bitcount==0)
return "";
$mid=($min+$max)/2;
if ($number>$mid)
return "1".$this->binEncode($number, $mid, $max,$bitcount-1);
else
return "0".$this->binEncode($number, $min, $mid,$bitcount-1);
}
private function binDecode($binary, $min, $max)
{
$mid=($min+$max)/2;
if (strlen($binary)==0)
return $mid;
$bit=substr($binary,0,1);
$binary=substr($binary,1);
if ($bit==1)
return $this->binDecode($binary, $mid, $max);
else
return $this->binDecode($binary, $min, $mid);
}
}
测试
<?php
require_once("db_config.php");
require_once('geohash.class.php');
$geohash=new Geohash;
//经纬度转换成Geohash
//获取附近的信息
$n_latitude = 34.236080797698;
$n_longitude = 109.0145193757;
echo "当前位置为：经度108.7455，纬度34.3608<br/><br/>
以下网点离我最近：";
//开始
$b_time = microtime(true);
//方案A，直接利用数据库存储函数，遍历排序
//方案B geohash求出附近，然后排序
//当前 geohash值
$n_geohash = $geohash->encode($n_latitude,$n_longitude);
//附近
$n = 3;
$like_geohash = substr($n_geohash, 0, $n);
$sql = 'select * from retailersinfotable where geohash like
"'.$like_geohash.'%"';
$query = mysql_query($sql);
if(mysql_num_rows($query))
{
while($row=mysql_fetch_array($query))
{
$data[] = array (
"latitude" =>$row["Latitude"],
"longitude"=>$row["Longitude"],
"name"
=>$row["RetailersName"],
);
}
}
//算出实际距离
foreach($data as $key=>$val)
{
$distance = getDistance($n_latitude,$n_longitude,$val['latitude'],$val['longitude']
);
$data[$key]['distance'] = $distance;
//排序列
$sortdistance[$key] = $distance;
}
//距离排序
array_multisort($sortdistance,SORT_ASC,$data);
//结束
$e_time = microtime(true);
echo "（计算耗时：" ;
echo $e_time - $b_time;
echo "秒）<br/>";
//var_dump($data);
foreach($data as $key=>$val)
{
echo "</br>";
echo $val['distance']. " 米-------".$val['name'];
}
/**
* @desc 根据两点间的经纬度计算距离
* @param float $latitude 纬度值
* @param float $longitude 经度值
*/
function getDistance($latitude1, $longitude1, $latitude2, $longitude2)
{
$earth_radius = 6371000; //approximate radius of earth in meters
$dLat = deg2rad($latitude2 - $latitude1);
$dLon = deg2rad($longitude2 - $longitude1);
/*
Using the
Haversine formula
http://en.wikipedia.org/wiki/Haversine_formula
http://www.codecodex.com/wiki/Calculate_Distance_Between_Two_Points_on_a_Globe
验证：百度地图 http://developer.baidu.com/map/jsdemo.htm#a6_1
calculate the distance
*/
$a = sin($dLat/2) * sin($dLat/2) + cos(deg2rad($latitude1)) * cos(deg2rad($latitude2)) * sin($dLon/2) * sin($dLon/2);
$c = 2 * asin(sqrt($a));
$d = $earth_radius * $c;
return round($d); //四舍五入
}
?>
5.备注
在纬度相等的情况下：
经度每隔0.00001度，距离相差约1米；
每隔0.0001度，距离相差约10米；
每隔0.001度，距离相差约100米；
每隔0.01度，距离相差约1000米；
每隔0.1度，距离相差约10000米。
在经度相等的情况下：
纬度每隔0.00001度，距离相差约1.1米；
每隔0.0001度，距离相差约11米；
每隔0.001度，距离相差约111米；
每隔0.01度，距离相差约1113米；
每隔0.1度，距离相差约11132米。
6.参考
查找附近点--Geohash方案讨论
Haversine formula球面距离公式
球面距离公式代码实现
球面距离公式验证
?>
