ISOんな都道府県コードなデータセットんな – 3166-2式

これとかで多重定義しまくっててウザいやつ。

ISO 3166 の件

2021-04-25 にみつけた pycountry

こういうものを扱う機会は個人的にめっちゃ多かったわりに、一度もちゃんとしたやつを使おうとしたことがないので、一度みとこうかと。こんなんやだぁ:

 1 # こうやって毎度自力るんじゃのーて、皆の衆による共有資材を使えんのけ?
 2 _PREFCODE = {
 3     '北海道': 1,
 4     '青森県': 2,
 5     '岩手県': 3,
 6     '宮城県': 4,
 7     '秋田県': 5,
 8     '山形県': 6,
 9     '福島県': 7,
10     '茨城県': 8,
11     '栃木県': 9,
12     '群馬県': 10,
13     '埼玉県': 11,
14     '千葉県': 12,
15     '東京都': 13,
16     '神奈川県': 14,
17     '新潟県': 15,
18     '富山県': 16,
19     '石川県': 17,
20     '福井県': 18,
21     '山梨県': 19,
22     '長野県': 20,
23     '岐阜県': 21,
24     '静岡県': 22,
25     '愛知県': 23,
26     '三重県': 24,
27     '滋賀県': 25,
28     '京都府': 26,
29     '大阪府': 27,
30     '兵庫県': 28,
31     '奈良県': 29,
32     '和歌山県': 30,
33     '鳥取県': 31,
34     '島根県': 32,
35     '岡山県': 33,
36     '広島県': 34,
37     '山口県': 35,
38     '徳島県': 36,
39     '香川県': 37,
40     '愛媛県': 38,
41     '高知県': 39,
42     '福岡県': 40,
43     '佐賀県': 41,
44     '長崎県': 42,
45     '熊本県': 43,
46     '大分県': 44,
47     '宮崎県': 45,
48     '鹿児島県': 46,
49     '沖縄県': 47,
50 }

てことな。

あるだろうあるに決まってる、とずっと理解しつつも調べてこなかったので、ISO の何番、JIS で何番の規格なのかについて「初めて知る」→ ISO-3166-2 ね、今の目的のものは。

そして、「あるだろうあるに決まってる、とずっと理解しつつも調べてこなかった」のは、まぁ当たり前だけど「それにまつわるインフラを手に入れることが即戦力になるとは限らなくて、むしろ見合わないに違いない」と思うから、だね。上の例だと、「所詮は47行のデータ」なので、管理の仕方の工夫さえなんとかすれば、「巨大かもしれないそのインフラは、きっといらんと思うんではないのかな」てことだわな。全世界をターゲットにするソフトウェアを開発するならあってもいいわけよ。けれども個人ユースでローカルな小さな見える化をしたいだけには、「きっと牛刀なのであろうなぁ」と。

案の定、て感じかなぁ、ひとまず「pycountry」。10MB、だったかな? まぁ当然のごとくデカい。ターゲットは 3166-2 だけでなくて:

  • 639-3 Languages
  • 3166 Countries
  • 3166-3 Deleted countries
  • 3166-2 Subdivisions of countries
  • 4217 Currencies
  • 15924 Scripts

このデータ、JSON なのだが、どこか本家から取り寄せて加工しての JSON なのかなぁ、と思ったら、Debian の pkg-isocodes によるものなんだって。

使うつもりの場合は、まぁなんなら「この JSON の存在だけでもありがてー」のかもしらんね、とは思う。python のパッケージとしてちゃんと使うつもりならたとえば:

1 >>> from pycountry import subdivisions
2 >>> subdivisions.lookup("JP-13")
3 Subdivision(code='JP-13', country_code='JP', name='Tokyo', parent_code=None, type='Prefecture')
4 >>> subdivisions.lookup("TokYo")
5 Subdivision(code='JP-13', country_code='JP', name='Tokyo', parent_code=None, type='Prefecture')
6 >>> [(t.name, t.code) for t in subdivisions.lookup("JP") if t.name.startswith("M")]
7 [('Miyazaki', 'JP-45'), ('Mie', 'JP-24'), ('Miyagi', 'JP-04')]

