はじめに
この記事は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 |
なんとトップ5のうち3つ(1位:栃木県、2位:岡山県、4位:滋賀県)をコートジボワールが独占しました。
輪郭が丸や四角に近く、目立った凹凸が無い都道府県の多くはコートジボワールが良いスコアを叩き出していそうです。
続いて形が特徴的な都道府県について、いくつか結果を見てみます。
- 新潟県(ジョージア)
- 確かに結構似ています
北海道(イギリス)
- 全く似ていません。さすがに北海道に形が似ている国は無さそうでしょうか
青森県(マリ)
- 下北半島の形が思い出せなくて諦めた人が描いた青森県のようです
- 東京都(マラウイ)
- これも結構似ていると思います。ちなみに高知県もマラウイでした
終わりに
GeoPandasとOpenCVを用いて各都道府県に似ている形状の国を紐づけることが出来ました。
リロングウェ(マラウイの首都)は東京都でいうと町田市のあたりになりそうです。
では、私は各都道府県を別の国で置き換えたエセ日本地図を完成させてきます。お読み頂きましてありがとうございました。
参考資料・出典
- 都道府県・国のイラストはillust image様より画像を利用させて頂きました
- 高崎市のイラストは素材Library.com様より画像を利用させて頂きました
- 日本地図のシルエットは日本地図様より画像を利用させて頂きました
- GeoPandasの使い方はこちらを参考にさせて頂きました