見出し画像

都道府県丸さランキング　製作ノート

柞刈湯葉

この記事は前回の「47都道府県丸さランキング」の製作ノートである。「ああいう検証記事はどういうテンションで作ってるのか」を知りたい非プログラマに向けて書いた。初歩的なことしかやってないので職業プログラマが見ればイライラする思うがご了承願いたい。

今回の目的は、県の丸さを「その県と同じ面積の円に、県の何%を収めることができるか」として定量化し、そのランキングを作ることである。

「丸さ」の例

地図を扱う、といっても画像データを扱うわけではない。地図はそもそも測量とかで得られた数値データから作るものであり、そして数値データの方がプログラム的には扱いやすい。なので「数値地図」というものを探す。国土交通省が出している国土数値情報を使う。

伊能忠敬の時代は日本地図を国外に持ち出すだけで大罪になった（シーボルト事件）が、いまや国が積極的に公開しているのだから大したものである。
データは「XML」という形式である。テキストエディタで開くとこうなっている。

<GM_PointArray.column>
    <jps:GM_Position.direct>
        <DirectPosition.coordinate>37.15421400 139.95701100</DirectPosition.coordinate>
        <DirectPosition.dimension>2</DirectPosition.dimension>
    </jps:GM_Position.direct>
</GM_PointArray.column>
<GM_PointArray.column>
    <jps:GM_Position.direct>
        <DirectPosition.coordinate>37.15420900 139.95749900</DirectPosition.coordinate>
        <DirectPosition.dimension>2</DirectPosition.dimension>
    </jps:GM_Position.direct>
</GM_PointArray.column>
<GM_PointArray.column>
    <jps:GM_Position.direct>
        <DirectPosition.coordinate>37.15420300 139.95766500</DirectPosition.coordinate>
        <DirectPosition.dimension>2</DirectPosition.dimension>
    </jps:GM_Position.direct>
</GM_PointArray.column>
<GM_PointArray.column>
    <jps:GM_Position.direct>
        <DirectPosition.coordinate>37.15419000 139.95782700</DirectPosition.coordinate>
        <DirectPosition.dimension>2</DirectPosition.dimension>

プログラミングに縁のない方はよく「プログラミングできる人はこういうのの意味が全部わかるんだろうな〜」と誤解するが全然そんなことはない。「なんか英語が < > に囲まれていて、その外側に数字が書いてあるな〜」くらいの認識にしかならない。

ただ今回は「これは栃木県の数値地図だ」という予備知識があるので、ちょいちょい書かれている 37.154 139.957 といった数字が北緯・東経っぽいと気づく。これらの数字は <DirectPosition.coordinate> という行に入っているので、全部取り出すことにする。

座標データをひとつひとつコピー＆ペーストするのは面倒なので、Python というプログラミング言語を使って抽出する。4行で書ける。Python というのはこういうダルい作業を自動化するために存在する。

import xml.etree.ElementTree as ET

tree = ET.parse('tochigi.xml')
root = tree.getroot()
border = [ point.text for point in root.findall('.//DirectPosition.coordinate') ]
['37.15423700 139.95680000',
 '37.13917000 140.03105800',
 '37.13862100 139.94537100',
 ...]

10万5010点の座標が得られた。手作業でやっていたら宇宙が熱的死するところだった。危ない危ない。

栃木県（平成22年）

こういう図の作成も Python で行える（matplotlib というライブラリを使う）。Python を使うとたいていのことはできる。

地元の方はお気づきだろうが、地図データが少し古いため合併でなくなった町がいくつか存在する。今回は県の形状がテーマなので問題ない。というか市町村の境界データは余計なので除去しなければならない。

数値地図のファイルを端から見ていくと、座標データはいくつかのグループに分かれている。「一番最初が県境データで、その後が市町村境界に違いない」とあたりを付けて id="c00001" という属性があるものだけ絞り込む。

border1 = [point.text for point in root.find('.//*[@id="c00001"]').findall('.//DirectPosition.coordinate')]
栃木県民向けシルエットクイズ

なんだこれ。栃木県の形ではない。地図と見比べると那須塩原市らしい。「黒磯駅があるところ」と言えば18きっぷ旅行者には通じやすいだろう。しかし今回欲しいのは栃木県境なのである。

ちゃんと仕様書を読めば「県境データはここに格納されています」と書いてあるかもしれないが、役所データの仕様書・Wi-Fi の利用規約・中学時代の小説というのは日本三大読みたくない文章である。他の方法を考えよう。

筆者はIQが53万あるのでうまい手が思いついた。那須塩原市の境界の北側は福島県との県境であり、西側は日光市との境界である。しかしこのデータには日光市側の境界線データも入っていて、それは那須塩原市の境界とぴったり一致する。

つまり、市町村境界データの重複している点を除けば、県境だけを取り出せるはずだ。やってみよう。

prefborder = [k for k, v in collections.Counter(border).items() if v == 1]
ユナイテッド栃木

できた。あまりに想定通りに行ったのでびっくりした。

さて、ここからが問題だ。県境データができたからといって、どこからどこまでが「県内」かを識別するのはそこまで簡単ではない。多角形の内外判定アルゴリズムを調べる。