ここまで書いたとこで、気付く人はもうとっくに気付くと思うけれど、ワタシの目下必要なもの、つまり「宮城県は4で、13は東京」という関係についてはまかなってくれない、てこと。この「とっても有意義で大事なインフラ」を入手することによってワタシが出来るようになったのは「Miyagi は4で、13は Tokyo」。そう、結局ワタシは「_PREFCODE」を自力り続ける必要がある。

pycountry についてはそうなんだけれど、うーん、元にした Debian の pkg-isocodes の方って、翻訳が管理されてるのよねぇ。ワタシが今やってるようなもの相手には明らかにツーマッチにはなるけれど、何かデカいソフトウェアでこの手のことが必要になった場合に、この PO に依存するのもいいかもしれない。pycountry で「Tokyo」まで引ければ「東京」を導ける。逆引きが面倒だろうね。

うーん、この先をどうしようか。「ワタシさえ良ければいい」のつもりなら、ワタシには全然必要ない検証になっちゃうのよねぇ。一応技術ネタとして「gettext」ネタになるので、役に立たないこともないとは思うのだけれども。(確かかなり前に python から gettext にアクセスするネタはやった気がする。)


少しばかり考えあぐねるも、…やるか。

とりあえず、最もバカなアプローチとしては「ja.po」を単にダウンロードしてきて「テキストファイルのまま、テキストファイルとして自力で解読する」のも必ずしも「絶対にナシ」とまでは言えない。難しいフォーマットではないから。継続行の扱いだのが「ダル」いだけで。ただまぁ「せっかく電池付き」なのだから、ちゃんとした方がお得はお得。メッセージのインストール場所に困り果てる Windows はともかくとして、Debian ユーザだったりしたならば、「そのまま使えちゃうぜ」の道も拓けないこともないわけだし、と。

と思ったのだが、ちゃんとするためにはそもそも「コンパイル」しないとね、てことなのだが、持ってたっけか msgfmt。あ、持ってた:

1 [me@host: ~]$ type -p msgfmt
2 /c/Program Files/emacs-27.2-x86_64/bin/msgfmt.exe

確か個別に GNU の本家からもダウンロードした記憶はあるんだけど、そうか、emacs が抱え込んでるのね。

ja.po を salsa.debian.org の WEB UI 経由で「ダウンロード」すると「iso_3166-2_ja.po」という名前で持ってくる(ハメになる)。「ちゃんと」するために、色々流儀に従う必要がある。たぶんこうだよな、合ってるか?:

1 [me@host: ~]$ mkdir -p _locale/ja_JP/LC_MESSAGES
2 [me@host: ~]$ msgfmt iso_3166-2_ja.po -o _locale/ja_JP/LC_MESSAGES/iso_3166-2.mo
3 [me@host: ~]$ # ↓「できてる?」確認
4 [me@host: ~]$ stat -c '%n: %s' _locale/ja_JP/LC_MESSAGES/iso_3166-2.mo
5 _locale/ja_JP/LC_MESSAGES/iso_3166-2.mo: 122879

「ja_JP」なのよね問題は。CPython Windows 版の locale ってとてもヘンで、正しいのかわからんのだがこんな妙なことになっとるんよね:

1 Python 3.9.2 (tags/v3.9.2:1a79785, Feb 19 2021, 13:44:55) [MSC v.1928 64 bit (AMD64)] on win32
2 Type "help", "copyright", "credits" or "license" for more information.
3 >>> import locale
4 >>> locale.setlocale(locale.LC_ALL, '')
5 'Japanese_Japan.932'

これがあるので Japanese_Japan.932/LC_MESSAGES/iso_3166-2.mo としなけりゃならんのか、となりそうなのだが、それは違う気がする。(というか「だから locale を扱うのはイヤ」なのだよ。東アジア人にはほとんど旨味がないからね、locale。) とりあえず ja_JP で行っとくことにする。

