使用OpenStreetMap和DBSCAN空间聚类捕捉最受瞩目的城市区域的通用框架
在本文中，我展示了一种快速且易于使用的方法，该方法能够根据从OpenStreetMap(OSM)收集的兴趣点(POI)来确定给定兴趣的热点。我首先收集 ChatGPT 上找到的几个类别的 POI 的原始数据，我假设它们代表某种时尚生活方式(例如咖啡馆、酒吧、市场、瑜伽工作室)；在将数据转换为便捷的 GeoDataFrame 后，我进行地理空间聚类，最后根据每个群集中不同的城市功能混合程度评估结果。
尽管我所称之为 “时髦”的主题和与之相关的 POI 类别有些主观性，但它们可以很容易地被其他主题和类别取代 —— 自动热点检测的方法保持不变。这种易于采纳的方法的优点包括：识别支持创新规划的本地创新中心，发现支持城市规划倡议的城市子中心，评估不同的商业市场机会，分析房地产投资机会或捕捉旅游热点。
所有图片均为作者创建。
1. 从OSM获取数据
首先，我获取目标城市的行政多边形。因为布达佩斯是我的家乡，所以为了便于 (实地) 验证，我使用了布达佩斯的数据。然而，由于我只使用了全球数据库的 OSM，因此这些步骤可以轻松地重现在 OSM 覆盖的世界其他地区。具体来说，我使用 OSMNx 包轻松地获取行政边界。
import osmnx as ox # 版本: 1.0.1city = 'Budapest'admin = ox.geocode_to_gdf(city)admin.plot()此代码块的结果：
现在，使用 OverPass API 下载落在布达佩斯行政边界范围内的 POI。在 amenity_mapping 列表中，我编制了一个与时尚生活方式相关联的 POI 类别列表。在此要注意的是，这是一个模糊且不基于专家的分类，使用本文介绍的方法，任何人都可以据此更新类别列表。此外，可以结合其他包含更精细多级分类的 POI 数据源，以更准确地描述给定主题。换句话说，这个列表可以根据您的需要进行任何改变 —— 从更好地涵盖时髦事物到调整此练习以适应其他主题分类(例如美食广场、购物区、旅游热点等)
注意：由于 OverPass 下载器返回了边界框内的所有结果，因此在此代码块的末尾，我使用 GeoPandas 的交叠功能过滤掉了边界之外的 POI。
import overpy # 版本: 0.6from shapely.geometry import Point # 版本: 1.7.1import geopandas as gpd # 版本: 0.9.0# 启动 APIapi = overpy.Overpass()# 获取包围边界框minx, miny, maxx, maxy = admin.to_crs(4326).bounds.T[0]bbox = ','.join([str(miny), str(minx), str(maxy), str(maxx)])# 定义感兴趣的 OSM 类别amenity_mapping = [ ("amenity", "cafe"), ("tourism", "gallery"), ("amenity", "pub"), ("amenity", "bar"), ("amenity", "marketplace"), ("sport", "yoga"), ("amenity", "studio"), ("shop", "music"), ("shop", "second_hand"), ("amenity", "foodtruck"), ("amenity", "music_venue"), ("shop", "books"),]# 遍历所有类别，调用 overpass api，并将结果添加到 poi_data 列表poi_data = []for idx, (amenity_cat, amenity) in enumerate(amenity_mapping): query = f"""node["{amenity_cat}"="{amenity}"]({bbox});out;""" result = api.query(query) print(amenity, len(result.nodes)) for node in result.nodes: data = {} name = node.tags.get('name', 'N/A') data['name'] = name data['amenity'] = amenity_cat + '__' + amenity data['geometry'] = Point(node.lon, node.lat) poi_data.append(data) # 将结果转换为 GeoDataFramegdf_poi = gpd.GeoDataFrame(poi_data)print(len(gdf_poi))gdf_poi = gpd.overlay(gdf_poi, admin[['geometry']])gdf_poi.crs = 4326print(len(gdf_poi))此代码块的结果是每个已下载的POI类别的频率分布：
2. 可视化POI数据
现在，可视化所有2101个POI：
import matplotlib.pyplot as pltf, ax = plt.subplots(1,1,figsize=(10,10))admin.plot(ax=ax, color = 'none', edgecolor = 'k', linewidth = 2)gdf_poi.plot(column = 'amenity', ax=ax, legend = True, alpha = 0.3)这段代码的结果：
这个图形很难解释 – 除了市中心非常拥挤，因此让我们尝试使用一个交互式可视化工具，Folium。
import foliumimport branca.colormap as cm# 获取城市的质心并设置地图x, y = admin.geometry.to_list()[0].centroid.xym = folium.Map(location=[y[0], x[0]], zoom_start=12, tiles='CartoDB Dark_Matter')colors = ['blue', 'green', 'red', 'purple', 'orange', 'pink', 'gray', 'cyan', 'magenta', 'yellow', 'lightblue', 'lime']# 转换gdf_poi的颜色设置amenity_colors = {}unique_amenities = gdf_poi['amenity'].unique()for i, amenity in enumerate(unique_amenities): amenity_colors[amenity] = colors[i % len(colors)]# 使用散点图可视化POIsfor idx, row in gdf_poi.iterrows(): amenity = row['amenity'] lat = row['geometry'].y lon = row['geometry'].x color = amenity_colors.get(amenity, 'gray') # 在颜色映射中找不到时，默认为灰色 folium.CircleMarker( location=[lat, lon], radius=3, color=color, fill=True, fill_color=color, fill_opacity=1.0, # 对于点标记没有透明度 popup=amenity, ).add_to(m)# 显示地图m这个地图的默认视图（你可以通过调整zoom_start=12参数来改变）：
然后，可以改变缩放参数并重新绘制地图，或使用鼠标放大：
或完全缩小：
3. 空间聚类
现在我手头有所有必要的POI了，我会使用DBSCAN算法，首先编写一个函数来对POI进行聚类。我只会调整DBSDCAN的eps参数，该参数实际上计量了聚类的特征尺寸，即要将哪些POI归为一组的距离。另外，我会将几何信息转换为本地CRS（EPSG:23700），以便使用国际单位进行工作。有关CRS转换的更多信息请见此处。
from sklearn.cluster import DBSCAN # 版本: 0.24.1from collections import Counter# 进行聚类def apply_dbscan_clustering(gdf_poi, eps): feature_matrix = gdf_poi['geometry'].apply(lambda geom: (geom.x, geom.y)).tolist() dbscan = DBSCAN(eps=eps, min_samples=1) # 可根据需要调整min_samples cluster_labels = dbscan.fit_predict(feature_matrix) gdf_poi['cluster_id'] = cluster_labels return gdf_poi# 转换为本地CRSgdf_poi_filt = gdf_poi.to_crs(23700) # 进行聚类eps_value = 50 clustered_gdf_poi = apply_dbscan_clustering(gdf_poi_filt, eps_value)# 打印带有聚类标识的GeoDataFrameprint('找到的聚类数量：', len(set(clustered_gdf_poi.cluster_id)))clustered_gdf_poi这个单元格的结果：
这里有1237个聚类——如果我们只关注舒适、时髦的热点地区，这个数量似乎有点过多。让我们先看一下聚类的大小分布，然后选择一个大小阈值——把只有两个POI的聚类称为热点可能并不合理。
clusters = clustered_gdf_poi.cluster_id.to_list()clusters_cnt = Counter(clusters).most_common()f, ax = plt.subplots(1,1,figsize=(8,4))ax.hist([cnt for c, cnt in clusters_cnt], bins = 20)ax.set_yscale('log')ax.set_xlabel('聚类大小', fontsize = 14)ax.set_ylabel('聚类数量', fontsize = 14)这个单元格的结果：
基于直方图中的间隙，让我们保留至少有10个POI的聚类！目前来说，这个工作假设足够简单了。然而，也可以通过更复杂的方式来实现，例如结合不同POI类型的数量或覆盖的地理区域。
to_keep = [c for c, cnt in Counter(clusters).most_common() if cnt>9]clustered_gdf_poi = clustered_gdf_poi[clustered_gdf_poi.cluster_id.isin(to_keep)]clustered_gdf_poi = clustered_gdf_poi.to_crs(4326)len(to_keep)这段代码显示满足筛选条件的聚类数量为15。
一旦我们有了这15个真正的时髦聚类，将它们放在地图上：
import foliumimport random# 获取城市的质心并设置地图的位置min_longitude, min_latitude, max_longitude, max_latitude = clustered_gdf_poi.total_boundsm = folium.Map(location=[(min_latitude+max_latitude)/2, (min_longitude+max_longitude)/2], zoom_start=14, tiles='CartoDB Dark_Matter')# 为每个聚类获取唯一的随机颜色unique_clusters = clustered_gdf_poi['cluster_id'].unique()cluster_colors = {cluster: "#{:02x}{:02x}{:02x}".format(random.randint(0, 255), random.randint(0, 255), random.randint(0, 255)) for cluster in unique_clusters}# 可视化POIfor idx, row in clustered_gdf_poi.iterrows(): lat = row['geometry'].y lon = row['geometry'].x cluster_id = row['cluster_id'] color = cluster_colors[cluster_id] # 创建一个点标记 folium.CircleMarker( location=[lat, lon], radius=3, color=color, fill=True, fill_color=color, fill_opacity=0.9, popup=row['amenity'], ).add_to(m)# 显示地图m
4.比较聚类
每个聚类都被视为一个时髦的聚类 – 然而，它们在某种程度上都必须是独特的，对吧？让我们通过比较它们所提供的POI类别组合的多样性来看看它们有多独特。
首先，注重多样性，并通过计算它们的熵来度量每个聚类中POI类别的多样性/变化。
import mathimport pandas as pddef get_entropy_score(tags): tag_counts = {} total_tags = len(tags) for tag in tags: if tag in tag_counts: tag_counts[tag] += 1 else: tag_counts[tag] = 1 tag_probabilities = [count / total_tags for count in tag_counts.values()] shannon_entropy = -sum(p * math.log(p) for p in tag_probabilities) return shannon_entropy# 创建一个字典，其中每个聚类都有自己的POI列表clusters_amenities = clustered_gdf_poi.groupby(by = 'cluster_id')['amenity'].apply(list).to_dict()# 计算并存储熵分数entropy_data = []for cluster, amenities in clusters_amenities.items(): E = get_entropy_score(amenities) entropy_data.append({'cluster' : cluster, 'size' :len(amenities), 'entropy' : E}) # 将熵分数添加到数据帧entropy_data = pd.DataFrame(entropy_data)entropy_data此代码块的结果：
接下来是对此表进行快速相关分析：
entropy_data.corr()
在计算了聚类ID、聚类大小和聚类熵之间的相关性后，大小和熵之间存在显著的相关性；然而，这远不能解释所有的多样性。显然，有些热点比其他热点更具多样性，而其他热点则更专注于某些特定领域。它们专门从事什么？通过将每个聚类的POI配置文件与聚类中每种POI类型的整体分布进行比较，并选择与平均水平相比最典型的三个POI类别来回答这个问题。
# 将poi配置文件打包到字典中clusters = sorted(list(set(clustered_gdf_poi.cluster_id)))amenity_profile_all = dict(Counter(clustered_gdf_poi.amenity).most_common())amenity_profile_all = {k : v / sum(amenity_profile_all.values()) for k, v in amenity_profile_all.items()}# 计算每个聚类的相对频率# 仅保留高于平均值（>1）的前三个候选项clusters_top_profile = {}for cluster in clusters: amenity_profile_cls = dict(Counter(clustered_gdf_poi[clustered_gdf_poi.cluster_id == cluster].amenity).most_common() ) amenity_profile_cls = {k : v / sum(amenity_profile_cls.values()) for k, v in amenity_profile_cls.items()} clusters_top_amenities = [] for a, cnt in amenity_profile_cls.items(): ratio = cnt / amenity_profile_all[a] if ratio>1: clusters_top_amenities.append((a, ratio)) clusters_top_amenities = sorted(clusters_top_amenities, key=lambda tup: tup[1], reverse=True) clusters_top_amenities = clusters_top_amenities[0:min([3,len(clusters_top_amenities)])] clusters_top_profile[cluster] = [c[0] for c in clusters_top_amenities] # 对于每个聚类，打印其顶级类别：for cluster, top_amenities in clusters_top_profile.items(): print(cluster, top_amenities)这段代码的结果：
顶级分类描述已经显示出一些趋势。例如，聚类17明显是用于饮酒，而聚类19也涉及音乐，可能是与派对有关。拥有书店、画廊和咖啡馆的聚类91显然是白天放松的地方，而聚类120中的音乐和画廊可能是任何酒吧爬行的良好热身。从分布中，我们还可以看出，在酒吧里喝酒总是合适的（或者，根据使用情况，我们应该考虑基于类别频率的进一步归一化）！
结论
作为一名当地居民，我可以确认，尽管使用了简单的方法，但这些聚类是完全合理的，并且很好地代表了所需的城市功能混合。当然，这是一个可以通过以下几种方式进行丰富和改进的快速试点项目：
- 依靠更详细的POI分类和选择
- 在进行聚类时考虑POI类别（语义聚类）
- 通过社交媒体评论和评分等方式丰富POI信息