「外側と直線で結ぶ際に境界線を奇数回通過すれば内側」らしい。なるほど〜。ただ「境界を通過」の判定が大変そうなのでやめる。もっと単純にマップ全体をラスタ化することにする。栃木県全体を緯度経度0.01°ごとの正方形マスに区切り、県境のあるマスを黄色に、そうでないマスを紫に塗る。

# 解像度を落とす
meshborder = np.unique((npprefborder * 100).astype(int), axis=1).T

# 原点に寄せる
meshborder[:,0] -= np.min(meshborder[:,0])-1
meshborder[:,1] -= np.min(meshborder[:,1])-1

pmap = np.zeros((np.max(meshborder[:,0])+2, np.max(meshborder[:,1])+2), dtype=int)
for bp in meshborder:
    pmap[bp[1], bp[0]] = 3

plt.imshow(pmap, origin={'lower', 'left'})
ラスタライズド栃木

丸さの判定ならこの解像度で十分だろう。あとは一番左上を「県外」とし、県外に隣接する県境でない点も「県外」とみなす。いわゆるバケツツールの要領で「県外」を決め、残った部分を「県内」にする。

# 左上隅は県外とみなす
pmap[0, 0] = 1
count = 0

xmax = pmap.shape[0]-1
ymax = pmap.shape[1]-1

# 県内0、県外1、県境3 として、隣マスとの差が1のマスがなくなるまで塗り続ける

while np.any(np.sum(np.abs(np.diff(pmap, axis=0)) == 1)) or np.any(np.sum(np.abs(np.diff(pmap, axis=1)) == 1)):
    count += 1
    it = np.nditer(pmap, flags=['multi_index'])
    while not it.finished:
        idx = it.multi_index
    
        if pmap[idx] == 1:
            if idx[0] < xmax and pmap[idx[0]+1, idx[1]] == 0:
                pmap[idx[0]+1, idx[1]] = 1
            if idx[0] > 0 and pmap[idx[0]-1, idx[1]] == 0:
                pmap[idx[0]-1, idx[1]] = 1
            if idx[1] < ymax and pmap[idx[0], idx[1]+1] == 0:
                pmap[idx[0], idx[1]+1] = 1 
            if idx[1] > 0 and pmap[idx[0], idx[1]-1] == 0:
                pmap[idx[0], idx[1]-1] = 1
        it.iternext()
    
plt.imshow(pmap, origin={'lower', 'left'})
アイソレーテッド栃木

このあたりでメルカトル図法の問題に気づいたが、これについては本編記事で触れたので省略。

栃木県は形状が単純なのでこの方法で問題ないが、海岸線のややこしい県では内側の海が「県内」と判定されてしまう。長崎県の大村湾などがその例である。

バケツツールの宿命的問題

大村湾と佐世保湾の海峡は170mしかないので、これを通過できるほどの解像度にすると計算量が爆発してしまう。こういう海は手作業で「県外」ポイントを置いていくことにする。あと、そもそも県境のデータがない富士山にも手動で県境をつけていく。ダルい。

if pref == 'aomori': # 十和田湖の県境
    pmap[143, 28:33] = 3
    pmap[142, 33:38] = 3
if pref == 'akita': # 十和田湖の県境
    pmap[145, 250:255] = 3
    pmap[146, 244:250] = 3
    pmap[33, 38] = 3
if pref == 'ishikawa': # 能登島南部に海判定を入れる
    pmap[94,161] = 1
    pmap[84,164] = 1
if pref == 'shizuoka': # 富士山境界
    pmap[170:177,131] = 3
# 退屈なことを Python がやってくれるとは限らない

どうにか県外・県内・県境の塗り分けができたので、あとは面積を計算して円を描いて内部の「県内」「県境」マスを数える。今から考えると「県境」は1/2マスで判定するべきだったがこの時は1マスで数えた。

内包率が最大になる中心点を見つけるアルゴリズムは、考えるのが面倒だったので端から走査することにした。これだと1県に1分くらいかかるが、他の県のデータをダウンロードしたり整形作業をすることで時間が潰せるのであまり問題なかった。

サークルド栃木

できた。ただ視覚的にダサいのでもう少しかっこよくする。matplotlib はグラフ描写ツールなのでデータごとに違う色を使うが、デザインの基本は無駄な色彩を減らすことである。日本語のテキストを入れるのはちょっと手間なので Impact フォントでローマ字にする。

完成

栃木県の画像に「栃木」と書いてるのは逆に新鮮である。あとは数値を Excel に並べてソートし、順位を決定した。

ただ、沖縄と東京に関しては島嶼部が遠すぎ、マップが広くなりすぎて計算量が爆発したので本土の部分だけを手動で切り出すことにした。このためまた XML ファイルと格闘することになるのだが、その話はもはや書くのもダルいので省略する。

以上。次回は標高データを使いたいので「日本でいちばん正規分布っぽい山は香川県の飯野山なのか？」をやろうと思います。

飯野山 by odysseygate CC BY-SA 2.0


文章で生計を立てる身ですのでサポートをいただけるとたいへん嬉しいです。メッセージが思いつかない方は好きな食べ物を書いてください。