MySQL5.7でPolygonの内外判定

# はじめに

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

# geometry型

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

https://dev.mysql.com/doc/refman/5.7/en/gis-data-formats.html

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


# Polygonの内外判定

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

https://dev.mysql.com/doc/refman/5.7/en/spatial-relation-functions-object-shapes.html

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

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

# 実際に試してみる

実際に動作して確かめてみましょう。今回ポリゴンのサンプルデータとして使うのは以下の記事でも用いた、[政府統計のe-Stat](https://www.e-stat.go.jp/gis/statmap-search?page=1&type=2&aggregateUnitForBoundary=A&toukeiCode=00200521&toukeiYear=2015&serveyId=A002005212015&datum=2000&coordsys=1&format=shape)からダウンロードできる東京都の境界データです。

https://qiita.com/s2hap/items/4b918a33c1b77c6e43f9

## 急いでる人向け

ちなみに今回用いたコードは全て[Github](https://github.com/mynkit/mysql-geometry)にあげているので、とりあえず動かすところまでやりたい方は[このリポジトリ](https://github.com/mynkit/mysql-geometry)をcloneしてREADMEの通りに実行してみてください。

## MySQLサーバ立ち上げ

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

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


```yaml: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サーバが立ち上がります

```sh
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にアクセス

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

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

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

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

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

```sql
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経由などでインストールしておきます。

```sh
pip install sqlalchemy pymysql
```

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


```python
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](https://www.e-stat.go.jp/gis/statmap-search?page=1&type=2&aggregateUnitForBoundary=A&toukeiCode=00200521&toukeiYear=2015&serveyId=A002005212015&datum=2000&coordsys=1&format=shape)から東京都全域の境界データをダウンロードしてから解凍し、どこかのフォルダに置いておきます。(そのディレクトの相対パスを`./data/A002005212015DDSWC13/`ということにしておきます。)

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

```sh
pip install geopandas
```

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

```python
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']]
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/437488/15e9b886-6751-42bf-5ea1-cc40fd78d8cd.png)


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関数で実行します。(先に断っておくと、おそらくこの処理は失敗します。)

```python
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したほうがうまくいく可能性が上がります。


```python
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)
```

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