GeoPandasは100万行の夢を見るか


<img src="https://qiita-image-store.s3.amazonaws.com/0/100567/531160a3-6a1b-92e2-aafc-4589d966e14b.png" width=60%>

**年末オープンジオかくし芸大会** でおなじみの、[FOSS4GAdventカレンダー2018](https://qiita.com/advent-calendar/2018/foss4g)の23日目です。

## はじめに

この記事は

- 地図や位置のデータをPythonで処理したい
- Python初学者
- QGISちょっとわかる
- PostGISやMySQLなどのサーバー系はちょっとしんどい
- 商用GISソフトは手元にない
- 暗中模索話だいすき

という方むけの検証記事です。どうぞよろしくおねがいいたします。

## きっかけ
　いろいろとアレで大量のcsvファイルをなだめたりすかしたりクリーニングしたりするために、PythonのライブラリであるPandasをよく使っています。Pandasについては[多くの知見](https://qiita.com/search?q=pandas)があるため割愛しますが、表形式のデータを取り扱うのに様々な命令があるため、結構便利ですし、今後もいろいろ使ってみたいなあと考えています。
※[こちらのパンダさん](https://qiita.com/ozo360):panda_face:とはちょっと関係がないかもしれませんが、どこかでなにかつながりがあるかもしれませんね。

さてPandasから地図(位置)のデータもさっくり処理できちゃうと、わたくしとしては業務的に便利かもしれないと思い調べてみたら、[GeoPandas](https://geopandas.readthedocs.io/en/latest/)というライブラリがありました。一応、[FOSS4Gの一員](http://wiki.geosteiniger.cl/mediawiki-1.22.7/index.php/FOSS4G_Software_List)でもあるようです。
(not integrated toolsってきになりますがまあいいか)

※2018年のGeopython conferenceでもGeoPandasのワークショップが開かれるなど、すごくマイナーなわけではなさそうです
http://2018.geopython.net/#w4

## ねらい
GeoPandasの導入や視覚化などについては、[id:sanvarieさんによる情報](https://www.gis-py.com/entry/2018/10/09/115021)や、Qiitaの[GeoPandasタグあたり](https://qiita.com/tags/geopandas)をご参考にしてみてください。

他のエントリでは、視覚化やデータ変換などについて紹介がありますが、わたしとしては

- pythonで空間結合がさっくりできちゃうと、なんかうれしい
- （あわよくば）100万行(レコード)とかもさくっと処理できるかしら

　というゴールで試してみます。

空間結合とは、大まかに言うとお店などのポイントデータに対して、メッシュや街区などのポリゴンデータから、属性を取ってくる処理です。串刺し処理ともいいますね(意外と引用可能な模式図とかないので、[例のチラシ](https://howmori.org/2018/07/26/foss4g-hokkaido-2018%E3%83%87%E3%82%B6%E3%82%A4%E3%83%B3/)裏の落書きですみません）
<img src="https://qiita-image-store.s3.amazonaws.com/0/100567/fc68deb5-de43-67de-a0b3-89c200776195.png" width=80%>
みんな大好きQGISでももちろんこのような処理は可能なのですが、対象レコードが万オーダーくらいになると、全然処理が終わらなかったり、ソフトそのものが落ちたりします。
また、空間処理をpythonのみで回せると、設計的な縛りにうまくとりこめたりするかもしれません。
ほかにも、明日までにAフォルダ内にあるポリゴン１０００ファイルに対し、Bフォルダファイルのポイント１０００ファイルについてすべての組み合わせで空間結合やんなきゃ…という場合とかは、このGeoPandasを活用すればできるとおもいます。

## GeoPandasをためしてみる

今回の実行環境スペックは

- ProductName:	Mac OS X
- ProductVersion:	10.13.6
- BuildVersion:	17G65
- Python 3.6.７
- GeoPandas 0.3.0
- QGIS 3.4.2

です。

あと、わたしのスペックは

- QGISよく使う
- PythonもPandasも理解はざっくり
- 好きな食べ物はじねんじょ

ぐらいのヒヨコですので、いいサジェスチョンがあればとても嬉しいです。

まず、[こちらのサイト](https://www.gis-py.com/entry/2018/10/09/115021)を参考にAnacondaに入れてみました。
なお、Anaconda上ではEnvironmentsタブでインストールされているかチェックをしたほうがよいです。
![©OpenStreetMap　Contributors/MIERUNE](https://qiita-image-store.s3.amazonaws.com/0/100567/a261e909-2588-f2f7-f32d-32adf1053e37.png)

## データ準備
まずは軽めのデータということで、QGISのベクタメニューから調査ツール→グリッド作成を使って、六角形グリッド60個(3857_hex500_n60)とランダムポイント180個(3857_point_n180)をGeoJSON形式で作成します。

※GeoPandasはSHPもダイレクトで読み込んでくれるんですが、ええまあ最近は[こんな流れ](http://switchfromshapefile.org/#shapefileisbad)もあったりしますので、GeoJSONでやってみましょう…）
※対応フォーマットはこの辺に情報があります
http://geopandas.org/io.html

※あと六角形ポリゴンなんて[大戦略](https://ja.wikipedia.org/wiki/%E5%A4%A7%E6%88%A6%E7%95%A5%E3%82%B7%E3%83%AA%E3%83%BC%E3%82%BA)ヲタ以外では普通あんま使わないとおもうのですが、生まれて初めて作りました。

![©OpenStreetMap Contributors/MIERUNE](https://qiita-image-store.s3.amazonaws.com/0/100567/76149159-65ae-39c7-7b4a-7f484acad413.png)

hogehogeというフォルダに保存したこれらのGeoJSONをGeoPandasで表示・確認する場合は

```python
import geopandas as gpd

# GeoJSON読み込み
pt = "/hogehoge/3857_point_n180.geojson"
hx = "/hogehoge/3857_hex500_n60.geojson"
df_pt = gpd.read_file(pt) #ランダムポイント
df_hx = gpd.read_file(hx) #HEXポリゴン

# GeoJSON描画
base = df_hx.plot(color='white', edgecolor='black')
df_pt.plot(ax=base, marker='o', color='red', markersize=5)
```
としてAnaconda上のJupyter notebookで
![image.png](https://qiita-image-store.s3.amazonaws.com/0/100567/4cce104c-2b0f-0d0f-bcdb-c45fd6680841.png)

このような可視化ができます、QGISと同じ分布で出ていますね。背景もなんかすれば出るのかもしれませんがまあいいや。

## いざ空間結合
では、空間結合をやってみましょう。GeoPandasでは
[**geopandas.sjoin**](http://geopandas.org/reference/geopandas.sjoin.html)を用いて処理することができます。

この処理はGeoPandasでは

`spj = geopandas.sjoin(df_pt, df_hx, op='within')`

の一行で済みます。
これの意味は 左から読むと

- **spj**というジオな箱(データフレーム)に結果を出力して頂戴
- **geopandas.sjoin**という機能で空間結合してね
- 刺すデータは**df_pt** というポイントデータ
- 刺される**df_hx** というポリゴンデータ
- なお、含まれる(within)という条件で空間結合ですよ

となります。

geopandasという名称は長いのでgpdと省略させて最初からコードをまとめると

```python
import geopandas as gpd

# GeoJSON読み込み
pt = "/hogehoge/3857_point_n180.geojson"
hx = "/hogehoge/3857_hex500_n60.geojson"
df_pt = gpd.read_file(pt) #ランダムポイント
df_hx = gpd.read_file(hx) #HEXポリゴン

# 空間結合
spj = gpd.sjoin(df_pt, df_hx, op='within')

# 先頭5行を確認
print(spj.head())

# GeoJSONで書き出し
output = "/hogehoge/rslt.geojson"
spj.to_file(output, driver="GeoJSON")


```
これを実行してみると、結果のデータフレーム出力では、ポイントデータにポリゴンデータの属性が結合されています

![image.png](https://qiita-image-store.s3.amazonaws.com/0/100567/276fedff-6cde-55a3-cad4-e00bcf4da44a.png)

また、出力されたGeoJSONをQGISで確認すると、空間結合を行った結果カラム（列）が追加されています。接尾語が長々とつくのはまあご愛嬌ということで。
![image.png](https://qiita-image-store.s3.amazonaws.com/0/100567/dbb4a81a-7ea1-9612-6dea-1c51d266fcf6.png)

これでGeoPandasで空間結合ができましたね

## 処理速いの？
それな。
夢の100万行処理にあたって速度は気にしたいので、timeモジュールを用いて、計測してみました。
出力も省略した測定用コードをこんな感じにしてみました。

```python
import time
import geopandas as gpd

# 開始時間取得
t1 = time.time() 

# GeoJSON読み込み
pt = "/hogehoge/3857_point_n180.geojson"
hx = "/hogehoge/3857_hex500_n60.geojson"
df_pt = gpd.read_file(pt) #ランダムポイント
df_hx = gpd.read_file(hx) #HEXポリゴン

# 空間結合
spj = gpd.sjoin(df_pt, df_hx, op='within')

# GeoJSONで書き出し
output = "/hogehoge/rslt.geojson"
spj.to_file(output, driver="GeoJSON")

# 処理後の時刻
t2 = time.time()

# 経過時間を表示
elapsed_time = t2-t1
print(f"経過時間：{elapsed_time}")
```
ということで結果は

`0.14秒`

でした。比較対象として同じデータをQGISの「場所で属性を結合」で行うと

`0.67秒`

だったので約80％の高速化がはかられています。
<img src="https://qiita-image-store.s3.amazonaws.com/0/100567/4c98bc90-161f-ec27-ca30-5ddd63485bd1.png" width=60%>

では、同様に6万個のポリゴンと10万行のポイントで空間結合をやってみると
GeoPandasでは
`56.80秒`
となりました

しかし、QGISですと**夜通し14時間経っても2％しか進まない状態**となりました…

![image.png](https://qiita-image-store.s3.amazonaws.com/0/100567/3c004780-98e6-747f-5a53-e69149cbf42a.png)

（いろいろ調整したら動くかもしれませんが）このペースで行くと約29日（700時間）ほどかかることになりますので、[OHP](https://www.google.com/search?q=OHP%E3%82%B7%E3%83%BC%E3%83%88&source=lnms&tbm=isch&sa=X&ved=0ahUKEwiOnZOewLPfAhUEv7wKHeLqCgoQ_AUIDigB&biw=1149&bih=655)とバイトさんで確認したほうが早いやんけ！」と憤る上司が現れてもおかしくない状況となります。

…おそらく29日かかる推定なので比較するべきではないと思いますが、参考までにグラフ化すると

<img src="https://qiita-image-store.s3.amazonaws.com/0/100567/213dc5c3-8dee-338d-e50b-cac9cb940421.png" width=60%>

となります。

## もっと速くならんのか

ということで、GeoPandasはQGISより早いとはいえ、10万行で56秒なら、夢の100万行だと概算10分以上はかかることなります。まあ待てなくもないのですが、もうすこしなんとかしたいところ。

もともとPythonは速くない言語ですので、ママチャリに新幹線なみの速度を求めるのは酷ですが、ちょっと調べてみました。

たとえば　**Pickle（酢漬け）** というモジュールによってデータ形式を変換して[高速化をはかる手法](https://note.nkmk.me/python-pandas-to-pickle-read-pickle/)がデータマイニング界隈などであります。今回の位置データをPickle化して試してみましたが、Pickle化を行うところでボトルネックが起きてしまうのか、6万個のポリゴンと10万行のポイントの空間結合で
`59.70秒`
と、ノーマルの56.80秒より遅くなっています。

<img src="https://qiita-image-store.s3.amazonaws.com/0/100567/ddcfe87f-e7fb-18f3-481d-c57a2b5dce0e.png" width=60%>

あと考えられることとしては、マルチプロッセンシングによる並列処理化を行う[Dask](https://orizuru.io/blog/machine-learning/pandas-dask/)なども使うという手もあるかもしれませんが、ためしていません。でも結局データ読み込み書き込みはそれなりにかかるので、ちょっと厳しいかなという予想をしています。

## 結局そういうことなのね

さらにいろいろ調べていたところ、10月にあったFOSS4G Belgium 2018でGeoPandasの発表があったようです。御仁がおっしゃるには[Cython](https://qiita.com/pashango2/items/45cb85390193d97523ca)使うと[PostGIS並みに速くなるよ！](https://jorisvandenbossche.github.io/talks/2018_FOSS4GBE_geopandas/#1)ということで、そうか、Cythonによるブーストか…という元も子もないオチに打ちひしがれました。
（ミニ四駆にガソリンエンジン載せるような感覚かな、と）
なお、GeoPandas+Cythonに関する[補足情報はここにもあります](https://jorisvandenbossche.github.io/blog/2017/09/19/geopandas-cython/)。
![by Joris Van den Bossche/FOSS4G Belgium　2018](https://qiita-image-store.s3.amazonaws.com/0/100567/d48f0cce-1d68-5e30-3df2-f8ff334e11ff.png)


まだまだ試すべき課題や壁が発見されてしまいましたが、今回使ったGeoJSONデータは[ここにおいています](https://www.dropbox.com/s/2icd66q9v3tbhp5/geopandas_sample_data.zip?dl=0)。ここまでお読みくださった方は、わたしの屍をコヤシにさらなるpython縛りのハイスピードシンプル空間データ処理をいっしょに考えてくださると嬉しいです。

※なお、QGISでドチャクソ遅い属性結合もGeoPandasでは[GISファイル直読みで属性結合できる](http://geopandas.org/mergingdata.html#attribute-joins)ので、そのうちためしてみたいですね。
以前はSHP→csv→Pandas→計算→属性結合でそこそこ手間がかかっていたので、これはありがたいかも。

## まとめ

今回の検証まとめとして
GeoPandasは…

- GISファイルが直読みできる
- QGISよりも早く空間結合処理ができる
- Python上でさまざまな処理に組み込める
- 100万行処理はなんとかできるけどもっと速度チューニングができるはず

ということがわかりました :panda_face: 

## この記事のライセンス
[![クリエイティブ・コモンズ・ライセンス](https://i.creativecommons.org/l/by/4.0/88x31.png)](http://creativecommons.org/licenses/by/4.0/deed.ja)
この記事は[CC BY 4.0（クリエイティブ・コモンズ 表示 4.0 国際 ライセンス）](http://creativecommons.org/licenses/by/4.0/deed.ja)の元で公開します、どうぞどうぞ。