PyLearn MAP 03: コロプレス図と公共施設の利便性

ここでは、PythonのFoliumを利用して、コロプレス図を書いたり、公共施設の利便性を求める方法を説明します。

0. データ取得

国土数値情報ダウンロードサービスにいって、「2. 政策区域」にある「行政区域」にいきます。そこで、千葉県の行政区域データを入手してください。N03-19_12_190101.geojsonがあれば当たりです。

1. 千葉県内すべて描画

In [1]:
import folium
import json

# 市町村の境界データの読み込み
f = open("N03-19_12_190101.geojson", 'r', encoding='utf-8_sig')
geojson = json.load(f)
f.close()

# 中心取得
import geocoder as ge
cent = ge.arcgis("千葉県")
cent = cent.latlng

# マップの作製
m1 = folium.Map(location=cent, zoom_start=8)
folium.Choropleth(geo_data=geojson).add_to(m1)
Out[1]:
<folium.features.Choropleth at 0x11f46c9e8>
In [2]:
m1
Out[2]:

市で区切った、チーバ君を可視化することができました。このように地区ごとで区切った図をコロプレス図といいます。

2. 習志野市の描画

さて、千葉すべてではなく、特定の市だけを可視化したい時があります。今回は例として、習志野市だけの可視化を試みます。このために、jsonファイルの何番目にターゲットとしたい市があるのか探します。例えば、100番目にある市の名前は、以下のコードでアクセスできます。

In [3]:
i=100
geojson["features"][i]["properties"]["N03_004"]
Out[3]:
'銚子市'

てきとうに数字を入れて、習志野市を探す、というのはNGです。for文でターゲットとなる市の要素を取得しましょう。

In [4]:
# 習志野市が格納されている区画のIDを取得
Target = "習志野市"
for i in range(len(geojson["features"])):
    txt = geojson["features"][i]["properties"]["N03_004"]
    if txt == Target:
        GetID = i
        break

print(Target, "は", str(GetID), "番です")

習志野市 は 617 番です

geo_dataにgeojsonと指定すると、千葉県内すべてになりますが、geojson["features"][GetID]とすると、習志野市だけを描画することができます。塗りつぶしに使うChoropleth関数は以下の引数を持ちます。

  • fill_color: 塗りつぶしの色
  • fill_opacity: 塗りつぶしの濃さ（0薄い〜1濃い）
  • line_opacity: 境界の色の濃さ（0薄い〜1濃い）
In [5]:
# マップの作製
import geocoder as ge
cent = ge.arcgis("習志野市")
cent = cent.latlng

m2 = folium.Map(location=cent, zoom_start=12)
folium.Choropleth(geo_data=geojson["features"][GetID],
                              fill_color='YlGn', # 塗りつぶしの色
                              fill_opacity=0.1, # 塗りつぶしの濃さ
                              line_opacity=0.9, # 境界の色の濃さ
                            ).add_to(m2)
Out[5]:
<folium.features.Choropleth at 0x11f5299e8>
In [6]:
m2
Out[6]:

習志野市だけを区切って取り出すことができました。

3. 人口情報の添加

続いて、区切りだせた習志野市に対し、人口情報を添加してみます。PyLearnMAP02の知見を使います。

人口データの入手先： 国土数値情報　ダウンロードサービスの下の方にある「500mメッシュ別将来推計人口（H29国政局推計）（shape形式版）」を使用しています。以下は、千葉県のデータを利用しています。MAP02では1kmメッシュデータを使用しましたが、今回は500mメッシュデータを利用します。

PyLearnMAP02のコードでは、千葉県全域の人口を取り出していましたが、今回は習志野市だけを取り出す必要があります。習志野市の市町村コードは、webで調べれば簡単に出てきます。

  • 習志野市の市町村コード: 12216

PyLearnMAP02で使用したデータから、習志野市が格納されている場所を探してみます。前回の復習で、src.record(i)[1]にi番目の地域の市町村コードが格納されているので、そこを走査すればokです。

In [7]:
import shapefile
src = shapefile.Reader('Mesh4_POP_12.shx',encoding='SHIFT-JIS')

TargetID = 12216 # 習志野市の市町村コード
Narashino = []
for i in range(0, len(src)):
    if src.record(i)[1] == TargetID:
        Narashino.append(i)

print("習志野市の格納場所: ")
print(Narashino)