なお、Unix 系であれば、この mo 置き場は慣習でほぼ完全に特定できて、「/usr/share/locale」(など)である。Windows にだって標準場所はあるのだが、たぶん「Microsoft 謹製アプリケーション」以外では使われてない、ような気がする。ちゃんと意識したことがないからわからんけれど。(少なくとも OSS はこれの存在をほとんど無視している。)そういうわけで、「mo にコンパイルできたのはいいが、Windows では置き場所に困る」という結論に(ほぼ)なる、ので、とにかく今は「カレントディレクトリに _locale ディレクトリ」という構造を前提にしておく。(つまり上のコマンドラインのまんま。)

のうえで。「gettext」は Python の標準添付ライブラリである。ここまで準備できたので、もう使える:

1 >>> import gettext
2 >>> tr = gettext.translation("iso_3166-2", "./_locale", ["ja_JP"])
3 >>> tr.gettext("Tokyo")
4 '東京'

てわけで、「pycountry + iso_3166-2.mo」:

1 >>> import gettext
2 >>> tr = gettext.translation("iso_3166-2", "./_locale", ["ja_JP"])
3 >>> # 「_("msgid")」とするのが gettext 流儀で、python としてはかなり異端。
4 >>> # が、まぁとりあえず従っとく。
5 >>> _ = tr.gettext
6 >>> #
7 >>> from pycountry import subdivisions
8 >>> [(t.name, _(t.name), t.code) for t in subdivisions.lookup("JP") if t.name.startswith("M")]
9 [('Miyazaki', '宮崎', 'JP-45'), ('Mie', '三重', 'JP-24'), ('Miyagi', '宮城', 'JP-04')]

うん、出来る、よかね。

繰り返すけれど、「そこまでやる必要あんの?」。てこと。少なくとも、今のワタシのちんまいニーズからすると、これは全然見合ってない。導入コストもランタイムコストも、それに保守コストも、高過ぎるであろうよ。もちろん、出来ると知っとくことは悪くはないし、大きなプロジェクトならここまでするのも価値はある、当然。

2021-06-02 追記: 節穴の件 + cartopy 関連のネタから芋づるで見つけた別のもの

ほんと、注意力散漫というか…。「locale」なんだけどさ、pycountry インストールすると、同梱されるじゃないの。なんで気付かなかったかね。ちゃんと mo の状態で配置されてくれるので、自分でコンパイルするとか書いたのは全部無駄骨。まぁ笑ってやってくれよ。

もう一つが追記の本題。「3166-2」と題してるネタからすると外れるネタなんだけれど、根っこの同じ 3166 のネタとして。

cartopy から WMTS を利用するネタで、「United States」範囲を描画しているんだけれど、こういうことをするときに「ちょうどいいあんばいの bbox (or extent)」を決めるのが存外ダルい。WEB マップも、基本的に「地図中心」の緯度/経度のみを意識させるインターフェイスになっていて、現在可視状態になっている範囲の bbox を教えてくれるものはない(と思う)。

まぁかなり前からこういうことは思っていたけれど、ただ、ワタシに関しては多くの場合は「日本全域もしくはその一部」だけを相手にし、そしてその場合は「一度わかってしまえば同じ範囲を使い続ける」し、たとえば気象庁提供の MSM GPV を可視化したいケースなら、基本的には得られる行列の範囲を使えばいい(pygrib を使えば直接その範囲が返ってくる)、ということで、今回みたいな「cartopy で色んな WMTS サーバ相手に遊び」たい際に発生する「あるときは USA、あるときは南極、あるときはウクライナ」の bbox を求めたくなるなんてことはなかった。

てわけで、ちょっと探してたら、「python」というキーワードを入れてないのに、country-bounding-boxes を見つけた。Natural Earth のデータを利用するらしい。これの shapefile から読み解くシカケなので、cartopy ネタの(5)でも触れた pyshp を使う。また、iso-3166 に基づくルックアップというインターフェイスのために、先に見つけた pycountry とは別の iso3166 に依存。

