LoginSignup
62
64

More than 5 years have passed since last update.

@kochizufan(Kohei OTSUKA)

一次元ハッシュコードによる空間半径検索

Posted at

注意

mysqlで空間検索をする場合、InnoDBでは、空間テーブルが作れないため、一次元ハッシュコードによる検索が必要になります。
MyISAMでも問題ない場合は、空間検索速度が速いので空間テーブルの利用を検討してください。

一次元ハッシュコード検索の基本

考え方の基本については、もっとも流行っている一次元ハッシュコードである、GeoHashについてのこちらの記事（ http://blog.masuidrive.jp/index.php/2010/01/13/geohash/ ）を確認してください。

  • 採用する一次元ハッシュコードの精度レベルいくらのメッシュが、求めようとする緯度で、メッシュの短辺が求めたい検索半径より大きくなるかを調べます。
  • 求めた精度レベルで、検索中心の経緯度のハッシュコードを求めます。
  • そのハッシュコードの、周辺8つ（自分自身も含めると9つ）のメッシュコードを求めます。
    9つ求める理由は、検索中心経緯度が中心メッシュのどこに位置しても、必ず検索半径が9つのメッシュの中に含まれるようにするためです。
  • 求めた9つのハッシュコードで、位置情報テーブルのハッシュコードカラムに前方一致検索をかけます。
  • 検索半径外のデータが表示されても問題なければ、この検索結果を表示すれば終了です。
    もし厳密に検索範囲内をチェックしたい場合は、この検索結果の全レコードに対して中心点からの距離を計算し、それでフィルタして正確な半径内データを洗い出します。

これが一次元ハッシュコードによる半径検索の基本的考え方です。

quadkey

世間では、GeoHashが流行っていますが、同じ一次元ハッシュコード検索をするならば、GeoHashより検索効率のよいquadkeyを使うべきと考えています。

quadkeyとは ⇒ http://msdn.microsoft.com/en-us/library/bb259689.aspx このURLの下半分、Tile Coordinates and Quadkeysあたり

quadkeyの方が優れている理由：

  • コードメッシュの形が、実距離に置いても正しく正方形
  • 地図ベースの技術なので当然、緯度が大きく違えば同じ精度レベルのメッシュでも辺の長さは異なりますが、基本的には1つのメッシュの中では、経度方向と緯度方向の長さが同じ、正方形になります。
    これは半径ベースの検索をする際、使いやすいと思います。
    GeoHashは元々長方形な上、メルカトル図法座標ではなく生経緯度から算出しているので、その割合も緯度ごとに大違いで、半径を含むメッシュサイズを求めるのがめんどいし、最適化されない。
  • メッシュ上位と下位の分割が、4分木でマイルド
    quadkeyは上の精度レベルから下の精度レベルへの分割が4分割なので、精度レベルの調整が効きやすく、無駄な範囲の検索が少なくなります。
    GeoHashは上の精度レベルから下の精度レベルへの分割が32分割なので、下の精度レベルでは足りないからと上のレベルに上げると一気に32倍の範囲を検索することになり、無駄が大きいです。
  • 隣接タイルの取得ロジックが簡単
  • メッシュのカバーする距離判定も簡単

といった点で優れています。

quadkeyの欠点は、GeoHashに比べてコードサイズが長くなってしまいデータ量が増大する（単純計算で同じ精度なら2.5倍）ことですが、きょうびBLOBデータも扱ったりする時代に8バイトコードが20バイトでもそこまで大差はないでしょうし、InnoDBであれば前方一致での検索にコード長さによるデメリットもないようなので（MyISAMだと若干ある模様）、気にする必要はないと考えています。

空間テーブル、GeoHash、quadkey間の比較

