はじめに

こんにちは。はじめまして。tarokamikazeです。

担当するプロダクトで急に「GoogleMapにアレをプロットしたいよね」とか「位置情報取って、近くにあるアレを表示しといて」という、ゆるふわな要件が来るときが一生に一度はあると思います。
そんなあなたに贈る、ゆるふわジオコーディング記事です。

基礎知識

住所

東京都千代田区千代田１−１ みたいな文字列です。
都道府県・市区町村区切りで「近く」を判定できるような要件だったらいいのですが、そんなに人生は甘くありません。ええ、甘くはなかったのです。

緯度経度

(緯度/latitude)35.685323, (経度/longitude)139.752768 みたいなfloatです。
住所から緯度経度への変換は、GooglePlaceAPI等を使えばさくっとできます。

GeoJSON

{"type":"Feature","geometry":{"type":"Point","coordinates":[139.752768,35.685323]},"properties":null} みたいなJSONです。
この例では一点だけを表現していますが、線とかエリアとかいろいろ表現できます。
何が便利って、 http://geojson.io/ にコピペするだけでカンタンに確認できることです！
位置情報関係の開発を行う際は必須と言っても過言ではないでしょう。

緯度経度検索のとっつきづらさ

たいていの要件について、「近く」というのは「特定の緯度経度からXメートル以内」とか「モバイルアプリで表示中の地図範囲」だったりします。
検索対象に緯度経度を付与しておけば、できそうな気がしますね。が、具体的にどうすればいいのでしょう。

皆様は、これらを検索するためのスキーマ設計やクエリをさくさく書けるでしょうか。
難しいんじゃないかと思います。

FYI: MySQLならば

MySQL MyISAMを使える環境の方にとっては、geometory型で一発なのですが、そんな幸せな方は多くはないでしょう。
Datastoreのようにクエリが貧弱なDBを用いている方は目も当てられません。

こちらの記事では、geometry型を使わないクエリの投げ方を紹介されています。

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;

文系エンジニアがメンテするには、難しすぎますね。

検索インデックスをつかおう(Quadkey)

世の中には「地図を一定区域に区切って一意の名前をつけちゃえ」という考え方があります。
geohashが最も有名ですが、こいつは切り方が長方形で、なんとも使いづらい。
ほぼ正方形できってくれて、文字列自体も扱いやすいQuadkey がおすすめです。

Quadkeyとは、地図を正方形に切っていって、0~3の4進法で表現しちゃおうというものです。
レベルが上がる（1正方形あたりの面積が1/4になる)ごとに、桁数が増えます。
レベル１６だと 1330021123102203 など。実際にどう切られるかはこちらで確認できます。

一次元ハッシュコードによる空間半径検索という記事であげられているメリットが、まさしくそのとおり！といった感じです。

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

コードメッシュの形が、実距離に置いても正しく正方形

地図ベースの技術なので当然、緯度が大きく違えば同じ精度レベルのメッシュでも辺の長さは異なりますが、基本的には1つのメッシュの中では、経度方向と緯度方向の長さが同じ、正方形になります。
これは半径ベースの検索をする際、使いやすいと思います。
GeoHashは元々長方形な上、メルカトル図法座標ではなく生経緯度から算出しているので、その割合も緯度ごとに大違いで、半径を含むメッシュサイズを求めるのがめんどいし、最適化されない。

メッシュ上位と下位の分割が、4分木でマイルド
quadkeyは上の精度レベルから下の精度レベルへの分割が4分割なので、精度レベルの調整が効きやすく、無駄な範囲の検索が少なくなります。
GeoHashは上の精度レベルから下の精度レベルへの分割が32分割なので、下の精度レベルでは足りないからと上のレベルに上げると一気に32倍の範囲を検索することになり、無駄が大きいです。

隣接タイルの取得ロジックが簡単

メッシュのカバーする距離判定も簡単

図は["Bing Maps Tile System" Microsoft](https://msdn.microsoft.com/ja-jp/library/bb259689.aspx) より引用

Datastoreでも前方部分一致のクエリが投げられる ので、このやり方ならGAE派でも使えますね。

実践(golang)

では実際のコードを書いてみましょう。
golangではorbという素晴らしいライブラリがあるので、こちらを活用します。
もともとはgo.geoという、Stravaでも使われたライブラリの後継です。

生成編

main.go

package main

import (
	"github.com/paulmach/orb"
	"github.com/paulmach/orb/geojson"
	"github.com/paulmach/orb/maptile"
	"strconv"
)

// quadkeyString returns quadkey string from tile.
func quadkeyString(t maptile.Tile) string {
	// see the original logic; https://github.com/paulmach/go.geo/blob/master/point.go#L149
	s := strconv.FormatInt(int64(t.Quadkey()), 4)

	// for zero padding
	zeros := "000000000000000000000000000000"
	return zeros[:((int(t.Z)+1)-len(s))/2] + s
}

func main() {
	name := "皇居"
	address := "東京都千代田区千代田１−１"

	// Creating a point
	lat := 35.685323
	lng := 139.752768
	p := orb.Point{lng, lat}

	// Change to a maptile.
	t := maptile.At(p, 17)
	println(quadkeyString(t)) // 13300211231022032

	// viewing as geojson
	c := geojson.NewFeatureCollection()
	gp := geojson.NewFeature(p)
	gp.Properties["name"] = name
	gp.Properties["address"] = address

	gt := geojson.NewFeature(t.Bound())
	gt.Properties["quadkey"] = quadkeyString(t)

	c = c.Append(gp).Append(gt)
	b, _ := c.MarshalJSON()
	println(string(b)) //{"type":"FeatureCollection","features":[{"type":"Feature","geometry":{"type":"Point","coordinates":[139.752768,35.685323]},"properties":{"address":"東京都千代田区千代田１−１","name":"皇居"}},{"type":"Feature","geometry":{"type":"Polygon","coordinates":[[[139.7515869140625,35.684071533140965],[139.75433349609375,35.684071533140965],[139.75433349609375,35.68630240145626],[139.7515869140625,35.68630240145626],[139.7515869140625,35.684071533140965]]]},"properties":{"quadkey":"13300211231022032"}}]}
}

保存したい場所マスタをDBにInsertするときに、緯度経度とQuadkey(13300211231022032みたいなやつ)も一緒に保存するとよいです。
可能な限り大きなレベル(23)で保存しておくと自在に前方一致検索できるので、カリカリにチューニングしなくていい場合はおすすめです。

次回予告

風邪でしんどいので、今回はここまで。
来週あたりに、GAE/Goで実践してみた記事をUpしたいと思います。