MySQL5.7でPolygonの内外判定

はじめに

ポリゴンの内外判定はPythonのshapelyなどを使うことで実現可能です。ただサーバサイドの言語に依存させたくないのと、MySQLでもポリゴンの内外判定ができるらしく、使い勝手もよさそうなので試してみました。

geometry型

MySQLではgeometry型という位置情報系のデータ型を用いることができます。

これによってポリゴンなどのデータもデータベース側に持つことができます

Polygonの内外判定

ST_Containsという関数を用いればポリゴンの内外判定ができます。

たとえばgeometriesテーブルという、geometry型のpolygonカラムを持つテーブルがあったとして、以下のようなクエリを叩けば、緯度経度35.67405, 139.77250の位置を含むレコードが返ってきます。

SELECT
    *
FROM
    geometries
WHERE ST_Contains(polygon, GeomFromText("POINT(139.77250 35.67405)")) = 1
;

実際に試してみる

実際に動作して確かめてみましょう。今回ポリゴンのサンプルデータとして使うのは以下の記事でも用いた、政府統計のe-Statからダウンロードできる東京都の境界データです。

急いでる人向け

ちなみに今回用いたコードは全てGithubにあげているので、とりあえず動かすところまでやりたい方はこのリポジトリをcloneしてREADMEの通りに実行してみてください。

MySQLサーバ立ち上げ

MySQLでデータベースを作成したいのですが、今回はdockerを使っていきます。

以下のようなdocker-composeファイルを用意します。

docker-compose.yml

version: "2"
services:
  mysql:
    container_name: mysql_host
    image: mysql:5.7.12
    ports:
      - 3306:3306
    environment:
      TZ: Asia/Tokyo
      MYSQL_ROOT_PASSWORD: rootpw
      MYSQL_DATABASE: sample_db
      MYSQL_USER: user
      MYSQL_PASSWORD: E9c9y6IwaDzNtLAfJn3qd3Sa20pOCz
      BIND-ADDRESS: 0.0.0.0
    tty: true
    command: mysqld --character-set-server=utf8mb4 --collation-server=utf8mb4_unicode_ci

詳しくは解説しないですが、以下のコマンドでmysqlサーバが立ち上がります

docker-compose up -d

docker-compose psして以下のように表示されればちゃんと立ち上がってます。

$ docker-compose ps
   Name                 Command               State           Ports         
----------------------------------------------------------------------------
mysql_host   docker-entrypoint.sh mysql ...   Up      0.0.0.0:3306->3306/tcp

mysqlにアクセス

ローカルからデータベースにログインしてみましょう。

# ローカルにMySQLがインストール済みなことが前提
mysql -h 127.0.0.1 -uroot -prootpw

でアクセスできるはずです。

データベース、geometriesテーブルの作成

ログイン後、以下のコマンドでsample_dbデータベースにgeometriesテーブルを作ります。

CREATE DATABASE IF NOT EXISTS sample_db;

USE sample_db;

DROP TABLE IF EXISTS `geometries`;

CREATE TABLE IF NOT EXISTS `geometries` (
  `id` int(11) NOT NULL AUTO_INCREMENT,
  `addressCode` varchar(255) NOT NULL,
  `address` varchar(255) NOT NULL,
  `polygon` geometry  NOT NULL,
  PRIMARY KEY (`id`)
) ENGINE=InnoDB DEFAULT CHARSET=utf8;

東京都の境界データをgeometriesテーブルに挿入

ここからは挿入するデータの量も多いのと、データ整形も必要なのでPythonで作業していきます。

まずPython経由でMySQLにアクセスする必要があるので、以下のパッケージをpip経由などでインストールしておきます。

pip install sqlalchemy pymysql

PythonからInsertしたいので、以下の関数を定義しておきます。

from sqlalchemy import create_engine

def execute(sql, user, password, host, port, db):
    con = f'mysql+pymysql://{user}:{password}@{host}:{port}/{db}'
    engine = create_engine(con)
    connect = engine.connect()
    connect.execute(sql)
    connect.close()