比較してみた結果を載せます。
条件は、

  • 日本の範囲に近い北緯20度～50度、東経120度～150度の領域に、ランダムに100万点のポイントを散りばめて、空間カラム、GeoHash、quadkeyの3値に直した形でレコードに持つテーブルを、MyISAMとInnoDBでそれぞれ準備。
    もちろん、システム非対応のためにインデックスを張れないInnoDBの空間カラム以外は、全てインデックスを張っています。
  • 走らせた検索は、東経135度、北緯35度を中心に半径10kmの検索。
  • 空間検索は空間インデックス/GeoHash/quadkey問わず、いきなり円内検索はできず、まず円を内包する矩形（quadkey,GeoHashでは近傍9メッシュ）で検索し、その後中心との距離を計算してフィルタしてやるしかないので、検索結果比較としては、以下を並べます。
    • 事前の矩形検索時の、検索速度の比較
    • 事前の矩形検索で、本来必要な検索結果数と比較してどれくらい余分に拾ってしまったか
    • その後の円内篩いまで実行した結果の、検索速度の比較（SQL上でやってます）
  • 検索速度は各々5回試行の平均、5回のうち最初2回は各sql交互に実行、3回は同じsqlを連続実行しています。
実施条件 ふるい前実行時間 ふるい前検索数 最終実行時間 最終検索数
空間検索:MyISAM 0.78ms 43 0.86ms 36
quadkey:MyISAM 1.96ms 269 3.22ms 36
GeoHash:MyISAM 1.52ms 664 7.04ms 36
quadkey:InnoDB 1.04ms 269 2.00ms 36
GeoHash:InnoDB 1.08ms 664 4.36ms 36

実験コード

上記実験の実行コードを載せます。

SQL検索条件取得コード

php上でのベンチマーク取得方法調査をサボったので、下記コードでWHERE条件を取得し、それをphpMyAdmin上で実行してテストしました。

SQL検索条件取得コード
<?php

require_once('./GeoHash.php');
require_once('./QuadKey.php');
require_once('./Spatial.php');

$gh = new GeoHash(35.0000000,135.0000000); //GeoHashのオブジェクトを東経135度、北緯35度を中心に生成
echo $gh->getWhere(10000) . "\n";          //上記で半径10kmの検索を行うWhere条件を表示

$qk = new QUadKey(35.0000000,135.0000000); //quadkeyのオブジェクトを東経135度、北緯35度を中心に生成
echo $qk->getWhere(10000) . "\n";          //上記で半径10kmの検索を行うWhere条件を表示

$sp = new Spatial(35.0000000,135.0000000); //空間検索のオブジェクトを東経135度、北緯35度を中心に生成
echo $sp->getWhere(10000) . "\n";          //上記で半径10kmの検索を行うWhere条件を表示
?>

実際に発行されたWHERE条件としては以下のような感じです。
いずれも、1行目が矩形条件、2行目が半径条件です。

空間検索：

空間検索
WHERE MBRWithin(`geometry`, GeomFromText('LineString(134.89033595285 34.910168471588, 135.10966404715 35.089831528412)', 4326) ) 
AND SQRT( POWER( ( X( `geometry` ) - 135 ) * 91187.58845252,2 ) + POWER( ( Y( `geometry` ) - 35 ) * 111319.49079327,2 ) ) < 10000

quadkey：

quadkey
WHERE ( `quadkey` LIKE '13211313131%' OR `quadkey` LIKE '13300202020%' OR `quadkey` LIKE '13300202021%' OR `quadkey` LIKE '13211313133%' OR `quadkey` LIKE '13300202022%' OR `quadkey` LIKE '13300202023%' OR `quadkey` LIKE '13211313311%' OR `quadkey` LIKE '13300202200%' OR `quadkey` LIKE '13300202201%' ) 
AND SQRT( POWER( ( X( `geometry` ) - 135 ) * 91187.58845252,2 ) + POWER( ( Y( `geometry` ) - 35 ) * 111319.49079327,2 ) ) < 10000

GeoHash：

GeoHash
WHERE ( `geohash` LIKE 'wyr8%' OR `geohash` LIKE 'wyrb%' OR `geohash` LIKE 'xn20%' OR `geohash` LIKE 'wypx%' OR `geohash` LIKE 'wypz%' OR `geohash` LIKE 'xn0p%' OR `geohash` LIKE 'wypw%' OR `geohash` LIKE 'wypy%' OR `geohash` LIKE 'xn0n%' ) 
AND SQRT( POWER( ( X( `geometry` ) - 135 ) * 91187.58845252,2 ) + POWER( ( Y( `geometry` ) - 35 ) * 111319.49079327,2 ) ) < 10000

