如何用Python/Geopandas裁剪交叉口中心区域的Geopandas几何图形？

以下是基于Python/Geopandas的具体实现步骤：

1. 加载数据并确认坐标系 #

首先加载你的道路数据集，确保其坐标系为EPSG:4326：

import geopandas as gpd

# 加载本地数据（支持shp、geojson等格式）
gdf = gpd.read_file("你的道路数据文件路径")
# 验证坐标系
print(gdf.crs)  # 预期输出：EPSG:4326

2. 转换为投影坐标系（关键步骤） #

EPSG:4326是地理坐标系（经纬度），直接计算距离或创建缓冲区会以度为单位，不符合实际需求。需转换为对应区域的投影坐标系（如UTM，单位为米）：

# 自动匹配数据所在区域的UTM坐标系
utm_crs = gdf.estimate_utm_crs()
gdf_utm = gdf.to_crs(utm_crs)

若自动匹配失效，可手动指定UTM代码（例如中国东部区域用EPSG:32650）。

3. 获取交叉口中心坐标 #

有两种方式获取中心：

• 自动计算整体中心（适合未知具体坐标的场景）

# 计算所有几何图形边界框的中心
center = gdf_utm.total_bounds[[0,1,2,3]].mean(axis=0)  # 返回 [x_center, y_center]

• 手动指定中心（已知交叉口经纬度时）

from shapely.geometry import Point

# 假设已知中心经纬度为 (lon, lat)
center_wgs84 = Point(lon, lat)
# 转换为UTM坐标系下的点
center_utm = gpd.GeoSeries([center_wgs84], crs="EPSG:4326").to_crs(utm_crs).iloc[0]
center = (center_utm.x, center_utm.y)

4. 创建裁剪范围缓冲区 #

围绕中心创建圆形或矩形裁剪范围，这里以50米半径的圆形为例：

from shapely.geometry import Point

# 创建中心点Shapely对象
center_point = Point(center)
# 创建缓冲区（单位为米，因使用UTM坐标系）
buffer = center_point.buffer(50)  # 可根据需求调整半径大小

若需要矩形范围，可直接构建多边形：

from shapely.geometry import Polygon

# 中心向四周扩展50米的矩形
min_x, min_y = center[0]-50, center[1]-50
max_x, max_y = center[0]+50, center[1]+50
buffer = Polygon([(min_x, min_y), (max_x, min_y), (max_x, max_y), (min_x, max_y)])

5. 执行裁剪操作 #

使用Geopandas内置方法完成裁剪：

# 方法1：用clip函数（简洁推荐）
clipped_gdf_utm = gpd.clip(gdf_utm, buffer)

# 方法2：用intersection操作（更灵活）
clipped_gdf_utm = gdf_utm[gdf_utm.geometry.intersects(buffer)]
clipped_gdf_utm["geometry"] = clipped_gdf_utm.geometry.intersection(buffer)

6. 转回原坐标系（可选） #

若需将结果恢复为EPSG:4326坐标系：

clipped_gdf = clipped_gdf_utm.to_crs("EPSG:4326")

7. 验证裁剪效果 #

通过绘图确认结果：

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(8,8))
clipped_gdf.plot(ax=ax, edgecolor='black', linewidth=2)
plt.show()

额外提示 #

• 可根据实际需求调整缓冲区半径，确保覆盖交叉口核心区域
• 裁剪后若出现空几何图形，可通过clipped_gdf = clipped_gdf[~clipped_gdf.geometry.is_empty]过滤

内容的提问来源于stack exchange，提问作者Paulo Henrique PH