2020年1月28日火曜日

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

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


環境


Windows10(1903)のWSL(Ubuntu 18.04)とJupyter Notebookを使用。



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


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


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


ダウンロードするシェープファイルは地図で見る統計(統計GIS)の境界データ。今回は以下の場所にあるファイルを使う。
境界データダウンロード＞小地域＞国勢調査＞2015年＞小地域
＞世界測地系緯度経度・Shapefile
都道府県ごとに全域のファイルがあるのでそれをダウンロードする。全47ファイルをダウンロードしたら、zip圧縮されているのですべて解凍して同じディレクトリに保存しておく。

また、以下の場所から定義書がダウンロードできるので、これをみるとどんな項目があるか確認できる。
境界データダウンロード＞小地域＞国勢調査＞2015年＞小地域
今回使うのは次の3項目。
  • PREF_NAME 都道府県名
  • AREA 面積(㎡)
  • JINKO 人口

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


GeoPandasで都道府県境図を作成


まずはGeoPandasとそのほか必要なパッケージなどをインストールする。


インストールされたバージョン。


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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
%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()


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


コロプレス図を作成する


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

コードは以下の通り。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
%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()

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

投稿者 時刻:
ラベル: , , ,

0 件のコメント:

コメントを投稿

登録: コメントの投稿 (Atom)