基于Python对比GeoJSON文件，计算骑行路线的自行车道占比

我是Python地理空间数据处理新手，现有两个GeoJSON文件：

• 包含某城市所有自行车道的参考文件
• 同一城市内起点到终点的LineString类型骑行路线文件

需要对比这两个文件，识别骑行路线中属于自行车道的路段，进而计算该路线在自行车道上的占比。目前仅知道Python可读取GeoJSON，但找不到无需地图对比两个GeoJSON对象的方法。

已完成操作 #

1. 加载骑行路线数据 #

with open(file_path, encoding = 'utf-8') as file:
        data = json.load(file)
route = pd\
            .json_normalize(data['journeys'], record_path=['sections', 'geojson', 'coordinates'])\
            .rename(columns={0: 'x', 1: 'y'})

route['index'] = 3

gdf = gpd.GeoDataFrame(route, geometry=gpd.points_from_xy(route.x, route.y), crs = "EPSG:4326") # type: ignore
route = gdf.groupby(['index'])['geometry'].apply(lambda x:  LineString(x.tolist()))
route = gpd.GeoDataFrame(route, geometry='geometry', crs="EPSG:4326") # type: ignore
route.reset_index(inplace=True)
route.plot(column='index')

2. 加载巴黎自行车道参考文件 #

with open(ref_path, encoding = 'utf-8') as file:
    data = json.load(file)

ref = pd.json_normalize(data)[['objectid', 'num_pave', 'lib_level', 'lib_classe', 'geo_shape.geometry.coordinates']]
#ref = gpd.GeoDataFrame(ref) # type: ignore
ref = ref.dropna()
ref = ref.explode('geo_shape.geometry.coordinates').explode('geo_shape.geometry.coordinates').reset_index()
ref['x'] = ref['geo_shape.geometry.coordinates'].apply(lambda x: x[0])
ref['y'] = ref['geo_shape.geometry.coordinates'].apply(lambda x: x[1])

ref['type'] = ref['x'].apply(lambda x: 1 if type(x) == list else 0)
ref.loc[ref['type'] == 1, 'x'] = ref.loc[ref['type'] == 1, 'x'].apply(lambda x: x[0])
ref.loc[ref['type'] == 1, 'y'] = ref.loc[ref['type'] == 1, 'y'].apply(lambda x: x[1])

ref = gpd.GeoDataFrame(ref, geometry = gpd.points_from_xy(ref['x'], ref['y']), crs = "EPSG:4326") # type: ignore
ref = ref.groupby(['objectid'])['geometry'].apply(lambda x:  LineString(x.tolist()))
ref = gpd.GeoDataFrame(ref, geometry='geometry', crs="EPSG:4326") # type: ignore
ref.reset_index(inplace=True)
ref.plot(column='objectid')

3. 尝试叠加但仅得到多点结果 #

使用gpd.overlay得到了Multipoint结果，但无法计算骑行路线与自行车道的重叠长度：

intersection = gpd.overlay(route, ref, how = 'intersection', keep_geom_type = False)

输出结果示例：

>>>  index  objectid                                           geometry
0       3       133    MULTIPOINT (2.38773 48.84134, 2.38793 48.84129)
1       3       136    MULTIPOINT (2.38804 48.84127, 2.38799 48.84128)
2       3       218    MULTIPOINT (2.39365 48.84099, 2.39378 48.84100)
...

解决方案 #

核心要点：坐标系转换与线相交计算 #

EPSG:4326是经纬度坐标系，用它计算长度会得到以度为单位的不准确结果，必须先转换为平面坐标系（巴黎推荐用EPSG:2154，法国Lambert-93投影）。同时，gpd.overlay对LineString相交的处理不够直接，应使用intersection方法直接计算线与线的交集，获取重叠线段后再统计长度。

步骤1：转换坐标系 #

# 将骑行路线转换为平面坐标系
route_proj = route.to_crs("EPSG:2154")
# 将自行车道参考数据转换为平面坐标系
ref_proj = ref.to_crs("EPSG:2154")

步骤2：计算骑行路线总长度 #

total_length = route_proj['geometry'].length.sum()

步骤3：计算重叠线段总长度 #

遍历每条自行车道，计算其与骑行路线的交集，仅保留线段类型的交集并累加长度：

overlap_length = 0.0
# 提取骑行路线的单条LineString（假设路线是单条连续线段）
route_line = route_proj.iloc[0]['geometry']

for _, bike_lane in ref_proj.iterrows():
    intersect = route_line.intersection(bike_lane['geometry'])
    # 只处理线段类几何，排除点或空几何
    if intersect.geom_type == 'LineString':
        overlap_length += intersect.length
    elif intersect.geom_type == 'MultiLineString':
        # 处理多线段情况
        for segment in intersect.geoms:
            overlap_length += segment.length

步骤4：计算占比 #

if total_length > 0:
    percentage = (overlap_length / total_length) * 100
    print(f"骑行路线在自行车道上的占比：{percentage:.2f}%")
else:
    print("骑行路线长度为0，无法计算占比")

说明：为什么之前得到Multipoint？ #

骑行路线与自行车道可能是近似重合但不完全共线，或者采样点不匹配，导致gpd.overlay仅识别到离散交点而非连续线段。直接使用intersection方法能更准确处理线与线的空间关系，包括部分重叠的线段。

内容的提问来源于stack exchange，提问作者CBO