quadkeyクラス

quadkeyクラス
<?php
/**
 * QuadKeyクラス
 *
 * GPS世界座標系の座標をQuadKeyに変換します。
 * 
 *
 * 2012/10/04 1.0      初期版
 * 2012/12/13 1.1      地球の端での境界条件を追加
 * 2012/12/19 1.2      変換ロジックを簡単化
 *
 * @package 
 * @access  public
 * @author  OHTSUKA Ko-hei
 * @create  2012/10/04
 * @version 1.2
 **/

class QuadKey {

    const M_ON_EQ       = 40075016.68557849;   //地球を半径6378.137kmの球として赤道での周(m)
    const DEG_TO_RAD    = 0.017453292519943295; //度をラジアンにするための定数

    const EarthRadius   = 6378137;
    const MinLatitude   = -85.05112878;
    const MaxLatitude   = 85.05112878;
    const MinLongitude  = -180;
    const MaxLongitude  = 180;

    private $codingMap=array();

    private $add = array( 'right'  =>  1,
                          'left'   => -1,
                          'top'    => -2,
                          'bottom' =>  2);

    public  $lat;
    public  $long;
    public  $hash;
    public  $level;

    public function QuadKey($alathash=null,$along=null,$alevel=null)
    {
        if ($alathash !== null) {
            if ($along !== null) {
                $this->setLatLng($alathash,$along,$alevel);
            } else {
                $this->setHash($alathash);
            }
        }
    }

    public function setLatLng($alat,$along,$alevel=null)
    {
        $this->lat   = $alat;
        $this->long  = $along;
        if ($alevel !== null) {
            $this->level = $alevel;
        } else {
            $this->level = 18;
        }

        $this->hash = $this->LatLngToQuadKey($this->lat,$this->long,$this->level);
    }

    public function setHash($ahash) {
        $this->hash = $ahash;
        $coord = $this->QuadKeyToLatLng($ahash);

        $this->lat   = $coord[0];
        $this->long  = $coord[1];

        $this->level = strlen($ahash);
    }

    public function getWhere($radius) {
        $this->setOptimizedLevel($radius);
        $list = $this->getNeighbors();

        if ($list === 'all_globe') {
            return 'WHERE 1 ';
        }

        //nullや重複があり得るようになったのでクリーニング
        $cleaning = array();

        foreach ($list as $qkey) {
            if ($qkey !== null) $cleaning[$qkey] = 1;
        }

        $where = join("%' OR `quadkey` LIKE '",array_keys($cleaning));
        $where = "WHERE ( `quadkey` LIKE '" . $where ."%' ) ";

        return $where;
    }

    /**
     * ┌─┬─┬─┐
     * │０│１│２│
     * ├─┼─┼─┤
     * │３│４│５│ 4=mine
     * ├─┼─┼─┤
     * │６│７│８│
     * └─┴─┴─┘
     * 
     */
    public function getNeighbors($srcHash=null) {
        if ($srcHash === null) {
            $srcHash = $this->hash;
        }
        // 地球全体とか大きすぎ
        if ($srcHash === '' || preg_match('/^[0-3]$/',$srcHash)) return 'all_globe';

        // 3x3の9マスを隣接として返す
        $matrix = array(NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);

        $matrix[4] = $srcHash;
        $matrix[1] = $this->calculateAdjacent($srcHash, 'top');
        $matrix[3] = $this->calculateAdjacent($srcHash, 'left');
        $matrix[5] = $this->calculateAdjacent($srcHash, 'right');
        $matrix[7] = $this->calculateAdjacent($srcHash, 'bottom');

        $matrix[0] = $this->calculateAdjacent($matrix[1], 'left');
        $matrix[6] = $this->calculateAdjacent($matrix[7], 'left');

        $matrix[2] = $this->calculateAdjacent($matrix[1], 'right');
        $matrix[8] = $this->calculateAdjacent($matrix[7], 'right');

        return $matrix;
    }

