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

目录

• 引言
• 1. 环境配置
• 2. 空间连接基础概念
• 2.1 什么是空间连接
• 2.2 空间连接的类型
• 2.3 常见空间关系
• 3. 实际案例-城市POI与行政区划分析
• 3.1 加载数据
• 3.2 点与多边形的空间连接
• 3.3 结果可视化
• 4. 高级空间连接操作
• 4.1 不同谓词的空间连接
• 4.2 空间连接与属性筛选结合
• 5. 地理计算核心操作
• 5.1 几何操作-缓冲区分析
• 5.2 叠加分析-区域覆盖范围计算
• 5.3 距离计算与最近邻分析
• 6. 高级地理空间分析案例
• 6.1 热点密度分析
• 6.2 空间聚类分析
• 6.3 网络可达性分析
• 7. 性能优化技巧
• 7.1 空间索引加速查询
• 7.2 并行处理大规模数据
• 8. 实际应用场景
• 8.1 交通可达性分析
• 8.2 商业选址分析
• 9. 总结与进阶建议

在地理信息系统(GIS)领域，空间连接和地理计算是处理地理数据的核心能力。GeoPandas作为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

空间连接是根据地理对象之间的空间关系合并数据集的操作。与传统数据库连接使用共同的键值不同，空间连接使用几何关系（如相交、包含、覆盖）作为连接条件。

GeoPandas提供了多种空间连接方法：

• 左连接(left join): 保留左侧DataFrame的所有记录
• 右连接(right join): 保留右侧DataFrame的所有记录
• 内连接(inner join): 仅保留两侧都满足空间关系的记录

GeoPandas支持以下空间关系类型：

• intersects: 两个几何对象有任何交集
• contains: 一个几何对象完全包含另一个
• within: 一个几何对象完全在另一个内部
• touches: 两个几何对象的边界相接触但内部不相交
• crosses: 两个几何对象有部分但不完全重叠
• overlaps: 两个几何对象有部分重叠且维度相同

让我们加载一些示例数据：

# 加载行政区划数据
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")

我们将POI点与行政区划进行空间连接，确定每个POI点位于哪个行政区内：

# 确保坐标系统一致
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之间的距离矩阵和最近设施：

# 计算所有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()

使用DBSCAN算法进行空间聚类：

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)

假设我们有道路网络数据，可以分析POI的网络可达性：

# 注意：这需要预先准备道路网络数据和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']} 的路径")

GeoPandas使用R树索引来加速空间查询：

# 确保已安装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}")

在这篇文章中，我们深入探讨了GeoPandas中的空间连接和地理计算功能。我们学习了如何：

• 执行不同类型的空间连接操作
• 进行几何操作如缓冲区分析
• 计算空间距离和最近邻
• 进行热点分析和空间聚类
• 优化地理空间处理性能
• 应用这些技术到实际场景如交通分析和商业选址

进阶学习路径

如果你想更深入地学习地理空间分析，可以考虑以下方向：

• 空间统计学：探索空间自相关和地理加权回归等高级统计方法
• 网络分析：使用NetworkX和OSMnx进行复杂的路网分析
• 栅格数据处理：学习使用rasterio进行栅格数据处理
• 时空数据分析：探索时间和空间维度的结合分析
• 3D GIS：使用Python进行三维空间数据可视化和分析

推荐资源

• 《Python地理空间分析指南》(Bonnici, Erik)
• GeoPandas官方文档: https://geopandas.org/
• Shapely文档: https://shapely.readthedocs.io/
• 空间数据科学实用指南: https://geographicdata.science/book/intro.html

通过掌握这些空间连接和地理计算技术，你将能够处理各种复杂的地理空间数据分析任务，从简单的空间关系查询到复杂的多因素地理建模。

以上就是Python通过Geopandas实现空间连接与地理计算的详细内容，更多关于Python Geopandas空间连接与地理计算的资料请关注风君子博客其它相关文章！

您可能感兴趣的文章:

• Python通过Geopandas实现地理空间数据可视化
• python中geopandas库安装出现各种问题的解决办法
• python库geopandas读取写入空间数据及绘图实例探索
• Python实战之利用Geopandas算出每个省面积
• python geopandas读取、创建shapefile文件的方法