MySQLの空間データ型の変換(2)~POINTの集合からLINESTRINGを作る~

この日記は、RDBMS-GIS(MySQL,PostgreSQLなど) Advent Calendar 2020 の24日目の記事です。

はじめに

 先日の日記で、LINESTRING や MULTIPOINT にある点の要素を、POINTのデータ群として得る方法のアイデアを紹介しました。
MySQLの空間データ型の変換(1)~MULTIPOINTやLINESTRINGからPOINTを得る~

 今回は、その続編として、POINT型のデータ群をつなげて LINESTRING にする方法のアイデアについて書きたいと思います。正直なところ「ちからワザ」です。とりあえずこのようなやり方で実現は可能だぞという、ひとつの思考実験的なものとしてお読みいただければと思います。
f:id:sakaik:20201218223043p:plain

実験データの用意

GEOMETRY型のカラムを持つ、以下のテーブルを作成し、データを投入します。とりあえず3件(3点)

CREATE TABLE t1 (id integer, p geometry, odr integer, cat varchar(1))
INSERT INTO t1 VALUES (1,  ST_GeomFromText('POINT(1 3)'), 10, "A");
INSERT INTO t1 VALUES (2,  ST_GeomFromText('POINT(2 2)'), 20, "A");
INSERT INTO t1 VALUES (3,  ST_GeomFromText('POINT(3 4)'), 30, "A");

この3つの点を、odr列の値の順につなげて LINESTRING 型の結果を得ることを目標とします。

思考実験

 処理のイメージとしては、

  • 必要なPOINT型の行を GROUP BY でひとかたまりとして扱って
  • SUM() 関数みたいな感じで、例えば GROUP_LINESTRING()みたいな集約関数を使って LINESTRING にする(その際ORDERも指定する)

なんてことができたらいいなと考えるわけですが、残念ながらそんなグルーピング関数はありません。

 仕方がないので、作戦変更し、POINT型の集合を得た後、自力で LINESTRING にすることを考えてみます。
LINESTRING にするためには、ST_GeomFromText('LINESTRING(1 1, 2 2, 3 3, 4 4)') のように、WKT文字列を構成してから ST_GeomFromText() で変換すれば良いわけです。

3点をLINESTRINGにしてみる

 ということで、先ほどの3つのPOINTをLINESTRINGに変換してみたのが以下のSQL
GROUP_CONCATを使って、抽出したPOINTのX,Yの集合をコンマ区切りでつなげます。GROUP_CONCAT の ORDER BY を使うことで、点の順序を指定できます。

SELECT ST_AsText(
    ST_GeomFromText(
    CONCAT("LINESTRING(", GROUP_CONCAT(ST_X(p), " ", ST_Y(p) ORDER BY odr), ")"))) a
 FROM t1
 WHERE cat="A";
+-------------------------+
| a                       |
+-------------------------+
| LINESTRING(1 3,2 2,3 4) |
+-------------------------+
1 row in set (0.00 sec)

 なんだかそれっぽいLINESTRINGが得られました。これ、単に文字列結合した結果が表示されているわけではなくて、一旦Geomに変換後、改めて表示用に ST_AsText() してるんですからね。
念のため、ORDER BY がちゃんと効いていることを確認するために、ORDER BY の DESC の動作も見ておきましょう。

mysql> SELECT ST_AsText(
    ->     ST_GeomFromText(
    ->     CONCAT("LINESTRING(", GROUP_CONCAT(ST_X(p), " ", ST_Y(p) ORDER BY odr DESC), ")"))) a
    ->  FROM t1
    ->  WHERE cat="A";
+-------------------------+
| a                       |
+-------------------------+
| LINESTRING(3 4,2 2,1 3) |
+-------------------------+

 期待通りに ORDER BY が働いているようです。

このクエリの制限と回避方法

 とりあえず力業で、やりたいことを実現させましたが、実はこのやり方は万能ではありません。
どこに落とし穴があるのか。それは GROUP_CONCAT() の文字数制限です。
デフォルトでは、GROUP_CONCATは 1024文字までを扱うことができます。この値は変更することができるので、実際に、先ほどのSQLのGROUP_CONCAT() が出力することになりそうな最大のサイズを見積もった上で、group_concat_max_len を設定すると良いでしょう。

mysql> show variables like 'group_con%';
+----------------------+-------+
| Variable_name        | Value |
+----------------------+-------+
| group_concat_max_len | 1024  |
+----------------------+-------+

もう少し点数の多い例でも確認

 3点だけではやや物足りないので、もう少し点数の多いPOINTを、LINESTRINGにしてみましょう。以下のSQLでデータを追加します。