習志野市の格納場所: 
[8521, 8550, 8551, 8554, 9835, 9836, 9837, 9838, 9839, 9840, 9841, 9842, 9844, 9845, 9874, 9875, 9876, 9877, 9878, 9879, 9880, 9881, 9882, 9883, 9884, 9885, 9886, 9887, 9888, 9891, 9892, 9893, 9895, 9913, 9914, 9916, 9917, 9918, 9919, 9920, 9921, 9922, 9923, 9924, 9925, 9926, 9927, 9928, 9929, 9930, 9931, 9932, 9933, 9934, 9935, 9936, 9939, 9954, 9962, 9965, 9969, 9970, 9973, 9974, 9975, 9976, 9977, 10013, 10014, 10015, 10016, 10017]

これで習志野市の経度緯度、人口データが格納されている場所がわかりました。続いて、習志野市だけのデータをappendさせていきます。これもMAP02を理解していればかけるはずです。

In [8]:
Place = []
MeshCode = []
Population = []

# (なぜか0番目に埼玉が混ざっているので1から)
for i in range(1, len(src)): # len(src)は千葉県内のメッシュの数 

    # 習志野市であれば、経度緯度、人口をappend
    for j in range(0, len(Narashino)):
        if Narashino[j] == i:

            # 座標の取得
            a = src.shape(i).points[0][0]
            b = src.shape(i).points[0][1]
            Place.append([float(b), float(a)]) # 経度と緯度を逆にし、Foliumの順に

            # メッシュコードの取得
            MeshCode.append(str(src.record(i)[0]))

            # 2010年の総人口を取得
            Population.append(src.record(i)[2])
    
# 習志野市の中心座標を取得
cent = ge.arcgis("習志野市")
cent = cent.latlng

習志野市のデータが取り出せたところで、地図を描いてみます。

In [9]:
import numpy as np

# 習志野市を中心とするマップの作成
m3 = folium.Map(location=cent, zoom_start=13)

# 習志野市の塗りつぶし
folium.Choropleth(geo_data=geojson["features"][GetID],
                              fill_color='YlGn', # 塗りつぶしの色
                              fill_opacity=0.1, # 塗りつぶしの濃さ
                              line_opacity=0.9, # 境界の色の濃さ
                            ).add_to(m3)



# 1kmメッシュ地点のマーカを設定
from folium.features import DivIcon

# 赤いボックス
def MakeBox_red(in_text):
    out_text = '''
            <div style="text-align: center;
            font-size: 5pt;
            color: white;
            padding-top: 0px;
            height: 10px;
            width: 23px;
            background: red;
            border: 0.5px solid;
            border-radius: 30%;">'''
    out_text = out_text+ in_text + '</div>'
    return out_text

# 青いボックス
def MakeBox_blue(in_text):
    out_text = '''
            <div style="text-align: center;
            font-size: 5pt;
            color: white;
            padding-top: 0px;
            height: 10px;
            width: 23px;
            background: blue;
            border: 0.5px solid;
            border-radius: 30%;">'''
    out_text = out_text+ in_text + '</div>'
    return out_text

MeanPop = np.mean(Population) # 500mメッシュ平均人口
print("500mメッシュ平均人口: ", np.round(MeanPop, 2), "人")
for i in range(0, len(Place)):
    if Population[i] > MeanPop: # 1万人オーバーなら、青
        html_text = str(int(Population[i]))
        folium.Marker(location = Place[i], icon=DivIcon(html=MakeBox_blue(html_text))).add_to(m3)
    else:
        html_text = str(int(Population[i]))
        folium.Marker(location = Place[i], icon=DivIcon(html=MakeBox_red(html_text))).add_to(m3)

m3

500mメッシュ平均人口:  2285.14 人
Out[9]:

500mメッシュの平均人口を求め、それ以上なら青、未満なら赤にしてみました。数字は人口（単位は人）です。海よりは人口が少ないですね。

4. 公共施設の利便性の算出

これだけじゃ面白くないので、設置されている公共施設がどのくらい利便性があるのか、計算してみます。例として、消防署を考えます。習志野市には、以下5点の消防署があるようです。（本部と中央は同じ場所なのでひとつ扱い）

  • 習志野市消防本部（千葉県習志野市鷺沼２丁目１−４３）
  • 習志野市東消防署（千葉県習志野市東習志野２丁目２−１５）
  • 習志野市東消防署藤崎出張所（千葉県習志野市藤崎６丁目２０−１１）
  • 習志野市中央消防署谷津奏の杜出張所（千葉県習志野市奏の杜２丁目１３−１）
  • 習志野中央消防署 秋津出張所（千葉県習志野市秋津３丁目７−１）
