なぜだか小学校一年生くらいに行きたい国の話題になって、「アメリカと外国っ!」と意気揚々としてた記憶が鮮明である。そしてなぜに「世界全国」という言葉がないのかが今でも気になる。全国なのだからワールドワイドなのではないの?
cartopy ネタとしての直接の (5) の続き、ではないのだけれども。
国土地理院配布の shapefile の最小単位がね、足りない場合もあるわけよ。考え始めてから気付いたんだけれど、「東京都大田区」という単位は「仙台市泉区」という単位よりも「大きい」のよ、同じでないの。して最小単位は「市町村 + 東京都の区」、ゆえに仙台市青葉区のポリゴンは得られない。
まぁ今のワタシが困ってるということじゃないし、将来的にも困る状況を想像しにくいのではあるけれど、皮算用つーか。
「東京２３区は特別だ」というのは、たとえば「２３区全体」と「八王子市」を同列に扱いたい場合がある、みたいな形で現れるが、その場合は別に問題はないわけね:
1 # -*- coding: utf-8 -*-
2 import shapely.geometry as sgeom
3 import shapely.ops as sops
4 from cartopy.io.shapereader import Reader as ShapeReader
5 import matplotlib.pyplot as plt
6 import cartopy.crs as ccrs
7
8
9 def collect(filename="polbnda_jpn"):
10 from itertools import groupby
11 # Source: Geospatial Information Authority of Japan website
12 # (https://www.gsi.go.jp/kankyochiri/gm_japan_e.html)
13 reader = ShapeReader(filename)
14 def keyfun(r):
15 attr = r.attributes
16 nam, laa = attr["nam"], attr["laa"]
17 shiku = laa.rpartition(" ")[-1]
18 if shiku == "Ku":
19 return nam, 0, shiku, "Ku"
20 else:
21 return nam, 1, shiku, laa
22
23 result = sorted(reader.records(), key=keyfun)
24 return groupby(result, keyfun)
25
26
27 def main():
28 geoms = []
29 for k, gr in [(k, list(gr)) for k, gr in collect()
30 if k[0] == "Tokyo To"]:
31 g = sops.cascaded_union([r.geometry for r in gr])
32 geoms.append((k, g))
33 # 描画対象が日本全体の場合の東京描画に困難さはないが、東京のみにスコープしたい
34 # 場合の「東京全部」という領域描画が難しい。日本でこの難しさがあるのは
35 # 東京と沖縄、鹿児島。こと東京に関しては、目的に応じて分割するしかなかろう。
36 # ここでは「最大面積部分 = 陸域部分を抽出」。
37 ugall = list(sops.cascaded_union([g for k, g in geoms])) # a list of POLYGON
38 ugall.sort(key=lambda g: -g.area)
39 bnds = ugall[0].bounds
40 ext = [bnds[0], bnds[2], bnds[1], bnds[3]]
41 #
42 fig = plt.figure()
43 fig.set_size_inches(16.53, 11.69)
44 ax1 = fig.add_subplot(projection=ccrs.PlateCarree())
45 ax1.set_extent(ext, ccrs.PlateCarree())
46 ax1.add_wms(
47 wms='https://ows.terrestris.de/osm/service',
48 layers=["OSM-WMS",])
49
50 ax1.add_geometries(
51 [g for k, g in geoms], ccrs.PlateCarree(),
52 facecolor='#ffff0066',
53 edgecolor='b', linewidth=0.5)
54 import sys, os
55 plt.savefig(os.path.basename(sys.argv[0]) + ".png", bbox_inches="tight")
56
57
58 main()
「足立区」単独で見ることもできるが２３区という一体で見ようと思えば見れる、つまりは小が大を兼ねているケースよね。
けれども「仙台市」は、行政区という意味ではそれ以上細かく分解できないので、「大が小を兼ねてくれない」ことがある。つまりたとえば「仙台市若林区」みたいな単位の境界が欲しいなら、ほかのものに頼るほかはない、てことになる。
StackExchange のを見てた。
OSM から overpass API で(Japan コミュニティにおけるお約束に基づいて)お取り寄せ可能、という情報については、ひとまず overpass API そのものについてはここに追記を書いた。が、Overpass QL は少し理解出来て、それで問い合わせできるようにはなったのだけれど、Overpass XML で書かれてる StackExchange 回答がわからず、overpass turbo にそのまま貼り付けても動かないし、自分なりに Overpass QL で admin_level=”9″とかの記述をしてみても、期待したポリゴンがほとんど取れない。なんか Point とか LineString がたくさん出てきたりはするんだけどね、全然うまくいかない。うーん、後回すか…。
直截的には、 国土交通省配布の GML で頑張る手がある。
GML が大好きな人なんかいるのか? いないことを願ってるよ。こと GML 3 になってからが特にヒドい。GML 2 だった頃はそもそも結構単純さを保っていたし、多くの OSS が頑張ってこれに取り組んでいた。ところが GML 3 になって、「様式美が大好きなお役所」がとてもとても、それはそれはとても頑張って「GML 3 化しました、ステキでしょ」とせっせと乗り換えしちゃったけれど、驚くほど OSS がこれに追従しなかった。こういう状態になったのはおよそ10年前だが、なんと、この状態は今でもほとんど変わってない。なんというかあれよ、結局のところ皆 OSM や GoogleMap に満足しちゃって、そっちにばかりリソースを費やしてるってことよ。
「GML 3 を扱えるお便利インフラ」はあるかもしれんし、ないかもしれん。あるならばそれはそれは「大層な」ライブラリであろう。その「大仰な」ライブラリなしでやるなら、こういうムゴいことになる:
1 # -*- coding: utf-8 -*-
2 import sys
3 import re
4 from xml.etree.ElementTree import ElementTree
5 from shapely.geometry import asPolygon, asLineString
6 import matplotlib.pyplot as plt
7 import cartopy.crs as ccrs
8
9
10 def read_from_mlitN03(fn):
11 """
12 fn: for example "N03-110331_04.xml"
13 """
14 tree = ElementTree()
15 tree.parse(fn)
16 ns = {
17 "ksj": "http://nlftp.mlit.go.jp/ksj/schemas/ksj-app",
18 "jps": "http://www.gsi.go.jp/GIS/jpgis/standardSchemas",
19 "ksjs": "http://nlftp.mlit.go.jp/ksj/schemas/ksj-app-cd",
20 }
21 aa01obj = tree.find("dataset/ksj:object/ksj:AA01/ksj:OBJ", ns)
22 points = {}
23 for p in aa01obj.findall("jps:GM_Point", ns):
24 pid = p.attrib["id"]
25 pos = p[0][0].find("DirectPosition.coordinate")
26 points[pid] = tuple(map(float, pos.text.split(" ")[::-1]))
27 #
28 curves = {}
29 for c in aa01obj.findall("jps:GM_Curve", ns):
30 cid = c.attrib["id"]
31 csg = c.findall("jps:GM_Curve.segment", ns)[0][0]
32 par = csg.findall("jps:GM_LineString.controlPoint", ns)[0][0]
33 curves[cid] = []
34 for pac in par:
35 p = pac[0]
36 if p.tag.endswith("GM_Position.indirect"):
37 ptref = list(p)[0].attrib["idref"]
38 pos = points[ptref]
39 else:
40 pos = list(p)[0]
41 pos = tuple(map(float, pos.text.split(" ")[::-1]))
42 curves[cid].append(pos)
43 #
44 surfas = {}
45 for s in aa01obj.findall("jps:GM_Surface", ns):
46 sid = s.attrib["id"]
47 sug = s.find("jps:GM_Surface.patch", ns)[0]
48 sbg = sug.find("jps:GM_Polygon.boundary", ns)[0]
49 t = sbg.find("jps:GM_SurfaceBoundary.exterior", ns)[0].find(
50 "jps:GM_CompositeCurve.generator", ns)
51 cid = re.sub(r"^_", "", t.attrib["idref"])
52 surfas[sid] = curves[cid]
53 #
54 result = {}
55 for ec01 in aa01obj.findall("ksj:EC01", ns):
56 eid = ec01.attrib["id"]
57 are = surfas[ec01.find("ksj:ARE", ns).attrib["idref"]]
58 result[eid] = dict(
59 prn=ec01.find("ksj:PRN", ns).text,
60 con=ec01.find("ksj:CON", ns).text,
61 cn2=ec01.find("ksj:CN2", ns).text,
62 are=asPolygon(are),
63 )
64 return result
65
66
67 def main():
68 # "N03-110331_04.xml" の "EC01_02141"は仙台市泉区
69 result = read_from_mlitN03("N03-110331_04.xml")["EC01_02141"]
70 g = result["are"]
71 ext = [g.bounds[0], g.bounds[2], g.bounds[1], g.bounds[3]]
72 #
73 fig = plt.figure()
74 fig.set_size_inches(16.53, 11.69)
75 ax1 = fig.add_subplot(projection=ccrs.PlateCarree())
76 ax1.set_extent(ext, ccrs.PlateCarree())
77 ax1.add_wms(
78 wms='https://ows.terrestris.de/osm/service',
79 layers=["OSM-WMS",])
80
81 ax1.add_geometries(
82 [g], ccrs.PlateCarree(),
83 facecolor='#ffff0066',
84 edgecolor='b', linewidth=0.5)
85 import sys, os
86 plt.savefig(os.path.basename(sys.argv[0]) + ".png", bbox_inches="tight")
87
88
89 if __name__ == '__main__':
90 main()
相当雑に扱ってるが、今ワタシがやりたいレベルのことは出来てる、一応。「GM_CompositeCurve.generator idref=”_c00001″」の意味なんか理解してないよ。けど「c00001」の GM_Curve を参照すれば、どうやらあってるっぽいし。まぁそれでいいんじゃないかな、うまくいってるうちは:
いい加減 GML がダメなものだと気付いて欲しいよ…。これ、どっかで前に書いたけどさ、「クラス図が美しい」以上のメリットなんかないんだよ。帯域の無駄遣い!
たぶん最善はこれ:
ある程度一括でほしい場合も標準地域コード インデックスがあるので、大変扱いやすい。geojson もしくは topojson でダウンロード出来る。geojson の方は lonlat なので難しくはないが、GML の構造を引きずってる? ちと気持ち悪い:
1 # -*- coding: utf-8 -*-
2 import io
3 import sys
4 import json
5 from shapely.geometry import asPolygon, asLineString
6 import matplotlib.pyplot as plt
7 import cartopy.crs as ccrs
8
9
10 def main():
11 # 『歴史的行政区域データセットβ版』（CODH作成）
12
13 result = json.load(io.open("04105.geojson", encoding="utf-8"))
14 g = asPolygon(result["features"][0]["geometry"]["coordinates"][0][0])
15 ext = [g.bounds[0], g.bounds[2], g.bounds[1], g.bounds[3]]
16 #
17 fig = plt.figure()
18 fig.set_size_inches(16.53, 11.69)
19 ax1 = fig.add_subplot(projection=ccrs.PlateCarree())
20 ax1.set_extent(ext, ccrs.PlateCarree())
21 ax1.add_wms(
22 wms='https://ows.terrestris.de/osm/service',
23 layers=["OSM-WMS",])
24
25 ax1.add_geometries(
26 [g], ccrs.PlateCarree(),
27 facecolor='#ffff0066',
28 edgecolor='b', linewidth=0.5)
29 import sys, os
30 plt.savefig(os.path.basename(sys.argv[0]) + ".png", bbox_inches="tight")
31
32
33 if __name__ == '__main__':
34 main()
結果の絵は同じ:
いくつか topojson でしか手に入らないものがあるのだが、それらは「仙台市」「さいたま市」などの「区」をまとめあげたものだけっぽいので、なくても困らないと思う。topojson の仕様はちょっと独特だが計算方法もコンパクトに書いてあるので、自力でやらなければならないとしてもどうにかなりそうだと思う。
案の定 python による実装がある。(optional に geopandas 依存だが、必須の依存として勝手に芋づるインストールすることはないのでそこはご安心を。)基本的に「geometry → topology」のためのライブラリだが、「topology → geojson」は可能で、一番簡単なのはたぶんこんな感じ:
1 # -*- coding: utf-8 -*-
2 import io
3 import sys
4 import json
5 from shapely.geometry import asPolygon, asLineString
6 import matplotlib.pyplot as plt
7 import cartopy.crs as ccrs
8 import topojson
9
10
11 def main():
12 # 『歴史的行政区域データセットβ版』（CODH作成）
13
14 # "04100.topojson" は「仙台市」全域。
15 result = json.load(io.open("04100.topojson", encoding="utf-8"))
16 # 「objectname="city"」は nii の "04100.topojson" の場合。これは個々に個々。
17 result = topojson.utils.serialize_as_geojson(result, objectname="city")
18 g = asPolygon(result["features"][0]["geometry"]["coordinates"][0][0])
19 ext = [g.bounds[0], g.bounds[2], g.bounds[1], g.bounds[3]]
20 #
21 fig = plt.figure()
22 fig.set_size_inches(16.53, 11.69)
23 ax1 = fig.add_subplot(projection=ccrs.PlateCarree())
24 ax1.set_extent(ext, ccrs.PlateCarree())
25 ax1.add_wms(
26 wms='https://ows.terrestris.de/osm/service',
27 layers=["OSM-WMS",])
28
29 ax1.add_geometries(
30 [g], ccrs.PlateCarree(),
31 facecolor='#ffff0066',
32 edgecolor='b', linewidth=0.5)
33 import sys, os
34 plt.savefig(os.path.basename(sys.argv[0]) + ".png", bbox_inches="tight")
35
36
37 if __name__ == '__main__':
38 main()
ちなみに、仙台市は「山奥の県境から太平洋まで全部仙台市、しかも隣が山形市」:
1 # -*- coding: utf-8 -*-
2 import io
3 import sys
4 import json
5 from shapely.geometry import asPolygon, asLineString
6 import shapely.ops as sops
7 import matplotlib.pyplot as plt
8 import cartopy.crs as ccrs
9 import topojson
10
11
12 def main():
13 # 『歴史的行政区域データセットβ版』（CODH作成）
14
15 # "04100.topojson" は「仙台市」全域。
16 res_s = json.load(io.open("04100.topojson", encoding="utf-8"))
17 # 「objectname="city"」は nii の "04100.topojson" の場合。これは個々に個々。
18 res_s = topojson.utils.serialize_as_geojson(res_s, objectname="city")
19 g_s = asPolygon(res_s["features"][0]["geometry"]["coordinates"][0][0])
20 # "06201.geojson" は山形市。
21 res_y = json.load(io.open("06201.geojson", encoding="utf-8"))
22 g_y = asPolygon(res_y["features"][0]["geometry"]["coordinates"][0][0])
23 #
24 g = sops.cascaded_union([g_s, g_y])
25 #
26 ext = [g.bounds[0], g.bounds[2], g.bounds[1], g.bounds[3]]
27 #
28 fig = plt.figure()
29 fig.set_size_inches(16.53, 11.69)
30 ax1 = fig.add_subplot(projection=ccrs.PlateCarree())
31 ax1.set_extent(ext, ccrs.PlateCarree())
32 ax1.add_wms(
33 wms='https://ows.terrestris.de/osm/service',
34 layers=["OSM-WMS",])
35 for g, c in [(g_s, '#ffff0066'), (g_y, '#00ff0066')]:
36 ax1.add_geometries(
37 [g], ccrs.PlateCarree(),
38 facecolor=c,
39 edgecolor='b', linewidth=0.5)
40 import sys, os
41 plt.savefig(os.path.basename(sys.argv[0]) + ".png", bbox_inches="tight")
42
43
44 if __name__ == '__main__':
45 main()
政令指定都市になるために合併を繰り返した、1990年になる直前くらいに、今のこのブチ抜き仙台市が完成した。ちなみにワタシは「宮城県泉市立」の中学に入学し、「宮城県仙台市立」の中学を卒業した。むろん同じ中学である。(泉市が吸収合併されたのが 1988 年。合併の是非について市民投票にかけられたので、全国的に一時有名になった。)
ところで、この Geoshapeリポジトリなのだけれど、非常に意欲的で、なおかつ、この成果がまさに「そうそう、それが欲しかったのよ」満載で、ちょっとというかかなりというか相当感動している。
GIS を扱う際に、「広がりの方ではなくて代表点」が必要になるケースは、やればやるほど「とても需要が多い」ことが理解出来てくると思っている。駅データを使いたいと思ったのも、動機はだいたい同じね。
「航空」に関しては、この問題には昔から一つのアプローチがずっと使われ続けている。なにかというと、「ARP」ね。Aerodrome Reference Point といって、「空港」というバカでかいエリアに対する「空港といえば」の代表点を求める統一的な手法。これは個々に個々で好き勝手に決めたり、あるいは「空港の位置」を必要とする事業者たちが自由裁量で考えるんではなくて、まずは「空港事業者」自身が公式に公表するオープンデータであり、しかもその決め方が決まっている。wikipedia が引用している ICAO による説明は「The aerodrome reference point shall be located near the initial or planned geometric centre of the aerodrome and shall normally remain where first established.」で、United States においては「”the approximate geometric center of all usable runway surfaces”, computed as a weighted average of the end of runway coordinates.」。(読み直すまでこの US 方式が全世界的なものだと誤解してた。)
いずれにしても、たとえば航空機が搭載するエアナビゲーションシステムが依拠する「空港の位置」は「空港コードと ARP の位置」というシンプルなペアを保持しておくだけで良く、世界中の空港をかき集めてもせいぜい万のオーダなので、かなり小さなデータとして保持しておける。また、「都市間距離」を示すのにも ARP が非常に便利である。便利、というか、巨大な GIS データに依存したくないなら真っ先に候補になるアプローチ。「北海道と東京の直線距離は?」という問いに、その正確な意味と正確な値が必要ない「概算」でいいなら、はっきりいって ARP を使う以上に楽な手はない。
航空以外でこうした代表点が簡単に見つかるかというと、これはなかなかなくて、「仕方なく centroid」とすることが多い。実際ここ数週間以内のワタシのネタでそれをやってる。けれどもこの centroid、かなり多くの目的向けには大抵「不適切」である。たとえば「北海道」という検索をすると北海道にパン・ズームするアプリケーションが必要だとして、北海道の centroid にパン・ズームすると、まぁ多くのユーザには不満だろう。だいたいは初期値としては札幌中心としてジャンプしてほしいと願うだろう。さらにこれが東京だとさらに不遇で、「東京都全域」の centroid は太平洋上である。
この代表点問題を解消しようとしているのがGeoshapeリポジトリのミッションの一つで、Geoshapeリポジトリ > 歴史的行政区域データセットβ版 > データセットから csv ダウンロード出来て、そこに代表点が与えられている:
1 entry_id,body,body_variants,suffix,prefname,prefname_variants,countyname,countyname_variants,ne_class,hypernym,latitude,longitude,address,code,valid_from,valid_to,source
2 01100A1972,札幌,札幌,市,北海道,北海道,札幌市,札幌市,市区町村/政令指定都市,北海道/札幌市,43.06197200,141.35437400,北海道札幌市,01100,1972-04-01,,1/札幌市役所/札幌市中央区北1条西2/P34-14_01.xml
3 01101A1972,中央,中央,区,北海道,北海道,札幌市,札幌市,市区町村,北海道/札幌市,43.05536000,141.34100500,北海道札幌市中央区,01101,1972-04-01,,1/札幌市中央区役所/札幌市中央区南3条西11/P34-14_01.xml
4 ... (省略) ...
5 03B0150017,和賀,和賀,村,岩手県,岩手県,和賀郡,和賀郡,市区町村,岩手県/和賀郡,39.30666100,140.97916500,岩手県和賀郡和賀村,03000,,,2/和賀庁舎/北上市和賀町横川目11-160/P34-14_03.xml/(19551001\03\03B0150017.wkt.txt)
6 04100A1989,仙台,仙台,市,宮城県,宮城県,仙台市,仙台市,市区町村/政令指定都市,宮城県/仙台市,38.26800800,140.86961700,宮城県仙台市,04100,1989-04-01,,1/仙台市役所/仙台市青葉区国分町3-7-1/P34-14_04.xml
7 04101A1989,青葉,青葉,区,宮城県,宮城県,仙台市,仙台市,市区町村,宮城県/仙台市,38.26908700,140.87037800,宮城県仙台市青葉区,04101,1989-04-01,,1/仙台市青葉区役所/仙台市青葉区1-5-1/P34-14_04.xml
8 04102A1989,宮城野,宮城野,区,宮城県,宮城県,仙台市,仙台市,市区町村,宮城県/仙台市,38.26618900,140.90984100,宮城県仙台市宮城野区,04102,1989-04-01,,1/仙台市宮城野区役所/仙台市宮城野区五輪2-12-35/P34-14_04.xml
9 04103A1989,若林,若林,区,宮城県,宮城県,仙台市,仙台市,市区町村,宮城県/仙台市,38.24418000,140.90074800,宮城県仙台市若林区,04103,1989-04-01,,1/仙台市若林区役所/仙台市若林区保春院前丁3-1/P34-14_04.xml
10 04104A1989,太白,太白,区,宮城県,宮城県,仙台市,仙台市,市区町村,宮城県/仙台市,38.22439500,140.87721200,宮城県仙台市太白区,04104,1989-04-01,,1/仙台市太白区役所/仙台市太白区長町南3-1-15/P34-14_04.xml
11 04105A1989,泉,泉,区,宮城県,宮城県,仙台市,仙台市,市区町村,宮城県/仙台市,38.32636100,140.88161500,宮城県仙台市泉区,04105,1989-04-01,,1/仙台市泉区役所/仙台市泉区泉中央2-1-1/P34-14_04.xml
12 04201A1968,仙台,仙台,市,宮城県,宮城県,,,市区町村,宮城県,38.26800800,140.86961700,宮城県仙台市,04201,1889-04-01,1989-04-01,1/仙台市役所/仙台市青葉区国分町3-7-1/P34-14_04.xml/(19851001\04\04201A1968.wkt.txt)
13 04202A1968,石巻,石巻,市,宮城県,宮城県,,,市区町村,宮城県,38.43445700,141.30290600,宮城県石巻市,04202,,,1/石巻市役所/石巻市穀町14-1/P34-14_04.xml
14 04203A1968,塩竈,塩竈,市,宮城県,宮城県,,,市区町村,宮城県,38.31436000,141.02203000,宮城県塩竈市,04203,,,1/塩竈市役所/塩竈市旭町1-1/P34-14_04.xml
15 04204A1968,古川,古川,市,宮城県,宮城県,,,市区町村,宮城県,38.57713200,140.95556500,宮城県古川市,04204,,2006-03-31,1/大崎市役所/大崎市古川七日町1-1/P34-14_04.xml/(20050101\04\04204A1968.wkt.txt)
16 04205A1968,気仙沼,気仙沼,市,宮城県,宮城県,,,市区町村,宮城県,38.90812700,141.57004400,宮城県気仙沼市,04205,,,1/気仙沼市役所/気仙沼市八日町1-1-1/P34-14_04.xml
17 04206A1968,白石,白石,市,宮城県,宮城県,,,市区町村,宮城県,38.00246700,140.61988800,宮城県白石市,04206,,,1/白石市役所/白石市大手町1-1/P34-14_04.xml
18 04207A1968,名取,名取,市,宮城県,宮城県,,,市区町村,宮城県,38.17150100,140.89184900,宮城県名取市,04207,,,1/名取市役所/名取市増田字柳田80/P34-14_04.xml
19 04208A1968,角田,角田,市,宮城県,宮城県,,,市区町村,宮城県,37.97701600,140.78154300,宮城県角田市,04208,,,1/角田市役所/角田市角田字大坊41/P34-14_04.xml
20 04209A1971,多賀城,多賀城,市,宮城県,宮城県,,,市区町村,宮城県,38.29380300,141.00436900,宮城県多賀城市,04209,1971-11-01,,1/多賀城市役所/多賀城市中央2-1-1/P34-14_04.xml
21 04210A1971,泉,泉,市,宮城県,宮城県,,,市区町村,宮城県,38.32636100,140.88161500,宮城県泉市,04210,1971-11-01,1988-03-01,1/仙台市泉区役所/仙台市泉区泉中央2-1-1/P34-14_04.xml/(19851001\04\04210A1971.wkt.txt)
22 04211A1971,岩沼,岩沼,市,宮城県,宮城県,,,市区町村,宮城県,38.10430300,140.86994900,宮城県岩沼市,04211,1971-11-01,,1/岩沼市役所/岩沼市桜1-6-20/P34-14_04.xml
23 ... (省略) ...
24 13100A1968,特別区部,特別区部,特別区部,東京都,東京都,特別区部,特別区部,市区町村/特別区部,東京都/特別区部,35.68956800,139.69171700,東京都特別区部,13100,1947-03-15,,1/東京都庁/東京都新宿区西新宿2-8-1/P28-13_13.xml
25 ... (省略) ...
「変遷」を記録するデータセットなので最新状態だけしか興味がない人にとっては利用に余分な一手間が必要となるものだが、別にそのくらいの手間は大したことではなかろう。たとえば:
1 # -*- coding: utf-8 -*-
2 import io
3 import csv
4 from datetime import date
5 from collections import namedtuple, defaultdict
6
7
8 def from_geoshape_city_csv():
9 # 『歴史的行政区域データセットβ版』（CODH作成）
10 # https://geonlp.ex.nii.ac.jp/dictionary/geoshape-city/
11
12 # 目的によって、どういうまとめ方が使いやすいかはいろいろ違うが、「現在」だけに
13 # 興味があるケース相手に一番楽なのは、「valid_to を過ぎているかどうか」で
14 # 二分すること、である。
15 now = date.today()
16 result_recent, result_past = defaultdict(list), defaultdict(list)
17 reader = csv.reader(io.open("geoshape-city.csv", encoding="utf-8"))
18 fields = next(reader)
19 Record = namedtuple("Record", fields)
20 def _map(rec):
21 import re
22 temp = []
23 for v in rec:
24 if re.match(r"\d+\.\d+", v):
25 temp.append(float(v))
26 elif re.match(r"\d+\-\d+\-\d+", v):
27 temp.append(date(*map(int, v.split("-"))))
28 else:
29 temp.append(v)
30 return Record(*temp)
31 for rec in (_map(rec) for rec in reader):
32 pc = rec.entry_id[:2] # iso-3166-2 都道府県コードに一致
33 if rec.valid_to and rec.valid_to < now:
34 result_past[pc].append(rec)
35 else:
36 result_recent[pc].append(rec)
37
38 return result_recent, result_past
39
40
41 #from_geoshape_city_csv()
代表点が情報として与えられることのメリットを実感してもらうには、やはり東京の例がいい:
1 # -*- coding: utf-8 -*-
2 import sys, os
3 import matplotlib.pyplot as plt
4 import cartopy.crs as ccrs
5
6
7 def main():
8 # windows ユーザ以外の方は fontproperties は自分でどうにかしなはれ。
9 fontprop = {"fname": "c:/Windows/Fonts/meiryo.ttc"}
10
11 #
12 centers = [
13 ("東京都全域の centroid",
14 (139.5623713172737, 35.07700753947662)),
15
16 ("東京都陸域の最大面積を持つポリゴンの centroid",
17 (139.4341927146594, 35.71538109419526)),
18
19 ("geoshape.ex.nii.ac.jp プロジェクト定義の東京都代表点 (東京都庁)",
20 (139.691717, 35.689568)),
21 ]
22 exp_km = 50
23 exp = 1. / (3600 * ((30 + 27) / 2)) * exp_km * 1000 # 約 exp_km km
24 for i, (tit, ct) in enumerate(centers):
25 ext = [ct[0] - exp, ct[0] + exp, ct[1] - exp, ct[1] + exp]
26 #
27 fig = plt.figure()
28 fig.set_size_inches(16.53, 11.69)
29 ax1 = fig.add_subplot(projection=ccrs.PlateCarree())
30 ax1.set_extent(ext, ccrs.PlateCarree())
31 ax1.add_wms(
32 wms='https://ows.terrestris.de/osm/service',
33 layers=["OSM-WMS",])
34 ax1.set_title(
35 "{} を中心とする一辺約 {} km 矩形".format(tit, exp_km * 2),
36 fontproperties=fontprop)
37 ax1.text(
38 ct[0], ct[1],
39 "lon={:.6f}\nlat= {:.6f}".format(*ct),
40 transform=ccrs.Geodetic()._as_mpl_transform(ax1),
41 fontproperties=fontprop,
42 color='black',
43 bbox=dict(facecolor='white', alpha=0.5, boxstyle='round'),
44 )
45 plt.savefig(
46 os.path.splitext(
47 os.path.basename(sys.argv[0]))[0] + ".{}.png".format(i + 1),
48 bbox_inches="tight")
49
50
51 if __name__ == '__main__':
52 main()
この絵がたとえば地図検索ソフトの検索結果だとして、しかも「東京都へ飛べ!」と指示した結果が一つ目の絵なのだとしたら、かなりのユーザはがっかりすることだろう。二つ目の絵は惜しいけれど、実は日本全国を描画範囲とするケースではこの位置がむしろ神奈川に属しているかのように見えてしまって、かなり不適切に感じると思う。やってみればわかる。これは、東京の陸域が意外に東西に長いことが理由。手持ちのポリゴンの粒度にも依存するが、たとえば polbnda_jpn だけに依存する場合は、新宿区の centroid を使う、ということはまぁ一応出来るが、日本全国あらゆる場所に対して同じようなことを考えたいならば、(たとえば polbnda_jpn には名古屋市中区のポリゴンはないので)結局のところ「代表点データベース」のようなものが必要となり、「ありがとう、geoshape プロジェクト」となるのだ。
にしても、「日本」に関してはこれでいいんだけれど、「世界」相手にこれに類することが必要になったケースはどうしたらいいもんだかね? 各国についてそれぞれ、かなぁ?