次に、政府統計のe-Statから東京都全域の境界データをダウンロードしてから解凍し、どこかのフォルダに置いておきます。(そのディレクトの相対パスを./data/A002005212015DDSWC13/ということにしておきます。)

shapefileをデータフレームとして扱いたいので、geopandasをインストールしておきます。

pip install geopandas

以下のようにshapefileを読み込み、軽くデータ整形しておきます。

import geopandas

gdf = geopandas.read_file('./data/A002005212015DDSWC13/h27ka13.shp')
gdf['addressCode'] = gdf.apply(
    lambda x: str(x['PREF']) + str(x['CITY']) + str(x['S_AREA']), axis=1)
gdf['address'] = gdf.apply(
    lambda x: x['PREF_NAME'] + x['CITY_NAME'] + (x['S_NAME'] if x['S_NAME'] is not None else ''), axis=1)
gdf[['addressCode', 'address', 'geometry']]

geometry型のカラムにデータを挿入するためには、

POLYGON ((140.34225 29.79356, 140.34149 29.793...))

となっているgeometry文字列を

GeomFromText('POLYGON ((140.34225 29.79356, 140.34149 29.793...))')

のように変換してあげる必要があります。

そのため、Insert文は以下のように生成でき(geo_text)、これを先ほど作ったexecute関数で実行します。(先に断っておくと、おそらくこの処理は失敗します。)

geo_text_list = [
    f"('{r['addressCode']}', '{r['address']}', GeomFromText('{r['geometry']}'))" for _, r in gdf.iterrows()]
geo_text = ','.join(_geo_text_list)
execute(
    f'''
    INSERT INTO
        geometries
    (addressCode, address, polygon)
    VALUES
    {geo_text}
    ;
    ''',
    user='root',
    password='rootpw',
    host='127.0.0.1',
    port=3306,
    db='sample_db'
)

ただこれだとInsertするデータが多すぎておそらく失敗します。そのため、いくつかの処理にわけてInsertしたほうがうまくいく可能性が上がります。

import numpy as np

geo_text_list = [
    f"('{r['addressCode']}', '{r['address']}', GeomFromText('{r['geometry']}'))" for _, r in gdf.iterrows()]
# geo_text_listをひとつの要素数が100個くらいになるように分割
split_num = 100
geo_text_list_split = [list(_) for _ in np.array_split(
    geo_text_list, len(geo_text_list)//split_num)]
print('insert into geometries.')
print(f'num of insert records: {len(geo_text_list)}')
for _geo_text_list in geo_text_list_split:
    geo_text = ','.join(_geo_text_list)
    execute(
        f'''
        INSERT INTO
            geometries
        (addressCode, address, polygon)
        VALUES
        {geo_text}
        ;
        ''',
        user='root',
        password='rootpw',
        host='127.0.0.1',
        port=3306,
        db='sample_db'
    )

おそらくこれで挿入がうまくいっているはずです。念のため、上のコード出力させたnum of insert recordsの値と、SELECT COUNT(*) FROM geometries;の出力結果が同じであることを確認しておいてください。

ST_Containsで内外判定

ここまでできたら、あとは冒頭のクエリをなげて試すだけです。

$ mysql -h 127.0.0.1 -uroot -prootpw sample_db
mysql> SELECT
    ->     id,
    ->     addressCode,
    ->     address
    -> FROM
    ->     geometries
    -> WHERE  ST_Contains(polygon, GeomFromText("POINT(139.77250 35.67405)")) = 1
    -> ;
+-------+-------------+-----------------------------------+
| id    | addressCode | address                           |
+-------+-------------+-----------------------------------+
| 17256 | 13102003001 | 東京都中央区銀座１丁目            |
+-------+-------------+-----------------------------------+
1 row in set, 1 warning (0.10 sec)

上のように返ってくれば成功です。