    /**
     * $dir = top / left / bottom / right
     */
    public function calculateAdjacent($srcHash, $dir=null) {
        if ($srcHash === null) return null;

        if ($dir === null) {
            $dir     = $srcHash;
            $srcHash = $this->hash;
        }

        //極域の境界条件
        if (preg_match('/^[01]+$/',$srcHash) && $dir == 'top') {
            return null;
        } else if (preg_match('/^[23]+$/',$srcHash) && $dir == 'bottom') {
            return null;
        }

        $tileXY = $this->QuadKeyToTileXY($srcHash);
        if ($dir == 'top') {
            $tileXY[1]--;
        } else if ($dir == 'bottom') {
            $tileXY[1]++;
        }
        if ($dir == 'left') {
            $tileXY[0]--;
        } else if ($dir == 'right') {
            $tileXY[0]++;
        }
        return $this->TileXYToQuadKey($tileXY[0],$tileXY[1],strlen($srcHash));
    }    

    /**
    * calcurate precise of given latlng
    */
    public function setOptimizedLevel($radius) {
        $this->level = $this->getOptimizedLevel($radius);

        $this->hash  = substr($this->hash,0,$this->level);
    }

    /**
    * calcurate optimized level of given radius
    */
    public function getOptimizedLevel($alat, $along=null, $radius=null) {
        if ($along === null) {
            $radius = $alat;
            $alat   = $this->lat;
            $along  = $this->long;
        }

        $norm_radius = $radius / cos($alat * self::DEG_TO_RAD);

        for ($lv=18;$lv>0;$lv--) {
            $lat_bit = self::M_ON_EQ / pow(2, $lv);

            if ($lat_bit >= $norm_radius) {
                return $lv;
            }
        }

        return 0;
    }

    private function _Clip($n, $minValue, $maxValue) {
        return min(max($n, $minValue), $maxValue);
    }

    public function MapSize($levelOfDetail) {
        return 256 << $levelOfDetail;
    }

    public function LatLngToPixelXY($lat, $lng, $levelOfDetail) {
        $lat = $this->_Clip($lat, self::MinLatitude, self::MaxLatitude);
        $lng = $this->_Clip($lng, self::MinLongitude, self::MaxLongitude);

        $x = ($lng + 180) / 360;
        $sinLat = sin($lat * M_PI / 180);
        $y = 0.5 - log((1 + $sinLat) / (1 - $sinLat)) / (4 * M_PI);

        $mapSize = $this->MapSize($levelOfDetail);
        $pixelX  = (int) $this->_Clip($x * $mapSize + 0.5, 0, $mapSize - 1);
        $pixelY  = (int) $this->_Clip($y * $mapSize + 0.5, 0, $mapSize - 1);

        return array($pixelX, $pixelY);
    }

    public function pixelXYToLatLng($pixelX, $pixelY, $levelOfDetail) {
        $mapSize = $this->MapSize($levelOfDetail);
        $x = ($this->_Clip($pixelX, 0, $mapSize - 1) / $mapSize) - 0.5;
        $y = 0.5 - ($this->_Clip($pixelY, 0, $mapSize - 1) / $mapSize);

        $lat = 90 - 360 * atan(exp(-$y * 2 * M_PI)) / M_PI;
        $lng = 360 * $x;

        return array($lat, $lng);
    }

    public function PixelXYToTileXY($pixelX, $pixelY) {
        $tileX = $pixelX / 256;
        $tileY = $pixelY / 256;
        return array((int)$tileX, (int)$tileY);
    }

    public function TileXYToPixelXY($tileX, $tileY) {
        $pixelX = $tileX * 256;
        $pixelY = $tileY * 256;
        return array($pixelX, $pixelY);
    }

    public function TileXYToQuadKey($tileX, $tileY, $levelOfDetail) {
        $xbin = substr(decbin($tileX), -$levelOfDetail);
        $ybin = substr(decbin($tileY), -$levelOfDetail);

        $quadKey = substr('00000000000000000000' . ($xbin + $ybin * 2) , -$levelOfDetail);

        return $quadKey;
    }

    public function QuadKeyToTileXY($quadKey) {
        $xbin  = preg_replace(array('/[13]/','/[02]/'),array('1','0'),$quadKey);
        $ybin  = preg_replace(array('/[01]/','/[23]/'),array('0','1'),$quadKey);

        $tileX = bindec($xbin);
        $tileY = bindec($ybin);

        return array($tileX, $tileY);
    }

    public function LatLngToTileXY($lat, $lng, $levelOfDetail) {
        list($px, $py) = $this->LatLngToPixelXY($lat, $lng, $levelOfDetail);
        return $this->PixelXYToTileXY($px, $py);
    }

    public function LatLngToQuadKey($lat, $lng, $levelOfDetail) {
        list($px, $py) = $this->LatLngToPixelXY($lat, $lng, $levelOfDetail);
        list($x, $y) = $this->PixelXYToTileXY($px, $py);

        $quadKey = $this->TileXYToQuadKey($x, $y, $levelOfDetail);
        return $quadKey;
    }

    public function QuadKeyToLatLng($quadKey) {
        list($x, $y) = $this->QuadKeyToTileXY($quadKey);
        list($pixelX, $pixelY) = $this->TileXYToPixelXY($x, $y);
        return $this->pixelXYToLatLng($pixelX, $pixelY, strlen($quadKey));
    }
}
?>

GeoHashクラス

経度180度をまたぐ場合や、北極南極付近で無限ループエラーが生じますが修正していません。

GeoHashクラス
<?php
/**
 * GeoHashクラス
 *
 * GPS世界座標系の座標をGeoHashに変換します。
 * 
 *
 * 2010/11/01 0.1      初期版
 * 2012/10/04 1.1      改変
 *
 * @package 
 * @access  public
 * @author  TakanoriYAMAZAKI
 * @modifier OHTSUKA Ko-hei
 * @create  2010/11/01
 * @modified 2012/10/04
 * @version 1.1
 **/

class GeoHash {

    const M_PER_DEGREE  = 111319.49079327357;   //緯度1度あたりのkm, 地球を半径378.137kmの球として
    const DEG_TO_RAD    = 0.017453292519943295; //度をラジアンにするための定数

    private $coding                                         ="0123456789bcdefghjkmnpqrstuvwxyz";
    private $codingMap=array();

    private $NEIGHBORS = array( 'right'  => array( 'even' => "bc01fg45238967deuvhjyznpkmstqrwx"),
                                'left'   => array( 'even' => "238967debc01fg45kmstqrwxuvhjyznp"),
                                'top'    => array( 'even' => "p0r21436x8zb9dcf5h7kjnmqesgutwvy"),
                                'bottom' => array( 'even' => "14365h7k9dcfesgujnmqp0r2twvyx8zb"));
    private $BORDERS   = array( 'right'  => array( 'even' => "bcfguvyz"),
                                'left'   => array( 'even' => "0145hjnp"),
                                'top'    => array( 'even' => "prxz"),
                                'bottom' => array( 'even' => "028b"));

    public  $lat;
    public  $long;
    public  $hash;
    public  $level;

    public function GeoHash($alathash=null,$along=null,$alevel=null)
    {
        //build map from encoding char to 0 padded bitfield
        for($i=0; $i<32; $i++)
        {
            $this->codingMap[substr($this->coding,$i,1)]=str_pad(decbin($i), 5, "0", STR_PAD_LEFT);
        }

        $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'];

        if ($alathash !== null) {
            if ($along !== null) {
                $this->setLatLng($alathash,$along,$alevel);
            } else {
                $this->setHash($alathash);
            }
        }
    }

    public function setLatLng($alat,$along,$alevel=null)
    {
        $this->lat   = $alat;
        $this->long  = $along;
        if ($alevel !== null) {
            $this->level = $alevel;
        } else {
            $this->level = 12;
        }

        $this->hash = substr($this->encode($this->lat,$this->long),0,$this->level);
    }

    public function setHash($ahash) {
        $this->hash = $ahash;
        $coord = $this->decode($ahash);

        $this->lat   = $coord[0];
        $this->long  = $coord[1];

        $this->level = strlen($ahash);
    }

    public function getWhere($radius) {
        $this->setOptimizedLevel($radius);
        $list = $this->getNeighbors();

        $where = join("%' OR `geohash` LIKE '",$list);
        $where = "WHERE ( `geohash` LIKE '" . $where ."%' ) ";

        return $where;
    }