In [10]:
FireName = [
    "千葉県習志野市鷺沼２丁目１−４３",
    "千葉県習志野市東習志野２丁目２−１５",
    "千葉県習志野市藤崎６丁目２０−１１",
    "千葉県習志野市奏の杜２丁目１３−１",
    "千葉県習志野市秋津３丁目７−１"]

FirePoint = []
for i in range(len(FireName)):
    FirePoint.append(ge.arcgis(FireName[i]).latlng)
    print(FireName[i], FirePoint[i])

千葉県習志野市鷺沼２丁目１−４３ [35.680536000000004, 140.0279445]
千葉県習志野市東習志野２丁目２−１５ [35.69166, 140.06911050000002]
千葉県習志野市藤崎６丁目２０−１１ [35.69019300000001, 140.0431545]
千葉県習志野市奏の杜２丁目１３−１ [35.68963500000001, 140.01124950000002]
千葉県習志野市秋津３丁目７−１ [35.66864700000001, 140.01572700000003]

消防署の経度・緯度が求められました。

In [11]:
# 習志野市を中心とするマップの作成
m4 = folium.Map(location=cent, zoom_start=13)

# 習志野市の塗りつぶし
folium.Choropleth(geo_data=geojson["features"][GetID],
                              fill_color='YlGn', # 塗りつぶしの色
                              fill_opacity=0.1, # 塗りつぶしの濃さ
                              line_opacity=0.9, # 境界の色の濃さ
                            ).add_to(m4)

# 500mメッシュ地点のマーカを設定
for i in range(0, len(Place)):
    if Population[i] > MeanPop: # 1万人オーバーなら、青
        html_text = str(int(Population[i]))
        folium.Marker(location = Place[i], icon=DivIcon(html=MakeBox_blue(html_text))).add_to(m4)
    else:
        html_text = str(int(Population[i]))
        folium.Marker(location = Place[i], icon=DivIcon(html=MakeBox_red(html_text))).add_to(m4)
        
        
# 消防署をプロット
for i in range(len(FireName)):
    folium.Marker(location = FirePoint[i]).add_to(m4)
    
m4
Out[11]:

消防署の位置を燻んだ青のマーカで表示しました。赤青ボックスの数字は、500mメッシュ地点の人口です。いい感じの位置にありそうなことがわかります。さて、この消防署がどのように市民をカバーしているのか、可視化してみます。このため、以下の仮定を置きます。

  • 火事が起きたとき、最も近い消防署が対応する。
  • ひとつのメッシュは、ひとつの消防署のみが対応する。

どちらも等価な条件ですが、念の為わかりやすく書いておきました。実際は厳密にはこうはならないと思いますが、現象を理解するには、単純化せねばなりません。そういうわけで、この仮定のもと、カバーする領域を可視化してみます。b

このためにまず、どの消防署が、500mメッシュ地点のどれをカバーするのか算出する関数を構築します。

In [12]:
from geopy import distance
import numpy as np

# FirePoint: リスト化された、消防署の経度・緯度
# Place: リスト化された、メッシュ地点の経度・緯度
def DecisionPoint(FirePoint, Place):

    # 以下の処理を全メッシュ地点に行う
    MinPoint = []
    for i in range(len(Place)):
        
        # メッシュ1地点と、消防署全地点の距離を求める
        temp_list = []
        for j in range(len(FirePoint)):
            temp_list.append(distance.distance(Place[i], FirePoint[j]).km)
            
        # 距離が最小の消防署のIDを取得
        MinPoint.append(np.argmin(temp_list))
        
    return MinPoint

この関数は、500mメッシュ地点に対して、すべての消防署との距離を算出し、最も距離が小さい消防署のIDを得る関数です。これにより、Populationに対する、使用される消防署がわかります。

In [13]:
MinPoint = DecisionPoint(FirePoint, Place)
for i in range(10): # たくさん出るので、10個のみ。len(MinPoint)で全部表示
    print("座標: ", Place[i], " / 消防署: ", MinPoint[i])

