【Python】新潟県に一番似ている国は○○？各都道府県に似た形の国を調べてみた

はじめに

この記事はNTTドコモアドベントカレンダーの4日目の記事です。

こんにちは、NTTドコモサービスイノベーション部の福島です。

こちらの日本地図、おかしなところがあります。気づきますでしょうか。正解は少し下にあります。

話は変わりますが、以前テレビを見ていると「日本地図の四国をオーストラリアに変えても気づかないのでは？」という検証をやっていました。

また、パスタが名物の群馬県高崎市の特集番組では、市の形もパスタの本場のイタリアと似ているというトリビアが紹介されていました。

名前は知っていても国の形までは知らない国って結構あるな、各都道府県の形に似ている国って知らないだけで実はあるんじゃないかな、と気になったので調べてみます。

冒頭の日本地図ですが、正解は九州の各県が別の国のシルエットに置き換えられていました。

実現方法

図形の形の類似度を数値化するライブラリが無いか調べてみると、最も有名な画像処理ライブラリであるOpenCVにMatch Shapesというピッタリな関数が用意されていました。

こちらは図形の拡大縮小や回転も考慮して比較を行ってくれるため、回転させた際に最も形が似ている都道府県と国をマッチさせてくれそうです。

国と都道府県の輪郭のデータさえ用意出来れば、この関数を使って各都道府県と各国の形の類似度が計算でき、各都道府県に対して最も形の似ている国が導けそうです。

国境・都道府県境のデータの取得

少し調べてみると、地理情報をデータとして扱うシステムはGIS（Geographic Information System）と呼ばれ、図形や座標の情報の格納にはシェープファイルというフォーマットが主流で使われているようです。後述する総務省のホームページから取得した各都道府県の輪郭のデータもこの形式でした。シェープファイルの概要についてはこちらやこちらの記事が詳しいです。

Pythonでこのシェープファイルを扱う場合、GeoPandasというライブラリを用いるのが主流のようです。データ分析界隈の皆さんならよくご存知のPandasの派生のようで、Pandasに慣れた人なら直感的に使えそうです。

国境データの取得

こちらのWorld Borders Datasetをダウンロードして利用しました。

各国の輪郭はGeoPandasのライブラリにデフォルトのデータセットとしても入っていたのですが、こちらはあまりにも頂点数が少なく簡易な形状だったのでWorld Borders Datasetを利用しています。

dataフォルダの配下にダウンロードしたshpファイルを格納し、下記で輪郭を取得しました。後述の都道府県データでも同じなのですが、複数の島からなる国が多数あるのですが、単一の図形同士の類似度を数値化するというアルゴリズムの制約から最も面積が大きい陸地1つを計算に採用しています。

import glob

import numpy as np 
import pandas as pd
import geopandas as gpd
from shapely import geometry
from tqdm import tqdm

shape_file = "data/TM_WORLD_BORDERS-0.3.shp"
gdf = gpd.read_file(shape_file)

dic_country = {}
for index, row in gdf.iterrows():
    country_name = row['ISO3']
    border = row['geometry']

    if isinstance(border, geometry.polygon.Polygon):
        points = border.exterior.coords[:-1] # 最後の点は最初の点と重なるので省いて良い
    elif isinstance(border, geometry.multipolygon.MultiPolygon): 
        area_list = [(b.area, b) for b in border.geoms]
        area_list = sorted(area_list, key=lambda x: x[0], reverse=True)
        largest_poly = area_list[0][1] # 複数の領域からなる国は面積が最大の領域を採用
        points = largest_poly.exterior.coords[:-1] # 最後の点は最初の点と重なるので省いて良い
    dic_country[country_name] = points

np.save('dic_country.npy', dic_country)

県境データの取得

都道府県の形状については、政府の統計ポータルサイトであるe-Statから取得しました。上記の国境データとほとんど同じ形式のため、輪郭を抽出するプログラムもほぼ同じです。1都道府県=1ファイルなので47ファイル分をダウンロードするのが少し大変でした。

import glob

import numpy as np 
import pandas as pd
import geopandas as gpd
from shapely import geometry
from tqdm import tqdm

dir_shape = ('./todoufuken') # e-Statから取得した境界データが格納されたフォルダ
shape_files = glob.glob(f'{dir_shape}/*.shp')

# 各都道府県のshpファイルの読み込み→都道府県名と境界の頂点群を保存
japan = None
for shape_file in tqdm(shape_files[1:3]):
    gdf = gpd.read_file(shape_file)
    pref = gdf[['PREF_NAME', 'geometry']].dissolve(by='PREF_NAME', aggfunc='sum')
    if japan is None:
        japan = pref
    else:
        # GeoDataFrameの結合
        japan = gpd.GeoDataFrame(pd.concat([japan, pref]), crs=japan.crs)