    /**
     * ┌─┬─┬─┐
     * │０│１│２│
     * ├─┼─┼─┤
     * │３│４│５│ 4=mine
     * ├─┼─┼─┤
     * │６│７│８│
     * └─┴─┴─┘
     * 
     */
    public function getNeighbors($srcHash=null) {
        if ($srcHash === null) {
            $srcHash = $this->hash;
        }

        // 3x3の9マスを隣接として返す
        $matrix = array(NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);

        $matrix[4] = $srcHash;
        $matrix[1] = $this->calculateAdjacent($srcHash, 'top');
        $matrix[3] = $this->calculateAdjacent($srcHash, 'left');
        $matrix[5] = $this->calculateAdjacent($srcHash, 'right');
        $matrix[7] = $this->calculateAdjacent($srcHash, 'bottom');

        $matrix[0] = $this->calculateAdjacent($matrix[1], 'left');
        $matrix[6] = $this->calculateAdjacent($matrix[7], 'left');

        $matrix[2] = $this->calculateAdjacent($matrix[1], 'right');
        $matrix[8] = $this->calculateAdjacent($matrix[7], 'right');

        return $matrix;
    }

    /**
     * $dir = top / left / bottom / right
     */
    public function calculateAdjacent($srcHash, $dir=null) {
        if ($dir === null) {
            $dir     = $srcHash;
            $srcHash = $this->hash;
        }

        $dir = strtolower($dir);
        $srcHash = strtolower($srcHash);

        $lastChr = substr($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)];
    }

    /**
    * Decode a geohash and return an array with decimal lat,long in it
    */
    public function decode($hash)
    {
        //decode hash into binary string
        $binary="";
        $hl=strlen($hash);
        for($i=0; $i<$hl; $i++)
        {
            $binary.=$this->codingMap[substr($hash,$i,1)];
        }

        //split the binary into lat and log binary strings
        $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);

        }

        //now concert to decimal
        $lat=$this->binDecode($blat,-90,90);
        $long=$this->binDecode($blong,-180,180);

        //figure out how precise the bit count makes this calculation
        $latErr=$this->calcError(strlen($blat),-90,90);
        $longErr=$this->calcError(strlen($blong),-180,180);

        //how many decimal places should we use? There's a little art to
        //this to ensure I get the same roundings as geohash.org
        $latPlaces=max(1, -round(log10($latErr))) - 1;
        $longPlaces=max(1, -round(log10($longErr))) - 1;

        //round it
        $lat=round($lat, $latPlaces);
        $long=round($long, $longPlaces);

        return array($lat,$long);
    }

    /**
    * Encode a hash from given lat and long
    */
    public function encode($lat,$long)
    {
        //how many bits does latitude need?    
        /*$plat=$this->precision($lat);
        $latbits=1;
        $err=45;
        while($err>$plat)
        {
            $latbits++;
            $err/=2;
        }

        //how many bits does longitude need?
        $plong=$this->precision($long);
        $longbits=1;
        $err=90;
        while($err>$plong)
        {
            $longbits++;
            $err/=2;
        }

        //bit counts need to be equal
        $bits=max($latbits,$longbits);*/

        //as the hash create bits in groups of 5, lets not
        //waste any bits - lets bulk it up to a multiple of 5
        //and favour the longitude for any odd bits
        $bits= 30;
        $longbits=$bits;
        $latbits=$bits;
        $addlong=1;
        while (($longbits+$latbits)%5 != 0)
        {
            $longbits+=$addlong;
            $latbits+=!$addlong;
            $addlong=!$addlong;
        }

        //encode each as binary string
        $blat=$this->binEncode($lat,-90,90, $latbits);
        $blong=$this->binEncode($long,-180,180,$longbits);

        //merge lat and long together
        $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;
        }

        //convert binary string to hash
        $hash="";
        for ($i=0; $i<strlen($binary); $i+=5)
        {
            $n=bindec(substr($binary,$i,5));
            $hash=$hash.$this->coding[$n];
        }

        return $hash;
    }

    /**
    * What's the maximum error for $bits bits covering a range $min to $max
    */
    private function calcError($bits,$min,$max)
    {
        $err=($max-$min)/2;
        while ($bits--)
            $err/=2;
        return $err;
    }

    /*
    * returns precision of number
    * precision of 42 is 0.5
    * precision of 42.4 is 0.05
    * precision of 42.41 is 0.005 etc
    */
    private function precision($number)
    {
        $precision=0;
        $pt=strpos($number,'.');
        if ($pt!==false)
        {
            $precision=-(strlen($number)-$pt-1);
        }

        return pow(10,$precision)/2;
    }

    /**
    * create binary encoding of number as detailed in http://en.wikipedia.org/wiki/Geohash#Example
    * removing the tail recursion is left an exercise for the reader
    */
    private function binEncode($number, $min, $max, $bitcount)
    {
        if ($bitcount==0)
            return "";

        #echo "$bitcount: $min $max<br>";

        //this is our mid point - we will produce a bit to say
        //whether $number is above or below this mid point
        $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);
    }

    /**
    * decodes binary encoding of number as detailed in http://en.wikipedia.org/wiki/Geohash#Example
    * removing the tail recursion is left an exercise for the reader
    */
    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);
    }

    /**
    * calcurate precise of given latlng
    */
    public function setOptimizedLevel($radius) {
        $this->level = $this->getOptimizedLevel($radius);

        $this->hash  = substr($this->hash,0,$this->level);
    }

    /**
    * calcurate optimized level of given radius
    */
    public function getOptimizedLevel($alat, $along=null, $radius=null) {
        if ($along === null) {
            $radius = $alat;
            $alat   = $this->lat;
            $along  = $this->long;
        }

        $lat_m_per_deg = self::M_PER_DEGREE;
        $lng_m_per_deg = cos($alat * self::DEG_TO_RAD) * self::M_PER_DEGREE;

        for ($lv=12;$lv>0;$lv--) {
            $lng_bit = ceil(  5 * $lv / 2.0 );
            $lat_bit = floor( 5 * $lv / 2.0 );
            $lng_grid_size = (360.0 / pow(2, $lng_bit)) * $lng_m_per_deg;
            $lat_grid_size = (360.0 / pow(2, $lat_bit)) * $lat_m_per_deg;

            if ($lng_grid_size >= $radius && $lat_grid_size >= $radius) {
                return $lv;
            }
        }

        return 0;
    }
}
?>

空間検索クラス

空間検索クラス
<?php
/**
 * Spatialクラス
 *
 * GPS世界座標系の座標をSpatial検索するための矩形をはきます。
 * 
 *
 * 2012/10/04 0.1      初期版
 *
 * @package 
 * @access  public
 * @author  OHTSUKA Ko-hei
 * @create  2012/10/04
 * @version 0.1
 **/

class Spatial {

    const M_PER_DEGREE  = 111319.49079327357;   //緯度1度あたりのkm, 地球を半径378.137kmの球として
    const DEG_TO_RAD    = 0.017453292519943295; //度をラジアンにするための定数

    public  $lat;
    public  $long;

    public function Spatial($alat,$along)
    {
        $this->lat   = $alat;
        $this->long  = $along;
    }

    public function getWhere($radius) {
        $lat_m_per_deg = self::M_PER_DEGREE;
        $lng_m_per_deg = cos($this->lat * self::DEG_TO_RAD) * self::M_PER_DEGREE;

        $d_lat = $radius / $lat_m_per_deg;
        $d_lng = $radius / $lng_m_per_deg;

        $where = "WHERE MBRWithin(`geometry`, GeomFromText('LineString(" . ($this->long - $d_lng ) . " " . ($this->lat - $d_lat ) . ", " . ($this->long + $d_lng ) . " " . ($this->lat + $d_lat ) . ")', 4326)) ";
        $where .= "AND SQRT(POWER((X(`geometry`) - " . $this->long . ") * " . $lng_m_per_deg . ",2) + POWER((Y(`geometry`) - " . $this->lat . ") * " . $lat_m_per_deg . ",2)) < " . $radius;

        return $where;
    }
}
?>
62
64
1

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
62
64