あった。
この話、どっから説明すりゃぁいいんだろう、て思う話ではあるのよね。少しだけ基礎的な話はしとこうか。
Python でもほかの言語でもなんでもいいんだけど、「リストと辞書(とか連想配列)どっちが適切かしら?」と決断する際の基準はなに? これは要するに「キー検索したいかどうか」、よね? けど「キー検索」って言うけど、これと何が違うのかな?
1 def find_by_key(iterable, key):
2 for it in iterable:
3 k, v = it # タプルで収められとるとして
4 if k == key:
5 return v
そう。辞書を使う決断の際に暗黙にあるのが、(上で示したコードのような)線形検索ではない「賢くて効率が良い」ことを期待するわけだね。最も愚直な辞書的実装の基礎となるアルゴリズムが「二分木」(バイナリツリー)ね。二分木はまぁ絵にするとかなりわかりやすいアルゴリズムなんだけど、まぁわからない人はそれこそ Wikipedia でも見ておくれ。
じゃぁ。「位置」と情報を紐付けたいとしたらどうする? 緯度, 経度、としておこうか。
「緯度経度での完全一致が高速に出来ますなんてステキだぁ!」で思考が終わるバカは置いといて(身近に実在しましたよ)、すぐにまさに「辞書じゃダメじゃん」てことが理解出来るはず。だって引きたい情報のピンポイントの緯度経度がわかってる状況なんか、万に一つとしてないからね。「10km 以内から欲しいぞ」とかなんでしょ、ふつー。
そう。「二分木」に相当する、二次元データを効率よく検索するためのインデクス構造が、R-Tree (等) ね。これも絵にすると割とわかりやすかったりするけど、二分木よりはさすがに凝ってる。かなり似てるけど。てわけで詳細は Wikipedia みて。かなりわかりやすく書かれてたはず。読んだの3、4年前だけど。
さてさて。「効率的に」は「暗黙の期待」だと言いました。そうなんだけどな、そもそも「位置から情報を引く」のってさ、「線形検索でさえうっとうしい」のは理解出来ます? たぶんこんなだよ:
1 def find_by_lomlat(iterable, lon, lat, torr):
2 for it in iterable:
3 data_lon, data_lat, data_value = it # タプルで収められとるとして
4 if ((data_lon - torr <= lon and lon < data_lon + torr) and \
5 (data_lat - torr <= lat and lat < data_lat + torr)):
6 # あるいは GeographicLib 的な高級なものを使って
7 return data_value
グリッドデータなんかなら numpy 使うともちっと「一撃で検索」でけたりもするけれど、そもそも対象にしたいのがいつでも密な行列にしたい(出来る)とは限らないわよね。要するに「世界中の空港一覧」なんてのに「全地球グリッド」なんか使えませんて。(作れないよそんなの。) のでね、いずれにしても(つまり「効率」考えなくとも) R-Tree 欲しいっす、となるわけ。楽でしょ、「緯度経度をキーに辞書的に」に近いノリで検索出来るんだから。
はい、基礎、おしまい。
でね。PostgreSQL + PostGIS などの地理情報を扱える RDB を使う、という選択肢がいまないとして、一応どメジャーなのがあるわけよ。libspatialindexね。C++ で書かれてて速いてのと、「インデクスファイル」まで関知してくれる、結構充実したものね。Python バインディングもある。(こやつの致命的な欠陥は、ドキュメントがヒドいこと。けどまぁなくてもなんとかなるでしょう。)
なんだけど、「pure Python 版」が欲しかったのです。なんでかってのは、android 機の kivyLauncher で使いたいとか、自分で環境制御出来ることが限られるレンタルサーバで使いたいとかね。
あったよ。pyrtree。が、Google Code にいる。つまり、プロジェクトとしてはもう死んでる(なのか忙しくて移行出来てないのかもしらんが)。pip とかでのインストールは出来ないんで、Google Code からソースコードダウンロードしてきて python setup.py install な。
必要なことが出来ればいいかぁ、てことで、「やろうと思ってること」が出来るかだけ確認してみる。こういうことだ:
- R-Tree にはコンテンツとして「情報へのキー」だけ突っ込む
- 位置以外の情報全ては(そのキーに紐付けた)ふつーの辞書
- さて、メモリ上に構築した R-Tree は永続化できるかしら? (つまり普通に pickle とか出来る?)
こんだけ、確認したいのは。保守されてないからといって、GIS の基礎を外してるとはひとまず思ってない。(たまにあるけどね、「モーフィング」の意味を知らずに「モーフィングツール」を公開しちゃう輩とかおるしの。)
まず、検索用のデータではなくて、「元データ」としてこんな csv があるとする:
1 code_query,code,name,location,lon,lat,elev,mag decl.,use
2 AGGH,AGGH,Honiara Henderson Airport,"Guadalcanal, Solomon Islands",160.053,-9.42833333,30.0 feet MSL. ,9° East ,Open to the Public
3 AYMO,AYMO,Momote Airport,"Manus, Papua New Guinea",147.42416667,-2.06,13.0 feet MSL. ,5° East ,Open to the Public
4 AYNZ,AYNZ,Nadzab Airport,"Morobe, Papua New Guinea",146.728,-6.56916667,231.0 feet MSL. ,6° East ,Open to the Public
5 ...
このままでは「検索しずれーでしょーよ」つぅことね。なので検索に適した形にしておきたい、と。位置情報インデクス以外については普通の辞書だけれどもこれは pickle 化するとする。位置情報インデクスも同じように pickle 化(素直に)出来る?
1 # -*- coding: utf-8 -*-
2 # --------------------------------------------------------------
3 # インデクスとデータを永続化するほうのスクリプト、のふり
4 # (元ネタデータが csv として。)
5 from pyrtree import RTree, Rect
6 import csv
7 rt = RTree()
8 reader = csv.reader(open("from_skyvector.csv"))
9 baseinfos = {}
10 for row in reader:
11 # ここは普通に csv てるだけだよ
12 _, code, name, location, lon, lat, elev, magdecl, use = row
13 if name == "name" or not name or not lon:
14 continue
15 # R-Tree にはキーだけ、他の情報は baseinfos に。
16 # この設計がたぶんフツーでしょう。インデクスにあれやこれや
17 # 詰め込んだら肥大するでしょーしな。
18 lon, lat = float(lon), float(lat)
19 baseinfos[code] = (name, location, lon, lat, elev, magdecl, use)
20
21 # この MBR のサイズは目安はあるかいねぇ? 対象によるか…
22 mbr = Rect(lon - 1/60., lat - 1/60., lon + 1/60., lat + 1/60.)
23 #mbr = Rect(lon - 1, lat - 1, lon + 1, lat + 1)
24 rt.insert(code, mbr)
25 #
26 # pickle が最善とも思えないけどまぁひとまず。
27 # (pyrtree 自身は最適なファイル構造には関知してない、と思う。)
28 import pickle
29 pickle.dump(rt, open("rt.pickle", "wb"), protocol=-1)
30 pickle.dump(baseinfos, open("bi.pickle", "wb"), protocol=-1)
うむ、素直に pickle 化出来た。「__reduce__」とか書かなくていい。
んじゃば検索してみろい:
1 # -*- coding: utf-8 -*-
2 # --------------------------------------------------------------
3 # インデクスとデータが永続化済みなことをあてにしてクエリをかける
4 # ほうのスクリプト、のふり
5 from pyrtree import Rect
6 import pickle
7 rt = pickle.load(open("rt.pickle", "rb"))
8 baseinfos = pickle.load(open("bi.pickle", "rb"))
9
10 # 検索したい対象:
11 query_range_dia = 60 * 1852. # NM to metre
12 lon, lat = 139.78116667, 35.55333333
13
14 # 「こんくらいの矩形内で検索せい」のために geographiclib のような
15 # 「厳密な」geodesic 計算をする必要はそんなにはないんだけれど、
16 # pyrtree 同様に pure python で済むのと、「ひどく遅い」ものでも
17 # ないので。(いや、遅いは遅いんだけれども、「pure python」が
18 # 今は大事なので。)
19 # - どうしても速度が気になるなら、簡単な近似計算は使えまする。
20 # - 緯度にもよるけれど、日本の場合はだいたい緯度一秒≒30m、
21 # 経度1秒≒27mくらい、で計算するとまぁまぁ合う。
22 from geographiclib.geodesic import Geodesic
23 direct = Geodesic.WGS84.Direct
24 dr1 = direct(lat, lon, 90, query_range_dia / 2) # 経度方向
25 dr2 = direct(lat, lon, 0, query_range_dia / 2) # 緯度方向
26 lon_d, lat_d = dr1["lon2"] - lon, dr2["lat2"] - lat
27 target_rc = Rect(lon - lon_d, lat - lat_d, lon + lon_d, lat + lat_d)
28
29 # query_point というのもある。ここでは query_rect を。
30 for r in rt.query_rect(target_rc):
31 if r.is_leaf():
32 code = r.leaf_obj()
33 print(code, baseinfos[code])
1 me@host: ~$ python aaa_q.py
2 ('RJTY', ('Yokota Airport', 'Tokyo, Japan', 139.34866667, 35.74866667, '463.0 feet MSL. ', '6\xc2\xb0 West ', 'Military'))
3 ('RJTT', ('Tokyo International. Airport', 'Tokyo, Japan', 139.78116667, 35.55333333, '21.0 feet MSL. ', '7\xc2\xb0 West', 'Open to the Public'))
4 ('RJTF', ('Chofu Airport', 'Tokyo, Japan', 139.528, 35.67166667, '139.0 feet MSL. ', '7\xc2\xb0 West ', 'Open to the Public'))
5 ('RJTK', ('Kisarazu Airport', 'Chiba, Japan', 139.90983333, 35.39833333, '10.0 feet MSL. ', '6\xc2\xb0 West ', 'Military'))
6 ('RJAA', ('Narita International. Airport', 'Chiba, Japan', 140.3855, 35.76533333, '135.0 feet MSL. ', '7\xc2\xb0 West ', 'Open to the Public'))
おっけー期待通り。うんうん、これならあたしのスマホに突っ込めるで。
ちなみに「速度」だけど、測ってません、し、測る気もなし。たぶん遅いんじゃないかな。速度気にするなら libspatialindex 使ってください。ただ当たり前だけど「インデクスを作る」方はそれなりの時間がかかる(体感で 3～5 秒くらいかな? csv は 5409 レコードで)けど、検索はさくっと返ってきます。のでね、性能がシビアに求められないなら実用十分なんじゃないかなと思う。それにしたってこんなん、バカなリニアサーチよりは速いに決まってる(と思う)。
