Python通过Geopandas实现空间连接与地理计算

# 安装必要的库
!pip install geopandas matplotlib contextily rtree shapely
# 导入库
import geopandas as gpd
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import contextily as ctx
from shapely.geometry import point, linestring, polygon

# 加载行政区划数据
districts = gpd.read_file("districts.geojson")
# 创建poi点数据
poi_data = {
    'name': ['咖啡店a', '餐厅b', '商场c', '学校d', '医院e', '公园f', '图书馆g', '超市h'],
    'type': ['餐饮', '餐饮', '购物', '教育', '医疗', '休闲', '文化', '购物'],
    'latitude': [39.91, 39.92, 39.93, 39.89, 39.91, 39.88, 39.90, 39.94],
    'longitude': [116.41, 116.42, 116.39, 116.38, 116.44, 116.43, 116.41, 116.40]
}
# 创建点几何
geometry = [point(xy) for xy in zip(poi_data['longitude'], poi_data['latitude'])]
pois = gpd.geodataframe(poi_data, geometry=geometry, crs="epsg:4326")

# 确保坐标系统一致
districts = districts.to_crs(pois.crs)
# 执行空间连接 - 点位于哪个区域
joined_data = gpd.sjoin(pois, districts, how="left", predicate="within")
print(f"连接后的数据包含 {len(joined_data)} 行，其中包括点位信息和所在区域")

# 绘制底图
fig, ax = plt.subplots(figsize=(12, 10))
# 绘制区域边界
districts.boundary.plot(ax=ax, color='black', linewidth=1)
# 使用分类变量绘制不同类型的poi
pois.plot(ax=ax, column='type', cmap='viridis', 
          legend=true, markersize=100, categorical=true)
# 添加标题
plt.title('poi分布与行政区划', fontsize=15)
# 添加背景底图
ctx.add_basemap(ax, crs=pois.crs.to_string(), source=ctx.providers.cartodb.positron)
plt.tight_layout()
plt.show()

# 使用不同空间关系进行连接
within_join = gpd.sjoin(pois, districts, how="inner", predicate="within")
intersects_join = gpd.sjoin(pois, districts, how="inner", predicate="intersects")
print(f"within连接结果: {len(within_join)} 行")
print(f"intersects连接结果: {len(intersects_join)} 行")
# 对于点数据，"within"和"intersects"通常产生相同结果
# 但对于线或面数据则不同

# 筛选特定类型的poi并确定它们所在的行政区
restaurant_pois = pois[pois['type'] == '餐饮']
restaurant_districts = gpd.sjoin(restaurant_pois, districts, how="left", predicate="within")
# 统计每个行政区内餐饮poi的数量
poi_count_by_district = restaurant_districts.groupby('district_name').size()
print("各行政区餐饮poi数量统计:")
print(poi_count_by_district)

# 为每个poi创建500米缓冲区
# 首先将数据转换为投影坐标系统以使用米为单位
pois_projected = pois.to_crs(epsg=3857)  # web mercator投影
pois_buffered = pois_projected.copy()
pois_buffered['geometry'] = pois_projected.buffer(500)  # 500米缓冲区
# 转回原始坐标系进行显示
pois_buffered = pois_buffered.to_crs(pois.crs)
# 绘制缓冲区
fig, ax = plt.subplots(figsize=(12, 10))
districts.boundary.plot(ax=ax, color='black', linewidth=1)
pois_buffered.plot(ax=ax, alpha=0.5, column='type', categorical=true)
pois.plot(ax=ax, color='red', markersize=20)
ctx.add_basemap(ax, crs=pois.crs.to_string(), source=ctx.providers.cartodb.positron)
plt.title('poi点位500米服务范围分析', fontsize=15)
plt.tight_layout()
plt.show()

# 将缓冲区与行政区划进行交集运算
districts_projected = districts.to_crs(epsg=3857)
pois_buffered_projected = pois_buffered.to_crs(epsg=3857)
# 选择一个特定区域进行示例
target_district = districts_projected.iloc[0]
# 计算该区域与各poi缓冲区的交集
intersections = []
for idx, poi_buffer in pois_buffered_projected.iterrows():
    if poi_buffer.geometry.intersects(target_district.geometry):
        intersection = poi_buffer.geometry.intersection(target_district.geometry)
        intersections.append({
            'poi_name': poi_buffer['name'],
            'poi_type': poi_buffer['type'],
            'geometry': intersection,
            'area': intersection.area  # 平方米
        })
