使用adm2 shapefile时地图着色多边形与预期不符的问题排查

问题描述 #

将风险数据按等级在地图上着色，读取shapefile与csv数据合并后，adm1层级的地图显示完全正常，但adm2层级出现严重偏差：着色的多边形与对应数据的位置错位（例如应为编号42的区域着色，实际却渲染了编号41的区域）。

使用的Python代码 #

#!/home/zmumba/anaconda3/bin/python

import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
from mpl_toolkits.basemap import Basemap  # optional ?

map_df = gpd.read_file("/home/zmumba/DA/Dekad_Progs/Shapefiles/Lesotho/geoBoundaries-LSO-ADM2-all/geoBoundaries-LSO-ADM2.shx")
risks_df=pd.read_csv("/home/zmumba/DA/Dekad_Progs/Output/H1dRrisks.csv")

map_df["risk"] = map_df.merge(risks_df, left_on="shapeName", right_on="District")["risk"]

colors = {1: "green", 2: "yellow", 3: "orange", 4: "red"}  # or a list
labels = {1: "no risk", 2: "low risk", 3: "medium risk", 4: "high risk"}
catego = map_df["risk"].astype(str).str.cat(map_df["risk"].map(labels), sep="- ")

fig, ax = plt.subplots(figsize=(5, 5))
plt.title(f'Risk of Heavy 24hr Rain: 20-24Nov', y=1.04)

map_df.plot(
column=catego,
categorical=True,
edgecolor="k",
linewidths=0.8,
alpha=0.7,
cmap=ListedColormap([c for r,c in colors.items() if r in map_df["risk"].unique()]),
legend=True,
legend_kwds={
"title": "Risk Level",
"shadow": True,
"loc": "lower right",
"fontsize": 10,
},
ax=ax,
)
ax.set_axis_off()
plt.savefig('H1dRriskmap.png', dpi=300)

plt.show()

补充信息 #

• 偏差示例：编号42对应的数据实际着色到41区域，编号43对应的数据着色到42区域
• 多边形质心坐标（来自shapefile本身，排除坐标源问题）：

  -29.3892289999999,  28.3056183629316  -> 41
  -29.2870792999999,  29.2715247780569  -> 42
  -29.2255776170000,  27.6546474532047  -> 43
  

• CSV数据格式：

  "District","risk"
  "name1",1
  "name2",1
  ...
  "name78",1
  

  risk值为1-4，对应不同风险等级
• 尝试R语言解决时遇到包缺失问题（如提示no package sf）

核心问题定位 #

问题出在合并数据的代码行：

map_df["risk"] = map_df.merge(risks_df, left_on="shapeName", right_on="District")["risk"]

• merge()默认使用内连接，如果adm2的shapefile中存在shapeName在csv中无匹配项，或csv存在多余项，合并后的DataFrame行数会与原map_df不一致
• 直接赋值时，合并后的risk列会按行索引对齐，但合并后行数减少会导致原map_df未匹配行被填充NaN，匹配行的顺序也可能被打乱，最终造成风险值与多边形错位

解决步骤 #

1. 改用左连接保留所有多边形
   替换合并代码，确保原shapefile的所有区域都被保留：

   map_df = map_df.merge(risks_df, left_on="shapeName", right_on="District", how="left")
   

   此方式会保留map_df所有行，匹配到的区域填充risk值，未匹配的留空，避免行错位。

2. 验证名称匹配一致性
   合并后检查匹配情况，排查名称拼写、大小写或空格差异：

   # 查看未匹配的行
   print(map_df[map_df["risk"].isna()])
   # 检查shapeName与District的唯一值数量
   print("shapeName数量:", map_df["shapeName"].nunique())
   print("District数量:", risks_df["District"].nunique())
   

3. 修复分类标签生成逻辑
   处理未匹配数据的NaN情况，避免报错：

   map_df["risk_label"] = map_df["risk"].map(labels).fillna("未匹配数据")
   catego = map_df["risk"].astype(str).str.cat(map_df["risk_label"], sep="- ").fillna("未匹配数据")
   

4. 可选：重置索引确保对齐
   若合并后索引混乱，可重置索引：

   map_df = map_df.reset_index(drop=True)
   

额外建议 #

• 读取shapefile时直接读取.shp主文件，而非.shx索引文件：gpd.read_file("/path/to/geoBoundaries-LSO-ADM2.shp")，避免读取不完整
• R语言的包缺失问题，可通过install.packages("sf")安装对应包（需确保网络正常或使用镜像源）

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