小粒つーか、まぁなんとも「大したことをしてない」んだけれど、正直「まさしくワタシがやりたかったことそのもの」は出来る:

 1 >>> from country_bounding_boxes import country_subunits_by_iso_code
 2 >>> 
 3 >>> for it in country_subunits_by_iso_code("USA"):
 4 ...     print(it.postal, it.brk_a3, it.brk_name, it.bbox)
 5 ... 
 6 AK USK Alaska (172.494824219, 51.3722167969, 179.779980469, 53.0129882813)
 7 AK USK Alaska (-178.19453125, 51.6036621094, -130.0140625, 71.4076660156)
 8 US USB U.S.A. (-124.709960938, 24.5423339844, -66.9870117187, 49.3696777344)
 9 HI USH Hawaii (-160.243457031, 18.9639160156, -154.804199219, 22.2231445312)
10 >>> for it in country_subunits_by_iso_code("JPN"):
11 ...     print(it.postal, it.brk_a3, it.brk_name, it.bbox)
12 ... 
13 BI JPB Bonin Is. (142.107128906, 26.6156738281, 142.202148438, 26.7264648437)
14 NN JPO Nanseishoto (123.679785156, 24.2660644531, 129.714648438, 28.5174804688)
15 SH JPS Shikoku (132.032617188, 32.7520019531, 134.738867188, 34.3583984375)
16 J JPN Japan (128.649121094, 30.24140625, 141.329199219, 45.4654785156)
17 HN JPH Honshu (130.889257813, 33.4869628906, 141.993164062, 41.5055664062)
18 HK JPK Hokkaido (139.820898438, 41.4232421875, 145.833007812, 45.5095214844)
19 KY JPY Kyushu (129.580078125, 31.0151367188, 132.00859375, 33.9277832031)
20 IZ JPI Izushoto (139.768945313, 33.0454589844, 139.873632813, 33.1292480469)

なので「アメリカ本土」をクリップするのに例えばこんなふうに出来る:

 1 # -*- coding: utf-8 -*-
 2 import matplotlib.pyplot as plt
 3 import cartopy.crs as ccrs
 4 from country_bounding_boxes import country_subunits_by_iso_code
 5 
 6 
 7 def main():
 8     url = 'https://mrdata.usgs.gov/mapcache/wmts'
 9     layers = ['mrds-PGE_OS', 'mrds-STN_D_G', 'mrds-DIT', 'mrdsc']
10 
11     fig = plt.figure()
12     fig.set_size_inches(16.53, 11.69)
13     ax = fig.add_subplot(projection=ccrs.GOOGLE_MERCATOR)
14     bbox = [
15         it.bbox for it in country_subunits_by_iso_code("USA")
16         if it.subunit == "United States"][0]
17     ax.set_extent([bbox[0], bbox[2], bbox[1], bbox[3]])
18     for layer in layers:
19         ax.add_wmts(url, layer)
20     ax.coastlines("110m")
21     import os, sys
22     plt.savefig(os.path.basename(sys.argv[0]) + ".png", bbox_inches="tight")
23     #plt.show()
24     plt.close(fig)
25 
26 
27 if __name__ == '__main__':
28     main()

この WMTS は本日時点での cartopy だとエラーになる。cartopy ネタ(6.5)参照。エラーに対する措置を済ませてあれば、結果はここうこうこうこう:
wmts_exam_03.py

むろんこいつは「Natural Earth のデータ管理方法」そのものに依存するので、サブ領域への分割がお気に召さないならこのプロジェクトは使えない。だいいち 3166-2 の subdivision に関してはガン無視だもの、こやつ。まぁ日本であれば、国土地理院配布の shapefile に基づいて、この parse.py 相当のことをすればいいんじゃないかな。

それと、そもそも country_subunits_by_iso_code をはじめとして、トップレベルにある「お便利関数」がどれもダサいというかモッサリしてて、場合によっては不自由に感じるんではないかと思う。たとえば show_all_bounding_boxes なんかかなり意味わからんし。結局「counties さえあればええんちゃうか」てことだと思うのよね:

1 >>> from country_bounding_boxes.generated import countries
2 >>> 
3 >>> for it in [it for it in countries if it.postal == "J"]:
4 ...     print(it.postal, it.brk_a3, it.brk_name, it.bbox)
5 ... 
6 J JAM Jamaica (-78.3395019531, 17.7149414062, -76.2107910156, 18.5222167969)
7 J JOR Jordan (34.95078125, 29.1904785156, 39.2927734375, 33.3722167969)
8 J JPN Japan (128.649121094, 30.24140625, 141.329199219, 45.4654785156)

2021-06-03 02:30追記:
日本の都道府県単位の bbox も欲しいんじゃい、のケースで「国土地理院配布の shapefile に基づいて、この parse.py 相当のことをすればいいんじゃないかな」と言った。実際に試みてみる:

結局 country_bounding_boxes の成果ってほとんどこれだけなのだよ
 1 #! py -3
 2 # -*- coding: utf-8 -*-
 3 import shapefile
 4 
 5 
 6 def extract_data(sh_fn):
 7     sf = shapefile.Reader(sh_fn)
 8     fields = [f[0] for f in sf.fields if isinstance(f, list)]
 9 
10     print("#! py -3")
11     print("# -*- coding: utf-8 -*-")
12     print("#")
13     print("# extracted from " + sh_fn)
14     print("")
15     print("from collections import namedtuple")
16     print("")
17     print("Part = namedtuple('Part', [")
18     print("    'bbox'," +
19           ','.join(["\n    '{}'".format(f) for f in fields]) + "])")
20 
21     shapes = sf.shapes()
22     records = sf.records()
23 
24     print("parts = [")
25     for i in range(0, len(records)):
26         rec = records[i]
27 
28         box = "({}, {}, {}, {})".format(*shapes[i].bbox)
29         fs = ','.join(['\n        {}={}'.format(k, repr(v))
30                        for (k, v) in zip(fields, rec)])
31         print('    Part(')
32         print('        bbox=' + box + ',' + fs + '),')
33 
34     print(']')
35 
36 if __name__ == "__main__":
37     extract_data("polbnda_jpn")

