通过经纬度区间实现附近查询
最近在做一个地图功能，其中有个功能就是根据当前经纬度获筛选附近的公司，第一次尝试，没有这方面的经验积累，网上找到了一种解决方案
PostgreSQL & PostGIS 这让一直受微软关照的我瞬间不知从何入手
后来找到了一种叫geohash 的解决方案，马上将数据库中的经纬度转换成geohash
附上geohash 算法 C#
public static class Geohash { #region Direction enum public enum Direction { Top = 0, Right = 1, Bottom = 2, Left = 3 } #endregion private const string Base32 = "0123456789bcdefghjkmnpqrstuvwxyz"; private static readonly int[] Bits = new[] { 16, 8, 4, 2, 1 }; private static readonly string[][] neighbors = { new[] { "p0r21436x8zb9dcf5h7kjnmqesgutwvy", // Top "bc01fg45238967deuvhjyznpkmstqrwx", // Right "14365h7k9dcfesgujnmqp0r2twvyx8zb", // Bottom "238967debc01fg45kmstqrwxuvhjyznp", // Left }, new[] { "bc01fg45238967deuvhjyznpkmstqrwx", // Top "p0r21436x8zb9dcf5h7kjnmqesgutwvy", // Right "238967debc01fg45kmstqrwxuvhjyznp", // Bottom "14365h7k9dcfesgujnmqp0r2twvyx8zb", // Left } }; private static readonly string[][] Borders = { new[] {"prxz", "bcfguvyz", "028b", "0145hjnp"}, new[] {"bcfguvyz", "prxz", "0145hjnp", "028b"} }; public static String CalculateAdjacent(String hash, Direction direction) { hash = hash.ToLower(); char lastChr = hash[hash.Length - 1]; int type = hash.Length % 2; var dir = (int)direction; string nHash = hash.Substring(0, hash.Length - 1); if (Borders[type][dir].IndexOf(lastChr) != -1) { nHash = CalculateAdjacent(nHash, (Direction)dir); } return nHash + Base32[neighbors[type][dir].IndexOf(lastChr)]; } public static List<string> Neighbors(String hash) { List<string> result = new List<string>(); result.Add(CalculateAdjacent(hash, Direction.Top)); result.Add(CalculateAdjacent(hash, Direction.Bottom)); result.Add(CalculateAdjacent(hash, Direction.Right)); result.Add(CalculateAdjacent(hash, Direction.Left)); result.Add(CalculateAdjacent(result[3], Direction.Top)); result.Add(CalculateAdjacent(result[2], Direction.Top)); result.Add(CalculateAdjacent(result[2], Direction.Bottom)); result.Add(CalculateAdjacent(result[3], Direction.Bottom)); return result; } public static void RefineInterval(ref double[] interval, int cd, int mask) { if ((cd & mask) != 0) { interval[0] = (interval[0] + interval[1]) / 2; } else { interval[1] = (interval[0] + interval[1]) / 2; } } public static double[] Decode(String geohash) { bool even = true; double[] lat = { -90.0, 90.0 }; double[] lon = { -180.0, 180.0 }; foreach (char c in geohash) { int cd = Base32.IndexOf(c); for (int j = 0; j < 5; j++) { int mask = Bits[j]; if (even) { RefineInterval(ref lon, cd, mask); } else { RefineInterval(ref lat, cd, mask); } even = !even; } } return new[] { (lat[0] + lat[1]) / 2, (lon[0] + lon[1]) / 2 }; } public static String Encode(double latitude, double longitude, int precision = 12) { bool even = true; int bit = 0; int ch = 0; string geohash = ""; double[] lat = { -90.0, 90.0 }; double[] lon = { -180.0, 180.0 }; if (precision < 1 || precision > 20) precision = 12; while (geohash.Length < precision) { double mid; if (even) { mid = (lon[0] + lon[1]) / 2; if (longitude > mid) { ch |= Bits[bit]; lon[0] = mid; } else lon[1] = mid; } else { mid = (lat[0] + lat[1]) / 2; if (latitude > mid) { ch |= Bits[bit]; lat[0] = mid; } else lat[1] = mid; } even = !even; if (bit < 4) bit++; else { geohash += Base32[ch]; bit = 0; ch = 0; } } return geohash; } }
PHP：
<?php class Geohash { private $bitss = array(16, 8, 4, 2, 1); private $neighbors = array(); private $borders = array(); private $coding = "0123456789bcdefghjkmnpqrstuvwxyz"; private $codingMap = array(); public function __construct() { $this->neighbors['right']['even'] = 'bc01fg45238967deuvhjyznpkmstqrwx'; $this->neighbors['left']['even'] = '238967debc01fg45kmstqrwxuvhjyznp'; $this->neighbors['top']['even'] = 'p0r21436x8zb9dcf5h7kjnmqesgutwvy'; $this->neighbors['bottom']['even'] = '14365h7k9dcfesgujnmqp0r2twvyx8zb'; $this->borders['right']['even'] = 'bcfguvyz'; $this->borders['left']['even'] = '0145hjnp'; $this->borders['top']['even'] = 'prxz'; $this->borders['bottom']['even'] = '028b'; $this->neighbors['bottom']['odd'] = $this->neighbors['left']['even']; $this->neighbors['top']['odd'] = $this->neighbors['right']['even']; $this->neighbors['left']['odd'] = $this->neighbors['bottom']['even']; $this->neighbors['right']['odd'] = $this->neighbors['top']['even']; $this->borders['bottom']['odd'] = $this->borders['left']['even']; $this->borders['top']['odd'] = $this->borders['right']['even']; $this->borders['left']['odd'] = $this->borders['bottom']['even']; $this->borders['right']['odd'] = $this->borders['top']['even']; 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); $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); } private function calculateAdjacent($srcHash, $dir) { $srcHash = strtolower($srcHash); $lastChr = $srcHash[strlen($srcHash) - 1]; $type = (strlen($srcHash) % 2) ? 'odd' : 'even'; $base = substr($srcHash, 0, strlen($srcHash) - 1); if (strpos($this->borders[$dir][$type], $lastChr) !== false) { $base = $this->calculateAdjacent($base, $dir); } return $base . $this->coding[strpos($this->neighbors[$dir][$type], $lastChr)]; } public function neighbors($srcHash) { //$geohashPrefix = substr($srcHash, 0, strlen($srcHash) - 1); $neighbors['top'] = $this->calculateAdjacent($srcHash, 'top'); $neighbors['bottom'] = $this->calculateAdjacent($srcHash, 'bottom'); $neighbors['right'] = $this->calculateAdjacent($srcHash, 'right'); $neighbors['left'] = $this->calculateAdjacent($srcHash, 'left'); $neighbors['topleft'] = $this->calculateAdjacent($neighbors['left'], 'top'); $neighbors['topright'] = $this->calculateAdjacent($neighbors['right'], 'top'); $neighbors['bottomright'] = $this->calculateAdjacent($neighbors['right'], 'bottom'); $neighbors['bottomleft'] = $this->calculateAdjacent($neighbors['left'], 'bottom'); return $neighbors; } 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); } }
demo
$length = 6; $geohash = new Geohash(); $hash = substr($geohash->encode('39.981998','116.306785'), 0, $length); $box = $geohash->neighbors($hash); $box[] = $hash; $in_str = "'".implode("','", $box)."'"; echo $hash.'<br/>'; echo $in_str;
这样就会生成包含自己的九宫格的范围区间
原理图
但是生成9个geohash 去in感觉性能太差了，然后不好限定范围，所以弃用了这个方案，选了另一个方案，根据中心点坐标算出需要范围的四个经纬度数据
原理图
理论上一定范围是要基于圆形的，
范围是个内切圆，这样圆心到四个角的范围就超出了指定范围，但是地图上显示的应该是方形的，所以我觉得应该也不用太纠结这点问题，现在附上源代码
C#
1 public class Location 2 { 3 public double Latitude { get; set; } 4 public double Longitude { get; set; } 5 public Location(double lat,double lng) 6 { 7 Latitude = lat; 8 Longitude = lng; 9 } 10 } 11 12 public class CoordDispose 13 { 14 /** 地球半径 */ 15 private static double EARTH_RADIUS = 6378137.0; 16 /** 范围距离 */ 17 private double distance; 18 /** 左上角 */ 19 private Location left_top = null; 20 /** 右上角 */ 21 private Location right_top = null; 22 /** 左下角 */ 23 private Location left_bottom = null; 24 /** 右下角 */ 25 private Location right_bottom = null; 26 27 private CoordDispose(double distance) 28 { 29 this.distance = distance; 30 } 31 32 private static double degrees(double d) 33 { 34 return d * (180 / Math.PI); 35 } 36 37 private void getRectangle(double lat, double lng) 38 { 39 40 // float dlng = 2 * asin(sin(distance / (2 * EARTH_RADIUS)) / cos(lat)); 41 // float dlng = degrees(dlng) // 弧度转换成角度 42 double dlng = 2 * Math.Asin(Math.Sin(distance / (2 * EARTH_RADIUS)) / Math.Cos(radians(lat))); 43 dlng = degrees(dlng); 44 45 // dlat = distance / EARTH_RADIUS 46 // dlng = degrees(dlat) # 弧度转换成角度 47 double dlat = distance / EARTH_RADIUS; 48 dlat = degrees(dlat); // # 弧度转换成角度 49 50 // left-top : (lat + dlat, lng - dlng) 51 // right-top : (lat + dlat, lng + dlng) 52 // left-bottom : (lat - dlat, lng - dlng) 53 // right-bottom: (lat - dlat, lng + dlng) 54 left_top = new Location(lat + dlat, lng - dlng); 55 right_top = new Location(lat + dlat, lng + dlng); 56 left_bottom = new Location(lat - dlat, lng - dlng); 57 right_bottom = new Location(lat - dlat, lng + dlng); 58 59 } 60 61 /// <summary> 62 /// 计算两个经纬度之间的直接距离 63 /// </summary> 64 public static double GetDistance(Location lt1, Location lt2) 65 { 66 double radLat1 = radians(lt1.Longitude); 67 double radLat2 = radians(lt2.Longitude); 68 double a = radLat1 - radLat2; 69 double b = radians(lt1.Latitude) - radians(lt2.Latitude); 70 double s = 2 * Math.Asin(Math.Sqrt(Math.Pow(Math.Sin(a / 2), 2) + Math.Cos(radLat1) * Math.Cos(radLat2) * Math.Pow(Math.Sin(b / 2), 2))); 71 s = s * EARTH_RADIUS; 72 s = Math.Round(s * 10000) / 10000; 73 return s; 74 } 75 /// <summary> 76 /// 计算两个经纬度之间的直接距离(google 算法) 77 /// </summary> 78 public static double GetDistanceGoogle(Location lt1, Location lt2) 79 { 80 double radLat1 = radians(lt1.Longitude); 81 double radLng1 = radians(lt1.Latitude); 82 double radLat2 = radians(lt2.Longitude); 83 double radLng2 = radians(lt2.Latitude); 84 double s = Math.Acos(Math.Cos(radLat1) * Math.Cos(radLat2) * Math.Cos(radLng1 - radLng2) + Math.Sin(radLat1) * Math.Sin(radLat2)); 85 s = s * EARTH_RADIUS; 86 s = Math.Round(s * 10000) / 10000; 87 return s; 88 } 89 90 public static double hav(double theta) 91 { 92 double s = Math.Sin(theta / 2); 93 return s * s; 94 } 95 96 private static double radians(double d) 97 { 98 return d * Math.PI / 180.0; 99 } 100 101 public static double getDistance(double lat0, double lng0, double lat1, 102 double lng1) 103 { 104 // from math import sin, asin, cos, radians, fabs, sqrt 105 106 // def hav(theta): 107 // s = sin(theta / 2) 108 // return s * s 109 110 // def get_distance_hav(lat0, lng0, lat1, lng1): 111 // "用haversine公式计算球面两点间的距离。" 112 // # 经纬度转换成弧度 113 // lat0 = radians(lat0) 114 // lat1 = radians(lat1) 115 // lng0 = radians(lng0) 116 // lng1 = radians(lng1) 117 118 // dlng = fabs(lng0 - lng1) 119 // dlat = fabs(lat0 - lat1) 120 // h = hav(dlat) + cos(lat0) * cos(lat1) * hav(dlng) 121 // distance = 2 * EARTH_RADIUS * asin(sqrt(h)) 122 123 // return distance 124 125 lat0 = radians(lat0); 126 lat1 = radians(lat1); 127 lng0 = radians(lng0); 128 lng1 = radians(lng1); 129 130 double dlng = Math.Abs(lng0 - lng1); 131 double dlat = Math.Abs(lat0 - lat1); 132 double h = hav(dlat) + Math.Cos(lat0) * Math.Cos(lat1) * hav(dlng); 133 double distance = 2 * EARTH_RADIUS * Math.Asin(Math.Sqrt(h)); 134 135 return distance; 136 } 137 138 public static Location[] getDegreeCoordinates(double lat, double lng, 139 double distance) 140 { 141 CoordDispose llu = new CoordDispose(distance); 142 llu.getRectangle(lat, lng); 143 Location[] locations = new Location[4]; 144 locations[0] = llu.left_top; 145 locations[1] = llu.right_top; 146 locations[2] = llu.left_bottom; 147 locations[3] = llu.right_bottom; 148 return locations; 149 } 150 }
PHP：
<?php class Location { public $Latitude; public $Longitude; public function __construct($lat,$lng) { $this->Latitude = $lat; $this->Longitude = $lng; } public function __tostring() { return $this->Longitude.','.$this->Latitude; } } class CoordDispose { //地球半径(米) private static $EARTH_RADIUS = 6378137.0; //范围距离 private $distance; //左上角 private $left_top; //右上角 private $right_top; //左下角 private $left_bottom; //右下角 private $right_bottom; private function __construct($dist) { $this->distance = $dist; } public function getRectangle($lat,$lng){ $dlng = 2 * asin(sin($this->distance / (2 * self::$EARTH_RADIUS)) / cos(self::radians($lat))); $dlng = self::degrees($dlng); $dlat = $this->distance / self::$EARTH_RADIUS; $dlat = self::degrees($dlat); // # 弧度转换成角度 $this->left_top = new Location($lat + $dlat, $lng - $dlng); $this->right_top = new Location($lat + $dlat, $lng + $dlng); $this->left_bottom = new Location($lat - $dlat, $lng - $dlng); $this->right_bottom = new Location($lat - $dlat, $lng + $dlng); } public static function hav($theta){ $s = sin($theta / 2); return $s * $s; } //角度数转换为弧度公式 private static function radians($d){ $pi = pi(); return $d * $pi / 180.0; } //弧度转换为角度数公式 private static function degrees($d){ $pi = pi(); return $d * (180.0 / $pi); } //计算两个经纬度之间的直接距离 public static function getDistance($lt1,$lt2){ $radLat1 = self::radians($lt1->Longitude); $radLat2 = self::radians($lt2->Longitude); $a = $radLat1 - $radLat2; $b = self::radians($lt1->Lat) - self::radians($lt2->Lat); $s = 2 * asin(sqrt(pow(Sin($a / 2), 2) + cos($radLat1) * cos($radLat2) * pow(sin($b / 2), 2))); $s = $s * self::$EARTH_RADIUS; $s = round($s * 10000) / 10000; return $s; } //计算两个经纬度之间的直接距离(google 算法) public static function getDistanceGoogle($lt1,$lt2){ $radLat1 = self::radians($lt1->Longitude); $radLng1 = self::radians($lt1->Lat); $radLat2 = self::radians($lt2->Longitude); $radLng2 = self::radians($lt2->Lat); $s = acos(cos($radLat1) * cos($radLat2) * cos($radLng1 - $radLng2) + sin($radLat1) * sin($radLat2)); $s = $s * self::$EARTH_RADIUS; $s = round($s * 10000) / 10000; return $s; } //以一个经纬度为中心计算出四个顶点 //$lt1 中心点坐标 //$distance 距离 public static function getDegreeCoordinates($lt1,$distance){ $llu = new self($distance); $llu->getRectangle($lt1->Latitude,$lt1->Longitude); $ret = array(); $ret[] = $llu->left_top; $ret[] = $llu->right_top; $ret[] = $llu->left_bottom; $ret[] = $llu->right_bottom; return $ret; } } ?>
demo
$d = new Location(116.412007, 39.947545); $v = CoordDispose::getDegreeCoordinates($d, 2000); var_dump($v);
将左上角 和右下角的经纬用做区间查询就可以获得所需范围内的数据了，当然用数据库只是用来测试的
分享只是为了帮助需要的人，图片和文字描述不是很好，还请将就着看了