# 创建交集geodataframe
intersections_gdf = gpd.geodataframe(intersections, crs=pois_buffered_projected.crs)
# 计算覆盖率
district_area = target_district.geometry.area
total_coverage_area = gpd.geoseries([geom['geometry'] for geom in intersections]).unary_union.area
coverage_ratio = total_coverage_area / district_area * 100
print(f"区域名称: {target_district['district_name']}")
print(f"区域总面积: {district_area/1e6:.2f} 平方公里")
print(f"poi服务覆盖面积: {total_coverage_area/1e6:.2f} 平方公里")
print(f"覆盖率: {coverage_ratio:.2f}%")

# 计算所有poi之间的距离矩阵（以米为单位）
def calculate_distance_matrix(gdf):
    n = len(gdf)
    dist_matrix = np.zeros((n, n))
    for i in range(n):
        for j in range(i+1, n):
            dist = gdf.iloc[i].geometry.distance(gdf.iloc[j].geometry)
            dist_matrix[i, j] = dist
            dist_matrix[j, i] = dist
    return pd.dataframe(dist_matrix, 
                        index=gdf.index, 
                        columns=gdf.index)
# 计算距离矩阵
distance_matrix = calculate_distance_matrix(pois_projected)
# 查找每个poi的最近邻
nearest_neighbors = {}
for idx in pois.index:
    distances = distance_matrix[idx].drop(idx)  # 排除自身
    nearest = distances.idxmin()
    nearest_neighbors[pois.loc[idx, 'name']] = {
        'nearest_poi': pois.loc[nearest, 'name'],
        'distance_meters': distances.min()
    }
print("每个poi的最近邻设施:")
for poi, info in nearest_neighbors.items():
    print(f"{poi} -> {info['nearest_poi']} (距离: {info['distance_meters']:.2f}米)")

# 使用kde进行热点分析
from scipy.stats import gaussian_kde
# 提取poi坐标
x = pois_projected.geometry.x
y = pois_projected.geometry.y
positions = np.vstack([x, y])
# 计算kde
kernel = gaussian_kde(positions)
# 创建网格进行评估
x_grid, y_grid = np.mgrid[x.min():x.max():100j, y.min():y.max():100j]
positions_grid = np.vstack([x_grid.ravel(), y_grid.ravel()])
# 计算密度
z = kernel(positions_grid).reshape(x_grid.shape)
# 可视化结果
fig, ax = plt.subplots(figsize=(12, 10))
districts_projected.boundary.plot(ax=ax, color='black', linewidth=1)
contour = ax.contourf(x_grid, y_grid, z, cmap='viridis', levels=20, alpha=0.75)
pois_projected.plot(ax=ax, color='red', markersize=20)
plt.colorbar(contour, ax=ax, label='密度')
plt.title('poi热点密度分析', fontsize=15)
plt.tight_layout()
plt.show()

from sklearn.cluster import dbscan
# 提取坐标
coords = np.vstack((pois_projected.geometry.x, pois_projected.geometry.y)).t
# 使用dbscan进行聚类
clustering = dbscan(eps=1000, min_samples=2).fit(coords)  # 1000米范围内至少2个点形成聚类
# 将聚类结果添加到原始数据
pois_projected['cluster'] = clustering.labels_
# 可视化聚类结果
fig, ax = plt.subplots(figsize=(12, 10))
districts_projected.boundary.plot(ax=ax, color='black', linewidth=1)
pois_projected.plot(ax=ax, column='cluster', categorical=true, 
                   legend=true, markersize=100)
plt.title('poi空间聚类分析 (dbscan)', fontsize=15)
plt.tight_layout()
plt.show()
# 统计每个聚类的poi数量
cluster_stats = pois_projected.groupby('cluster').size()
print("各聚类poi数量:")
print(cluster_stats)

# 注意：这需要预先准备道路网络数据和networkx库
import networkx as nx
# 加载道路网络数据
roads = gpd.read_file("roads.geojson")
# 创建网络图
g = nx.graph()
# 从道路数据中构建网络
for idx, road in roads.iterrows():
    coords = list(road.geometry.coords)
    for i in range(len(coords) - 1):
        g.add_edge(coords[i], coords[i+1], 
                  weight=road.geometry.length,
                  road_id=road['id'])
# 找到距离每个poi最近的道路节点
nearest_nodes = {}
for idx, poi in pois_projected.iterrows():
    min_dist = float('inf')
    nearest_node = none
    for node in g.nodes():
        dist = poi.geometry.distance(point(node))
        if dist < min_dist:
            min_dist = dist
            nearest_node = node
    nearest_nodes[idx] = nearest_node
# 计算两个poi之间的最短路径
origin_idx = 0  # 起点poi索引
dest_idx = 3    # 终点poi索引
origin_node = nearest_nodes[origin_idx]
dest_node = nearest_nodes[dest_idx]
try:
    # 计算最短路径
    shortest_path = nx.shortest_path(g, origin_node, dest_node, weight='weight')
    path_length = nx.shortest_path_length(g, origin_node, dest_node, weight='weight')
    print(f"从 {pois.loc[origin_idx, 'name']} 到 {pois.loc[dest_idx, 'name']} 的最短路径长度: {path_length:.2f}米")
    # 可视化路径
    path_coords = shortest_path
    path_line = linestring(path_coords)
    path_gdf = gpd.geodataframe({'geometry': [path_line]}, crs=pois_projected.crs)
    fig, ax = plt.subplots(figsize=(12, 10))
    roads.plot(ax=ax, color='gray', linewidth=1)
    path_gdf.plot(ax=ax, color='red', linewidth=3)
    pois_projected.iloc[[origin_idx, dest_idx]].plot(ax=ax, color='blue', markersize=100)
    plt.title(f'从 {pois.loc[origin_idx, "name"]} 到 {pois.loc[dest_idx, "name"]} 的最短路径', fontsize=15)
    plt.tight_layout()
    plt.show()
except nx.networkxnopath:
    print(f"无法找到从 {pois.loc[origin_idx, 'name']} 到 {pois.loc[dest_idx, 'name']} 的路径")

# 确保已安装rtree库
# pip install rtree
# 创建空间索引
districts = districts.set_index('district_id', drop=false)
districts_sindex = districts.sindex
# 使用空间索引加速查询
def fast_sjoin_points_to_polygons(points_gdf, polygons_gdf):
    # 创建结果列表
    results = []
    # 遍历点
    for idx, point in points_gdf.iterrows():
        # 使用空间索引查找可能相交的多边形
        possible_matches_idx = list(polygons_gdf.sindex.intersection(point.geometry.bounds))
        possible_matches = polygons_gdf.iloc[possible_matches_idx]
        # 精确检查点是否在多边形内
        precise_matches = possible_matches[possible_matches.contains(point.geometry)]
        # 如果找到匹配项
        if len(precise_matches) > 0:
            for match_idx, match in precise_matches.iterrows():
                result = point.to_dict()
                result.update({
                    'district_id': match['district_id'],
                    'district_name': match['district_name']
                })
                results.append(result)
    if results:
        return gpd.geodataframe(results, crs=points_gdf.crs)
    else:
        return gpd.geodataframe([], crs=points_gdf.crs)
# 使用优化的函数进行空间连接
import time
start_time = time.time()
optimized_result = fast_sjoin_points_to_polygons(pois, districts)
end_time = time.time()
print(f"优化后的空间连接耗时: {end_time - start_time:.4f}秒")
# 对比标准sjoin
start_time = time.time()
standard_result = gpd.sjoin(pois, districts, how="left", predicate="within")
end_time = time.time()
print(f"标准空间连接耗时: {end_time - start_time:.4f}秒")

from joblib import parallel, delayed
# 并行处理缓冲区计算
def parallel_buffer(gdf, buffer_distance, n_jobs=-1):
    # 定义单个几何体缓冲区计算函数
    def buffer_geom(geom):
        return geom.buffer(buffer_distance)
    # 并行执行
    buffered_geoms = parallel(n_jobs=n_jobs)(
        delayed(buffer_geom)(geom) for geom in gdf.geometry
    )
    # 创建新的geodataframe
    result = gdf.copy()
    result['geometry'] = buffered_geoms
    return result
# 使用并行处理创建缓冲区
start_time = time.time()
pois_buffered_parallel = parallel_buffer(pois_projected, 500)
end_time = time.time()
print(f"并行缓冲区计算耗时: {end_time - start_time:.4f}秒")
# 对比标准处理
start_time = time.time()
pois_buffered_standard = pois_projected.copy()
pois_buffered_standard['geometry'] = pois_projected.buffer(500)
end_time = time.time()
print(f"标准缓冲区计算耗时: {end_time - start_time:.4f}秒")