INSERT INTO t1 VALUES (10,  ST_GeomFromText('POINT(1 2)'), 21, "T");
INSERT INTO t1 VALUES (11,  ST_GeomFromText('POINT(1 3)'), 25, "T");
INSERT INTO t1 VALUES (12,  ST_GeomFromText('POINT(1 4)'), 31, "T");
INSERT INTO t1 VALUES (13,  ST_GeomFromText('POINT(3 0)'), 12, "T");
INSERT INTO t1 VALUES (14,  ST_GeomFromText('POINT(3 2)'), 16, "T");
INSERT INTO t1 VALUES (15,  ST_GeomFromText('POINT(3 3)'), 24, "T");
INSERT INTO t1 VALUES (16,  ST_GeomFromText('POINT(3 4)'), 29, "T");
INSERT INTO t1 VALUES (17,  ST_GeomFromText('POINT(4 6)'), 33, "T");
INSERT INTO t1 VALUES (18,  ST_GeomFromText('POINT(5 2)'), 55, "T");
INSERT INTO t1 VALUES (19,  ST_GeomFromText('POINT(5 3)'), 50, "T");
INSERT INTO t1 VALUES (20, ST_GeomFromText('POINT(5 4)'), 36, "T");
INSERT INTO t1 VALUES (21, ST_GeomFromText('POINT(5 0)'), 60, "T");
INSERT INTO t1 VALUES (22, ST_GeomFromText('POINT(7 2)'), 53, "T");
INSERT INTO t1 VALUES (23, ST_GeomFromText('POINT(7 3)'), 39, "T");
INSERT INTO t1 VALUES (24, ST_GeomFromText('POINT(7 4)'), 34, "T");

 先ほどうまくいったクエリを少し改良して、、、、、

mysql> SELECT ST_AsText(
    ->     ST_GeomFromText(
    ->     CONCAT("LINESTRING(", GROUP_CONCAT(ST_X(p), " ", ST_Y(p) ORDER BY odr), ")"))) a
    ->  FROM t1
    ->  WHERE cat="T";
+-------------------------------------------------------------------------+
| a                                                                       |
+-------------------------------------------------------------------------+
| LINESTRING(3 0,3 2,1 2,3 3,1 3,3 4,1 4,4 6,7 4,5 4,7 3,5 3,7 2,5 2,5 0) |
+-------------------------------------------------------------------------+

 LINESTRING が得られました。目的達成です! でも文字だけの出力結果を見ても何だかよく分からないので、やっぱりビジュアルに表示してみたいですね。MySQL Workbench を使って見てみましょう。
その際、クエリ結果を人間が見るわけじゃないので、ST_AsText()は外します。

SELECT ST_GeomFromText(
    CONCAT("LINESTRING(", GROUP_CONCAT(ST_X(p), " ", ST_Y(p) ORDER BY odr), ")")) a
 FROM t1
WHERE cat="T";


f:id:sakaik:20201224002629p:plain

 メリークリスマス! Merry Xmas!!

カジュアルにmysqlを落とす(Windows)

 先日の日記を書くために色々試している最中、あまりにも潔く mysql クライアントコマンドがさくっと落ちるので、発生条件を調べてみました。原因の特定には至っていません。

起こっていること

 Windows 10 のMySQLにて、mysql> プロンプトが出ているところで、単純なSELECT文を叩くと、無言で mysql があっさり落ちる。

mysql> SELECT * FROM n03_20_200101 WHERE RECID<20;

D:\>

使用データ

 先日の日記に書いた方法で取り込んだ「全国」のデータ(千葉県ではなく全国)の入ったテーブルを使って再現させることができます。ふつうに blob に放り込んだ値を HEX() 取ったりしても再現できるかもしれません。あ、blob である必要すらなくて、varchar とか text とかでも良さそう。

sakaik.hateblo.jp

発生条件

 どうやら、結果セットに含まれるカラム値の サイズまたは行のサイズが大きすぎるときに落ちる気がする。

mysql> SELECT RECID, LENGTH(HEX(SHAPE)),LENGTH(ST_AsText(SHAPE)) FROM n03_20_200101 WHERE RECID<20;
+-------+--------------------+--------------------------+
| RECID | LENGTH(HEX(SHAPE)) | LENGTH(ST_AsText(SHAPE)) |
+-------+--------------------+--------------------------+
|     0 |              57826 |                    55812 |
|     1 |              57634 |                    55409 |
|     2 |              27394 |                    26233 |
|     3 |              25666 |                    24674 |
|     4 |              47362 |                    45461 |
|     5 |             162210 |                   156154 |
|     6 |              52898 |                    50970 |
|     7 |              22178 |                    21460 |
|     8 |              35874 |                    34414 |
|     9 |              56194 |                    54322 |
|    10 |               1154 |                     1079 |
|    11 |                386 |                      344 |
|    12 |                642 |                      582 |
(略)

