まちのクラスタリングをやりたい(2)
全国都市クラスタリング可視化マップ作ってみました。
「全国の都市が、都市相互でどのような近隣関係を持つのか」を可視化することがテーマ。図には「都市名：人口」がマークされています。
札幌の都市圏は、道東・道北の都市とつながりつつ、青森方向にも繋がりがあることが分かります。
青森・函館・札幌という拠点都市以外にはそれほど都市は存在しませんね。
続いて、首都圏。
この東京圏の図を解像度高く見やすくしたものは、また別に作ったので、機会があれば。
まあこんな感じで、甲信越～北陸～名古屋～近畿と作っていきました。
以上のような成果物ができましたとさ。
…作り方は以下のように進めていきました。
①人口情報と、経度緯度情報の結合 ward4,ward5 → ward.6 とりあえず人口5万人以上の都市について、市区町村の名前をキーにして、経度緯度データを結合する
df = open('ward3.txt', 'r', encoding='shift_jis')
string = df.readlines()
a=[]
for x in string:
a.append(x.split('\t'))
print(a)
df2 = open('ward4.txt', 'r', encoding='shift_jis')
string2 = df2.readlines()
a2=[]
for x in string2:
a2.append(x.split('\t'))
print(a2)
for i in a2:
k = 0
for j in a:
if(i[0] == j[0]):
if(i[1] == j[1] or i[1] == j[2]):
i.append([j[3],j[4],j[5]])
k = 1
if(k == 0):
print("NO NAME")
print(i[2])
pprint.pprint(a2)
②複数の区をもつ市は、複数の経度緯度のデータを持つようになる。その場合、すべての区の重心を求めて代表点とする ward6 → List Object
df = open('ward6.txt', 'r', encoding='shift_jis')
string = df.readlines()
a=[]
for x in string:
a.append(x.split('\t'))
b = []
m = ''
n = ''
for i in a:
if(i[0]) == str('['):
c = 1
lng = float(i[5])
lat = float(i[6])
tdfk = i[1]
city = i[2]
popl = i[3]
if(i[0] == str('')):
c += 1
lng += float(i[5])
lat += float(i[6])
if(i[8] == str(']]\n')):
m = str(round(lng / c, 6))
n = str(round(lat / c, 6))
b.append([tdfk,city,popl,m,n])
pprint.pprint(b)
とりあえず出力されたものがこちら
[['[神奈川県]', '[横浜市]', 3724844, '35.453767', '139.5838'],
['[大阪府]', '[大阪市]', 2691185, '34.671937', '135.50801'],
['[愛知県]', '[名古屋市]', 2295638, '35.14962', '136.926371'],
['[北海道]', '[札幌市]', 1952356, '43.052385', '141.364933'],
['[福岡県]', '[福岡市]', 1538681, '33.585779', '130.38472'],
['[兵庫県]', '[神戸市]', 1537272, '34.685926', '135.148421'],
['[神奈川県]', '[川崎市]', 1475213, '35.580315', '139.614373'],
['[京都府]', '[京都市]', 1475183, '34.999513', '135.754851'],
['[埼玉県]', '[さいたま市]', 1263979, '35.896661', '139.63885'],
['[広島県]', '[広島市]', 1194034, '34.407758', '132.463318'],
['[宮城県]', '[仙台市]', 1082159, '38.266032', '140.887951'],
…
人口100万人以上までコピペしたけど、こういった調子で人口5万人前後の都市（最後尾は滋賀県高島市）まで続く。
（ついでに言うと、自分が扱っているデータは、1次データでないし計測年月も確認してないから、数値として当てにしない方がいいです）
③データの階層クラスリングを実行 （重みづけ有のWard法）
INTA = []
for i in b:
INTA.append([[i[0]+i[1],i[2],float(i[3]),float(i[4])]])
def DSfor(C):
x = 0.0
y = 0.0
c = 0.0
for i in C:
x = x + i[2]
y = y + i[3]
c = c + 1.0
xmean = x/c
ymean = y/c
Ax = 0.0
Ay = 0.0
for i in C:
Ax = Ax + (i[2] - xmean)**2.0
Ay = Ay + (i[3] - ymean)**2.0
return (56.25) *Ax + (126.5625) *Ay
def J(g,h,U):
C = U[g]+U[h]
RevS = []
DSRev = 0.0
d = 0
for i in U:
if(not (d == g or d == h)):
RevS.append(d)
DSRev = DSRev + DSfor(i)
d = d + 1
DSsum = DSRev + DSfor(C)
return DSsum
def updateU(U):
minDS = J(0,1,U)
ug1 = 0
ug2 = 1
for i in range(len(U)):
for j in range(len(U)):
if(not i == j):
if(J(i,j,U) < minDS):
minDS = J(i,j,U)
ug1 = i
ug2 = j
print(ug1,ug2,minDS)
V = []
for i2 in range(len(U)):
if(not (i2 == ug1 or i2 == ug2)):
V.append(U[i2])
V.append(U[ug1]+U[ug2])
return V
for i in range(519):
c = updateU(c)
pprint.pprint(c[-1])
すると、こういうデータになる
…
50 80 407.2156215596323
[['[福岡県][宗像市]', 96516, 33.805431, 130.540718],
['[福岡県][福津市]', 58781, 33.786866, 130.469886],
['[福岡県][古賀市]', 57959, 33.728787, 130.469987],
['[福岡県][福岡市]', 1538681, 33.585779, 130.38472],
['[福岡県][春日市]', 110743, 33.532571, 130.470281],
['[福岡県][大野城市]', 99525, 33.536307, 130.478641],
['[福岡県][筑紫野市]', 101081, 33.496342, 130.515657],
['[福岡県][太宰府市]', 72168, 33.512799, 130.523906]]
40 141 413.8777469799072
[['[香川県][丸亀市]', 110010, 34.289506, 133.797706],
['[香川県][坂出市]', 53164, 34.316319, 133.86051],
['[岡山県][倉敷市]', 477118, 34.585013, 133.772084],
['[岡山県][総社市]', 66855, 34.672807, 133.746528],
['[岡山県][岡山市]', 719474, 34.644351, 133.953009],
['[岡山県][玉野市]', 60736, 34.491943, 133.945987]]
37 48 420.5537239345434
[['[富山県][南砺市]', 51327, 36.557513, 136.875336],
['[石川県][金沢市]', 465699, 36.561325, 136.656205],
['[石川県][白山市]', 109287, 36.514427, 136.565892],
['[石川県][野々市市]', 55099, 36.519704, 136.609971]]
…
重み、クラスターのまとまりからなるセット。
このデータを編集して作成した次のtxtデータをもとに、クラスターの座標点（N>=3）に対する凸包の多角形を求め、ひもづける。
import copy
df = open('ward7.txt', 'r', encoding='shift_jis')
string = df.readlines()
a=[]
for x in string:
a.append(x.split('\t'))
def det(p,q):
return (p[0]*q[1]-p[1]*q[0])
def sub(p,q):
return (p[0]-q[0],p[1]-q[1])
def get_convex_hull(set_points):
points = copy.deepcopy(set_points)
n = len(points)
points.sort() #素直にsortすると各配列の先頭要素をkeyとしてsortします
size_convex_hull = 0
ch = []
if n == 2:
ch.append((points[0][0],points[0][1]))
ch.append((points[1][0],points[1][1]))
ch.append((0,0))
else:
for L in range(n):
while size_convex_hull > 1:
v_cur = sub(ch[-1], ch[-2])
v_new = sub(points[L], ch[-2])
if det(v_cur, v_new) > 0:
break
size_convex_hull -= 1
ch.pop()
ch.append((points[L][0],points[L][1]))
size_convex_hull += 1
t = size_convex_hull
for L in range(n-2,-1,-1):
while size_convex_hull > t:
v_cur = sub(ch[-1], ch[-2])
v_new = sub(points[L], ch[-2])
if det(v_cur, v_new) > 0:
break
size_convex_hull -= 1
ch.pop()
ch.append((points[L][0],points[L][1]))
size_convex_hull += 1
return ch[:-1]
k=[]
j=[]
for i in a:
if(not i[0] == ''):
if(not j==[]):
j.append(get_convex_hull(j[1]))
k.append(j)
j=[]
j.append(i[0])
j.append([])
else:
i[4] = float(i[4])
i[5] = float(i[5])
i[6] = float(i[6].split(']')[0])
j[1].append((i[5],i[6],i[4],i[2]))
pprint.pprint(k)
④QGISにもちこむために、Well known text(WKT)形式にジオメトリ定義を変更したｃｓｖファイルを作成しようとするが、挫折
⑤ここで、Matplotlibで可視化する方針に変えて、テキストデータをきちんと構造的に読み取るプログラム作る
#!/usr/bin/python
# coding: UTF-8
import pprint
import os
os.chdir('C:\\Users\\leafh\\Desktop')
df = open('ward8.txt', 'r', encoding='shift_jis')
string = ""
for i in df.readlines():
string += i
df.close()
a = "".join(string)
b = []
for i in a:
b.append(i)
c=[]
iL=[]
def read_list(stl,start_sym = '[', end_sym =']', param = "JOINING", excp_syms = ('(',")"), remv_syms = ('\n'," ")):
i = 0
while True:
if not type(stl) == list:
print("ERROR? : This type is not the list-type")
break #case1
if i > len(stl)-1:
break #case2
elif stl[i] == str(start_sym):
c.append(i)
i += 1 #case3-A
elif stl[i] == str(end_sym):
if c == []:
print("ERROR? : start symbol-end symbol is not matched.")
print(i,iL)
i += 1
else:
nl,_ = (c[-1],0)
nh,_ = (i,0)
STL = list(stl[nl+1:nh])
stl = stl[:nl] + [STL] + stl[nh+1:]
c.pop()
iL.append(nl)
stl[nl] = read_list(stl[nl],start_sym,end_sym,param,excp_syms,remv_syms)
i,_ = iL[-1],0
iL.pop()
i += 1 #case3-B
elif type(stl[i]) == str:
if (param == False):
i += 1
else:
if stl[i] in excp_syms:
i += 1 #case4-A
elif stl[i] in remv_syms:
stl.pop(i) #case4-B: count doesnot proceed by removing
elif (param == "NOT JOINING"):
i += 1
elif i > len(stl)-2:
i += 1
elif not(type(stl[i+1]) == str):
i += 1
else:
if ((stl[i+1] == str(start_sym)) or (stl[i+1] == str(end_sym))):
i += 1 #case4-A
else:
if stl[i+1] in excp_syms:
i += 1 #case4-A
elif stl[i+1] in remv_syms:
stl.pop(i+1) #case4-B: count doesnot proceed by removing
else:
stl[i] = stl[i]+stl[i+1]
stl.pop(i+1) #case4-C: count doesnot proceed by join&pop
elif type(stl[i]) == list:
iL.append(i)
stl[i] = read_list(stl[i],start_sym,end_sym,param,excp_syms,remv_syms)
i,_ = iL[-1],0
iL.pop()
i += 1
else:
i += 1
return stl
b2 = read_list(b,start_sym ='(', end_sym =')', param = "JOINING", excp_syms = ['[',"]",',',"'"],remv_syms = ['\n'," "])
b3 = read_list(b2,start_sym ='[', end_sym =']', param = "JOINING", excp_syms = [',',"'"],remv_syms = ['\n'," "])
b4 = read_list(b3,start_sym ='*', end_sym ='*', param = "NOT JOINING", excp_syms = [],remv_syms = [',',"'"])
pprint.pprint(b4)
これで、たとえば第406項を取ると、
In : b4[0][406] ... Out : ['602.5981908', [['36.042126', '139.399959', '91437.0', ['埼玉県'], ['東松山市']], ['35.95717', '139.402905', '101679.0', ['埼玉県'], ['坂戸市']], ['35.934515', '139.393098', '70255.0', ['埼玉県'], ['鶴ヶ島市']], ['35.852942', '139.412213', '152405.0', ['埼玉県'], ['狭山市']], ['35.835766', '139.391058', '148390.0', ['埼玉県'], ['入間市']], ['35.855731', '139.327734', '80715.0', ['埼玉県'], ['飯能市']], ['35.907796', '139.339026', '56520.0', ['埼玉県'], ['日高市']], ['35.666339', '139.315806', '577513.0', ['東京都'], ['八王子市']], ['35.705755', '139.353548', '111539.0', ['東京都'], ['昭島市']], ['35.787996', '139.27583', '139232.0', ['東京都'], ['青梅市']], ['35.72892', '139.294119', '80954.0', ['東京都'], ['あきる野市']], ['35.738451', '139.326932', '58395.0', ['東京都'], ['福生市']], ['35.767222', '139.310944', '55833.0', ['東京都'], ['羽村市']], ['35.799672', '139.46861', '340386.0', ['埼玉県'], ['所沢市']], ['35.754597', '139.468489', '149956.0', ['東京都'], ['東村山市']], ['35.728577', '139.477456', '190005.0', ['東京都'], ['小平市']], ['35.710335', '139.463191', '122742.0', ['東京都'], ['国分寺市']], ['35.637004', '139.446307', '146631.0', ['東京都'], ['多摩市']], ['35.683885', '139.44138', '73655.0', ['東京都'], ['国立市']], ['35.671337', '139.394926', '186283.0', ['東京都'], ['日野市']], ['35.754861', '139.387399', '71229.0', ['東京都'], ['武蔵村山市']], ['35.714014', '139.407843', '176295.0', ['東京都'], ['立川市']], ['35.745369', '139.426593', '85157.0', ['東京都'], ['東大和市']]], [['35.637004', '139.446307'], ['35.666339', '139.315806'], ['35.72892', '139.294119'], ['35.787996', '139.27583'], ['36.042126', '139.399959'], ['35.799672', '139.46861'], ['35.728577', '139.477456']]]
となる。この各項は各クラスターグループの情報を表し、
b4[0][406][0]がクラスターの偏差二乗和値、
b4[0][406][1]がクラスターに含まれるグループ、
b4[0][406][2]がクラスターの凸包多角形の頂点座標の一順を記述するリスト。
⑥Matplotlibで可視化
plt.figure(figsize=(35,35), dpi = 450)
ax = plt.gca()
ax.set_xlim(137,141)
ax.set_ylim(34,38)
plys=[]
for p in b5[0]:
if len(p[2])>2:
poly = Polygon(xy = p[2], fc='orange', ec='k', alpha = 0.15, linewidth =0.05)
plys.append(poly)
p[1]
for p in plys:
ax.add_patch(p)
plt.text(135,40,"North Center(135,40)",fontsize=1)
for spot_id in range(len(Xbox)):
if Pbox[spot_id] < 200000:
plt.text(Xbox[spot_id],Ybox[spot_id], Nbox[spot_id]+":"+str(int(Pbox[spot_id])), fontsize=2, alpha =0.3)
elif Pbox[spot_id] < 500000:
plt.text(Xbox[spot_id],Ybox[spot_id], Nbox[spot_id]+":"+str(int(Pbox[spot_id])), fontsize=5, alpha =0.8)
else:
plt.text(Xbox[spot_id],Ybox[spot_id], Nbox[spot_id]+":"+str(int(Pbox[spot_id])), fontsize=7)
plt.text(135,35,"The Center(135,35)",fontsize=1)
plt.axis('scaled')
ax.set_aspect('equal')
plt.savefig('JapanTokyo.png')
部分拡大して取り込みたければ、plt.axis('scaled') の部分をなくしてあげると良い。
⑦改善点
横浜市、大阪市などは、「区」が重要になってくるのだけれどもそれを省いて今回はクラスタリングしています。
・本来だったら、クラスターに制限をかけて、ある一定範囲の面積内にするとか、人口をカウントしたうえで、大きい人口の場所が小さい人口の場所を取り込むようにした方がいい気もする
・愛媛の宇和島市、岐阜の高山市のような、5~10万人規模の都市で「近隣都市」と呼べるような周囲の都市が存在しない都市があり、あるていどクラスタリングが進んだ後の地域間のクラスタリングに影響を与えている
・静岡浜松は愛知との結びつきの方が強いイメージだが、関東圏としてクラスタリングされている
・山陰地方の扱いが雑い。山陰の中での結びつきが軽視されている。やはり2次元座標だと、山間部かどうかを考慮していないから？
今回初めて、何かしらのデータを可視化することに挑戦してみました。データ自体の信頼性とかは二の次で、一通りのフローを達成できたので満足です。
それにしても難しいね。
