「立体地形図を手作りしてみようかな、っと(8.5) 兼 六芒星」では 地理院地図WEBとか Google Map でチマチマ調べてプロットしたんだけれども。
OpenStreetMap のネタはそのうちまとまったものを書いてもいいかなとは思っている。というか何年か前に仕事で調べて結構な資料書いたんだけど、一般向けにしつつ、ついでに最新情報を仕込みたい。OpenStreetMap は「良いもの」でもあり「良くないもの」でもある、全体的に「わんぱくな」プロジェクトなんだけれど、そういったことも書けたらいいかなと思うし。とはいえこれは「膨大な」ものであって、まぁそんなにすぐにはやらない。
というかね。「立体地形図を手作りしてみようかな、っと(8.5) 兼 六芒星」で思い出したわけですよ。あぁそういや、狙った獲物だけ情報ダウンロード出来るわ、と。
OpenStreetMap プロジェクト、というのは、一つには「Google Map Alternative」なわけね。なおかつ、「WikiPedia いんすぱいあど」な、「誰でも地図作りに参加出来る」。これがこのプロジェクトの一番大きなミッションなわけな。
このプロジェクトの至上命題はもう2つあって、一つには「誰もが使える(ソフトウェアの)インフラ整備」。Mapnik とか、mod_tile とか。OpenStreetMap との直接の関係は不明だけど、OpenLayers との親和性も高い。OpenLayers は OSGeo のプロジェクト。OpenStreetMap が OSGeo と関係あるのかないのか、よくわかんない。もう一つが、「自由な GIS データ利用」。ダウンロードしてデータだけ使える、ということ。
このね、「ダウンロードしてデータだけ使う」のデータダウンロード手段の一つに、「XAPI」(ざっぴー)があるのです。HTTP GET で、「XPath いんすぱいあど」な API で取得する。無謀なデータ取得方法には「アジア全域」「日本全域」のような全データ取得があるけれど、圧縮した状態でうん十うん百「ギカ」なんてサイズで、そんじょそこらの「ぴぃしぃ」では取り回し出来ないほどの大規模データなので、現実的じゃない。(間違って持ってこれても、osm2pgsql がメモリエラーで落ちたり、落ちなくても丸一日かかっても処理が終わらなかったり。)それに較べると XAPI での「狙い撃ちデータ取得」は、大変現実的。
てな話。
この話、「思い出した、今でも同じこと出来るかな?」程度のことで、今回のエントリではどうデータ活用するかのとこまで行くつもりがない。単に「ダウンロードしとこっかな」でおしまい。
利用に必要な情報は3つ。Xapiと、Map Features、それと、OSM ファイルの構造。OSM ファイルの構造説明は…今回はやめとく。GML と違って簡単な構造で、OSGeo の Simple Feature Model ほどには単純過ぎないもの、で、別に説明が難解なものでもないんだけれど、今回は面倒なんで。簡単に言えば OSM ファイルの構造を意識したクエリを書きます、Xapi では。
Xapi を「XPath 的なもの」とみなすには相当制約ありすぎるけれど、単一条件での絞り込みくらいは出来る。本当の制約は、「bbox queries have a maximum of 10 square degrees」という明示的な制約と、結果セットが大きくなると Bad Request となる暗黙の制約の2つ。ので、「0.5度刻みくらいでチマチマと」取りに行く、というのが正解。
てわけで、「寺社仏閣系と、鉄道の駅くらいは欲しいよねぇ」をまとめて wget するスクリプト：
1 # -*- coding: utf-8 -*-
2
3 #
4 # build download list
5 #
6 import math
7 lon_step = 0.5
8 lat_step = 0.5
9
10 #http://www.overpass-api.de/api/xapi?*[religion=*][bbox=139.0608,35.6227,139.5608,35.7227]
11 #http://www.overpass-api.de/api/xapi?*[building=train_station][bbox=139.0608,35.6227,139.5608,35.7227]
12
13 bottom_left = (138.5, 33.5)
14 top_right = (142.5, 39.5)
15
16 _BASEURI = "http://www.overpass-api.de/api/xapi?"
17 with open("._downloadlist", "w") as fo:
18 for x in range(int(math.ceil((top_right[0] - bottom_left[0]) / lon_step))):
19 for y in range(int(math.ceil((top_right[1] - bottom_left[1]) / lat_step))):
20 min_long = bottom_left[0] + x * lon_step
21 min_lat = bottom_left[1] + y * lat_step
22 max_long = bottom_left[0] + (x + 1) * lon_step
23 max_lat = bottom_left[1] + (y + 1) * lat_step
24
25 # use XAPI
26 fo.write(
27 _BASEURI + "*[religion=*][bbox=%.4f,%.4f,%.4f,%.4f]\n" % (
28 min_long, min_lat, max_long, max_lat))
29 fo.write(
30 _BASEURI + "*[building=train_station][bbox=%.4f,%.4f,%.4f,%.4f]\n" % (
31 min_long, min_lat, max_long, max_lat))
32
33 #
34 # call wget with created list
35 #
36 import subprocess
37 rc = subprocess.call(" ".join(["wget", "-c", "-i", "._downloadlist"]), shell=True)
ちょいちょいバカだけど、いいんです。はなから真面目にやる気ないんだし。「building=train_station」はわかりやすいでしょう。「religion=*」はね、「religion=shinto」「religion=buddhist」なんてのが入ってますわよ。キリスト教系のもひっかかってくるよ。
religion=* の方はこんな感じ：
1 <?xml version="1.0" encoding="UTF-8"?>
2 <osm version="0.6" generator="Overpass API">
3 <note>The data included in this document is from www.openstreetmap.org. The data is made available under ODbL.</note>
4 <meta osm_base="2015-10-10T18:48:02Z"/>
5
6 <node id="442152305" lat="35.4371048" lon="139.3671385">
7 <tag k="amenity" v="place_of_worship"/>
8 <tag k="denomination" v="catholic"/>
9 <tag k="name" v="カトリック厚木教会"/>
10 <tag k="religion" v="christian"/>
11 </node>
12 <node id="448176637" lat="35.4853823" lon="139.3950393">
13 <tag k="amenity" v="place_of_worship"/>
14 <tag k="name" v="水上山 龍源院"/>
15 <tag k="religion" v="buddhist"/>
16 </node>
17 <node id="564634835" lat="35.4613393" lon="139.3985197">
18 <tag k="amenity" v="place_of_worship"/>
19 <tag k="name" v="弥生神社"/>
20 <tag k="religion" v="shinto"/>
21 <tag k="website" v="www11.plala.or.jp/yayoi1/"/>
22 </node>
23 <node id="564637579" lat="35.4679175" lon="139.4020130">
24 <tag k="amenity" v="place_of_worship"/>
25 <tag k="denomination:ja" v="曹洞宗"/>
26 <tag k="name" v="常泉院"/>
27 <tag k="religion" v="buddhist"/>
28 </node>
29 <node id="599912426" lat="35.4463402" lon="139.3755433">
30 <tag k="amenity" v="place_of_worship"/>
31 <tag k="name" v="稲荷山 安養院"/>
32 <tag k="religion" v="shinto"/>
33 </node>
34 <node id="663578611" lat="35.4377708" lon="139.4085008"/>
35 <node id="663578612" lat="35.4377494" lon="139.4086266"/>
36 <node id="663578613" lat="35.4378636" lon="139.4085247"/>
37 <node id="663578614" lat="35.4378421" lon="139.4086505"/>
38 <node id="676016142" lat="35.3933806" lon="139.4620820">
39 <tag k="amenity" v="place_of_worship"/>
40 <tag k="name" v="円行八幡宮"/>
41 <tag k="name:ja_rm" v="Engyō Hachimangū"/>
42 <tag k="religion" v="shinto"/>
43 </node>
44
45 <!-- ... -->
46 <way id="52061810">
47 <nd ref="663578613"/>
48 <nd ref="663578614"/>
49 <nd ref="663578612"/>
50 <nd ref="663578611"/>
51 <nd ref="663578613"/>
52 <tag k="amenity" v="place_of_worship"/>
53 <tag k="name" v="五社神社"/>
54 <tag k="religion" v="shinto"/>
55 <tag k="wheelchair" v="limited"/>
56 </way>
57 <way id="55910762">
58 <nd ref="702017371"/>
59 <nd ref="702017374"/>
60 <nd ref="702017376"/>
61 <nd ref="702017379"/>
62 <nd ref="702017381"/>
63 <nd ref="702017371"/>
64 <tag k="amenity" v="place_of_worship"/>
65 <tag k="name" v="鈴鹿明神社"/>
66 <tag k="note" v="estimated"/>
67 <tag k="religion" v="shinto"/>
68 </way>
69 <!-- ... -->
70
71 </osm>
「駅」のほうも似たようなもん：
1 <?xml version="1.0" encoding="UTF-8"?>
2 <osm version="0.6" generator="Overpass API">
3 <note>The data included in this document is from www.openstreetmap.org. The data is made available under ODbL.</note>
4 <meta osm_base="2015-10-10T18:50:01Z"/>
5
6 <node id="257362961" lat="35.9401590" lon="139.6097084"/>
7 <node id="257362962" lat="35.9401166" lon="139.6093643"/>
8 <node id="257362963" lat="35.9403205" lon="139.6093125"/>
9 <node id="257362964" lat="35.9403665" lon="139.6096631"/>
10 <node id="264813350" lat="35.6446804" lon="139.5014334"/>
11 <node id="428413445" lat="35.6256778" lon="139.7243784"/>
12 <!-- ... -->
13 <node id="1730435566" lat="35.6786902" lon="139.6934965">
14 <tag k="amenity" v="vending_machine"/>
15 <tag k="vending" v="public_transport_tickets"/>
16 </node>
17 <node id="1730435569" lat="35.6786924" lon="139.6936946"/>
18 <node id="1730435575" lat="35.6787152" lon="139.6935148">
19 <tag k="access" v="yes"/>
20 <tag k="entrance" v="main"/>
21 </node>
22 <node id="1730435578" lat="35.6789679" lon="139.6938729"/>
23 <!-- ... -->
24 <node id="2217952046" lat="35.6616200" lon="139.6671654">
25 <tag k="access" v="yes"/>
26 <tag k="entrance" v="main"/>
27 <tag k="name" v="下北沢駅入口"/>
28 <tag k="source" v="image,2013-03-23;Bing"/>
29 </node>
30 <node id="2237786589" lat="35.7060898" lon="139.6837092"/>
31 <!-- ... -->
32 <way id="23753050">
33 <nd ref="257362961"/>
34 <nd ref="257362962"/>
35 <nd ref="1400371223"/>
36 <nd ref="257362963"/>
37 <nd ref="257362964"/>
38 <nd ref="257362961"/>
39 <tag k="building" v="train_station"/>
40 <tag k="name:ja" v="宮原駅"/>
41 <tag k="note" v="estimated"/>
42 <tag k="source" v="Bing,2007-04"/>
43 </way>
44 <!-- ... -->
45 <relation id="3067866">
46 <member type="way" ref="155588179" role="outer"/>
47 <member type="way" ref="228930110" role="inner"/>
48 <member type="way" ref="228930111" role="inner"/>
49 <member type="way" ref="228930112" role="inner"/>
50 <tag k="building" v="train_station"/>
51 <tag k="type" v="multipolygon"/>
52 </relation>
53 <relation id="3808609">
54 <member type="way" ref="98964842" role="outer"/>
55 <member type="way" ref="286088576" role="inner"/>
56 <member type="way" ref="27223744" role="inner"/>
57 <tag k="building" v="train_station"/>
58 <tag k="building:levels" v="3"/>
59 <tag k="name" v="大宮駅"/>
60 <tag k="type" v="multipolygon"/>
61 </relation>
62 <!-- ... -->
63 <relation id="4856156">
64 <member type="way" ref="342276644" role="outer"/>
65 <member type="way" ref="342276673" role="outer"/>
66 <member type="way" ref="342276629" role="outer"/>
67 <member type="way" ref="342276670" role="outer"/>
68 <member type="way" ref="342276622" role="outer"/>
69 <member type="way" ref="342276645" role="outer"/>
70 <member type="way" ref="342276648" role="outer"/>
71 <member type="way" ref="342276672" role="outer"/>
72 <member type="way" ref="342276633" role="outer"/>
73 <member type="way" ref="342276639" role="outer"/>
74 <member type="way" ref="342276618" role="outer"/>
75 <member type="way" ref="342276667" role="outer"/>
76 <member type="way" ref="342276630" role="outer"/>
77 <member type="way" ref="342276657" role="outer"/>
78 <member type="way" ref="342276678" role="outer"/>
79 <member type="way" ref="342276631" role="outer"/>
80 <member type="way" ref="342276666" role="outer"/>
81 <member type="way" ref="342276638" role="outer"/>
82 <member type="way" ref="342276647" role="outer"/>
83 <member type="way" ref="342276609" role="outer"/>
84 <member type="way" ref="342276607" role="outer"/>
85 <member type="way" ref="342276660" role="outer"/>
86 <member type="way" ref="342276612" role="outer"/>
87 <member type="way" ref="342276650" role="outer"/>
88 <member type="way" ref="342276621" role="outer"/>
89 <member type="way" ref="342276668" role="outer"/>
90 <member type="way" ref="342276643" role="outer"/>
91 <member type="way" ref="342276654" role="outer"/>
92 <member type="way" ref="342276626" role="outer"/>
93 <member type="way" ref="145407027" role="outer"/>
94 <member type="way" ref="342276625" role="outer"/>
95 <member type="way" ref="342276653" role="outer"/>
96 <member type="way" ref="342276642" role="outer"/>
97 <member type="way" ref="342276655" role="outer"/>
98 <member type="way" ref="342276632" role="outer"/>
99 <tag k="addr:block_number" v="9"/>
100 <tag k="addr:city" v="千代田区"/>
101 <tag k="addr:country" v="JP"/>
102 <tag k="addr:full" v="東京都千代田区丸の内1-9-1"/>
103 <tag k="addr:housenumber" v="1"/>
104 <tag k="addr:neighbourhood" v="1丁目"/>
105 <tag k="addr:postcode" v="100-0005"/>
106 <tag k="addr:province" v="東京都"/>
107 <tag k="addr:quarter" v="丸の内"/>
108 <tag k="alt_name:en" v="Tokyo Station (Marunouchi)"/>
109 <tag k="alt_name:ja" v="東京駅 (丸の内)"/>
110 <tag k="building" v="train_station"/>
111 <tag k="name" v="東京駅"/>
112 <tag k="name:en" v="Tokyo Station"/>
113 <tag k="name:ja" v="東京駅"/>
114 <tag k="type" v="multipolygon"/>
115 </relation>
116
117 </osm>
ほとんど OSM の構造、わかったでしょ? node, way, relation。node で「点」、way や relation でポリゴンを作ってるというわけ。建造物の細かな構造もこうやって管理されてるわけ。
さーて、持ってきたはいいけど、どう使おうかなぁ…。Postgis に突っ込んでしまうのが一番「王道」ではあるんだけれど、素のままもしくは独自管理の何か(libspatialindex とか使って)という手もないではないしな。とはいっても、そこそこ複雑な構造ではあるから、独自するなら「ちゃんとプログラムしないといけない」しな。だいたいにして、ワタシの今の目的だと、「何か＝中心点」で十分なんだけど、そう使うにはポリゴンやらからわざわざ抽出せねばならん。
うーん、ゆっくり考えよっと。とりあえずはテキストエディタ(わたしの場合は emacs)で開いて検索する、というだけでも、「地理院地図WEBとか Google Map でチマチマ調べ」るよりはずっと快適なんだよね。だってさ、「地図上の場所でプロットされているの図」が欲しいんじゃなくて、緯度経度だけ知りたいんだから。GUI は全くいらんのです。というかその地図 WEB GUI が激重だからこそ、なかなかにストレスなわけで。テキストエディタでの検索ならさ、ネットワーク不調でも問題ないわけね。しばらくはそういう使い方だわな、ワタシの場合。
ところで、「WikiPedia いんすぱいあど」には厳重に注意ね。ライセンス問題は、「皆が行儀良ければ問題がない」ことになっているけれど、WikiPedia がそうではないのは誰でも知っていること(なにせ精神病的なパペットが暗躍する場所だ)。OpenStreetMap のデータにそんなものが紛れ込んでいるはずがない、なんて考えるのはお人好し過ぎです。
あぁ、あとねぇ、今のうちにいっておくと、「bbox が重複しないようにデータ取得」しても、結果的にデータ重複します。そりゃそうなのだ。bbox とのインターセクションで持ってくるんだから、ポリゴンは両方の bbox でヒットしてしまうのは当たり前。でな、これが osm2pgsql で問題起こします。というか、検証した何年か前(多分2年前)は問題起こしてた。そのときは手作業で(自作スクリプトで)重複を取り除く加工をしてから osm2pgsql にかける、という手間のかかることをしてた。ま、こんな情報でも誰かの役に立つかもしれんから、一応。
2021-05-24追記:
6年経過しての追記なので「さぞ進化した追記なのであろう」って期待はせぬべし。ワタシもまさに6年ぶりにやってみた、てだけなのだから。
xapi は6年前と変わらず使える。今回の追記で改めて書かねばと思ったポイントは二つの観点から。
一つ目。6年前は「制約」として
- 「bbox queries have a maximum of 10 square degrees」という明示的な制約
- 結果セットが大きくなると Bad Request となる暗黙の制約
の二つを説明したけれど、改めてやってみると実はもう一つあった。連続で短い時間内にクエリをかけると「429: Too Many Requests」を返す。6年前時点でも起こってたのかなぁ? まぁいわゆる DOS 攻撃に対する防御策だとは思うんだけれど、とにかくこれのおかげで、わしらクライアントとしては「お静かに」使う必要がある。
そんななので、たとえば(139.5, 35.5)～(140.5, 36.0)に対して [building=*] を 1/24度刻みで問い合わせたい、みたいなヘビィなことをするには、こんな感じ:
汎用スクリプトではなくて直値埋め込みだよ
1 # -*- coding: utf-8 -*-
2
3 import os
4 import math
5 import time
6 from urllib.request import quote
7
8 lon_step = 1.0/24
9 lat_step = 1.0/24
10
11 bottom_left = (139.5, 35.5)
12 top_right = (140.5, 36.0)
13
14 _BASEURI = "http://www.overpass-api.de/api/xapi?"
15 _downloadlist = []
16 for x in range(int(math.ceil((top_right[0] - bottom_left[0]) / lon_step))):
17 for y in range(int(math.ceil((top_right[1] - bottom_left[1]) / lat_step))):
18 min_long = bottom_left[0] + x * lon_step
19 min_lat = bottom_left[1] + y * lat_step
20 max_long = bottom_left[0] + (x + 1) * lon_step
21 max_lat = bottom_left[1] + (y + 1) * lat_step
22
23 # use XAPI
24 _downloadlist.append(
25 _BASEURI + "*[building=*][bbox=%.4f,%.4f,%.4f,%.4f]" % (
26 min_long, min_lat, max_long, max_lat))
27 import subprocess
28 time.sleep(10)
29 while len(_downloadlist):
30 uri = _downloadlist[0]
31 t = uri.rpartition("/")[-1].replace("*", "%2A").replace("?", "@")
32 suc = True
33 if os.path.exists(t):
34 print(t)
35 else:
36 try:
37 subprocess.check_call(["wget", "-c", uri])
38 time.sleep(int(60*4))
39 except Exception as e:
40 print(e)
41 time.sleep(int(60*5))
42 suc = False
43 if suc:
44 _downloadlist.pop(0)
みたいな感じに「そーっと」お取り寄せるしかない。むろんこれはとんでもない長時間かかる。「suc」が全部 True になる超強運だとしても、2880分(48時間)かかる計算。まぁこの sleep があってもなくてもどのみち時間がかかるクエリではあるんで、許容不可能ってほどでもないのだけれどもね。
もう一つの観点は、xapi のもう少し「狙い撃ちしたい」場合の話。上でやったように「いろいろな寺社仏閣データが欲しい」てことでなくて、例えば「東京ドームのポリゴンが欲しい」みたいに狙い撃ちしたいケースね。こういう場合、位置は正確にわかってるわけだから、スクリプトはこんなんでいいと思う:
特に凝ったことはしない
1 # -*- coding: utf-8 -*-
2 import subprocess
3
4
5 _BASEURI = "http://www.overpass-api.de/api/xapi?"
6
7
8 def main(args):
9 # bbox を、与えられた中心緯度経度と「幅」から求めるのだが、
10 # 精密な地理演算を pygrib や geographiclib で行う必要までは
11 # ないと思う。ここでは緯度方向での近似値「1秒≒30m」を使う。
12 clon, clat, width = args.clon, args.clat, args.width
13 degpermetre = 1 / (30.0 * 3600)
14 ext = (width / 2) * degpermetre
15 bbox = (clon - ext, clat - ext, clon + ext, clat + ext)
16 fq = ""
17 if args.feature_query:
18 fq = "[{}]".format(args.feature_query)
19 uri = _BASEURI + "*{fq}[bbox={bb}]".format(
20 fq=fq, bb=",".join(["{:.4f}".format(p) for p in bbox]))
21 subprocess.check_call(["wget", "-c", uri])
22
23
24 if __name__ == '__main__':
25 import argparse
26 ap = argparse.ArgumentParser()
27 ap.add_argument("clon", type=float)
28 ap.add_argument("clat", type=float)
29 ap.add_argument("width", type=float) # in metre
30 ap.add_argument("--feature_query", help="for example: 'religion=*'")
31 args = ap.parse_args()
32 main(args)
今問題にしたいのは本当に東京ドームを狙い撃ちで取得したい場合に「feature_query」として何を使えばいいのか、て話で、これがさ、まぁ Map_features を読んで「これだろう」と思ってもうまくいかんのよね。正解を偶然でもヒットさせてみて初めて正解がわかる:
抜粋
1 <tag k="building" v="yes"/>
2 <tag k="leisure" v="stadium"/>
3 <tag k="lit" v="yes"/>
4 <tag k="name" v="東京ドーム"/>
5 <tag k="name:en" v="Tokyo Dome"/>
6 <tag k="name:es" v="Tokyo Dome"/>
7 <tag k="name:ja" v="東京ドーム"/>
8 <tag k="name:ru" v="Стадион Токио Доум"/>
9 <tag k="name:zh" v="東京巨蛋"/>
10 <tag k="sport" v="baseball;multi"/>
11 <tag k="wikidata" v="Q733748"/>
12 <tag k="wikipedia" v="ja:東京ドーム"/>
「leisure=stadium」なんて想像もつかんし、「building=yes」ってなんだよ、と。(building=* なんてすさまじいデータ量になるが、それをやってみればめでたくヒットする、てことだよ。)
この正解を探り当ててみて初めて、以下が最も手早いのだとわかる:
上のスクリプトが zzz.py という名前だとして
1 [me@host: ~]$ py -3 zzz.py 139.75130 35.70516 2000 --feature_query='name:en=Tokyo Dome'
正解がわかってから改めて正解を取り寄せることは意味があるのか、てのは、まぁあるにはある。だって:
上のスクリプトが zzz.py という名前だとして
1 [me@host: ~]$ py -3 zzz.py 139.75130 35.70516 2000 --feature_query='building=*'
だと、東京ドーム以外の「building」が大量に入ってくるでしょ。ゲット出来た xml から目的の東京ドームのものを抽出する手間だって結構バカにならん。
ところで、最初のも先のも「wget を呼び出しちまえ」という雑な決断が迷惑な場合もある。とするなら例えばこんな:
ローカルキャッシュしてくれるサードパーティライブラリを使うともっとより良いが、ここでは依存物なしで
1 # -*- coding: utf-8 -*-
2 # python 3.x 用。いい加減 2.7 は捨てようね、特殊な事情が
3 # ある人以外は。
4 import io
5 import os
6 import hashlib
7 import tempfile
8 import urllib.request
9
10
11 _BASEURI = "http://www.overpass-api.de/api/xapi?"
12
13
14 def main(args):
15 # bbox を、与えられた中心緯度経度と「幅」から求めるのだが、
16 # 精密な地理演算を pygrib や geographiclib で行う必要までは
17 # ないと思う。ここでは緯度方向での近似値「1秒≒30m」を使う。
18 clon, clat, width = args.clon, args.clat, args.width
19 degpermetre = 1 / (30.0 * 3600)
20 ext = (width / 2) * degpermetre
21 bbox = (clon - ext, clat - ext, clon + ext, clat + ext)
22 fq = ""
23 if args.feature_query:
24 fq = "[{}]".format(urllib.request.quote(args.feature_query))
25 uri = _BASEURI + "*{fq}[bbox={bb}]".format(
26 fq=fq, bb=",".join(["{:.4f}".format(p) for p in bbox]))
27 resf = os.path.join(
28 tempfile.gettempdir(),
29 "xapi_" + hashlib.md5(uri.encode("utf-8")).hexdigest())
30 if not os.path.exists(resf):
31 cont, msg = urllib.request.urlretrieve(uri)
32 os.rename(cont, resf)
33 print(io.open(resf, encoding="utf-8").read())
34
35
36
37 if __name__ == '__main__':
38 import argparse
39 ap = argparse.ArgumentParser()
40 ap.add_argument("clon", type=float)
41 ap.add_argument("clat", type=float)
42 ap.add_argument("width", type=float) # in metre
43 ap.add_argument("--feature_query", help="for example: 'religion=*'")
44 args = ap.parse_args()
45 main(args)
こんな具合。
ここまで出来れば、「得られた xml (oml) から欲しいポリゴンをば」の初めの一歩れる:
1 # -*- coding: utf-8 -*-
2 # python 3.x 用。いい加減 2.7 は捨てようね、特殊な事情が
3 # ある人以外は。
4 import io
5 import os
6 import hashlib
7 import tempfile
8 import urllib.request
9 import xml.etree.ElementTree as ET
10
11
12 _BASEURI = "http://www.overpass-api.de/api/xapi?"
13
14
15 def main(args):
16 # bbox を、与えられた中心緯度経度と「幅」から求めるのだが、
17 # 精密な地理演算を pygrib や geographiclib で行う必要までは
18 # ないと思う。ここでは緯度方向での近似値「1秒≒30m」を使う。
19 # (軽度方向で解釈する場合より大きな bbox になる。)
20 clon, clat, width = args.clon, args.clat, args.width
21 degpermetre = 1 / (30.0 * 3600)
22 ext = (width / 2) * degpermetre
23 bbox = (clon - ext, clat - ext, clon + ext, clat + ext)
24 fq = ""
25 if args.feature_query:
26 fq = "[{}]".format(urllib.request.quote(args.feature_query))
27 uri = _BASEURI + "*{fq}[bbox={bb}]".format(
28 fq=fq, bb=",".join(["{:.4f}".format(p) for p in bbox]))
29 resf = os.path.join(
30 tempfile.gettempdir(),
31 "xapi_" + hashlib.md5(uri.encode("utf-8")).hexdigest())
32 if not os.path.exists(resf):
33 cont, msg = urllib.request.urlretrieve(uri)
34 os.rename(cont, resf)
35 #print(io.open(resf, encoding="utf-8").read())
36 parsed = ET.parse(resf)
37 nodes = {} # id: (lon, lat)
38 ways = {}
39 for tle in parsed.getroot():
40 if tle.tag == "node":
41 id = tle.attrib["id"]
42 lon = float(tle.attrib["lon"])
43 lat = float(tle.attrib["lat"])
44 nodes[id] = (lon, lat)
45 for tle in parsed.getroot():
46 if tle.tag == "way":
47 id = tle.attrib["id"]
48 ways[id] = {"pts": []}
49 for che in tle:
50 if che.tag == "nd":
51 refid = che.attrib["ref"]
52 ways[id]["pts"].append(nodes[refid])
53 elif che.tag == "tag":
54 ways[id][che.attrib["k"]] = che.attrib["v"]
55 print(ways)
56
57
58
59 if __name__ == '__main__':
60 import argparse
61 ap = argparse.ArgumentParser()
62 ap.add_argument("clon", type=float)
63 ap.add_argument("clat", type=float)
64 ap.add_argument("width", type=float) # in metre
65 ap.add_argument("--feature_query", help="for example: 'religion=*'")
66 args = ap.parse_args()
67 main(args)
無論、ここまで出来てるなら、おそらく以下がほとんどゴールであろう:
1 # -*- coding: utf-8 -*-
2 # python 3.x 用。いい加減 2.7 は捨てようね、特殊な事情が
3 # ある人以外は。
4 import io
5 import os
6 import hashlib
7 import tempfile
8 import urllib.request
9 import xml.etree.ElementTree as ET
10
11 from shapely.geometry import asPolygon
12
13
14 _BASEURI = "http://www.overpass-api.de/api/xapi?"
15
16
17 def main(args):
18 # bbox を、与えられた中心緯度経度と「幅」から求めるのだが、
19 # 精密な地理演算を pygrib や geographiclib で行う必要までは
20 # ないと思う。ここでは緯度方向での近似値「1秒≒30m」を使う。
21 # (軽度方向で解釈する場合より大きな bbox になる。)
22 clon, clat, width = args.clon, args.clat, args.width
23 degpermetre = 1 / (30.0 * 3600)
24 ext = (width / 2) * degpermetre
25 bbox = (clon - ext, clat - ext, clon + ext, clat + ext)
26 fq = ""
27 if args.feature_query:
28 fq = "[{}]".format(urllib.request.quote(args.feature_query))
29 uri = _BASEURI + "*{fq}[bbox={bb}]".format(
30 fq=fq, bb=",".join(["{:.4f}".format(p) for p in bbox]))
31 resf = os.path.join(
32 tempfile.gettempdir(),
33 "xapi_" + hashlib.md5(uri.encode("utf-8")).hexdigest())
34 if not os.path.exists(resf):
35 cont, msg = urllib.request.urlretrieve(uri)
36 os.rename(cont, resf)
37 #print(io.open(resf, encoding="utf-8").read())
38 parsed = ET.parse(resf)
39 nodes = {} # id: (lon, lat)
40 ways = {}
41 for tle in parsed.getroot():
42 if tle.tag == "node":
43 id = tle.attrib["id"]
44 lon = float(tle.attrib["lon"])
45 lat = float(tle.attrib["lat"])
46 nodes[id] = (lon, lat)
47 for tle in parsed.getroot():
48 if tle.tag == "way":
49 id = tle.attrib["id"]
50 ways[id] = {}
51 pts = []
52 for che in tle:
53 if che.tag == "nd":
54 refid = che.attrib["ref"]
55 pts.append(nodes[refid])
56 elif che.tag == "tag":
57 ways[id][che.attrib["k"]] = che.attrib["v"]
58 ways[id]["geometry"] = asPolygon(pts)
59 #print(ways)
60
61
62 if __name__ == '__main__':
63 import argparse
64 ap = argparse.ArgumentParser()
65 ap.add_argument("clon", type=float)
66 ap.add_argument("clat", type=float)
67 ap.add_argument("width", type=float) # in metre
68 ap.add_argument("--feature_query", help="for example: 'religion=*'")
69 args = ap.parse_args()
70 main(args)
ところで全然関係ないんだけど、pyproj 2.3.0 以降の機能を使うと、m2 での面積が簡単に求まる:
1 # -*- coding: utf-8 -*-
2 import pyproj
3 import shapely.wkt
4
5
6 # This polygon data was extracted from www.openstreetmap.org,
7 # the data is made available under ODbL.
8 _TOKYODOME_POLYGON = shapely.wkt.loads(
9 "POLYGON ((139.7513598 35.7064612, \
10 139.7515235 35.7064901, \
11 139.7517468 35.7064952, \
12 139.7521794 35.7064623, \
13 139.7523926 35.7064262, \
14 139.7525588 35.7063778, \
15 139.7526593 35.7063386, \
16 139.7527745 35.706248, \
17 139.7528709 35.7061439, \
18 139.7529162 35.7061581, \
19 139.7530749 35.7060648, \
20 139.753143 35.7058579, \
21 139.7531508 35.705673, \
22 139.7531475 35.7054376, \
23 139.7531262 35.705252, \
24 139.7530789 35.7050554, \
25 139.7529928 35.7048943, \
26 139.7528655 35.7047525, \
27 139.7526958 35.7046419, \
28 139.7525916 35.7045943, \
29 139.7524667 35.7045587, \
30 139.7523944 35.7045434, \
31 139.7523461 35.7045366, \
32 139.7522891 35.704533, \
33 139.7522608 35.7045617, \
34 139.7521553 35.7045742, \
35 139.7518103 35.7046149, \
36 139.7515323 35.7046478, \
37 139.7514454 35.70461, \
38 139.7513705 35.7046338, \
39 139.7512855 35.7046646, \
40 139.7510745 35.7047411, \
41 139.7509984 35.7048903, \
42 139.7510216 35.7049229, \
43 139.7509335 35.7050014, \
44 139.750846 35.7051343, \
45 139.7508066 35.7052455, \
46 139.7507985 35.7052842, \
47 139.7507889 35.70533, \
48 139.75078 35.7054681, \
49 139.7508092 35.7057256, \
50 139.7508613 35.705993, \
51 139.750925 35.7061552, \
52 139.7509995 35.70625, \
53 139.7511073 35.7063438, \
54 139.751219 35.7064118, \
55 139.7513598 35.7064612))")
56
57 wgs84 = pyproj.Geod(ellps="WGS84")
58 print(abs(wgs84.geometry_area_perimeter(_TOKYODOME_POLYGON)[0]))
これによれば、「東京ドームの面積は 38,687 m2」ということになるのだけれど、wikipedia 情報の「グラウンド面積：13,000 m2」「建築面積は46,755 m2」のどちらとも一致しないのよね。なので「ワタシのことどれくらい愛してるの」との問いに「東京ドーム３つ分」などと答えたい場合にとても困る。テレビとかでこの計算してるときって、どの値を使ってるんだろうね? (一応正方形換算での一辺が、13,000 m2 だと 114m、38,687 m2 だと 197m、46,755 m2 だと 216m となり、まぁ微々たる差だとも言える、けれども。)
2021-05-31追記:
少し色々お取り寄せてみようとやってたのだが、ちょっと shapely の遅延処理によって、扱いにくいエラーが起こるんだよね。おかしなポリゴンとかを、その場では検出出来なくて、時限爆弾的に例外となる。これだとちょっと困るので、例えばこうかなと:
1 # -*- coding: utf-8 -*-
2 # python 3.x 用。いい加減 2.7 は捨てようね、特殊な事情が
3 # ある人以外は。
4 import io
5 import os
6 import hashlib
7 import tempfile
8 import urllib.request
9 import xml.etree.ElementTree as ET
10
11 from shapely.geometry import asPolygon, asLineString
12
13
14 _BASEURI = "http://www.overpass-api.de/api/xapi?"
15
16
17 def main(args):
18 # bbox を、与えられた中心緯度経度と「幅」から求めるのだが、
19 # 精密な地理演算を pygrib や geographiclib で行う必要までは
20 # ないと思う。ここでは緯度方向での近似値「1秒≒30m」を使う。
21 # (軽度方向で解釈する場合より大きな bbox になる。)
22 clon, clat, width = args.clon, args.clat, args.width
23 degpermetre = 1 / (30.0 * 3600)
24 ext = (width / 2) * degpermetre
25 bbox = (clon - ext, clat - ext, clon + ext, clat + ext)
26 fq = ""
27 if args.feature_query:
28 fq = "[{}]".format(urllib.request.quote(args.feature_query))
29 uri = _BASEURI + "*{fq}[bbox={bb}]".format(
30 fq=fq, bb=",".join(["{:.4f}".format(p) for p in bbox]))
31 resf = os.path.join(
32 tempfile.gettempdir(),
33 "xapi_" + hashlib.md5(uri.encode("utf-8")).hexdigest())
34 if not os.path.exists(resf):
35 cont, msg = urllib.request.urlretrieve(uri)
36 os.rename(cont, resf)
37 #print(io.open(resf, encoding="utf-8").read())
38 parsed = ET.parse(resf)
39 nodes = {} # id: (lon, lat)
40 ways = {}
41 for tle in parsed.getroot():
42 if tle.tag == "node":
43 id = tle.attrib["id"]
44 lon = float(tle.attrib["lon"])
45 lat = float(tle.attrib["lat"])
46 nodes[id] = (lon, lat)
47 for tle in parsed.getroot():
48 if tle.tag == "way":
49 thisway = {}
50 pts = []
51 for che in tle:
52 if che.tag == "nd":
53 refid = che.attrib["ref"]
54 pts.append(nodes[refid])
55 elif che.tag == "tag":
56 thisway[che.attrib["k"]] = che.attrib["v"]
57 if pts[0] == pts[-1]:
58 poly = asPolygon(pts)
59 else:
60 poly = asLineString(pts)
61 try: # check if valid polygon
62 poly.type # cause Exception from shapely if it is not valid
63 thisway["geometry"] = poly
64 ways[tle.attrib["id"]] = thisway
65 except ValueError:
66 pass
67 def _str(s):
68 import sys
69 if s:
70 return s.encode(
71 sys.getfilesystemencoding(), errors="replace").decode(
72 sys.getfilesystemencoding(), errors="replace")
73 for k in ways.keys():
74 print(k, _str(ways[k].get("name")), ways[k]["geometry"])
75
76
77 if __name__ == '__main__':
78 import argparse
79 ap = argparse.ArgumentParser()
80 ap.add_argument("clon", type=float)
81 ap.add_argument("clat", type=float)
82 ap.add_argument("width", type=float) # in metre
83 ap.add_argument("--feature_query", help="for example: 'religion=*'")
84 args = ap.parse_args()
85 main(args)
なおこのコードは、このネタ内で実際に使ってみている。
2021-06-04追記:
最初にこのネタを書いた際に気付いていたのかどうかとか、前回の追記時点で気付いていたのかとかは聞かないで。自分でもわからない。が、とにかく「xapi のほかに、overpass API がある」。ここ最近のこうしたものには付き物の「オンラインお試しぁー」がある。ご自分でもお試しぁーてみた:
して、この overpass API への python ラッパーが、案の定ある。いくつかある。overpass-api-python-wrapper、overpy、osmpy。
API が over http のよくあるやつなわけなので、python ラッパーは結局はどれも同じ構造になってて、まぁぶっちゃけ「与えられたリクエストを素直にブン投げる」だけなのだが、だとするなら出来るだけ依存物が少ないものが良い、として選ぶなら、osmpy が真っ先に候補から外れる。なんで panda 依存なんだ。意味がわからない。ソースをぱっと見た感じだと、overpy で十分そうだったので、そちらを入れてみた。とするならば:
1 # -*- coding: utf-8 -*-
2 from shapely.geometry import asPolygon
3 import overpy
4 #
5 api = overpy.Overpass()
6 #
7 result = api.query("""
8 [timeout:60][bbox:40.7647452, -74.0066133, 40.8367796, -73.9138199];
9 way["name"="Central Park"];
10 out;
11 >; /* recurse down */
12 out;
13 """)
14 pts = []
15 for way in result.ways: # この結果の場合は way は一つのみ。
16 #print(way)
17 for node in way.nodes:
18 #print(node)
19 pts.append((node.lon, node.lat))
20 poly = asPolygon(pts)
21 print(poly)
1 POLYGON ((-73.9577049 40.8003107, -73.95768270000001 40.8003048, -73.9576536 40.8002923, -73.9559873 40.7995974, -73.95559249999999 40.7994309, -73.95512600000001 40.799237, -73.95270429999999 40.7982164, -73.952353 40.7980661, -73.95103829999999 40.79751, -73.95097389999999 40.7974831, -73.9502608 40.7971833, -73.94978829999999 40.7969865, -73.9497673 40.7969775, -73.9497623 40.796969, -73.9497571 40.7969586, -73.94975909999999 40.7969467, -73.9497626 40.7969378, -73.9497731 40.7969178, -73.9497839 40.7968924, -73.94979379999999 40.7968656, -73.9498019 40.7968361, -73.9498014 40.7967913, -73.9497688 40.7967264, -73.9497575 40.7967123, -73.9497422 40.7966959, -73.94972439999999 40.7966815, -73.9497027 40.7966668, -73.9496802 40.7966549, -73.9496618 40.7966473, -73.94964109999999 40.7966385, -73.94962990000001 40.7966288, -73.9496189 40.7966161, -73.9496116 40.796606, -73.9496061 40.7965935, -73.9496631 40.7965179, -73.9498487 40.7962602, -73.9499364 40.7961426, -73.95030319999999 40.7956392, -73.95074940000001 40.7950221, -73.9508105 40.7949343, -73.9512675 40.7943243, -73.95168630000001 40.7937493, -73.9518348 40.7935351, -73.9522388 40.7929767, -73.9530523 40.7918681, -73.9531393 40.7917411, -73.95352560000001 40.7912074, -73.9552971 40.7888635, -73.9554172 40.7886997, -73.9558772 40.788053, -73.95757330000001 40.7857067, -73.9586307 40.7843102, -73.9588854 40.7839459, -73.9605506 40.781666, -73.9610479 40.78096, -73.9611242 40.7808832, -73.96148100000001 40.7804146, -73.9615543 40.7803339, -73.96146400000001 40.7802617, -73.96179789999999 40.7798051, -73.9618727 40.7797018, -73.9622595 40.7791667, -73.9623347 40.7790627, -73.96271950000001 40.7785304, -73.9627898 40.7784331, -73.9631581 40.7779237, -73.9633211 40.7776951, -73.96332580000001 40.7776925, -73.9633312 40.7776914, -73.963339 40.7776932, -73.9633886 40.7777161, -73.9634089 40.7776896, -73.96347369999999 40.7775944, -73.963494 40.7775657, -73.96351110000001 40.7775416, -73.9637321 40.777239, -73.9640526 40.776801, -73.9656177 40.7746546, -73.9662265 40.7738129, -73.968529 40.770653, -73.9699119 40.7687628, -73.96997639999999 40.7686749, -73.9704097 40.768135, -73.9704918 40.7680224, -73.9708816 40.7674883, -73.97095299999999 40.7673907, -73.9715326 40.7665973, -73.9722517 40.7656121, -73.9724768 40.765304, -73.9725696 40.7651788, -73.9726515 40.7650644, -73.97303170000001 40.7652192, -73.9731694 40.7652316, -73.9732963 40.7651947, -73.9736176 40.7647452, -73.9809915 40.7678561, -73.9811787 40.7678877, -73.98119579999999 40.7678954, -73.9812094 40.7679176, -73.9812166 40.7679435, -73.9812024 40.7679925, -73.98119579999999 40.7680234, -73.9811947 40.7680529, -73.98119579999999 40.7680759, -73.9811983 40.7681234, -73.98120539999999 40.7681647, -73.98122499999999 40.7682236, -73.9812424 40.7682578, -73.9812666 40.7682979, -73.98129539999999 40.7683388, -73.9813532 40.7684005, -73.9814075 40.7684558, -73.98137850000001 40.7684879, -73.98129710000001 40.7685782, -73.98126379999999 40.7686118, -73.9809437 40.7690783, -73.9808882 40.7691598, -73.98083080000001 40.7692432, -73.9803189 40.7699787, -73.9797918 40.7706802, -73.9797693 40.7707071, -73.97921119999999 40.7714634, -73.97906039999999 40.7716807, -73.9789928 40.7717669, -73.97864939999999 40.772256, -73.9781445 40.7729298, -73.97807330000001 40.7730567, -73.97719619999999 40.7742602, -73.97717419999999 40.7742872, -73.9762831 40.7755038, -73.97586870000001 40.7760769, -73.97574059999999 40.7762343, -73.9748348 40.7774861, -73.9741143 40.7784888, -73.97371200000001 40.7790432, -73.9734706 40.7793708, -73.97340149999999 40.7794669, -73.973142 40.7798279, -73.97265470000001 40.7804911, -73.9725365 40.7806502, -73.9716566 40.781843, -73.9715029 40.7820647, -73.9714695 40.7821005, -73.9704966 40.7834445, -73.9702619 40.783752, -73.9700242 40.7840726, -73.96970760000001 40.7845135, -73.96965899999999 40.784584, -73.96925640000001 40.7851403, -73.96922290000001 40.7852014, -73.96911299999999 40.7853303, -73.9690437 40.7854255, -73.96823980000001 40.7865399, -73.9673372 40.7877836, -73.9672728 40.7878736, -73.9661212 40.7894415, -73.96565099999999 40.7900902, -73.96460449999999 40.7915356, -73.96451570000001 40.7916521, -73.9639878 40.7923871, -73.9626287 40.7942007, -73.9613682 40.795929, -73.9613456 40.7959607, -73.95996169999999 40.7978697, -73.95977600000001 40.7981253, -73.95906069999999 40.7991133, -73.9589862 40.7992085, -73.9589604 40.7992378, -73.9589165 40.7992931, -73.9588184 40.7994129, -73.9587263 40.7995306, -73.9584444 40.7999135, -73.9582594 40.8001649, -73.9582454 40.8001744, -73.9582322 40.8001812, -73.95822250000001 40.8001852, -73.9582116 40.8001889, -73.9582004 40.8001916, -73.9581901 40.8001932, -73.9581702 40.8001953, -73.958139 40.8001996, -73.9581067 40.8002035, -73.9580742 40.8002065, -73.958039 40.8002111, -73.9579531 40.8002352, -73.9579135 40.8002443, -73.9578853 40.8002523, -73.9578609 40.8002612, -73.9578357 40.8002711, -73.9578121 40.8002816, -73.9577602 40.8003077, -73.9577424 40.8003117, -73.95772650000001 40.8003135, -73.9577049 40.8003107))
xapi 経験で喰らった制約は、当然こちらでも同じなのだが、一応 API としてタイムアウトの制御が出来たり、python ラッパーが結構リトライを頑張ったりと、xapi だけ使ってるよりはいいんではないかと思う。
ただし。今一度おさらいをしておく。「全データが必要だ」ということならば本当に全データのアーカイブを丸ごとダウンロード出来る(または大陸別)。本当にそれが必要ならそちらを使う必要がある。現実には無理だと思うけどね、少なくとも個人が自分の PC でどうにか出来るようなデータサイズではない。よほどのモンスターマシンが必要、それをするには。だとしても(何度でも言うが)本当にそれが必要ならそちらを使う必要がある。なんで何度も言うかというと、「xapi であれ、overpass API であれ、これは「個別データをちまちまと少しだけ、狙い撃ちで」」ダウンロードするためのものであって、決してこれを使って「あらゆるデータを取っちゃるぜ」と考えてはいけないし、実際問題として「それは絶対に不可能」。有限時間内に処理を終えられるとは思えない、というかあなたが生きている間に処理が終わることは非常に考えにくい。
さらにいえば、on-the-fly のノリで使うようなものでもない。ここで例にしたような「bbox を小さく絞っての name=”Central Park”」くらいならそういう使い方も出来るけれど、ちょっとでも緩い検索をしようもんなら、すぐに「サーバから結果が返ってこない」だったり「サーバが不平を言ってきたり」する事象にエンカウントすることになる。つまり普通は「一回だけお取り寄せにいって自PCに取り込んでおいて、取り込んだそれを使う」という使い方を主体に考えるべき。