という情報を元に、RECID=5 の行のみを抽出してみる。

mysql> SELECT RECID, LENGTH(SHAPE), ST_AsText(SHAPE) FROM n03_20_200101 WHERE RECID=5;

D:\>

うん、落ちた。AsTextとかよくわかんないという人は、代わりに HEX(SHAPE)でも良いです。

閾値はどこに

 試行錯誤で範囲を絞ったので、以下のクエリで閾値捜索の最後のツメを。

SELECT RECID, LENGTH(SHAPE) a, LENGTH(ST_AsText(SHAPE)) b FROM n03_20_200101 
  WHERE LENGTH(ST_AsText(SHAPE)) >101296 AND LENGTH(ST_AsText(SHAPE))<101719 ORDER BY b ;
+--------+-------+--------+
| RECID  | a     | b      |
+--------+-------+--------+
|  99545 | 54305 | 101297 |
|  31451 | 53921 | 101297 |
|  59657 | 54225 | 101307 |OK
|  33266 | 53905 | 101361 |OK
|  48032 | 54049 | 101376 |NG
|  32203 | 54245 | 101437 |
|  11066 | 53301 | 101528 |
| 115488 | 49233 | 101718 |
+--------+-------+--------+
8 rows in set (28.52 sec)
SELECT RECID, LENGTH(SHAPE), ST_AsText(SHAPE) FROM n03_20_200101 WHERE RECID=33266; --OK
SELECT RECID, LENGTH(SHAPE), ST_AsText(SHAPE) FROM n03_20_200101 WHERE RECID=48032; --NG

1文字ずつくわえて見ると RECID=33266 は1文字でも加えたら、もうアウト。

SELECT RECID, LENGTH(SHAPE), CONCAT(ST_AsText(SHAPE),"X") FROM n03_20_200101 WHERE RECID=59657; --OK
SELECT RECID, LENGTH(SHAPE), CONCAT(ST_AsText(SHAPE),"X") FROM n03_20_200101 WHERE RECID=33266; --NG

 もしかして、出力カラムのサイズではなくレコードサイズなのでは?と 余計な取得列を少しずつ削って見ると、アタリ。他の列をすべて取得やめてみた結果が以下のもの。

SELECT CONCAT(ST_AsText(SHAPE),"XXXXX") FROM n03_20_200101 WHERE RECID=33266;  //OK
SELECT CONCAT(ST_AsText(SHAPE),"XXXXXX") FROM n03_20_200101 WHERE RECID=33266; //NG

 101361+5= 101366 byte まではOK、 101367バイトがNGと推測できそうです。

環境

 Windows 10 上に、MySQL Installer を使ってインストールしたもの。サーバは utf8mb4、
クライアントは cp932 でも utfmb4 でも発生。

mysql> status
--------------
mysql  Ver 8.0.22 for Win64 on x86_64 (MySQL Community Server - GPL)

Connection id:          30
Current database:       shptest
Current user:           root@localhost
SSL:                    Cipher in use is TLS_AES_256_GCM_SHA384
Using delimiter:        ;
Server version:         8.0.22 MySQL Community Server - GPL
Protocol version:       10
Connection:             localhost via TCP/IP
Server characterset:    utf8mb4
Db     characterset:    utf8mb4
Client characterset:    cp932
Conn.  characterset:    cp932
TCP port:               3306
Binary data as:         Hexadecimal
Uptime:                 3 hours 35 min 55 sec

Threads: 4  Questions: 117  Slow queries: 8  Opens: 190  Flush tables: 3  Open tables: 111  Queries per second avg: 0.009
mysql> status
--------------
C:\Program Files\MySQL\MySQL Server 8.0\bin\mysql.exe  Ver 8.0.22 for Win64 on x86_64 (MySQL Community Server - GPL)

Connection id:          29
Current database:       shptest
Current user:           root@localhost
SSL:                    Cipher in use is TLS_AES_256_GCM_SHA384
Using delimiter:        ;
Server version:         8.0.22 MySQL Community Server - GPL
Protocol version:       10
Connection:             localhost via TCP/IP
Server characterset:    utf8mb4
Db     characterset:    utf8mb4
Client characterset:    utf8mb4
Conn.  characterset:    utf8mb4
TCP port:               3306
Binary data as:         Hexadecimal
Uptime:                 3 hours 35 min 21 sec

Threads: 3  Questions: 113  Slow queries: 8  Opens: 190  Flush tables: 3  Open tables: 111  Queries per second avg: 0.008

