0 導入庫
import osimport pandas as pd
pd.set_option('display.max_rows',5)import osmnx as oximport geopandas as gpd
from shapely.geometry import Point1 讀取數據
假設我們有 如下的數據：
1.1 新加坡室外基站位置數據
cell_station=pd.read_csv('outdoor_LTE.csv')
cell_station
1.2 新加坡路網openstreetmap數據
G=ox.graph_from_place('Singapore,Singapore',network_type='drive')
ox.plot_graph(G)
1.2.1 從openstreetmap數據中提取路網數據
road_network=ox.utils_graph.graph_to_gdfs(G,nodes=False)
road_network
1.3 出行軌跡數據
traj=pd.read_csv('processed_dart_outdoor_3d.csv')
traj
其中latitude和longitude 是用戶位置，帶cell的是對應的基站位置，new_installation_id是用戶id，timestamp_5s是時刻
1.3.1 出行軌跡轉GeoDataFrame
points = [Point(xy) for xy in zip(traj.longitude, traj.latitude)]
points_gdf = gpd.GeoDataFrame(traj, geometry=points)
points_gdf
2 判斷用戶點在不在路網上
2.1 為每條道路創建非常小的緩沖區
在赤道附近，經緯度坐標系統中的一個度大約等于地球表面上的111公里，所以這里的buffer相當于1m
road_network_buffered = road_network.geometry.buffer(0.00001)
#將路網線幾何對象緩沖一定距離（例如，1米），創建一個新的GeoDataFrame
road_network_buffered
'''
u v key
25451929 6749812859 0 POLYGON ((103.87103 1.29515, 103.87066 1.29508...
25455287 1637003462 0 POLYGON ((103.87412 1.29550, 103.87413 1.29550......
10732302222 259401350 0 POLYGON ((103.90657 1.30628, 103.90657 1.30628...
10806629050 2325064861 0 POLYGON ((103.90709 1.30698, 103.90709 1.30698...
Length: 45583, dtype: geometry
'''緩沖區轉化為geoDataFrame
road_network_buffered_gdf = gpd.GeoDataFrame(geometry=road_network_buffered)
road_network_buffered_gdf
2.2 判斷每個點在不在路網的buffer上
points_in_road_network = gpd.sjoin(points_gdf,road_network_buffered_gdf, how="inner", op='within')
points_in_road_network-
gpd.sjoin()函數：執行空間連接操作。它將兩個GeoDataFrame基于空間關系合并。【基于點（points_gdf）是否在多邊形（road_network_buffered_gdf）內部】 -
how="inner"：指定連接類型為內連接。這意味著結果中只會包含在points_gdf中的點，并且這些點必須位于road_network_buffered_gdf內部。不在緩沖區內的點將被排除在外。 -
op='within'：指定空間操作類型為“within”，即查找outdoor_traj中哪些點位于road_network_buffered_gdf的緩沖區多邊形內部。
- 但是sjoin會存在一個問題：如果一個points_gdf中的點同時在兩條路段的buffer中，結果中會分別出現這個點+一條路段buffer 的兩個結果
- ——>一個時刻一個用戶id，只保留一條即可?
points_in_road_network_in_road=points_in_road_network.drop_duplicates(subset=['new_installation_id','timestamp_5s'])
points_in_road_network_in_road
3 不在路網的點和最近路段的距離
3.1 找到不在路網的用戶點
traj_remain=traj.iloc[traj.index.difference(points_in_road_network_in_road.index)]
traj_remain
同樣，生成對應的GeoDataFrame
geometry = [Point(xy) for xy in zip(traj_remain['longitude'], traj_remain['latitude'])]
traj_remain_gdf = gpd.GeoDataFrame(traj_remain, geometry=geometry)3.2 將經緯度坐標轉化為墨卡托坐標?
轉換成墨卡托坐標之后，兩個點之間的距離單位就是米了
# 轉換坐標系到UTM【橫軸墨卡托】
utm_projection = "EPSG:32648"
# 新加坡對應的EPSG代碼# 設置原始CRS為WGS 84 (EPSG:4326)
traj_remain_gdf.set_crs("EPSG:4326", inplace=True)
#這是GPS數據常用的坐標系統，其EPSG代碼為4326road_network_utm = road_network.to_crs(utm_projection)
traj_remain_utm = traj_remain_gdf.to_crs(utm_projection)
3.3 獲取距離
from shapely.ops import nearest_points
import pandas as pd# 創建一個空列表來存儲距離
distances = []# 計算距離
for point in traj_remain_utm.geometry:#遍歷每一個用戶點nearest_geom_index = list(road_network_utm.sindex.nearest(point, 1))[1]nearest_geom = road_network_utm.geometry.iloc[nearest_geom_index]# 獲取最近的路段（使用空間索引）distance = point.distance(nearest_geom)distances.append(distance.values[0])# 計算并存儲距離traj_remain_utm['distance_to_nearest_road'] = distances
# 將距離列表添加到outdoor_traj_not_in_network_utm DataFrametraj_remain_utm['distance_to_nearest_road'].describe()
'''
count 330825.000000
mean 29.847753
std 65.107624
min 1.106306
25% 3.725888
50% 9.576145
75% 44.843000
max 4582.239106
Name: distance_to_nearest_road, dtype: float64
'''
