Pandas+Matplotlibで市町村の人口を可視化する

%matplotlib inline
from itertools import chain
import pandas as pd
import numpy as np
import shapefile
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon
from matplotlib.collections import PatchCollection
import seaborn as sns

def list_to_dict(lst):
    return dict(zip(range(0, len(lst)), lst))

def list_flatten(lst):
    # 2重リストを1重にする
    return list(chain.from_iterable(lst))

def list_combine(lst, *args):
    for l in args:
        lst.extend(l)
    return tuple(lst)

# 点列と分割点のリストからPolygonのリストを作る
def make_polygons(row):
    row = row.tolist()
    axis = row[0]
    parts = list(row[1])
    parts.append(len(axis))
    axis_list = [axis[parts[i]:parts[i + 1]] for i in range(len(parts) - 1)]
    Polygons_list = [Polygon(x) for x in axis_list]
    return Polygons_list

# データの読み込み
sf = shapefile.Reader("japan_rvis_h25/japan_rvis_h25")

sf.shapes()[:10] # これら１つ１つがポリゴンデータ

[<shapefile._Shape at 0x108418cf8>,
 <shapefile._Shape at 0x108418f28>,
 <shapefile._Shape at 0x108418f60>,
 <shapefile._Shape at 0x108418b00>,
 <shapefile._Shape at 0x108418048>,
 <shapefile._Shape at 0x1084180f0>,
 <shapefile._Shape at 0x108418208>,
 <shapefile._Shape at 0x108418d68>,
 <shapefile._Shape at 0x1084182e8>,
 <shapefile._Shape at 0x108418390>]

field_name = [x[0] for x in sf.fields[1:]] # 属性データの項目
records = [[x.decode("shift-jis") if isinstance(x, bytes) else x for x in y]
                   for y in sf.iterRecords()]  #属性データのSHIFT-JISをデコード

sr_points = pd.Series(
    [i.shape.points for i in sf.shapeRecords()]).to_frame().rename(
        columns={0: "AXIS"})  # ポリゴンデータ
sr_parts = pd.Series([s.parts for s in sf.shapes()]).to_frame().rename(
    columns={0: "PARTS"})  # ポリゴンの分割点

field_name # 属性データの項目

['KEN', 'SICHO', 'GUN', 'SEIREI', 'SIKUCHOSON', 'JCODE', 'P_NUM', 'H_NUM']

#　ポリゴンデータ（1行が1つの市区町村に対応）
sr_points.head(5)

# 分割点のリスト
sr_parts.head(15)

# Dataframeでまとめる
df = pd.DataFrame(records).rename(
            columns = list_to_dict(field_name)
        )[["KEN", "SEIREI", "SIKUCHOSON","JCODE","P_NUM"]].assign(
            # 空白の数値データを0で埋める
            P_NUM = lambda df: df["P_NUM"].apply(lambda x: 0 if isinstance(x, str) else x),
        ).applymap(
            # 文字列のスペースを消去
            lambda x: x.strip() if isinstance(x, str) else x
        ).applymap(
            # 文字数ゼロをNaNにする
            lambda x: np.nan if x == "" else x
        ).join(
            # ポリゴンと分割点を統合
            [sr_points, sr_parts]
        ).assign(
            # MatplotlibのPolygonオブジェクトのインスタンスの列をつくる
            POLYGON = lambda df: df[["AXIS", "PARTS"]].apply(lambda x: make_polygons(x), axis=1)
        ).drop(
            # いらなくなったポリゴンと分割点を除去
            ["AXIS", "PARTS"], axis=1
        )

# データの確認。札幌市が区ごとに別のデータになっているのが分かる。
df.head(10)

# 政令指定都市一覧
df["SEIREI"].dropna().unique()

array(['札幌市', '仙台市', 'さいたま市', '千葉市', '横浜市', '川崎市', '相模原市', '新潟市', '静岡市',
       '浜松市', '名古屋市', '京都市', '大阪市', '堺市', '神戸市', '岡山市', '広島市', '北九州市',
       '福岡市', '熊本市'], dtype=object)

df_agg = df.drop(
                df.ix[df["SEIREI"].notnull()].index  # 政令指定都市である地域を削除
            ).append(
                # 政令指定都市を統合したデータを追加
                df.groupby(["KEN","SEIREI"]).agg(
                    {"P_NUM": np.sum ,
                     "POLYGON": list_combine,
                     "JCODE": lambda sr: sr.iloc[0]}
            ).reset_index().assign(
                # ポリゴンデータも統合
                POLYGON = lambda df: df["POLYGON"].apply(list).apply(list_flatten)
            )
        ).drop("SEIREI", axis=1).sort_values("JCODE").reset_index(drop=True)

# 札幌市をまとめた。（SIKUCHOSONがNaNになってしまってはいるが）
df_agg.head(10)

def make_majority_index(df):
    half_total_P_NUM = df.query("KEN == @prefecture")["P_NUM"].sum() / 2  # 県の人口の半分
    max_index = df.query("KEN == @prefecture")["P_NUM"].idxmax()  # 最大の人口を抱える市区町村のインデックス
    # 人口の多い地区を上から順番に、人口の半分を超えるまで選択する
    majority_index = df.query("KEN == @prefecture").sort_values(
        "P_NUM", ascending=False).assign(
            cum_P_NUM=lambda df: df["P_NUM"].cumsum()).query(
                "cum_P_NUM < @half_total_P_NUM").index
    if majority_index.size == 0:
        # 1つの地域だけで半分を超えた場合はその最大地域を採用
        majority_index = [max_index]
    return majority_index

prefecture = "北海道"
majority_index = make_majority_index(df_agg)  # 人口が多い方(分割で赤く塗る方)の市区町村のインデックスの列

# 県人口の半数を占める市町村
df_agg.query("index in @majority_index")

# 上で選択したmajority_indexを用いてPolygonの塗り分けを選ぶ
Polygon_red = df_agg.query("index in @majority_index")["POLYGON"]
Polygon_blue = df_agg.query("KEN == @prefecture and not index in @majority_index")["POLYGON"]

# プロット
with sns.axes_style('white'):
    fig = plt.figure(figsize=(15,10))
    ax = fig.add_subplot(111)

    ax.add_collection(PatchCollection(list_flatten(Polygon_red.tolist()), facecolor="red", edgecolor='k', linewidths=.2))
    ax.add_collection(PatchCollection(list_flatten(Polygon_blue.tolist()), facecolor="blue", edgecolor='k', linewidths=.2))
    ax.tick_params(axis='x', which='both', top='off', bottom='off', labelbottom='off')
    ax.tick_params(axis='y', which='both', left='off', right='off', labelleft='off')
    ax.set_aspect('auto')
    ax.autoscale_view()
    plt.show()