結局分からず仕舞い

 現時点では Linux での動作を試していないので、これが Windows固有の現象なのか OS問わず発生するものなのか、分かりません。また、サーバのエラーログには記録なしで、一切の調査の足がかりがないのが、なかなか辛いですね。 まぁ、許容サイズを超えてエラーになるのはまだしも、ストンと落ちてしまう(接続が切れてしまう)のはいただけないなぁというのが私の感覚です。

 mysql params でざっと見てみた感じでは、それっぽい制限値はなさそうで、もにゃもにゃした気分です。
variable / MySQL Parameters


 これ、なんなんですかね。

その晩の追記

 その後いろいろ試してみたら、とりあえず上で書いた「閾値」候補は、一旦とりさげ。
例えばこれ、

SELECT RECID, LENGTH(SHAPE), ST_AsText(SHAPE) FROM n03_20_200101 WHERE RECID=33266; --OK
SELECT RECID, LENGTH(SHAPE), ST_AsText(SHAPE) FROM n03_20_200101 WHERE RECID=48032; --NG

実行順序を逆にしたら、NGだったやつは通って、OKだったやつがダメになりました。もう少し再現性ある内容を書く前に、yoku0825さんが試してくれたものを紹介すると、これだけで落ちるらしい。

mysql> SELECT repeat('a', 101367);
(結果が表示される)

mysql> SELECT 'a';

D:\>

 「次のクエリをダメにするクエリ」みたいなのがあるっぽいと考えたほうがいいのかも?


もちょっと追記。「次のクエリをダメにするクエリ」という観点で試してみたら、こんな落とし方もあります、

mysql> SELECT repeat('a', 100000);
(結果が表示される)

mysql> STATUS

D:\>

shapefileをMySQLに取り込む!shp2sqlの紹介

この日記は、RDBMS-GIS(MySQL,PostgreSQLなど) Advent Calendar 2020 の7日目ぶん(あいていたので後から埋めてます)の記事です。

この日記は

 世の中の地理情報データ(位置の情報を含んだデータ)には、shapefile というファイル形式で公開されているものが非常にたくさんあります。この日記では、shapefileのデータをMySQL に取り込む方法として、私が開発している shp2sql というツールを紹介します。

shp2sql の入手

 shp2sql は単一の Python スクリプトです。github にて公開しています。以下のURLから shp2sql.py を取得してください。pythonの動作環境が必要です。
github.com

サンプルデータの入手

 ここでは、国土交通省が公開している「国土数値情報」を利用します。日本の国土に関する様々なデータが公開されていますが、ここでは以下のURLから、「2.政策区域」-「行政地域」-「行政区域(ポリゴン)」と辿ったページより、千葉・令和2年(N03-20200101_12_GML.zip)を取得して説明を続けます。ダウンロードしたzipファイルを展開しておきます。

https://nlftp.mlit.go.jp/ksj/index.html

MySQL側の用意

 MySQL側に、shapefileを取り込むためのデータベース(スキーマ)を作成しておきます。既存のスキーマにデータを取り込みたい場合は、作成不要です。 また、取り込み先のテーブルは作成しておく必要はありません。

CREATE DATABASE shptest;

shp2sql の実行

 展開した(ここでは千葉県の)shapefile があるフォルダに移動します。話をシンプルにするために ここでは、shp2sql.py も同じフォルダに存在していることにします。
 ここで使うshapefileは、拡張子が shp のものの他に、dbf, shx の計3つが必要になります。


 引数なしで shp2sql を実行すると、usage が表示されます。このへんのオプションは、まだまだ今後変更される可能性があります。-h を付けるともう少し詳細に表示されます。

D:\work\> shp2sql.py
usage: shp2sql.py [-h] [-c CHARSET] [-s SRID] [-d] [-x] [-T TABLENAME] input_filename
shp2sql.py: error: the following arguments are required: input_filename


 shapefileMySQLに取り込むには、基本的には、取り込みたい shapefile 名を指定して実行すれば良いのですが、文字コードutf-8 以外の場合は、その文字コードを教えてあげる必要があります。上記ダウンロードした ファイル群は cp932 ですので、-c オプションで与えてあげます。shp2sql は、shapefile を読み込んでSQL文を吐いてくれるだけのツールなので、その結果を mysqlクライアントに渡します。
 先ほどダウンロードした千葉県データの場合(N03-20_12_200101.shp)を取り込む実行例です。

D:\work\>shp2sql.py N03-20_12_200101.shx -ccp932 | mysql -uroot -p shptest


これだけで、shptest データベース内に N03-20_12_200101 テーブルが作られ、メタデータ並びに位置情報データが登録されます。

MySQL Workbench で確認

MySQL Workbench で、全件を見てみるとこんな感じ。SHAPE列に「BLOB」型のデータが入っていることがわかります。

f:id:sakaik:20201220223458p:plain


ひとつの行を選択して、右側から「Form Editor」を選ぶと、こんなふうに1行のデータが表示され、GEOMETRY系のカラムについては、その形も表示されます。
f:id:sakaik:20201220223705p:plain

右側から「Spatial View」で、千葉県全域のデータが入っていることも確認できます。
f:id:sakaik:20201220223410p:plain

もちろん、SELECT文にWHERE句を付けて、自分の見たいデータだけを抽出して確認もできます。みなさんは基本的な教養として我孫子市のコードが 12222 ということは知っていると思うので、それで検索してみた結果です。
f:id:sakaik:20201220223944p:plain
 ほら、我孫子市の形をしているでしょう。 (大丈夫、地元民でも、利根川手賀沼が描かれていないと「このへんこんな形をしていたかなぁ」とか迷ったりしますのでw)

制限事項とか今後の予定とか

 shp2sql.py は、現在開発途上です。まずは動作するものをと目指している段階です。以下の制限事項ならびに、改善を予定しています。

  • 一部のデータ型には未対応:POLYGON/LINESTRING/POINT のような基本的なデータ型には対応していますが、まだ全てのデータ型には対応していません。未対応のデータフォーマットファイルを見つけたら、ぜひお知らせください。調査します。
  • 遅いです:そんなに大きなデータも扱わないので、、と考えていたのですが、(千葉県だけでなく)全国のデータを一括で取り込もうとしていつまでも終わらず苦労しました。相当遅いです。single transaction 化や、INS/UPD ではなく一発のINSERTで実行するようにする等、改善予定です。

おわりに

 一度動く環境を作ってしまえば気楽に shapefileMySQL に取り込むことができる shp2sql.py。ぜひお手元に環境を作成して、試してみてください。取り込みができなかった(エラーとなった)ファイルがありましたら、そのファイルの入手方法とともにお知らせください。特に未対応ファイルトーマットについては優先順位を少しあげて対応を考えたいと思います。(全体としての優先順位は、「まず、ありとあらゆるshapefileMySQLに取り込めるようにすること」それから「速度改善」「便利機能」「きめ細かなエラー処理」です)


2020/12/21追記

 滅茶苦茶処理が遅かった話について。 上記全国のファイルのデータ登録を試みたら、半日以上経っても終わらなかったのですが、原因判明&手元では改善済み。
 解析の都合で、現在のツールでは、INSERT & UPDATE の処理を行っています。この 大量のUPDATEをする際の対象テーブルに、インデックス(プライマリキー)が張られていないというオチでした。PK張ったら、数秒で終わるようになりました。なんだ、色々改良点を考えたけど、その辺りは今のままでもいいのかもしれない(笑)。

MySQLでデフォルトデータベースをナシにする方法

タイトルを見てピンと来ない人も多いと思いますが、MySQLで、「use DB名」してデフォルトデータベースを一度指定すると、もう、何も指定していない状態には戻れないようですというお話です(ようです、というのは、私が方法を知らないだけかもしれないということ)。

MySQLに、デフォルトデータベース(以下デフォルトDB)を指定せずに接続した直後は、

mysql> SELECT DATABASE();
+------------+
| DATABASE() |
+------------+
| NULL       |
+------------+

 となっています。デフォルトDBを指定していないのだから当然ですよね。
デフォルトDBを指定(以下の例では test データベース)すると、DATABASE()の結果でも変更を確認することができます。

mysql> use test
Reading table information for completion of table and column names
You can turn off this feature to get a quicker startup with -A

Database changed
mysql> SELECT DATABASE();
+------------+
| DATABASE() |
+------------+
| test       |
+------------+

 さて、一度指定したこの デフォルトDB。なかったことにする(イメージとしては「スキーマの外に出る」という感じ)ことはできないのでしょうか。

unuse
とか
use null
とかのようなことができればよいのですが、残念ながらそういうコマンドはなさそうです(ご存知の方がいたら教えてください!)。


ということで、編み出した裏技(というほどでもない)。

以下のようにすると、デフォルトDBが なし になります。

CREATE DATABASE z;
use z;
DROP DATABASE z; 

 適当な名前の(既存ではない)データベースを作成し、それをデフォルトDBとしてから、そのデータベースをDROPする、ということをやっています。以下のように1行にまとめて実行することもできます。

