GeoPandas裁剪哥德堡shapefile后绘图报IndexError数组维度错误

操作流程 #

1. 提取哥德堡行政边界 #

# 提取哥德堡边界
boundary = adm_2.loc[(adm_2["NAME_2"] == "Göteborg")].copy()

print(boundary)
# 输出
#     ID_0    ISO NAME_0  ID_1    NAME_1  ID_2    NAME_2  TYPE_2  ENGTYPE_2   NL_NAME_2   VARNAME_2   geometry
# 246 222 SWE Sweden  21  Västra Götaland 247 Göteborg    Kommuner    Municipality    None    None    MULTIPOLYGON (((12.09368 57.67838, 12.09466 57...

# 裁剪矢量文件
roads = geopandas.clip(roads_sweden, boundary)
buildings = geopandas.clip(buildings_sweden, boundary)

# print(roads.shape) # 输出(9747, 7)
# print(roads_sweden.shape) # 输出(148403, 7)

# print(buildings.shape) # 输出(1163, 4)
# print(buildings_sweden.shape) # 输出(15253, 4)

2. 问题描述 #

裁剪后的数据结构和原瑞典全量道路、建筑shapefile完全一致，但绘图时抛出错误：IndexError: too many indices for array: array is 1-dimensional, but 2 were indexed，单独调用roads.plot()也会报错，原全量数据可正常绘图，尝试重置索引后问题仍未解决。

绘图代码

fig, ax = plt.subplots(figsize=(12, 10), facecolor='white', edgecolor='k')
boundary.plot(ax=ax, facecolor= "white", edgecolor="black", linewidth=.5)
roads.plot(ax=ax)
buildings.plot(ax=ax)

数据预览

print(roads_sweden.head())
# 输出
# osm_id  name    ref type    oneway  maxspeed    geometry
# 0   1240    Klensmedsvägen None    unclassified    0   0   LINESTRING (17.99027 59.29686, 17.99182 59.296...
# 1   1241    Hyvelvägen None    residential 0   0   LINESTRING (17.99273 59.29666, 17.99212 59.295...
# 2   1242    Spikvägen  None    residential 0   0   LINESTRING (17.99353 59.29640, 17.99299 59.295...
# 3   1243    Bultvägen  None    residential 0   0   LINESTRING (17.99439 59.29617, 17.99385 59.295...
# 4   1245    Tångvägen None    residential 0   0   LINESTRING (17.99691 59.29512, 17.99637 59.294...

print(roads.head())
# 输出
# osm_id  name    ref type    oneway  maxspeed    geometry
# 0   3846609 None    None    motorway_link   1   0   LINESTRING (12.09100 57.67369, 12.09090 57.673...
# 1   4040303 E6  None    motorway    0   0   LINESTRING (12.00496 57.84283, 12.00525 57.841...
# 2   4040436 None    None    motorway_link   1   0   LINESTRING (12.00643 57.79800, 12.00560 57.796...
# 3   4040439 None    E 20    motorway    1   0   LINESTRING (11.99463 57.71530, 11.99492 57.715...
# 4   4040441 None    E 20    motorway    1   0   LINESTRING (11.99520 57.71580, 11.99434 57.71536)

解决方案 #

该错误通常由裁剪后产生的空几何对象、坐标参考系（CRS）不匹配导致，按以下步骤排查即可解决：

1. 统一所有数据的坐标参考系
   裁剪前先确认行政边界和待裁剪的矢量数据CRS一致，不一致时先做投影转换：

# 确认CRS一致
assert boundary.crs == roads_sweden.crs == buildings_sweden.crs, "CRS不匹配，请先统一投影"
# 不一致时执行转换，以边界的CRS为准
roads_sweden = roads_sweden.to_crs(boundary.crs)
buildings_sweden = buildings_sweden.to_crs(boundary.crs)

2. 删除裁剪后的空几何对象
   边界裁剪过程中可能会生成无效的空几何，直接删除即可：

roads = roads[~roads.is_empty].reset_index(drop=True)
buildings = buildings[~buildings.is_empty].reset_index(drop=True)

3. 替换裁剪方法（可选）
   如果clip方法仍报错，可改用overlay方法实现交集裁剪，兼容性更好：

roads = geopandas.overlay(roads_sweden, boundary, how='intersection')
buildings = geopandas.overlay(buildings_sweden, boundary, how='intersection')

4. 绘图时指定固定颜色（可选）
   部分版本的geopandas会在默认颜色映射时出错，绘图时显式指定颜色可避免该问题：

fig, ax = plt.subplots(figsize=(12, 10), facecolor='white', edgecolor='k')
boundary.plot(ax=ax, facecolor= "white", edgecolor="black", linewidth=.5)
roads.plot(ax=ax, color='#999999', linewidth=0.3)
buildings.plot(ax=ax, color='#cccccc', edgecolor='#666666', linewidth=0.2)
plt.show()

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