【MATLAB】地域メッシュ統計を使った市町村ポリゴンデータの描画

MATLABのMapping toolboxを使用してshapefileを読み込み，任意の市町村のポリゴンデータを地図上に描画する方法を紹介します。

地域メッシュ統計とは

地域メッシュ統計とは，緯度・経度に基づき地域をメッシュ区域にわけて編成されたデータです。詳細は総務省統計局のこちらのページにまとめられています。地域メッシュの歴史などをより詳しく知りたい方はこちらを参照してください。

地域メッシュの区分

地域メッシュは，メッシュサイズに応じて以下の5区分に分けられています。メッシュ一辺の最小単位は約250m。日本のメッシュの区分図はこちらのpdfを参照してください。

市町村別メッシュコード

総務省統計局がこちらのページで都道府県別にメッシュコードをcsv形式で公開をしています。このメッシュコードは後述する地域メッシュデータをダウンロードする際に使用します。csvファイルは，都道府県市区町村コード（5桁），市区町村名，基準地域メッシュコード（3次地域メッシュ）で構成されています。愛知県を例に具体的なファイルの中身を見ていきます。A列が都道府県市区町村コード（5桁），B列が市区町村名，C列が基準地域メッシュ・コード（上表の第3次地域区画）です。

地域メッシュデータのダウンロード

メッシュデータは，政府統計のオープンデータサイト e-Stat で公開されています。こちらがメッシュデータの公開先ページです。データはShapefile形式，KML形式，GML形式から選択でき，座標データも緯度経度か平面直角座標系かを選択することができます。MATLABで扱うのであれば，Shapefile が簡単だと思います。「世界測地系緯度経度」であれば，緯度経度でプロットを作成することができますが，世界測地系座標系だとmに変換された絶対距離でデータが描画されます。今回は，「世界測地系緯度経度・Shapefile」を使用します。データは都道府県単位で検索することができるので，5次メッシュ・愛知県で絞り込みます。

そうすると，M5136からM5337の6つのデータが表示されます。

今回は豊田市のメッシュデータを描画したいと思います。01 でダウンロードした市町村別メッシュコードのcsvファイルから「豊田市」を選択します。C列の基準メッシュ・コードの左の4桁と一致するもの，つまりM5237のファイルをダウンロードします。ファイルはzip形式になっていてこれを解凍すると4つのファイルを取得することができます。ただし，このままではまだデータは「豊田市」に絞り込まれておらず，追加の処理が必要です。その処理について説明していきます。

市町村単位のshpファイルの描画

引き続き，豊田市を例に市町村単位のメッシュデータをMATLABを使用して描画する方法を紹介します。

用意するデータ

• ダウンロードしたshapeファイル一式

• 豊田市に絞り込んだ市町村メッシュコードファイル：先ほどのcsvファイルを豊田市にフィルタリングしたもののカラム名を英語にして別ファイルに保存しています。ファイル名は，toyota.xlsx です。ファイルの中身はこんな感じ。

shapefileの読み込み

MATLABのmapping toolboxを使用します。このtoolboxを導入すると使用できる readgeotable関数がとても優秀です。readgeotable関数を実行すると，.shpファイルと関連ファイルを読み込み，ポリゴンと関連情報をテーブルとして返してくれます。テーブルの各カラムはこちらの定義書に定義されています。

% shape fileの読み込み</mark>
shp = readgeotable("MESH05237.shp");

実行結果

ちなみに読み込んだこの shp は，geoplot関数で簡単に描画できます。試しに描画してみたのが左のマップです。なんだか黒い四角が描画されました。これを拡大していくと右のマップになります。ちゃんと四角形で領域が区切られたポリゴンファイルが描画されているのがわかります。5次メッシュなので，この小さな四角形の一辺が約250mということになります。

figure()
geoplot(shp)

実行結果

豊田市のポリゴンデータのみ抽出する

まずは，豊田市に絞り込んだ市町村別メッシュコードファイルを読み込みます。

toyota = readtable("toyota.xlsx");

実行結果

shp のうち，toyota.mesh_code と shp.KEY_CODE の左8文字とが一致するものが，豊田市のポリゴンということになります。ただし，このtoyota.mesh_code と shp.KEY_CODEの左8文字は，3次メッシュコードなので，解像度は1km四方ということになります。そのため厳密には，市町村の境界付近では，この3次メッシュの中に複数の市町村が含まれる可能性があり，完璧に豊田市だけを抜き出しているわけではないことは注意が必要です。あくまで，「豊田市と一部その周辺市町村」のポリゴンデータであるということになります。

toyota.mesh_code と shp.KEY_CODEをそれぞれkeyにして，toyotaテーブルとshpテーブルとを連結して，豊田市のポリゴンデータを得ます。

shp.KEY_CODE_double = str2double(extractBefore(shp.KEY_CODE,9)); % KEY_CODE の左から8文字のみ抽出し，doubleに変換
toyota_shp = innerjoin(toyota,shp,'LeftKeys','mesh_code','RightKeys','KEY_CODE_double'); 

実行結果

もともと shp に入っていた以外のカラムは，描画する際に不要なのでテーブルから削除します。

toyota_shp = removevars(toyota_shp, {'pref_code','city','mesh_code'});

実行結果

これで豊田市（と一部その周辺市町村）のポリゴンファイルを抽出することができました。

試しにプロットしてyahoo地図と比較してみます。

figure()
geoplot(toyota_shp)

実行結果

無事に抽出できていることがわかります。実際の豊田市の地図と比較しても，うまく抽出できていることが確認できます。

コードサマリ

% shape fileの読み込み
shp = readgeotable("MESH05237.shp");

% shape fileの描画
figure()
geoplot(shp)

% 豊田市の地域メッシュコードデータ読み込み
toyota = readtable("toyota.xlsx");

% 豊田市のポリゴンデータの抽出
shp.KEY_CODE_double = str2double(extractBefore(shp.KEY_CODE,9)); % KEY_CODE の左から8文字のみ抽出し，doubleに変換
toyota_shp = innerjoin(toyota,shp,'LeftKeys','mesh_code','RightKeys','KEY_CODE_double'); 
toyota_shp = removevars(toyota_shp, {'pref_code','city','mesh_code'}); % 不要なカラムを削除

% 豊田市のポリゴンデータの描画
figure()
geoplot(toyota_shp)

まとめ

豊田市を例に市町村別にshapefileのデータを描画する方法を紹介しました。ほかの市町村でも基本的なやり方は一緒です。

1. 総務省統計局の都道府県別のメッシュコードファイルから，抽出したい市区町村のファイルを作る
2. 該当の市区町村が含まれるshape file をダウンロードする
3. 1.で作成した市区町村のメッシュコード情報と2.で取得したshape fileから，該当市区町村のポリゴンデータのみを抽出する