mysql> CREATE DATABASE z; use z; DROP DATABASE z;
Query OK, 1 row affected (0.01 sec)
Database changed
Query OK, 0 rows affected (0.01 sec)
mysql> SELECT DATABASE();
+------------+
| DATABASE() |
+------------+
| NULL       |
+------------+


 なお、この方法は「やってみたらできた」ということを書いたものですので、本手順による悪影響等についての検証はしていません。大切な環境、大切な場面でお使いになる場合は、ご自身で十分に調査、検証をなさってください。


 まぁ、デフォルトDBの指定を なし にしたい、なんてことは、たぶんないとは思うんですけどね。
ちょっぴり ネタ気味の話題でした。

MySQLの空間データ型の変換(1)~MULTIPOINTやLINESTRINGからPOINTを得る~

この日記は、RDBMS-GIS(MySQL,PostgreSQLなど) Advent Calendar 2020 の 17日目のエントリーです。

はじめに

 MySQLで取り扱うことができる空間情報の型には、点、線、ポリゴン があります。それぞれ にそれらの集合を扱える型が存在してるので、都合6つとなります。これら相互の変換について考えてみたいと思います。

f:id:sakaik:20201217232839p:plain

なお、緑色線は本日記および今後の日記にて紹介を予定しているものです(スマートなやり方ではなく力ワザ(ちからわざ)のものも含む)。赤い点線は、今のところ試したことがないもの。

今回使うデータ

 今回のお試しの範囲では、テーブルにデータを入れておく必要もないのですが、実際のSQL利用イメージに近づけるためにこのようにしています。以下のテーブルおよびデータを登録しておきます。

mysql> CREATE TABLE s1 (id integer, sp GEOMETRY);
mysql> INSERT INTO s1 VALUES (1, ST_GeomFromText('MULTIPOINT(1 1, 2 2, 3 2, 4 5, 6 6)'));
mysql> INSERT INTO s1 VALUES (2, ST_GeomFromText('LINESTRING(1 0, 2 3, 4 8, 3 3, 2 7)'));

 こんな感じで格納されています。

mysql> SELECT id, ST_AsText(sp) FROM s1;
+------+-------------------------------------------+
| id   | ST_AsText(sp)                             |
+------+-------------------------------------------+
|    1 | MULTIPOINT((1 1),(2 2),(3 2),(4 5),(6 6)) |
|    2 | LINESTRING(1 0,2 3,4 8,3 3,2 7)           |
+------+-------------------------------------------+

今回の目標

  • 1つの MULTIPOINT として格納されているものを、それぞれの点にバラして POINT として得ます。
  • 1つの LINESTRING として格納されているものを、点にバラして POINT として得ます。

LINESTRING から点を抽出(点に分解)

 まず先に、専用の関数があってわかりやすい LINESTRING から考えてみます。
線を構成するポイントの数を得る ST_NumPoints() 、そして、N番目のポイントを得る ST_PointN() の動作を確認してみましょう。

mysql> SELECT id, ST_AsText(sp), ST_NumPoints(sp) n, ST_AsText(ST_PointN(sp, 2)) p
    ->   FROM s1
    ->  WHERE ID=2;
+------+---------------------------------+------+------------+
| id   | ST_AsText(sp)                   | n    | p          |
+------+---------------------------------+------+------------+
|    2 | LINESTRING(1 0,2 3,4 8,3 3,2 7) |    5 | POINT(2 3) |
+------+---------------------------------+------+------------+
1 row in set (0.00 sec)

 5個のポイントで構成され、2番目は POINT(2 3) であるという結果が得られました。
次に、このひとつの LINESTRING データをPOINTにバラしたい、つまり5件のデータにしたいということなので、このまま5件に分解してみます。CTEがあって良かった! MySQL 8.0 バンザイ!

WITH RECURSIVE seq(n) AS (SELECT 1 UNION SELECT n+1 FROM seq WHERE n<(SELECT ST_NumPoints(sp) FROM s1 WHERE id=2))
SELECT id, seq.n, ST_AsText(sp) FROM s1 JOIN seq ON (s1.id=2);

+------+------+---------------------------------+
| id   | n    | ST_AsText(sp)                   |
+------+------+---------------------------------+
|    2 |    1 | LINESTRING(1 0,2 3,4 8,3 3,2 7) |
|    2 |    2 | LINESTRING(1 0,2 3,4 8,3 3,2 7) |
|    2 |    3 | LINESTRING(1 0,2 3,4 8,3 3,2 7) |
|    2 |    4 | LINESTRING(1 0,2 3,4 8,3 3,2 7) |
|    2 |    5 | LINESTRING(1 0,2 3,4 8,3 3,2 7) |
+------+------+---------------------------------+
5 rows in set (0.00 sec)

 いい感じです。
