Pythonでシェープファイルから都道府県境界図を作成する

PythonでGIS用シェープファイルを可視化するでは、PythonライブラリのGeoPandasとMatplotlibでシェープファイルを読みこんで日本地図を表示した。今回は、シェープファイルをダウンロードして、何かと使うことのありそうな都道府県境界図をPythonで作成してみる。

環境

Windows10(1903)のWSL(Ubuntu 18.04)とJupyter Notebookを使用。
$ lsb_release -dr Description: Ubuntu 18.04.3 LTS Release: 18.04 $ python3 -V Python 3.6.9 $ pip3 show pandas Name: pandas Version: 0.25.3 ... $ pip3 show matplotlib Name: matplotlib Version: 3.1.2 ... $ jupyter --version 4.4.0

境界図のシェープファイル

都道府県境界図を作成するためにシェープファイルをダウンロードする。国土数値情報　ダウンロードサービスや地図で見る統計(統計GIS)からダウンロードできる。国土数値情報からダウンロードできるシェープファイルだと全国のデータが1ファイルにまとまっているので作成が簡単だと思ったのだが、都道府県境図にしようとしても市町村区境の図になってしまった。そこで、今回は地図で見る統計(統計GIS)から都道府県ごとのシェープファイルをダウンロードして、それを連結して都道府県境図を作成する。

シェープファイルのダウンロード

ダウンロードするシェープファイルは地図で見る統計(統計GIS)の境界データ。今回は以下の場所にあるファイルを使う。

境界データダウンロード＞小地域＞国勢調査＞2015年＞小地域
＞世界測地系緯度経度・Shapefile

都道府県ごとに全域のファイルがあるのでそれをダウンロードする。全47ファイルをダウンロードしたら、zip圧縮されているのですべて解凍して同じディレクトリに保存しておく。

また、以下の場所から定義書がダウンロードできるので、これをみるとどんな項目があるか確認できる。

境界データダウンロード＞小地域＞国勢調査＞2015年＞小地域

今回使うのは次の3項目。

• PREF_NAME 都道府県名
• AREA 面積(㎡)
• JINKO 人口

ちなみに、北海道全域のファイルをそのまま表示すると以下のような国勢調査の小地域の境界図となる。

GeoPandasで都道府県境図を作成

まずはGeoPandasとそのほか必要なパッケージなどをインストールする。
$ sudo apt install libspatialindex-dev $ pip3 install geopandas descartes rtree mapclassify

インストールされたバージョン。
$ pip3 list | grep -e geopandas -e descartes -e Rtree -e mapclassify descartes (1.1.0) geopandas (0.6.2) mapclassify (2.2.0) Rtree (0.9.3)

準備ができたので、GeoPnadasで都道府県境界図を作成する。手順としては、ダウンロードした都道府県ごとのシェープファイルは小地域境のデータなので、これを都道府県境のデータにする。GISではこういう処理をディゾルブと言うらしく、GeoPandasではそのためのメソッドが用意されている（Aggregation with dissolve）。ここでは都道府県境界図を作成するのでPREF_NAME（都道府県名）でディゾルブする。また、人口などの数値はaggfuncで指定した集計がされる。市区町村境界図を作成する場合はCITY_NAME（区町村名）でディゾルブすればよい。ディゾルブしたデータを連結していき、最終的に日本地図とする。これをシェープファイルとして保存。以下はそのPythonコード。

%matplotlib inline
import matplotlib.pyplot as plt

import os
from pathlib import Path
import glob

import pandas as pd
import geopandas as gpd

# シェープファイル保存ディレクトリ
SHAPE_DIR = os.path.join(str(Path.home()), 'estat/sgis_shape')

def create_outline_map():
    print('*** Creating outline map ***')
    print('Loading from {}'.format(SHAPE_DIR))

    jpn = None

    shpfiles = glob.glob(os.path.join(SHAPE_DIR, '*.shp'))
    for shp_path in shpfiles:
        gdf = gpd.read_file(shp_path)

        # byでディゾルブする単位を指定する
        # byで指定したカラムがindexになる
        # AREAとJINKOは後でコロプレス図で使用
        pref = gdf[['PREF_NAME', 'AREA', 'JINKO', 'geometry']].dissolve(by='PREF_NAME', aggfunc='sum')

        if jpn is None:
            jpn = pref
        else:
            # GeoDataFrameの結合
            jpn = gpd.GeoDataFrame(pd.concat([jpn, pref]), crs=jpn.crs)

    # シェープファイルにはindex情報は保存されないので都道府県名をカラムにしておく
    jpn['PREF_NAME'] = jpn.index

    print()
    print(jpn.info())
    print(jpn.head())
    
    # shapeファイル保存
    jpn.to_file('jpnmap.shp', encoding = 'shift-jis')

    # 境界図のプロット
    fig, ax = plt.subplots(figsize = (20,16)) 
    jpn.plot(ax=ax, facecolor='white', edgecolor='black', linewidth=0.5)
    ax.set_axis_off()
    plt.axis('equal')
    plt.show()

if __name__ == '__main__':
    create_outline_map()

以下のような都道府県境界図が作成された。

コロプレス図を作成する

ダウンロードしたファイルには面積と人口のデータも含まれている。せっかく境界図のシェープファイルを作成したので、人口密度を計算してコロプレス図を作成する。

コードは以下の通り。

%matplotlib inline
import matplotlib.pyplot as plt

import geopandas as gpd

def plot_choropleth():
    # 作成した都道府県境界図のシェープファイルの読み込み
    jpn = gpd.read_file('jpnmap.shp', encoding = 'shift-jis')

    # AREA（面積）の単位は平方メートルなので平方キロメートルに直して人口密度を計算
    jpn['DENSITY'] = jpn['JINKO'] / (jpn['AREA'] / (10**6))

    print()
    print(jpn.info())
    print(jpn.head())

    fig, ax = plt.subplots(figsize = (12,8))

    # 人口密度のコロプレス図
    jpn.plot(column = 'DENSITY', edgecolor = "black",
                    scheme='quantiles', cmap='YlOrRd', ax=ax, 
                    legend = True)

    ax.set_axis_off()
    plt.axis('equal');
    plt.show()

if __name__ == '__main__':
    plot_choropleth()

次のようなコロプレス図が作成された。