座標:  [35.65416666666667, 140.01875]  / 消防署:  4
座標:  [35.6625, 140.0125]  / 消防署:  4
座標:  [35.6625, 140.01875]  / 消防署:  4
座標:  [35.6625, 140.025]  / 消防署:  4
座標:  [35.66666666666667, 140.00625]  / 消防署:  4
座標:  [35.670833333333334, 140.0]  / 消防署:  4
座標:  [35.670833333333334, 140.00625]  / 消防署:  4
座標:  [35.66666666666667, 140.0125]  / 消防署:  4
座標:  [35.66666666666667, 140.01875]  / 消防署:  4
座標:  [35.670833333333334, 140.0125]  / 消防署:  4

こんな感じで、500mメッシュ地点に対する、対応する消防署のIDが出てきます。これを使って、消防署のカバー領域を可視化してみます。

In [14]:
# 習志野市を中心とするマップの作成
m5 = folium.Map(location=cent, zoom_start=13)

# 習志野市の塗りつぶし
folium.Choropleth(geo_data=geojson["features"][GetID],
                              fill_color='YlGn', # 塗りつぶしの色
                              fill_opacity=0.1, # 塗りつぶしの濃さ
                              line_opacity=0.9, # 境界の色の濃さ
                            ).add_to(m5)

# 500mメッシュ地点のマーカを設定
for i in range(0, len(Place)):
    if Population[i] > MeanPop: # 1万人オーバーなら、青
        html_text = str(int(Population[i]))
        folium.Marker(location = Place[i], icon=DivIcon(html=MakeBox_blue(html_text))).add_to(m5)
    else:
        html_text = str(int(Population[i]))
        folium.Marker(location = Place[i], icon=DivIcon(html=MakeBox_red(html_text))).add_to(m5)
        
        
# 消防署をプロット
for i in range(len(FireName)):
    folium.Marker(location = FirePoint[i]).add_to(m5)
    
# 消防署と500mメッシュのラインを表示    
for i in range(len(Place)):
    folium.PolyLine(locations=[FirePoint[MinPoint[i]], Place[i]], weight=2, color='black').add_to(m5)

m5
Out[14]:

消防署ごとのカバー領域が可視化できました。くすんだ青マーカは消防署、赤青ボックスの数字は、500mメッシュ地点の人口です。これは大変みやすく良い結果ですが、もう少し理系っぽくなるために、統計量のようなものを出して、消防署がどれだけ地域に貢献しているのか、定量化してみます。まずは、各消防署が何人の人口をカバーしているのか求めてみます。

In [15]:
CoverPop = []
for i in range(len(FirePoint)):
    s = 0
    for j in range(len(Population)):
        if MinPoint[j] == i: # そのメッシュ地点に対応している消防署のみをカウント
            s = s + Population[j]
        
    CoverPop.append(s) # 消防署iが対応する総人数をappend
    
for i in range(len(FirePoint)):
    print("消防署", i, " ( ", FireName[i], ") ", np.round(CoverPop[i], 0), "人")

消防署 0  (  千葉県習志野市鷺沼２丁目１−４３ )  34927.0 人
消防署 1  (  千葉県習志野市東習志野２丁目２−１５ )  29317.0 人
消防署 2  (  千葉県習志野市藤崎６丁目２０−１１ )  41398.0 人
消防署 3  (  千葉県習志野市奏の杜２丁目１３−１ )  28811.0 人
消防署 4  (  千葉県習志野市秋津３丁目７−１ )  30077.0 人

クリックしたらこの情報が反映されるようにします。凡例なども入れてみます。説明すると長いので興味ある人は解読をしてください。

In [47]:
# 習志野市を中心とするマップの作成
m6 = folium.Map(location=cent, zoom_start=13)

# 習志野市の塗りつぶし
folium.Choropleth(geo_data=geojson["features"][GetID],
                              fill_color='YlGn', # 塗りつぶしの色
                              fill_opacity=0.1, # 塗りつぶしの濃さ
                              line_opacity=0.9, # 境界の色の濃さ
                            ).add_to(m6)

# 500mメッシュ地点のマーカを設定
for i in range(0, len(Place)):
    if Population[i] > MeanPop: # 1万人オーバーなら、青
        html_text = str(int(Population[i]))
        folium.Marker(location = Place[i], icon=DivIcon(html=MakeBox_blue(html_text))).add_to(m6)
    else:
        html_text = str(int(Population[i]))
        folium.Marker(location = Place[i], icon=DivIcon(html=MakeBox_red(html_text))).add_to(m6)
        
        