ついでに、何番目の要素を取得したいかの数字も一緒に出力するようにしてみました。
この数字を使えば、最初に確認した ST_PointN() を使って、ほしいPOINTを得られますね。

WITH RECURSIVE seq(n) AS (SELECT 1 UNION SELECT n+1 FROM seq WHERE n<(SELECT ST_NumPoints(sp) FROM s1 WHERE id=2))
SELECT id, seq.n, ST_AsText(ST_PointN(sp, seq.n)) p FROM s1 JOIN seq ON (s1.id=2)

+------+------+------------+
| id   | n    | p          |
+------+------+------------+
|    2 |    1 | POINT(1 0) |
|    2 |    2 | POINT(2 3) |
|    2 |    3 | POINT(4 8) |
|    2 |    4 | POINT(3 3) |
|    2 |    5 | POINT(2 7) |
+------+------+------------+
5 rows in set (0.00 sec)

 はい。これで 1本の LINESTRING を構成する要素を POINT として得ることができました。いわゆるタテモチ・ヨコモチ変換の空間情報版みたいですね。

MULTIPOINT を POINT に

 こちら、LINESTRINGのような関数が存在していないので困りました。が、GEOMETRY として考えることで、LINESTRINGと同様のことができると分かったので、紹介します。

 対象データは、冒頭で登録した MULTIPOINT のデータです。

mysql> SELECT id, ST_AsText(sp) sp FROM s1 WHERE id=1;
+------+-------------------------------------------+
| id   | sp                                        |
+------+-------------------------------------------+
|    1 | MULTIPOINT((1 1),(2 2),(3 2),(4 5),(6 6)) |
+------+-------------------------------------------+

 ここでは、ST_NumGeometries() と ST_GeometryN() を使用します。

mysql> WITH RECURSIVE seq(n) AS (SELECT 1 UNION SELECT n+1 FROM seq WHERE n<(SELECT ST_NumGeometries(sp) FROM s1 WHERE id=1))
    -> SELECT id, seq.n, ST_AsText(ST_GeometryN(sp, seq.n)) p FROM s1 JOIN seq ON (s1.id=1)
    -> ;
+------+------+------------+
| id   | n    | p          |
+------+------+------------+
|    1 |    1 | POINT(1 1) |
|    1 |    2 | POINT(2 2) |
|    1 |    3 | POINT(3 2) |
|    1 |    4 | POINT(4 5) |
|    1 |    5 | POINT(6 6) |
+------+------+------------+

 はい。欲しい結果が得られました。
なお、先ほどの LINESTRING に対して ST_NumGeometries()を使おうとすると、うまくいきません。
ST_NumGeometries() は、 LINESTRING に対しては NULLを返すようです。

mysql> SELECT id,ST_NumGeometries(sp) FROM s1;
+------+----------------------+
| id   | ST_NumGeometries(sp) |
+------+----------------------+
|    1 |                    5 |
|    2 |                 NULL |
+------+----------------------+

おわりに

 とりあえず LINESTRING 型および MULTIPOINT 型から、それらを構成する POINT のデータ群へと変換することが実現できました。
パフォーマンス的なことは検証していませんし、抽出条件に相当するものが2か所に出ている点はもう少しきれいに書けるのではないかと思いますが、まずは最悪の方法だとしても、方法がないわけではないことを示せたことは一歩前進かなと思います。
 これを読んで、「なんてひどいSQLだ」と感じた方は(というか感じてほしいです)、ぜひスマートな書き方に挑戦していただき、ブログ等で紹介してほしいです。

 これ以外の変換については、来週のエントリで紹介予定です。カワザです(笑)。

MySQLのGIS機能で円の範囲に含まれる点を検索する試行

この日記は、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 でも、MySQLGIS機能活用事例として発表をされました。
 資料は、Oracleさんの以下の資料ページからダウンロードすることができます(要ログイン)。たぶん未来にこの日記をご覧になる方は、その後、どんどん資料が追加されていて今回の資料が潜ってしまっていると思うので「db tech showcase ONLINE 2020」をキーワードに探してみてください。

MySQL :: 資料ダウンロード

MySQL :: バス停検索サービスで周辺バス停を取得する重いSQLの改善



あー、面白かった。
やっぱり、実際のサービスに役に立つSQLを考えるのは楽しいですね。
全国の MySQL 8.0 で位置情報を扱っていて/扱おうとしていてちょっと困っちゃっている方。なんなら相談に乗りますよ。 MySQL 5.7 で頑張っている方は、とりあえず話は 8.0 に上げてからだ!(笑)


