この日記は、RDBMS-GIS(MySQL,PostgreSQLなど) Advent Calendar 2020 の16日目の記事です。
何をしたいか
先日の日記で雑に作成した大量のPOINTデータを使って、ある点から一定距離内にある点を検索します。普通にやると結構時間がかかりますが、空間インデックスがきちんと使用されて高速に検索できるところがポイントです。POINTの話だけに。
先日のデータ
先日の日記で書いたスクリプトを使って、テーブル sp1 に、約29万件のデータを作成しました。
mysql> select COUNT(*) FROM sp1; +----------+ | COUNT(*) | +----------+ | 286572 | +----------+
テーブル構造は、シンプルにこんな感じ。数値id と POINT型(SRID=4326)の sp の2カラムだけです。
mysql> SELECT * FROM sp1 LIMIT 3; +----+------------------------------------------------------+ | id | sp | +----+------------------------------------------------------+ | 1 | 0xE61000000101000000FAEDEBC039736140B98D06F016C04140 | | 2 | 0xE61000000101000000C1A8A44E40736140B98D06F016C04140 | | 3 | 0xE6100000010100000088635DDC46736140B98D06F016C04140 | +----+------------------------------------------------------+ 3 rows in set (0.00 sec) mysql> SELECT id, ST_AsText(sp) FROM sp1 LIMIT 3; +----+-------------------------+ | id | ST_AsText(sp) | +----+-------------------------+ | 1 | POINT(35.5007 139.6008) | | 2 | POINT(35.5007 139.6016) | | 3 | POINT(35.5007 139.6024) | +----+-------------------------+ 3 rows in set (0.01 sec)
とりあえず距離を表示してみる
ここからは、ある点P（POINT(35.685 139.762)）を起点にして、ここからの距離を見ていくことにします。
とりあえず、ただテーブル上の点と、それらの点Pからの距離だけを表示する例。全件表示しても意味がないので、LIMITかけています。
SELECT id, ST_AsText(sp) sp, ST_Distance_Sphere(sp, ST_GeomFromText('POINT(35.685 139.762)', 4326)) d FROM sp1 LIMIT 5; +----+-----------------------------------+--------------------+ | id | sp | d | +----+-----------------------------------+--------------------+ | 1 | POINT(35.5007 139.6008) | 25148.125196022018 | | 2 | POINT(35.5007 139.6016) | 25106.26817499381 | | 3 | POINT(35.5007 139.6024) | 25064.550018612503 | | 4 | POINT(35.5007 139.6032) | 25022.97142142526 | | 5 | POINT(35.5007 139.60399999999998) | 24981.533080267982 | +----+-----------------------------------+--------------------+ 5 rows in set (0.00 sec)
点Pから200m以内にある点を検索
まずは単純に、この ST_Distance_Sphere() を使って考えてみます。
SELECT id, ST_AsText(sp) p ,
ST_Distance_Sphere(sp, ST_GeomFromText('POINT(35.685 139.762)', 4326)) d
FROM sp1
WHERE ST_Distance_Sphere(sp, ST_GeomFromText('POINT(35.685 139.762)', 4326)) <200
ORDER BY d;
+--------+----------------------------------------------+--------------------+
| id | p | d |
+--------+----------------------------------------------+--------------------+
| 131966 | POINT(35.68480000000051 139.76239999999962) | 42.42302952190574 |
: : : :
| 131963 | POINT(35.68480000000051 139.75999999999962) | 181.99742870837645 |
| 132469 | POINT(35.68550000000051 139.7639999999996) | 188.99547708775145 |
| 132464 | POINT(35.68550000000051 139.75999999999962) | 188.99547715803516 |
+--------+----------------------------------------------+--------------------+
22 rows in set (2.59 sec) 22件が該当し、検索時間は 2.59秒。結構かかっています。
explainの結果は省略しますが、インデックスが使われていないようです。
ちなみに、参考情報ですが、ST_AsText() の結果ではなく、緯度と経度それぞれの数値を得たい場合は、ST_Latitude()、ST_Longitude() を使えば良いです。
SELECT id, ST_AsText(sp) p , ST_Latitude(sp) lat, ST_Longitude(sp) lng, ST_Distance_Sphere(sp, ST_GeomFromText('POINT(35.685 139.762)', 4326)) d FROM sp1 WHERE ST_Distance_Sphere(sp, ST_GeomFromText('POINT(35.685 139.762)', 4326)) <200 ORDER BY d; +--------+----------------------------------------------+--------------------+--------------------+--------------------+ | id | p | lat | lng | d | +--------+----------------------------------------------+--------------------+--------------------+--------------------+ | 131966 | POINT(35.68480000000051 139.76239999999962) | 35.68480000000051 | 139.76239999999962 | 42.42302952190574 | | 131965 | POINT(35.68480000000051 139.76159999999962) | 35.68480000000051 | 139.76159999999962 | 42.42302958257237 | | 132467 | POINT(35.68550000000051 139.76239999999962) | 35.68550000000051 | 139.76239999999962 | 66.303956142051 | ：
インデックスを使わせるために
ST_Within() を使うことを考えます。 点Pを中心とする半径200メートルのPOLYGONを生成し、それとテーブル上の各点との包含を調べるという作戦です。
ST_Buffer() は、ST_Buffer(中心点, 半径) の形で使用します。
ここで、
中心点： ST_GeomFromText('POINT(35.685 139.762)', 4326)
半径： 200 * 180.0 / PI() / 6378137.0
となります。中心点は、説明不要ですね。
半径は、これ、中心点を緯度経度で与えていますので、半径も度で与える必要があるんです。なので地球半径の値や円周率を使って半径200mに相当する度を算出しています。
実際にやってみると、、、
mysql> SELECT ST_AsText( -> ST_Buffer( -> ST_GeomFromText('POINT(35.685 139.762)', 4326), -> 200 * 180.0 / PI() / 6378137.0 -> )) a; ERROR 3618 (22S00): st_buffer(POINT, ...) has not been implemented for geographic spatial reference systems.
なんと。。ST_Bufferは、地理座標系に対応していないのだと。これは事件です！ 地球上のある一点を中心とする縁を、MySQLはまだ描けないということです。参った。
仕方がないので、誤差が出ることを承知で、強引に、デカルト座標（内部的には SRID=0）で生成させた円を、SRID 4326 のものだと見做して処理を続行してみます。
まず SRIDなし（内部的にはゼロ）の中心点から円を描かせます 。 SRID=0は axis order が逆になっていることに注意。
mysql> SELECT ST_AsText( -> ST_Buffer( -> ST_GeomFromText('POINT(139.762 35.685)'), -> 200 * 180.0 / PI() / 6378137.0 -> )) a\G *************************** 1. row *************************** a: POLYGON((139.76379663056824 35.685,139.76376210881565 35.6853505052361,139.7636598702095 35.68568754075255,139.76349384372202 35.68599815446345,139.76327040965808 35.68627040965809,139.76299815446345 35.686493843722026,139.76268754075255 35.68665987020948,139.7623505052361 35.68676210881566,139.762 35.68679663056824,139.7616494947639 35.68676210881566,139.76131245924745 35.68665987020948,139.76100184553655 35.686493843722026,139.76072959034192 35.68627040965809,139.76050615627798 35.68599815446345,139.7603401297905 35.68568754075255,139.76023789118435 35.6853505052361,139.76020336943176 35.685,139.76023789118435 35.6846494947639,139.7603401297905 35.68431245924746,139.76050615627798 35.68400184553656,139.76072959034192 35.68372959034191,139.76100184553655 35.68350615627798,139.76131245924745 35.68334012979052,139.7616494947639 35.68323789118435,139.762 35.683203369431766,139.7623505052361 35.68323789118435,139.76268754075255 35.68334012979052,139.76299815446345 35.68350615627798,139.76327040965808 35.68372959034191,139.76349384372202 35.68400184553656,139.7636598702095 35.68431245924746,139.76376210881565 35.6846494947639,139.76379663056824 35.685)) 1 row in set (0.01 sec)
できたデカルト上の円の数値をそのまま、SRID=4326 のものだと見做してみます。しつこく繰り返しますが、正確さの観点からは、全くナンセンスで意味のない読み替えをしています。
mysql> SELECT ST_AsText( -> ST_SRID( -> ST_Buffer( -> ST_GeomFromText('POINT(139.762 35.685)'), -> 200 * 180.0 / PI() / 6378137.0 -> ), 4326)) a\G *************************** 1. row *************************** a: POLYGON((35.685 139.76379663056824,35.6853505052361 139.76376210881565,35.68568754075255 139.7636598702095,35.68599815446345 139.76349384372202,35.68627040965809 139.76327040965808,35.686493843722026 139.76299815446345,35.68665987020948 139.76268754075255,35.68676210881566 139.7623505052361,35.68679663056824 139.762,35.68676210881566 139.7616494947639,35.68665987020948 139.76131245924745,35.686493843722026 139.76100184553655,35.68627040965809 139.76072959034192,35.68599815446345 139.76050615627798,35.68568754075255 139.7603401297905,35.6853505052361 139.76023789118435,35.685 139.76020336943176,35.6846494947639 139.76023789118435,35.68431245924746 139.7603401297905,35.68400184553656 139.76050615627798,35.68372959034191 139.76072959034192,35.68350615627798 139.76100184553655,35.68334012979052 139.76131245924745,35.68323789118435 139.7616494947639,35.683203369431766 139.762,35.68323789118435 139.7623505052361,35.68334012979052 139.76268754075255,35.68350615627798 139.76299815446345,35.68372959034191 139.76327040965808,35.68400184553656 139.76349384372202,35.68431245924746 139.7636598702095,35.6846494947639 139.76376210881565,35.685 139.76379663056824)) 1 row in set (0.00 sec)
一応、変換の結果 内部的にSRIDが 4326へと変更されたことを確認しておきましょうか。↑の ST_AsText のところを ST_SRID に変えただけです。
mysql> SELECT ST_SRID( -> ST_SRID( -> ST_Buffer( -> ST_GeomFromText('POINT(139.762 35.685)'), -> 200 * 180.0 / PI() / 6378137.0 -> ), 4326)) a; +------+ | a | +------+ | 4326 | +------+
SRID 4326 だと認識されていますね。
改めて、点Pから200m内の点を抽出する
ということでデカルト座標で生成した円ポリゴンを、強引に SRID 4326 の数値であるかのように読み替えたものに含まれる、テーブル上の点(sp)を抽出するSQLです。 ST_Within() を使っています。
SELECT id, ST_AsText(sp) p , ST_Distance_Sphere(ST_GeomFromText('POINT(35.685 139.762)', 4326), sp) d FROM sp1 WHERE ST_Within( sp, ST_SRID( ST_Buffer( ST_GeomFromText('POINT(139.762 35.685)'), 200 * 180.0 / PI() / 6378137.0 ), 4326) ) ORDER BY d; +--------+----------------------------------------------+--------------------+ | id | p | d | +--------+----------------------------------------------+--------------------+ | 131966 | POINT(35.68480000000051 139.76239999999962) | 42.42302952190574 | | 131965 | POINT(35.68480000000051 139.76159999999962) | 42.42302958257237 | : : | 130964 | POINT(35.6834000000005 139.76239999999962) | 181.54307375208234 | | 130963 | POINT(35.6834000000005 139.76159999999962) | 181.54307376625917 | +--------+----------------------------------------------+--------------------+ 18 rows in set (0.00 sec)
でも、あれっ？最初に ST_Distance_Sphere() で抽出したときには 確か22件でしたよね。今回18件というのはそれと相違します。
そう。何度も何度も強調した「無理矢理SRID 4326 だと思い込ませた円」というところがポイントなのです。今回のケース、平面座標上で作った円（の座標の値）をそのまま局面上の座標だと見做してしまうと、実は円が小さすぎるのです。そのため、本来検索にマッチすべき行が漏れてしまうと言うわけです。
決め技は「大きめに取ってきてから絞る」
もうゴールは近いです。円が小さいとわかってるなら、円を大きくすればいいじゃない。はい。強引なのです。
とりあえず 1.2 倍くらいしてみて、どーせ距離も算出しているのだから、後から欲しい距離(今回は200m)内のものだけを選べば良いわけです。
いきますよ。じゃん！
mysql> SELECT id, ST_AsText(sp) p , -> ST_Distance_Sphere(ST_GeomFromText('POINT(35.685 139.762)', 4326), sp) d -> FROM sp1 -> WHERE -> ST_Within( -> sp, -> ST_SRID( -> ST_Buffer( -> ST_GeomFromText('POINT(139.762 35.685)'), -> 200*1.2 * 180.0 / PI() / 6378137.0 -> ), 4326) -> ) -> HAVING d<200 -> ORDER BY d; +--------+----------------------------------------------+--------------------+ | id | p | d | +--------+----------------------------------------------+--------------------+ | 131966 | POINT(35.68480000000051 139.76239999999962) | 42.42302952190574 | : : | 131963 | POINT(35.68480000000051 139.75999999999962) | 181.99742870837645 | | 132469 | POINT(35.68550000000051 139.7639999999996) | 188.99547708775145 | | 132464 | POINT(35.68550000000051 139.75999999999962) | 188.99547715803516 | +--------+----------------------------------------------+--------------------+ 22 rows in set (0.01 sec)
22件が得られました。最後に出ている 181mくらいより外側のものが、先ほどは漏れてしまったわけです。このSQLのポイントは HAVING を使用しているということ。各行の距離 d を算出したあとで、その値に応じて絞るので、このようになります（GROUP BY 以外で HAVING を使ったことがない人も多いんじゃないかな）。
誤差はどれくらい
200mの範囲だとえいやっと 1.2倍すれば十分に欲しい円が含まれるようになりました。半径が大きくなるほどに、強引に「見做した」差は広がっていきますので、たとえば 20000m の円にすると、1.5倍くらいしないと欲しい円が含まれるようになりません。 使用したい半径がある程度小さいものを対象にする場合に、今回の手法が使えると考えた方が良いでしょう。
時間計測
最初にやった、ST_Distance_Sphere() での雑な結果。 先頭が絞った距離(メートル)、後半が mysql での実行結果行（マッチ行数と処理時間）です。
<200: 22 rows in set (2.59 sec) <2000: 2236 rows in set (2.74 sec) <4000: 8936 rows in set (3.10 sec) <8000: 35750 rows in set (4.55 sec) <16000: 140796 rows in set (10.22 sec) <20000: 205393 rows in set (13.75 sec)
今回の最終結果のSQLだと、このようになります。括弧内は、1.2倍の代わりに何倍にしたかを示しています（記述がないものは 1.2倍）。
<200 : 22 rows in set (0.01 sec) <2000 : 2212 rows in set (0.25 sec) <4000 : 8846 rows in set (0.98 sec) <4000(*1.3) : 8936 rows in set (1.11 sec) <8000 :35368 rows in set (3.89 sec) <8000(*1.4) :35750 rows in set (5.02 sec) <16000 :140043 rows in set (15.51 sec) <16000(*1.5):140796 rows in set (19.13 sec) <20000(*1.5):205393 rows in set (22.03 sec)
途中で、ST_Distance_Sphere() を使用した場合との処理時間の逆転が発生しています。想像ですが、範囲が広いというよりも、結果セットとして該当するものが多くなることにより、なんらかの処理に要するコストのほうが勝ってしまうことが起きているのかな、と雑に考えています。もうちょっと調べてみたいところですね（誰かに期待）。
今回の実験のきっかけ
「最寄りバス停を地図や住所から無料で探せる - バス停検索」さんに検索の遅さ相談をいただいて、色々調べたものでした。バス停検索 buste.in さんは、先日開催された db tech showcase ONLINE 2020 でも、MySQLのGIS機能活用事例として発表をされました。
資料は、Oracleさんの以下の資料ページからダウンロードすることができます（要ログイン）。たぶん未来にこの日記をご覧になる方は、その後、どんどん資料が追加されていて今回の資料が潜ってしまっていると思うので「db tech showcase ONLINE 2020」をキーワードに探してみてください。
MySQL :: バス停検索サービスで周辺バス停を取得する重いSQLの改善
あー、面白かった。
やっぱり、実際のサービスに役に立つSQLを考えるのは楽しいですね。
全国の MySQL 8.0 で位置情報を扱っていて／扱おうとしていてちょっと困っちゃっている方。なんなら相談に乗りますよ。 MySQL 5.7 で頑張っている方は、とりあえず話は 8.0 に上げてからだ！(笑)
- - - -追記2020/12/17- - - -
発表された福田さん自身の slideshare でも、資料が公開されました（ログイン不要）。
- - - - 追記ここまで - - - -