polbnda_jpn.shp を手にしていて、それ相手に実行すればこういう結果が得られるわけだ:

  1 #! py -3
  2 # -*- coding: utf-8 -*-
  3 #
  4 # extracted from polbnda_jpn
  5 
  6 from collections import namedtuple
  7 
  8 Part = namedtuple('Part', [
  9     'bbox',
 10     'f_code',
 11     'coc',
 12     'nam',
 13     'laa',
 14     'pop',
 15     'ypc',
 16     'adm_code',
 17     'salb',
 18     'soc'])
 19 parts = [
 20     Part(
 21         bbox=(140.991607666016, 42.781600952148416, 141.50590515136702, 43.1903991699219),
 22         f_code='FA001',
 23         coc='JPN',
 24         nam='Hokkai Do',
 25         laa='Sapporo Shi',
 26         pop=1930496,
 27         ypc=2014,
 28         adm_code='01100',
 29         salb='UNK',
 30         soc='JPN'),
 31     Part(
 32         bbox=(140.69290161132804, 41.710933685302734, 141.18690490722702, 42.01013183593753),
 33         f_code='FA001',
 34         coc='JPN',
 35         nam='Hokkai Do',
 36         laa='Hakodate Shi',
 37         pop=274485,
 38         ypc=2014,
 39         adm_code='01202',
 40         salb='UNK',
 41         soc='JPN'),
 42     Part(
 43         bbox=(140.84489440918003, 43.06146621704096, 141.292892456055, 43.239933013915994),
 44         f_code='FA001',
 45         coc='JPN',
 46         nam='Hokkai Do',
 47         laa='Otaru Shi',
 48         pop=127224,
 49         ypc=2014,
 50         adm_code='01203',
 51         salb='UNK',
 52         soc='JPN'),
 53     # ... (snip) ...
 54     Part(
 55         bbox=(122.93349314581921, 24.43718285160969, 123.04370882684186, 24.473849721726832),
 56         f_code='FA001',
 57         coc='JPN',
 58         nam='Okinawa Ken',
 59         laa='Yonaguni Cho',
 60         pop=1551,
 61         ypc=2014,
 62         adm_code='47382',
 63         salb='UNK',
 64         soc='JPN'),
 65     Part(
 66         bbox=(139.82659912109403, 35.61246490478524, 139.83929443359398, 35.63126754760743),
 67         f_code='FA001',
 68         coc='JPN',
 69         nam='Tokyo To',
 70         laa='UNK',
 71         pop=-99999999,
 72         ypc=0,
 73         adm_code='UNK',
 74         salb='UNK',
 75         soc='JPN'),
 76     Part(
 77         bbox=(139.79010009765602, 35.59385299682624, 139.819152832031, 35.61315155029298),
 78         f_code='FA001',
 79         coc='JPN',
 80         nam='Tokyo To',
 81         laa='UNK',
 82         pop=-99999999,
 83         ypc=0,
 84         adm_code='UNK',
 85         salb='UNK',
 86         soc='JPN'),
 87     Part(
 88         bbox=(139.80560302734403, 35.58953475952153, 139.82679748535196, 35.607334136962905),
 89         f_code='FA001',
 90         coc='JPN',
 91         nam='Tokyo To',
 92         laa='UNK',
 93         pop=-99999999,
 94         ypc=0,
 95         adm_code='UNK',
 96         salb='UNK',
 97         soc='JPN'),
 98     Part(
 99         bbox=(139.79342651367202, 35.585582733154304, 139.81179809570304, 35.59846496582033),
100         f_code='FA001',
101         coc='JPN',
102         nam='Tokyo To',
103         laa='UNK',
104         pop=-99999999,
105         ypc=0,
106         adm_code='UNK',
107         salb='UNK',
108         soc='JPN'),
109     Part(
110         bbox=(140.0500258933917, 31.43783834431019, 140.05171961897773, 31.440171031688116),
111         f_code='FA001',
112         coc='JPN',
113         nam='Tokyo To',
114         laa='UNK',
115         pop=-99999999,
116         ypc=0,
117         adm_code='UNK',
118         salb='UNK',
119         soc='JPN'),
120     Part(
121         bbox=(140.286499023438, 30.472333908081076, 140.31419372558602, 30.495466232299787),
122         f_code='FA001',
123         coc='JPN',
124         nam='Tokyo To',
125         laa='UNK',
126         pop=-99999999,
127         ypc=0,
128         adm_code='UNK',
129         salb='UNK',
130         soc='JPN'),
131     Part(
132         bbox=(140.34124834007002, 29.79278603571808, 140.34274182145802, 29.793919386842585),
133         f_code='FA001',
134         coc='JPN',
135         nam='Tokyo To',
136         laa='UNK',
137         pop=-99999999,
138         ypc=0,
139         adm_code='UNK',
140         salb='UNK',
141         soc='JPN'),
142     Part(
143         bbox=(136.80044555664097, 34.996093749999986, 136.82168579101602, 35.017635345458956),
144         f_code='FA001',
145         coc='JPN',
146         nam='Aichi Ken',
147         laa='UNK',
148         pop=-99999999,
149         ypc=0,
150         adm_code='UNK',
151         salb='UNK',
152         soc='JPN'),
153 ]

cartopy ネタの(5)と違って、要するに「geometry を bbox という要約情報に変換する以外のことはしてない」わけで、有意義かと言われると「そうでもない」。ただ、「bbox だけが欲しいケースって結構多いんだよな」ということでこれを正当化する気になれるかどうかだけが問題だ。そう、多いんだよね、確かに。


2021-06-03 08:10追記:
いやぁ。やっぱりさぁ、「geometry を bbox という要約情報に変換する以外のことはしてない」だとか、あるいは「cascade_union 済のものを予め」みたいなのって可搬性がなくて、「その場しのぎ」なのだよねぇ。全情報を別フォーマットに変換して置いておく、以外は、結局何かを失ってしまう。結局はこういうことがしたいんだと思うのよね:

 1 # -*- coding: utf-8 -*-
 2 from collections import OrderedDict
 3 import shapely.geometry as sgeom
 4 import shapely.ops as sops
 5 from cartopy.io.shapereader import Reader as ShapeReader
 6 import matplotlib.pyplot as plt
 7 import cartopy.crs as ccrs
 8 from cartopy.io.img_tiles import Stamen
 9 
10 
11 def collect(cond=lambda a, g: True, attrfilter=lambda f: True):
12     # Source: Geospatial Information Authority of Japan website
13     #         (https://www.gsi.go.jp/kankyochiri/gm_japan_e.html)
14     filename = "polbnda_jpn"
15     reader = ShapeReader(filename)
16     result = []
17     for r in reader.records():
18         attr = r.attributes
19         if cond(attr, r.geometry):
20             result.append((
21                 {k: v
22                  for k, v in attr.items()
23                  if attrfilter(k)}, r.geometry))
24     return result
25 
26 
27 def main():
28     # 東京23区だけ抽出してみる。
29     res = collect(
30         cond=lambda a, g: a["nam"] == "Tokyo To" and \
31         a["laa"].endswith(" Ku"),
32         attrfilter=lambda f: f in ("laa",))
33     geom = sops.cascaded_union([g for a, g in res])
34     if geom.geom_type == 'MultiPolygon':
35         bnds, geoms = geom.bounds, list(geom)
36     else:
37         bnds, geoms = geom.bounds, [geom]
38     ext = [bnds[0], bnds[2], bnds[1], bnds[3]]
39     tiler = Stamen(style="terrain")
40     fig = plt.figure()
41     fig.set_size_inches(16.53, 11.69)
42     ax1 = fig.add_subplot(projection=tiler.crs)
43     ax1.set_extent(ext, ccrs.PlateCarree())
44     ax1.gridlines()
45     ax1.add_image(tiler, 12)
46     #
47     ax1.add_geometries(
48         geoms, ccrs.PlateCarree(),
49         facecolor='#ffff0033',
50         edgecolor='b', linewidth=2)
51     import sys, os
52     plt.savefig(os.path.basename(sys.argv[0]) + ".png", bbox_inches="tight")
53 
54 
55 main()

xxxxx2.py
この例の場合は「東京23区」というのが興味の対象で、これが可能なのは、当然各区のポリゴンが用意されているからであって、「cascade_union 済」としたい場合の決断を「都道府県単位」としてしまったら最後、「東京23区」を取り出すことは出来ない。また、bbox について、cartopy.io.shapereader.Reader、shapefile.Reader のいずれからであっても各 shape レコードの bbox (bounds) を取得できるけれど、当然「東京23区」として cascade_union したポリゴンの bounds はそれらとは違って、(コードでやって見せている通り)改めて求める必要がある。あ、でも country_bounding_boxes のアプローチは「東京23区の bbox」を各 bbox の cascade_union で求めることが出来るので、「bbox だけが欲しい」というニーズには答え続けられるのか…。

なんというか「よく使うグルーピング単位についての演算済ポリゴンがあると便利だよなぁ、bounds だけ欲しいケースだってやっぱないとは言えんし」と思う一方で、ニーズのバリエーションが思ったより多くて、整理しづらいよなぁと思う。結局は、国土地理院や Natural Earth が整理している最小単位はまずは維持した方が、結局はお得なのかもしれんなぁと思った、て話。もう一つ例:

 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()

xxxxx3.py
これは「東京23区」というグルーピング以外は、残り全部もともとの polbnda_jpn による「laa」グルーピングに従う、としている。こんなことも結構よくやりたくなると思うけれど、「prededined として皆が使えるとハッピーになる」というほどの可搬性はなく、ゆえに「欲しい人が欲しいときに似たやり方で都度ご自由に」が正解だ、てことになるってこった。にしても、何か汎用ユーティリティとして切り出せんもんかねぇ? 毎度やるのはそれなりに鬱陶しいぞこれ。