# 消防署をプロット
for i in range(len(FireName)):
    folium.Marker(location = FirePoint[i],
                          popup = "<p>Fire station: "+ str(i)+"</p><p>Cover: " + str(int(CoverPop[i])) + "  persons</p>").add_to(m6)
    
# 消防署と500mメッシュのラインを表示    
for i in range(len(Place)):
    folium.PolyLine(locations=[FirePoint[MinPoint[i]], Place[i]], weight=2, color='black').add_to(m6)


# 凡例はhtmlで入れるしかない模様
legend_html = '''
     <div style="position: fixed; 
      top: 15px; left: 50px; width: 210px; height: 67px; Background: white;
      border:2px solid grey; z-index:9999; font-size:14px; padding-left: 5px">
      Blue box: High population area.<br>
      Red box: Low population area.<br>
      Marker: Fire station
     </div>
     '''
m6.get_root().html.add_child(folium.Element(legend_html))

m6
Out[47]:

消防署2がかなり頑張っていることがわかります。藤崎出張所ですね。消防署0が本部で一番大きい？ので、1位に来て欲しかったところですが、そうでもないようです。どこもだいたい3万人くらい〜それ以上は対応しているので、まあまあいい配置であることがわかります。

続いて、人口一人当たり、この配置がどれだけ嬉しいのか求めてみます。この指標として、人口一人当たりの、消防署までのユークリッド距離を使用します。市民の視点から見れば消防署が近くにあるほどありがたいため、一人当たりどれだけの距離に消防署があるのかを求めるという視点は、結構重要です。なお本来はユークリッド距離よりも現実的に歩く距離などを求めたほうがいいかもしれませんが、それを求めるのはかなり大変なので、ユークリッド距離で代用します。この求め方は以下の通り。

$$e = \frac{1}{P}\sum_{i \in I}\sum_{j \in J}p_i c_{ij} x_{ij}$$

$p_i$はメッシュ地点$i$の人口、$c_{ij}$は消防署$j$からメッシュ地点$i$へのユークリッド距離、$x_{ij}$はメッシュ地点$i$が消防署$j$ののカバー領域なら1、そうでなければ0を取る変数、$P$は習志野市の総人口です。$I$はすべてのメッシュ地点を包括する集合、$J$はすべての消防署を包括する集合です。また、今回はすべてのメッシュ地点は1箇所のみの消防署から供給を受けるという仮定があるので、次の制約条件が課せられています。

$$\sum_{j \in J}x_{ij} = 1, \forall i \in I$$

何故この式がすべてのメッシュ地点は1箇所のみの消防署から供給を受けるという仮定を意味しているのかについては、自分で翻訳してみてください。

In [16]:
P = np.sum(Population) # 習志野市の総人口

PCX_each_i = []
for i in range(len(Place)):
    # メッシュ地点iとそれをカバーする消防署との距離を算出
    # x_ij=1の項のみを考慮して距離を算出しているので、sum c_ij * x_ij と等価
    CX_i = distance.distance(Place[i], FirePoint[MinPoint[i]]).km

    # CXに対し、人口をかける。総人口で割る
    PCX_i = Population[i] * CX_i / P
    
    PCX_each_i.append(PCX_i)
    #print(PCX_each_i[i])

Eval = np.sum(PCX_each_i)
print("習志野市の消防署に対する1人当たりの距離は", np.round(Eval, 3), "kmです。")

習志野市の消防署に対する1人当たりの距離は 0.828 kmです。

$e$を求めることができました。一人当たりだいたい消防署まで 830m くらいの距離があるようです。厳密な値ではありませんが、時速60kmで走るとすると、だいたい1分くらいで到着できる位置に消防署が建てられている、と考えることもできます。

さて、これは最適でしょうか。ここでいう最適とは「$e$が最も低い」と定義します。当然のことながら、最適であるはずがありません。もっと良い消防署5箇所の配置があるはずです。

このような問題を、施設配置問題 (Facility Location Problem) と呼び、数理最適化、オペレーションズリサーチの分野で幅広く研究されています。

地図上のどこにどのような施設を何箇所配置するのが最適なのか、やってみたい人は、ここらのプログラミング技術を取得する必要があります。

消防署に限らず、中学校だとか、図書館だとか、救急車の配置だとか、病院だとか、とにかくいろいろな施設の最適配置が研究されています。興味がある人はgoogle scholarなどで「オペレーションズリサーチ、最適配置」などのキーワードで調べてみると良いでしょう。

最後に保存しておしまい。

In [48]:
m6.save("FireStation.html")