japan['PREF_NAME'] = japan.index

# 複数の領域（離島など）からなる都道府県に対する処理
dic_border = {}
for index, row in japan.iterrows():
    pref_name = row['PREF_NAME']
    border = row['geometry']
    print(border)
    input()
    if isinstance(border, geometry.polygon.Polygon):
        points = border.exterior.coords[:-1] # 最後の点は最初の点と重なるので省いて良い
    elif isinstance(border, geometry.multipolygon.MultiPolygon): 
        area_list = [(b.area, b) for b in border.geoms]
        area_list = sorted(area_list, key=lambda x: x[0], reverse=True)
        largest_poly = area_list[0][1] # 複数の領域からなる都道府県は面積が最大の領域を採用
        points = largest_poly.exterior.coords[:-1] # 最後の点は最初の点と重なるので省いて良い
    dic_border[pref_name] = points
    print(points[:10])
    input()

assert len(dic_border)==47 # 47都道府県
np.save('dic_prefecture.npy', dic_border)

類似度の計算

ここまでで国境の形状、都道府県境の形状がnumpyファイルとして保存されました。

あとはこれを読み込んで類似度の総当りをする下記のプログラムを動かすだけです。

import numpy as np
import cv2

dic_prefecture = np.load('dic_prefecture.npy', allow_pickle=True).item()
dic_country = np.load('dic_country.npy', allow_pickle=True).item()

for pref_name, pref_points in dic_prefecture.items():
    pref_cnts = []
    for p in pref_points:
        pref_cnts.append(list(map(int, map(lambda x: x*100, p)))) # 緯度経度をそのままintにすると丸め誤差が大きすぎるので100倍する
    pref_cnts = np.array(pref_cnts)
    pref_cnts -= np.min(pref_cnts, axis=0) # 原点基準に移動。デバッグ時に図形を表示する時用

    min_score = 10000.
    min_country = ""
    for country_name, country_points in dic_country.items():
        country_cnts = []
        for p in country_points:
            country_cnts.append([list(map(int, map(lambda x: x*100, p)))])
        country_cnts = np.array(country_cnts)
        country_cnts -= np.min(country_cnts, axis=0) # 原点基準に移動。デバッグ時に図形を表示する時用

        retval = cv2.matchShapes(pref_cnts, country_cnts, 1, 0.0)
        if retval < min_score:
            min_score = retval
            min_country = country_name
            
    print(pref_name, min_country, min_score)

結果の表はこちらです。

特に類似度が高かった（距離が小さかった）ペアのトップ5は太字にしています。

都道府県	最も形が似ている国・地域	輪郭の距離（小さいほど似ている）
北海道	イギリス	0.119
青森県	マリ	0.146
岩手県	ニウエ	0.031
宮城県	モルドバ共和国	0.051
秋田県	セルビア	0.058
山形県	ウガンダ	0.026
福島県	ブルキナファソ	0.058
茨城県	モルドバ共和国	0.110
栃木県	コートジボワール	0.009
群馬県	ケニア	0.055
埼玉県	ウォリス・フツナ	0.085
千葉県	ナミビア	0.103
東京都	マラウイ	0.141
神奈川県	ペルー	0.073
新潟県	ジョージア	0.327
富山県	マルティニーク	0.082
石川県	モルディブ	0.940
福井県	フィリピン	0.221
山梨県	ブルンジ	0.038
長野県	北マリアナ諸島	0.053
岐阜県	イラク	0.073
静岡県	西サハラ	0.129
愛知県	ボスニア・ヘルツェゴビナ	0.064
三重県	パキスタン	0.157
滋賀県	コートジボワール	0.016
京都府	モンゴル	0.060
大阪府	中華人民共和国	0.109
兵庫県	イラク	0.034
奈良県	サウジアラビア	0.036
和歌山県	トケラウ	0.069
鳥取県	パナマ	0.461
島根県	イスラエル	0.614
岡山県	コートジボワール	0.010
広島県	セントルシア	0.026
山口県	トルクメニスタン	0.108
徳島県	バハマ	0.054
香川県	レバノン	0.062
愛媛県	トンガ	0.066
高知県	マラウイ	0.533
福岡県	アイルランド	0.067
佐賀県	ベネズエラ・ボリバル共和国	0.062
長崎県	フォークランド（マルビナス）諸島	0.119
熊本県	ジブチ	0.024
大分県	リビア	0.087
宮崎県	中央アフリカ共和国	0.015
鹿児島県	バングラデシュ	0.122
沖縄県	ラトビア	0.248