https://buste.in/buste.in

  • - - -追記2020/12/17- - - -

発表された福田さん自身の slideshare でも、資料が公開されました(ログイン不要)。

www2.slideshare.net

  • - - - 追記ここまで - - - -

 

MySQLに超・雑に大量のPOINTデータ(位置データ)を作成する方法

これは、RDBMS-GIS(MySQL,PostgreSQLなど) Advent Calendar 2020 の11日目のエントリだったはずのものです。実際には3日ほど遅れてしまいました。

 MySQLで、空間情報を扱う動作の確認をしたいときに、ちょっとばかり多めのデータが欲しくなることがあります。
実際の世の中のデータを取り込む形で実現しても良いのですが、ある一定の範囲内に存在する点の分量を自由に増減できない点は不便です。

 ということで、とってもザツに大量データを作成する機会があったので、本日記ではその方法を披露。

とりあえずテーブルを作る

 とりあえず ID めいたものと、POINT型のカラムを持つテーブルを作成します。
今回は、点(POINT)しか入れるつもりがないので、GEOMETRY型ではなく POINT 型にしてみました。
この後、空間インデックスを張るので、SRID の指定が必須です。とりあえずよく見る WGS84 にて。

CREATE TABLE sp1 (
  id INTEGER NOT NULL AUTO_INCREMENT PRIMARY KEY,
  sp POINT SRID 4326 NOT NULL
);

空間インデックスを張る

 用途に依るのだけど、とりあえず今回の後続の目的では高速に検索することを試したかったので(というか、普通、DBにデータをつっこんだら高速に検索したいものです)、POINT型のカラムに空間インデックスを張っておきます。

CREATE SPATIAL INDEX idx_sp1_sp ON sp1(sp);

データを生成するPythonスクリプトを作成

 以下のファイルを作成しておきます。今回は雑な作業なので z という名前のファイルとして作成しました。
作成方法も、 $ cat > z としてからペタっと貼り付けて ctrl-D するという徹底したザツさで。

lat_st = 35.5
lat_en = 35.9
lat_ivl= 0.0007
lng_st = 139.6
lng_en = 140
lng_ivl= 0.0008
ins_str = "INSERT INTO sp1 (sp) VALUES (ST_GeomFromText('POINT(%s %s)', 4326));"

lat=lat_st
while lat<lat_en:
  lat=lat+lat_ivl
  lng=lng_st
  while lng<lng_en:
    lng=lng+lng_ivl
    print (ins_str % (lat, lng))

 少し解説すると、北緯 lat_stからlat_enの範囲(上記では 35.5~35.9)、東経lng_stからlng_enの範囲(同 139.6~140)に、 緯度は lat_ivlおき(同 0.0007)、経度は lng_ivlおき(同 0.0008) の点を生成します。 正確に言うと、生成するというか、そういう点データをINSERTするSQL文を作成します。
 指定範囲に、格子状に点を敷き詰めるイメージですね。バラつきが要求される場合には使えませんが、とりあえず指定範囲にいっぱい点があれば良い、という程度のものには非常に便利です。

実行

 MySQL の test データベースにテーブルが作成されている場合の例。

$ python z | mysql -uroot -p test

または

$ python3 z | mysql -uroot -p test


 結果例(抜粋):以下のようなINSERT文が作成され、上記コマンド例では そのまま mysql に喰わされます。

:
INSERT INTO sp1 (sp) VALUES (ST_GeomFromText('POINT(35.63580000000037 139.73519999999968)', 4326));
INSERT INTO sp1 (sp) VALUES (ST_GeomFromText('POINT(35.63580000000037 139.73599999999968)', 4326));
INSERT INTO sp1 (sp) VALUES (ST_GeomFromText('POINT(35.63580000000037 139.73679999999968)', 4326));
INSERT INTO sp1 (sp) VALUES (ST_GeomFromText('POINT(35.63580000000037 139.73759999999967)', 4326));
INSERT INTO sp1 (sp) VALUES (ST_GeomFromText('POINT(35.63580000000037 139.73839999999967)', 4326));
INSERT INTO sp1 (sp) VALUES (ST_GeomFromText('POINT(35.63580000000037 139.73919999999967)', 4326));
INSERT INTO sp1 (sp) VALUES (ST_GeomFromText('POINT(35.63580000000037 139.73999999999967)', 4326));
INSERT INTO sp1 (sp) VALUES (ST_GeomFromText('POINT(35.63580000000037 139.74079999999967)', 4326));
:


 もうちょっとちゃんとやるなら、BEGIN/COMMIT で括るとか、マルチINSERT化するとかで高速化の余地があります。


 次回の日記では、ここで生成したデータを使って、遊びます。