あった。
R-Tree な話
2016-09-11に書いた当初のやつ
この話、どっから説明すりゃぁいいんだろう、て思う話ではあるのよね。少しだけ基礎的な話はしとこうか。
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 レコードで)けど、検索はさくっと返ってきます。のでね、性能がシビアに求められないなら実用十分なんじゃないかなと思う。それにしたってこんなん、バカなリニアサーチよりは速いに決まってる(と思う)。
2021-06-07の追記
もともとのノリがかなりいい加減で、目下の自分の目的に叶えばいいや、てことで書いたネタだったんだけれど、なぜだか常に一定数需要があるみたいで。今になって自分でも「R*-Tree ネタ」が必要になっているのだけれども、そういうことならここへの追記の方がいいかなと思って、新たなページを起こすことなく、ここに書き連ねることにした。
だとすると、いくつかの観点からの追記が必要で。
一つ目は、2016年に書いたヤツで紹介した pyrtree の「python 3 対応」について。2016年時点だとワタシ自身もメインは python 2.7 だったので、python 3.x で動作しないことに気付いてなかった。最低でも絶対インポートの問題。python 3.8+ の場合はもう一つ、time.clock の問題がある。そして、まぁ自分で措置したい場合の一番のネックは _NodeCursor の __gt__ ね。ワタシはまだ措置方法がわかってない。
まぁそもそもすでにメンテナンスされてないプロジェクトだからさ、文句言わないこと。
あとね、書いた時に自覚してたのかがもうわからないんだけどね、「libspatialindex」を紹介してるけれど、これ、見返してみても python バインディングがあるわけではないんだよね。文章から察するに、python バインディングが当然のごとく存在している、ということを前提にしてたっぽいけど、実際はこのプロジェクト自身はそれを一切提供してない。だとして、2016年時点で libspatialindex への python バインディングそのものが存在していなかったのかといえばそんなことはなくて、rtree 0.7.0 は 2011年12月30日のリリース。まぁ今言っても詮無いことだけどね、書くならここまで書くべきであったと反省している。
ここまでは「2016年に書いたこと」に直接関係する追記だったが、以降は 2021 年の今だからこそ、の追記。
まず、「pure Python 実装での R-Tree」は、つい最近選択肢が増えた。rtreelib。Python 3.6+ を要求するが、普通はこれは問題ないだろう。一つには Python 3.5 がすでにその命を終えていること、もう一つは、2.7 でそれが必要なら、紹介済みの pyrtree を使えばいいわけだから。作者自身が「まだ始まったばかりなのよさ許してちょんまげ」言うてる通り、目指す機能セットのフル実装にまだ至ってない。けれども、ワタシがここで書いた程度のニーズになら既に耐えられると思う。挿入と検索は出来るんだから。なんとなく発展性の意味で楽しみなプロジェクトなのではないかという気がする。
rtreelib を使って、2016年にやったこととほぼ同じことを、今回は geoshape プロジェクトの成果を使おうとするノリで:
1 # -*- coding: utf-8 -*-
2 # 「インデクス作成」のふりをしでかすほうのスクリプト。
3 import io
4 import csv
5 import pickle
6 from datetime import date
7 from rtreelib import RTree, Rect
8 #from rtreelib import RStarTree, Rect # R*-Treeを使いたいならこっち
9
10
11 def _load_from_geoshape_city_csv():
12 # 『歴史的行政区域データセットβ版』（CODH作成）
13 # https://geonlp.ex.nii.ac.jp/dictionary/geoshape-city/
14 def _map(rec, fields):
15 import re
16 temp = []
17 for v in rec:
18 if re.match(r"\d+\.\d+", v):
19 temp.append(float(v))
20 elif re.match(r"\d+\-\d+\-\d+", v):
21 temp.append(date(*map(int, v.split("-"))))
22 else:
23 temp.append(v)
24 return {k: v for k, v in zip(fields, temp)}
25 reader = csv.reader(
26 io.open("geoshape-city.csv", encoding="utf-8"))
27 fields = next(reader)
28 return (_map(rec, fields) for rec in reader)
29
30
31 def _build_geoshape_city_index():
32 # pyrtree でやったのとほぼ同じノリで、「キー項目だけを R-Tree に」。
33 # 同じく、「今回もめんどっちぃので pickle やがろう」。
34 body = {} # key: entry_id
35 rt = RTree()
36 for rec in _load_from_geoshape_city_csv():
37 # 「代表点」
38 refpt = rec["longitude"], rec["latitude"]
39 # MBR については、先の例に似た「10NM くらい矩形」とするが、
40 # めんどうなので 緯度1秒=30m、経度1秒=27m の近似計算で。
41 mbr = Rect(
42 refpt[0] - 5*1852*(1./3600/27),
43 refpt[1] - 5*1852*(1./3600/27),
44 refpt[0] + 5*1852*(1./3600/30),
45 refpt[1] + 5*1852*(1./3600/30))
46 rt.insert(rec["entry_id"], mbr)
47 body[rec["entry_id"]] = rec
48 #
49 pickle.dump(body, io.open("geoshape-city-body.pickle", "wb"))
50 pickle.dump(rt, io.open("geoshape-city-rt.pickle", "wb"))
51
52
53 if __name__ == '__main__':
54 _build_geoshape_city_index()
1 # -*- coding: utf-8 -*-
2 # 「検索」のふりをしでかすほうのスクリプト。
3 import io
4 import pickle
5
6
7 def main():
8 body = pickle.load(io.open("geoshape-city-body.pickle", "rb"))
9 rt = pickle.load(io.open("geoshape-city-rt.pickle", "rb"))
10 #
11 for found in rt.query((140.88, 38.32)):
12 entry_id = found.data
13 data = body[entry_id]
14 if not data["valid_to"]:
15 print(data["address"])
16
17
18
19 if __name__ == '__main__':
20 main()
1 宮城県富谷市
2 宮城県仙台市宮城野区
3 宮城県仙台市青葉区
4 宮城県仙台市
5 宮城県仙台市泉区
6 宮城県宮城郡泉村
7 宮城県仙台市若林区
いいね、はっきりいって pyrtree と使い勝手は完全に同じ。
せっかくなので、「rtree」もお試してみる。libspatialindex が C++ なことから、本来であれば Windows 版が大変なことになりがちだが、これは spatialindex*.dll が同梱された wheel の形で配布されるので、(少なくとも Python 3.6+ 以降向けとしては)問題ない。使い勝手はまぁおおむね同じといえば同じ:
1 # -*- coding: utf-8 -*-
2 # 「インデクス作成」のふりをしでかすほうのスクリプト(rtree パッケージ版)。
3 import io
4 import csv
5 import pickle
6 from datetime import date
7 from rtree import index
8
9
10 def _load_from_geoshape_city_csv():
11 # 『歴史的行政区域データセットβ版』（CODH作成）
12 # https://geonlp.ex.nii.ac.jp/dictionary/geoshape-city/
13 def _map(rec, fields):
14 import re
15 temp = []
16 for v in rec:
17 if re.match(r"\d+\.\d+", v):
18 temp.append(float(v))
19 elif re.match(r"\d+\-\d+\-\d+", v):
20 temp.append(date(*map(int, v.split("-"))))
21 else:
22 temp.append(v)
23 return {k: v for k, v in zip(fields, temp)}
24 reader = csv.reader(
25 io.open("geoshape-city.csv", encoding="utf-8"))
26 fields = next(reader)
27 return (_map(rec, fields) for rec in reader)
28
29
30 def _build_geoshape_city_index():
31 # pyrtree でやったのとほぼ同じノリで、「キー項目だけを R-Tree に」。
32 # 「今回もめんどっちぃので pickle やがろう」は原則踏襲だが、
33 # rtree パッケージではもとから永続化ファイルをサポートしているので、
34 # それを利用。(index.Index にファイル名を与えるだけ。)
35 body = {} # key: entry_id
36 rt = index.Index("geoshape-city-rt-2")
37 for rec in _load_from_geoshape_city_csv():
38 # 「代表点」
39 refpt = rec["longitude"], rec["latitude"]
40 # MBR については、先の例に似た「10NM くらい矩形」とするが、
41 # めんどうなので 緯度1秒=30m、経度1秒=27m の近似計算で。
42 mbr = (
43 refpt[0] - 5*1852*(1./3600/27),
44 refpt[1] - 5*1852*(1./3600/27),
45 refpt[0] + 5*1852*(1./3600/30),
46 refpt[1] + 5*1852*(1./3600/30))
47 # なんだかんだで pyrtree, rtreelib の RTree にぶち込める
48 # データには制約がないのだが、libspatialindex 由来の rtree
49 # パッケージの場合、insert に与える「id」を一つ発明する
50 # 必要があり、なおかつそれは long integer。「id」といいつつ
51 # ユニークである必要はない。なんだろね、「好きに使え」のノリなのに、
52 # 与えるのを強制されるこの感じ、妙なインターフェイスだ。
53 # pyrtree, rtreelib の RTree でやったノリでの「キー値」は、
54 # obj= のほうに渡せる。こちらは pickle 可能なオブジェクトで
55 # ある必要がある。
56 entry_id = rec["entry_id"]
57 rt.insert(hash(entry_id), mbr, obj=entry_id)
58 body[rec["entry_id"]] = rec
59 #
60 pickle.dump(body, io.open("geoshape-city-body-2.pickle", "wb"))
61
62
63 if __name__ == '__main__':
64 _build_geoshape_city_index()
1 # -*- coding: utf-8 -*-
2 # 「検索」のふりをしでかすほうのスクリプト(rtree パッケージ版)。
3 import io
4 import pickle
5 from rtree import index
6
7
8 def main():
9 body = pickle.load(io.open("geoshape-city-body-2.pickle", "rb"))
10 rt = index.Index("geoshape-city-rt-2")
11 #
12 # objects=False とすると、「仕方なく渡した id」だけの返却になる。
13 # これでいい場合もあろうし、そうでない場合もあろう。管理が「数値キー」の
14 # データセットの場合だと、それでもいいんだろうなぁ。
15 # (無論「hash」が一意である保証が持てる場合も「それでいい」。)
16 for found in rt.intersection((140.88, 38.32), objects=True):
17 entry_id = found.object
18 data = body[entry_id]
19 if not data["valid_to"]:
20 print(data["address"])
21
22
23
24 if __name__ == '__main__':
25 main()
1 宮城県仙台市若林区
2 宮城県仙台市宮城野区
3 宮城県仙台市
4 宮城県仙台市青葉区
5 宮城県宮城郡泉村
6 宮城県仙台市泉区
7 宮城県富谷市
pyrtree 版と同じノリで「インデクスと実体を分けて考える」としたけれど、一応 rtree のノリとしては「obj=」に実体を突っ込むことを念頭に置いているんだとは思う。ただ、libspatialindex の考える「id」が long int であることの制約を考えると、たぶんワタシがやったやり方の方が可搬性があるのだろうと思う。
pure python でないものとしてはもう一つ、python_prtree があり、これも新しい。C++ で書かれているものの、Windows 版のものも python 3.x に関しては wheel の形で提供されているので、インストールには困らない。(名前の通りこれはオリジナルの R-Tree そのものではなくて「Priority R-Tree」。意味はワタシにきかないで。)
python_prtree については、まず、ドキュメントを読むだけで使いにくさがわかることになっている。何がマズいって、rtree パッケージでいうところの「obj=」に相当するものがなく、そして rtree パッケージと同じく id が long int である必要があるので、「long int の id の発明と、その id とそれに対応するオブジェクトへの対応管理表の「再」発明」が必要である、てことね。それだけでなくて、まぁ「C 拡張」ならではというかさ、問題が起こった場合のエラー報告がまぁはっきりいって零点で、全然先に進めなくて。なんというか「今じゃない」て感じだなぁ。一応やるだけはやってみたが、結構ムゴい:
1 # -*- coding: utf-8 -*-
2 # 「インデクス作成」のふりをしでかすほうのスクリプト(python_prtree パッケージ版)。
3 import io
4 import csv
5 import pickle
6 from datetime import date
7 from python_prtree import PRTree2D
8
9
10 def _load_from_geoshape_city_csv():
11 # 『歴史的行政区域データセットβ版』（CODH作成）
12 # https://geonlp.ex.nii.ac.jp/dictionary/geoshape-city/
13 def _map(rec, fields):
14 import re
15 temp = []
16 for v in rec:
17 if re.match(r"\d+\.\d+", v):
18 temp.append(float(v))
19 elif re.match(r"\d+\-\d+\-\d+", v):
20 temp.append(date(*map(int, v.split("-"))))
21 else:
22 temp.append(v)
23 return {k: v for k, v in zip(fields, temp)}
24 reader = csv.reader(
25 io.open("geoshape-city.csv", encoding="utf-8"))
26 fields = next(reader)
27 return (_map(rec, fields) for rec in reader)
28
29
30 def _build_geoshape_city_index():
31 # pyrtree でやったのとほぼ同じノリで、「キー項目だけを R-Tree に」。
32 # 「今回もめんどっちぃので pickle やがろう」は原則踏襲だが、
33 # python_prtree パッケージではもとから永続化ファイルをサポートしているので、
34 # それを利用。(save メソッドがある。)
35 body = {} # key: hash(entry_id)
36 # PRTree2D は「insert」で逐次構築できそうに思うが、「パフォーマンスチューニン
37 # グされてない」という説明以上の問題を抱えているらしく、「なんか起こったぜ」と
38 # 叱られるのみで一向に動作しない。結局コンストラクタに全部渡す方法しか今は動作
39 # しない模様。
40 mbrs = {}
41 for rec in _load_from_geoshape_city_csv():
42 # 「代表点」
43 refpt = rec["longitude"], rec["latitude"]
44 # MBR については、先の例に似た「10NM くらい矩形」とするが、
45 # めんどうなので 緯度1秒=30m、経度1秒=27m の近似計算で。
46 # なお、MBR は xmin, xmax, ymin, ymax(, zmin, zmax) で
47 # 与えなければならないのがほかと違う。
48 mbr = (
49 refpt[0] - 5*1852*(1./3600/27),
50 refpt[0] + 5*1852*(1./3600/30),
51 refpt[1] - 5*1852*(1./3600/27),
52 refpt[1] + 5*1852*(1./3600/30))
53 # pyrtree, rtreelib, rtree の3つと比較すると一番マズい。
54 # rtree と同じで「long int の id」を必要とするが、
55 # PRTree2D は要はその id と mbr との紐づけしか担わない。
56 entry_id = rec["entry_id"]
57 id = hash(entry_id) # まぁこれがちゃんと一意なら問題はないのだけどね。
58 #rt.insert(id, mbr) # こうしたいんだよほんとは。Unknown Exception てなにゃの
59 body[id] = rec
60 mbrs[id] = mbr
61 #
62 ids = list(body.keys())
63 rcs = [list(mbrs[id]) for id in ids]
64 rt = PRTree2D(ids, rcs)
65 rt.save("geoshape-city-rt-3.bin")
66 pickle.dump(body, io.open("geoshape-city-body-3.pickle", "wb"))
67
68
69 if __name__ == '__main__':
70 _build_geoshape_city_index()
1 # -*- coding: utf-8 -*-
2 # 「検索」のふりをしでかすほうのスクリプト(python_prtree パッケージ版)。
3 import io
4 import pickle
5 from python_prtree import PRTree2D
6
7
8 def main():
9 body = pickle.load(io.open("geoshape-city-body-3.pickle", "rb"))
10 rt = PRTree2D("geoshape-city-rt-3.bin")
11 #
12 # rtree パッケージ版と同じく long int の id での授受だが、 rtree
13 # パッケージ版と違って、「対応するオブジェクト」と紐づけるキーと
14 # id が一致している必要がある、というのがキツい制約。しかも
15 # bb での のクエリしかできない。(bb は xmin, xmax, ymin, ymax
16 # (, zmin, zmax)。)
17 for entry_id_hash in rt.query((140.88, 140.88, 38.32, 38.32)):
18 data = body[entry_id_hash]
19 if not data["valid_to"]:
20 print(data["address"])
21
22
23
24 if __name__ == '__main__':
25 main()
1 宮城県仙台市宮城野区
2 宮城県仙台市若林区
3 宮城県仙台市青葉区
4 宮城県富谷市
5 宮城県仙台市泉区
6 宮城県仙台市
7 宮城県宮城郡泉村
openmp による並列クエリ(batch_query)はおいしいかもと思う。今後に期待、という感じかなぁ?