# 假设有人口数据
population = gpd.read_file("population_grid.geojson")
# 选择交通设施poi
transport_pois = pois[pois['type'].isin(['公交站', '地铁站', '火车站'])]
# 创建交通设施缓冲区（500米步行距离）
transport_pois_projected = transport_pois.to_crs(epsg=3857)
transport_buffers = transport_pois_projected.copy()
transport_buffers['geometry'] = transport_pois_projected.buffer(500)
# 与人口网格进行空间连接
population_projected = population.to_crs(epsg=3857)
population_covered = gpd.sjoin(population_projected, transport_buffers.to_crs(population_projected.crs), 
                              how="inner", predicate="intersects")
# 计算覆盖的总人口
total_population = population_projected['population'].sum()
covered_population = population_covered['population'].sum()
coverage_percentage = (covered_population / total_population) * 100
print(f"总人口: {total_population}")
print(f"公共交通500米覆盖人口: {covered_population}")
print(f"覆盖率: {coverage_percentage:.2f}%")

# 假设我们有以下数据层
# - 人口密度
# - 竞争对手位置
# - 交通可达性
# - 商业中心
# 创建一个网格评价系统
from shapely.geometry import box
# 创建研究区域网格
xmin, ymin, xmax, ymax = districts_projected.total_bounds
cell_size = 500  # 500米网格
rows = int((ymax - ymin) / cell_size)
cols = int((xmax - xmin) / cell_size)
grid_cells = []
for i in range(cols):
    for j in range(rows):
        x1 = xmin + i * cell_size
        y1 = ymin + j * cell_size
        x2 = xmin + (i + 1) * cell_size
        y2 = ymin + (j + 1) * cell_size
        grid_cells.append(box(x1, y1, x2, y2))
grid = gpd.geodataframe({'geometry': grid_cells}, crs=districts_projected.crs)
# 计算各因素得分
# 1. 人口密度得分
grid['pop_score'] = np.random.randint(1, 10, size=len(grid))  # 模拟数据
# 2. 竞争对手距离得分
competitors = pois_projected[pois_projected['type'] == '餐饮']
grid['comp_dist'] = grid.geometry.apply(
    lambda g: min([g.distance(comp) for comp in competitors.geometry]) if len(competitors) > 0 else 0
)
grid['comp_score'] = grid['comp_dist'] / grid['comp_dist'].max() * 10
# 3. 交通可达性得分
transit = pois_projected[pois_projected['type'].isin(['公交站', '地铁站'])]
grid['transit_dist'] = grid.geometry.apply(
    lambda g: min([g.distance(t) for t in transit.geometry]) if len(transit) > 0 else float('inf')
)
grid['transit_score'] = 10 - (grid['transit_dist'] / grid['transit_dist'].max() * 10)
# 4. 商业中心距离得分
business_centers = pois_projected[pois_projected['type'] == '购物']
grid['biz_dist'] = grid.geometry.apply(
    lambda g: min([g.distance(bc) for bc in business_centers.geometry]) if len(business_centers) > 0 else float('inf')
)
grid['biz_score'] = 10 - (grid['biz_dist'] / grid['biz_dist'].max() * 10)
# 计算综合得分
grid['total_score'] = (
    grid['pop_score'] * 0.4 + 
    grid['comp_score'] * 0.2 + 
    grid['transit_score'] * 0.2 + 
    grid['biz_score'] * 0.2
)
# 可视化结果
fig, ax = plt.subplots(figsize=(12, 10))
districts_projected.boundary.plot(ax=ax, color='black', linewidth=1)
grid.plot(ax=ax, column='total_score', cmap='rdylgn', alpha=0.7, 
          legend=true, legend_kwds={'label': '选址适宜性得分'})
pois_projected.plot(ax=ax, color='blue', markersize=10)
plt.title('商业选址适宜性评价', fontsize=15)
plt.tight_layout()
plt.show()
# 找出得分最高的位置
best_locations = grid.sort_values('total_score', ascending=false).head(3)
print("最适宜开设新店铺的3个位置:")
for idx, loc in best_locations.iterrows():
    centroid = loc.geometry.centroid
    print(f"位置 {idx}: 坐标 ({centroid.x}, {centroid.y}), 得分: {loc['total_score']:.2f}")