国土交通省の国土数値情報 ダウンロードサービスで公開されている国土数値情報 500mメッシュ別将来推計人口（H29国政局推計）のデータは Shape 形式 （文字エンコードは Shift-JIS） ですので、GeoJSON 形式に変換します。 変換には PostGIS を使います。
PostGIS 環境は Docker で用意します。 Docker のイメージは mdillon/postgis を使います。 PostgreSQL のバージョンは 10.3、PostGIS のバージョンは 2.4 です。 事前に Docker を利用できるようにしておいてください。Jupyter を実行しているユーザーを docker グループに追加しておくと良いでしょう。
ノートブック内で Python から PostgreSQL にアクセスするには Psycopg を使います。
In [1]:
!docker run --name postgis -d -p 5432:5432 mdillon/postgis:10
Python の PostgreSQL ライブラリである psycopg2 をインストールします。
In [2]:
!pip install psycopg2-binary
デーモンとして起動したコンテナに接続し、PostGIS や各種ライブラリのバージョンを確認します。
version() と postgis_full_version() 関数で確認できます。
In [3]:
import pprint
import psycopg2
conn = psycopg2.connect("host=127.0.0.1 user=postgres")
cur = conn.cursor()
cur.execute('SELECT version(), postgis_full_version()')
pprint.pprint(cur.fetchone())
「国土数値情報 500mメッシュ別将来推計人口（H29国政局推計）」を利用します。 サイトで地域を選択し、利用条件を読んで利用用途のアンケートに記入してからファイルをダウンロードします。
| 項目 | 内容 |
|---|---|
| ファイル名 | m500-17_GML.zip |
| ファイル容量 | 108.58MB |
| 年度 | 平成29年 |
| 測地系 | 世界測地系 |
| 地域 | 全国 |
ZIP ファイルには Shape ファイルが含まれますので、拡張子が ".shp" のものと、その他に ".dbf" などいくつかのファイルがあります。 ZIP ファイルは /tmp フォルダに置いてあるものとします。 展開したファイルを合計すると760MB程度になります。
In [4]:
import zipfile
with zipfile.ZipFile('/tmp/m500-17_GML.zip') as myzip:
print(myzip.namelist())
myzip.extractall()
PostGIS パッケージに含まれる shp2pgsql を使って Shape ファイルを PostGIS 向けの SQL に変換します。
SQL でのテーブル名は mesh4_pop_00 になります。
テーブル名を指定するには shp2pgsql の第二引数に与えます。
In [5]:
!docker run --rm -v `pwd`:/tmp mdillon/postgis:10 \
shp2pgsql -W Shift_JIS /tmp/Mesh4_POP_00.shp 2>/dev/null | gzip > /tmp/mesh4_pop_00.sql.gz
SQLを実行します。Shape ファイルからの変換結果は膨大な SQL 文になりますので、個別にクエリを実行するように調整します。 カーソル変数 cur は上で定義済みとします。
In [6]:
import gzip
import io
buffer = io.StringIO()
with gzip.open('/tmp/mesh4_pop_00.sql.gz', 'rt') as fp:
for line in fp:
s = str(line).strip()
buffer.write(s)
if s.endswith(';'):
cur.execute(buffer.getvalue())
buffer = io.StringIO()
テーブルに登録されたデータ件数などを集計します。
In [7]:
cur.execute('''
SELECT COUNT(*),
COUNT(DISTINCT mesh_id), MIN(LENGTH(mesh_id::VARCHAR)), MAX(LENGTH(mesh_id::VARCHAR)),
COUNT(DISTINCT SUBSTR(LPAD(mesh_id::VARCHAR, 9, '0'), 1, 4)),
COUNT(DISTINCT city_code), MIN(LENGTH(city_code::VARCHAR)), MAX(LENGTH(city_code::VARCHAR)),
COUNT(DISTINCT SUBSTR(LPAD(city_code::VARCHAR, 5, '0'), 1, 2)),
SUM(pop2010), SUM(pop2020), SUM(pop2025), SUM(pop2030),
SUM(pop2035), SUM(pop2040), SUM(pop2045), SUM(pop2050)
FROM mesh4_pop_00''')
cur.fetchone()
Out[7]:
レコード数は 477,172 件、メッシュIDの上位４桁の一次メッシュは151種類、住所コードの上位２桁の都道府県コードは47種類あることを確認できます。Shapeファイルからの変換時に住所コードは数値として処理されているため、先頭のゼロが欠落して4桁になっている場合も見受けられます。また、人口は2010年以降は減少していく推計であることも確認できます。
In [8]:
columns = ['gid', 'mesh_id', 'city_code']
for p in ('', 'a', 'b', 'c', 'd'):
columns.append('pop2010{}'.format(p))
columns.extend(['pop{}{}'.format(y, p) for y in range(2020, 2051, 5)])
columns.extend(['index{}{}'.format(y, p) for y in range(2020, 2051, 5)])
columns.append('geom')
print('{} columns'.format(len(columns)))
In [9]:
cur.execute("SELECT * FROM mesh4_pop_00 LIMIT 1")
r = cur.fetchone()
print(len(r))
dict(zip(columns, r))
Out[9]:
都道府県コードごとのデータ件数と人口を確認します。 クエリ結果に合わせたカラム名と組み合わせて pandas の DataFrame にまとめます。
In [10]:
import pandas as pd
cur.execute('''
SELECT substr(lpad(city_code::varchar, 5, '0'), 1, 2) AS prefcode,
COUNT(*) AS cnt,
sum(pop2010) AS pop2010,
sum(pop2020) AS pop2020,
sum(pop2025) AS pop2025,
sum(pop2030) AS pop2030,
sum(pop2035) AS pop2035,
sum(pop2040) AS pop2040,
sum(pop2045) AS pop2045,
sum(pop2050) AS pop2050
FROM mesh4_pop_00
GROUP BY 1 ORDER BY 1
''')
cols = ['pref', 'cnt', 'pop2010']
cols.extend(['pop{}'.format(y) for y in range(2020, 2051, 5)])
df = pd.DataFrame(cur.fetchall(), columns=cols)
df
Out[10]:
グラフを描画してデータを確認します。
In [11]:
%matplotlib inline
In [12]:
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_style("whitegrid")
都道府県ごとのレコード数をグラフに描画します。
In [13]:
fig = plt.figure(figsize=(12, 3))
ax = fig.add_subplot(1, 1, 1)
ax.bar(df.pref, df.cnt)
Out[13]:
都道府県コードごとの人口をグラフに描画します。
subplots で8行分の領域を作成し、それぞれに2010年の人口と2020年以降の推計人口を表示します。
In [14]:
fig, ax = plt.subplots(8, 1, figsize=(12, 16), sharex=True, sharey=True)
ax[0].set_title('Population 2010')
ax[1].set_title('Population Estimation 2020')
ax[2].set_title('Population Estimation 2025')
ax[3].set_title('Population Estimation 2030')
ax[4].set_title('Population Estimation 2035')
ax[5].set_title('Population Estimation 2040')
ax[6].set_title('Population Estimation 2045')
ax[7].set_title('Population Estimation 2050')
ax[7].set_xlabel('Prefecture code')
ax[0].bar(df.pref, df.pop2010)
ax[1].bar(df.pref, df.pop2020)
ax[2].bar(df.pref, df.pop2025)
ax[3].bar(df.pref, df.pop2030)
ax[4].bar(df.pref, df.pop2035)
ax[5].bar(df.pref, df.pop2040)
ax[6].bar(df.pref, df.pop2045)
ax[7].bar(df.pref, df.pop2050)
Out[14]:
In [15]:
import json
sql = """
SELECT row_to_json(featurecollection)
FROM (
SELECT
'FeatureCollection' AS type,
array_to_json(array_agg(feature)) AS features
FROM (
SELECT
'Feature' AS type,
ST_AsGeoJSON(geom)::json AS geometry,
row_to_json((
SELECT p FROM (
SELECT
pop2010 AS population,
pop2020 AS pop2020,
pop2025 AS pop2025,
pop2030 AS pop2030,
pop2035 AS pop2035,
pop2040 AS pop2040,
pop2045 AS pop2045,
pop2050 AS pop2050,
mesh_id::varchar AS mesh,
lpad(city_code::varchar, 5, '0') AS city_code
) AS p
)) AS properties
FROM mesh4_pop_00
WHERE substr(lpad(city_code::varchar, 5, '0'), 1, 2) = %s
) AS feature
) AS featurecollection
"""
for p in range(47):
pref = '{:02d}'.format(p+1)
cur.execute(sql, (pref,))
geo = cur.fetchone()
with open('pref-{}-population.geojson'.format(pref), 'w') as fp:
json.dump(geo[0], fp)
GeoJSON ファイルを ZIP ファイルにまとめます。
In [16]:
import zipfile
with zipfile.ZipFile('/tmp/prefecture-population-geojson.zip', 'w') as myzip:
for p in range(47):
pref = '{:02d}'.format(p+1)
myzip.write('pref-{}-population.geojson'.format(pref))
テーブルの行ごとに GeoJSON の Feature に変換してから一次メッシュコードごとにデータを集約する SQL を定義して実行します。一次メッシュコードごとの GeoJSON ファイルをディスクに保存します。
In [17]:
import json
sql = """
SELECT mesh1,
row_to_json((SELECT c FROM (SELECT
'FeatureCollection' AS type,
features
) AS c)) AS featurecollection
FROM (
SELECT
mesh1,
array_agg(feature) AS features
FROM (
SELECT
substr(lpad(mesh_id::VARCHAR, 9, '0'), 1, 4) AS mesh1,
row_to_json((SELECT f FROM (SELECT
'Feature' AS type,
ST_AsGeoJSON(geom)::JSON AS geometry,
row_to_json((SELECT p FROM (SELECT
pop2010 AS population,
pop2020 AS pop2020,
pop2025 AS pop2025,
pop2030 AS pop2030,
pop2035 AS pop2035,
pop2040 AS pop2040,
pop2045 AS pop2045,
pop2050 AS pop2050,
mesh_id::VARCHAR AS mesh,
lpad(city_code::VARCHAR, 5, '0') AS city_code
) AS p)) AS properties
) AS f)) AS feature
FROM mesh4_pop_00
) AS t1
GROUP BY mesh1
) AS t2
ORDER BY mesh1
"""
cur.execute(sql)
for record in cur:
with open('mesh1-{}-population.geojson'.format(record[0]), 'w') as fp:
json.dump(record[1], fp)
GeoJSON ファイルを ZIP ファイルにまとめます。
In [18]:
import glob
import zipfile
with zipfile.ZipFile('/tmp/mesh1-population-geojson.zip', 'w') as myzip:
for f in glob.glob('mesh1-*-population.geojson'):
myzip.write(f)
形状データ無しのCSVファイルを出力します。
In [19]:
import csv
q = 'SELECT {} FROM mesh4_pop_00 ORDER BY 1'.format(','.join(columns[:-1]))
cur.execute(q)
with open('/tmp/mesh4-population.csv', 'w') as fp:
writer = csv.writer(fp, quoting=csv.QUOTE_ALL)
writer.writerow(columns[:-1])
for record in cur:
# Format each cell as string
r = [record[0], record[1]] # gid, mesh_id
r.append('{:05d}'.format(record[2])) # city_code
for v in record[3:]:
r.append('{:.02f}'.format(v))
writer.writerow(r)
In [ ]: