In [1]:
import osmnx as ox
import geopandas as gpd
import pandas as pd
from shapely.geometry import Point
from geopandas import GeoSeries
from shapely.ops import cascaded_union
import networkx as nx
import numpy as np
import shapely.geometry
from tqdm import tqdm
import matplotlib.pyplot as plt
In [2]:
import pysal as ps
from scipy.sparse import csr_matrix
In [3]:
# 台北的 shapefile
tw = gpd.read_file('./../../../../宸訊/lisa/data/taiwan_county/taiwan_county.shp')
tpe = tw.query("COUNTY == '台北市' or COUNTY == '新北市'")
tpe['geometry'].plot()
Out[3]:
In [4]:
tpe = tpe.to_crs(epsg = 4326)
tpe = tpe.reset_index(drop = True)
In [5]:
# 合併台北、新北市
tpe = gpd.GeoSeries(cascaded_union(tpe['geometry']))
tpe.plot()
Out[5]:
In [6]:
# 用 open street map 取得路網
G = ox.graph_from_polygon(tpe[0], network_type='drive')
ox.plot_graph(G)
Out[6]:
In [7]:
proj_G = ox.project_graph(G, to_crs = {'init':'epsg:3826'}) # 投影
nodes_proj, edges_proj = ox.graph_to_gdfs(proj_G, nodes=True, edges=True) # 節點+邊
edges = edges_proj[['geometry', 'length', 'name', 'osmid', 'u', 'v']]
edges['geometry'].plot()
Out[7]:
In [8]:
edges.head()
Out[8]:
bus stop¶
In [9]:
stop = pd.read_csv('./busstion.csv', encoding = "big5")
In [10]:
stop['Freq'] = stop['Freq'].apply(lambda t: t if t != np.inf else 1/18)
In [11]:
stop = stop[['BusInfo.Id', 'BusInfo.longitude', 'BusInfo.latitude', 'Freq']]
In [12]:
stop = stop.groupby('BusInfo.Id').mean().reset_index()
In [13]:
bus_geom = stop[['BusInfo.longitude', 'BusInfo.latitude']].apply(Point, 1)
In [14]:
stop['geometry'] = bus_geom
In [15]:
stop = gpd.GeoDataFrame(stop, crs = {'init':'epsg:4326'}, geometry = bus_geom)
In [16]:
stop = stop.to_crs(epsg = 3826)
In [65]:
stop.plot(markersize = 0.5, figsize = (8, 8))
Out[65]:
In [18]:
# 找每個公車站在 Graph 上的最近 node
neighbor = ox.utils.get_nearest_nodes(proj_G, stop['geometry'].x, stop['geometry'].y)
neighbor = list(neighbor)
len(neighbor)
Out[18]:
In [19]:
stop['osmid'] = neighbor
nodes = nodes_proj[['lat', 'lon', 'geometry', 'osmid']]
In [20]:
len(nodes)
Out[20]:
In [21]:
len(stop)
Out[21]:
In [22]:
stops = pd.merge(stop, nodes, on = 'osmid', how = 'left')
In [23]:
len(stops)
Out[23]:
In [24]:
stops = stops.groupby('osmid').mean()
In [25]:
stops = stops.reset_index()
In [26]:
stops.head()
Out[26]:
In [27]:
# 再轉一次
stops['lat'] = stops['lat'].astype(float)
stops['lon'] = stops['lon'].astype(float)
geometry = stops[['lon', 'lat']].apply(Point, 1)
stops = gpd.GeoDataFrame(stops, crs = {'init':'epsg:4326'}, geometry = geometry)
stops = stops.to_crs(epsg = 3826)
In [28]:
# 只是確定一下位置有沒有問題
ax = nodes.plot()
stops.plot(ax = ax)
Out[28]:
In [29]:
# 取 geometry 站點位置
stops_loc = stops[['osmid', 'geometry']]
stops_loc = pd.DataFrame(stops_loc)
stops_loc = gpd.GeoDataFrame(stops_loc, geometry = stops_loc['geometry'], crs = {'init':'epsg=3826'})
stops_loc['x'] = stops_loc['geometry'].x
stops_loc['y'] = stops_loc['geometry'].y
stops_loc = stops_loc[['osmid', 'x', 'y']]
len(stops_loc)
Out[29]:
In [30]:
# 計算兩點間最短路徑
def Distance(inp):
start = inp[0]
end = inp[1]
try:
route = nx.shortest_path(proj_G, start, end, weight='length')
# route 的距離
nodes_ = nodes.loc[route, 'geometry'].tolist()
length = shapely.geometry.LineString(nodes_).length
return length
except:
print(start, end)
In [31]:
bandwidth = 1000
In [32]:
# buffer
print('buffer')
stop_copy = stops.copy()
stop_copy['geometry'] = stops.geometry.buffer(bandwidth)
stop_copy = gpd.sjoin(stop_copy, stops, op='contains')
stop_copy = stop_copy[['osmid_left', 'osmid_right']]
# 最短路徑，建立sparce matrix
stop_copy = stop_copy.query('osmid_left != osmid_right') # 扣掉自己
# 加回位置
stop_copy = pd.merge(stop_copy, stops_loc, left_on = 'osmid_left', right_on = 'osmid')
stop_copy = stop_copy.rename(columns = {'x':'center_x', 'y':'center_y'})
stop_copy = pd.merge(stop_copy, stops_loc, left_on = 'osmid_right', right_on = 'osmid')
stop_copy = stop_copy.rename(columns = {'x':'neighbor_x', 'y':'neighbor_y'})
del stop_copy['osmid_x']
del stop_copy['osmid_y']
In [42]:
len(stop_copy)
Out[42]:
In [43]:
# 算出兩兩節點之間最短路徑
stop_copy['distance'] = stop_copy[['osmid_left', 'osmid_right']].apply(Distance, 1)
In [35]:
# 建立 kernel weight
stop_neighbor = stop_copy.query('distance < {}'.format(bandwidth))
stop_neighbor['weight'] = stop_neighbor['distance'].apply(lambda t: 1 - t/bandwidth) # kernel function = trianguler
# 製作 sparse matrix id 與 站牌 id 對應
stops_num = len(stop_neighbor['osmid_left'].drop_duplicates())
ids = list(range(stops_num))
osmid = stop_neighbor['osmid_left'].unique()
osmid_id = pd.DataFrame({
'index': ids,
'osmid': osmid
})
stop_sparse = pd.merge(stop_neighbor, osmid_id, left_on = 'osmid_left', right_on = 'osmid')
stop_sparse = stop_sparse.rename(columns = {'index':'index_left'})
stop_sparse = pd.merge(stop_sparse, osmid_id, left_on = 'osmid_right', right_on = 'osmid')
stop_sparse = stop_sparse.rename(columns = {'index':'index_right'})
del stop_sparse['osmid_x']
del stop_sparse['osmid_y']
# sparse matrix
row = np.array(stop_sparse['index_left'].astype(int))
col = np.array(stop_sparse['index_right'].astype(int))
weight = np.array(stop_sparse['weight'].astype(float))
# 建立權重
sparse = csr_matrix((weight, (row, col)), shape=(max(row)+1, max(row)+1))
wsp = ps.lib.weights.WSP(sparse)
w = ps.lib.weights.WSP2W(wsp)
# 對應的 frequency
stop_freq = pd.merge(osmid_id, stops[['osmid', 'Freq']].drop_duplicates(), on = 'osmid')
value = np.array(stop_freq['Freq'])
# global moran's I
gmI = ps.explore.esda.Moran(value, w, two_tailed=False).I
print(gmI)
In [36]:
lm = ps.explore.esda.Moran_Local(value, w)
In [37]:
lm.p_sim
Out[37]:
In [38]:
lm.q
Out[38]:
In [39]:
stop_freq['p'] = lm.p_sim
In [40]:
stop_freq['q'] = lm.q
In [41]:
stop_plot = pd.merge(stop_freq, stops_loc, on = 'osmid')
In [44]:
HH = stop_plot.query('p < 0.05 and q == 1').copy()
LL = stop_plot.query('p < 0.05 and q == 3').copy()
In [49]:
HH['geometry'] = HH[['x', 'y']].apply(shapely.geometry.Point, axis=1)
LL['geometry'] = LL[['x', 'y']].apply(shapely.geometry.Point, axis=1)
In [75]:
HH = gpd.GeoDataFrame(HH, crs={'init': 'epsg:3826'})
LL = gpd.GeoDataFrame(LL, crs={'init': 'epsg:3826'})
# HH = HH.to_crs(epsg=4326)
In [274]:
HH.to_file('HH.geojson', driver='GeoJSON')
In [275]:
# hh = stop_plot.query('q == 1')
In [47]:
%matplotlib ipympl
In [76]:
ax = edges.plot(alpha = 0.1, figsize = (8, 8), color = 'gray')
HH.plot(ax = ax, c = 'red', markersize = 0.7)
LL.plot(ax = ax, c = "blue", markersize = 0.7)
Out[76]:
In [77]:
ax = edges.plot(alpha = 0.05, figsize = (10, 10), color = 'gray')
plt.scatter(HH['x'], HH['y'], c = 'red', s = 0.8)
plt.scatter(LL['x'], LL['y'], s = 0.8)
Out[77]: