シェープファイルの読み込みと可視化
日本政府が公開している位置情報付きのデータのShapefile (シェープファイル)を扱って何か解析したいと思い立って、今回トライしたことを記事に書きたいと思います。
最初はPython でやろうとしたのですが、シェープファイルの読み込み用とプロット用にそれぞれ別のモジュールのインストールが必要で、しかも何種類かあり、前提モジュールも結構あって悩みました。
私は企業に属している身なので、ツールを選定するにも、メンテナンスがちゃんとされているか(個人ではなくチームでメンテナンスをしているか、脆弱性が発覚したときにパッチの提供があるか)とか、ライセンス的に問題ないか(法務部門が規約を読んでOK を出せるか、社員数で有償版を使わないといけなくならないか)とか、考えて選ぶようにしてます。
今回はシェープファイルを扱うのに必要なPython のモジュールを調査して実際に試してみたのですが、1時間ぐらい掛けてShapefile 読み込みまではできたのですが、その後可視化するところでまたモジュール選定から立ち戻ることになり、一旦白紙にしました。今回はサクッとMATLAB でやったので、またトライしたら記事にしたいと思います。
e-Stat から国勢調査の統計データをダウンロード
さて、今回扱いたいデータは、政府統計の総合窓口「e-Stat」 からダウンロードできる国勢調査の人口データを含むシェープファイルです。
それではe-Stat からデータをダウンロードします。詳細は後述しますが、e-Stat のデータの直リンクは以下のとおりです。
【e-Stat】統計GISデータダウンロード 2015年国勢調査 世界測地系緯度経度・Shapefile 東京都のリンク
新しいデータのほうが良いのですが、国勢調査の2015年では「小地域 (町丁・字等別)」というのがあるのですが、2020年のものでは現時点(2022年4月)ではまだありません。人口集中地区のみです。
なので2015年のデータをダウンロードします。
さらに投影座標の種類とファイルフォーマットによってデータが分かれています。
迷いますよね。ここでは「世界測地系緯度経度」にしてファイルフォーマットは「Shapefile」を選びます。
続いて対象地域の都道府県を選択しますが、ここでは「東京都」とし、「13000 東京都全域」の「世界測地系緯度経度・Shapefile」の水色ボタンをクリックすると、データがZip 形式でダウンロードされます。
ダウンロードしたZip ファイルを解凍すると、Shapefile 一式 (.dbf、.prj、.shp、.shx)が含まれています。
MATLAB のreadgeotable でシェープファイルを読込
MATLAB のMapping Toolbox を使います。
以前の記事でApple Watch のヘルスケアデータをエクスポートしてXML ファイルとワークアウトが格納されたGPX ファイル をMATLAB のreadgeotable とgeoplot でたった数行で可視化しましたが、最近のMATLAB は位置情報の扱いが本当に楽になりました。
今回はバージョンR2022a を使います。
まずはreadgeotable でシェープファイルを読み取ります。R2021b からMapping Toolbox に登場した関数で、シェープファイル、GPX ファイル、GeoJSON などの位置情報データを簡単に読み込めます。
t = readgeotable('h27ka13.shp');Code language: Matlab (matlab)シェープファイルを読み取るときは、.shp の拡張子を読み取れば自動的に付随するデータも読み取ります。古いバージョンから使えるshaperead という関数もあるのですが、シェープファイルを構造体として読み取るので、その後扱いやすいテーブル型に変換しないといけなかったりして、読み取りに少々時間が掛かってしまいますので、readgeotable がオススメです。
【MATLABドキュメント】readgeotable
こちらが実行後のワークスペースです。
東京都全域のデータが6010行×36列のテーブルとして読み取れました。東京都の区市町村の丁目毎にデータが分かれています。PREF_NAME やCITY_NAME などの日本語文字も文字化けせずに読めていますね。
白地図に東京都の丁目毎の人口をプロット
続いて地図にプロットしていきます。オープンソースだと別のモジュールが必要な場合がありますが、MATLAB では同じくMapping Toolbox で行けます。
ここで使うのはgeoplot という関数。以前の記事でGPX ファイルをプロットするときにも使ったのですが、あのときはMATLAB 本体のgeoplot (ドキュメントページはこちら)を使いました。今回使うのは、Mapping Toolbox のgeoplot (ドキュメントページはこちら)です。
実はR2022a からMapping Toolbox にMATLAB 本体のgeoplot を拡張する機能が追加されたのです。これを使うと、白地図の上にジオポリゴンを重ねることができます。
この3行でできます。
geoplot(t, ColorVariable="JINKO");
cb = colorbar;
cb.Label.String = "人口";Code language: Matlab (matlab)こちらが実行結果のFigure です。見やすくするため島しょ部を除いて表示しています。
geoplot のColorVariable のオプションでどのカラムをカラーバーの色として使用するか選択できます。
シェープファイルを簡単に読み込めて可視化できて、本当に便利になりましたね。
No responses yet