Python实现点归属多边形匹配：关联命案坐标与社区SHP数据

一、先搞定空间连接全NaN的核心问题 #

空间连接返回全NaN，90%以上是坐标系不匹配或者点数据没转成正确的几何对象，按以下步骤排查修正：

1. 读取社区SHP文件 #

import geopandas as gpd
import pandas as pd
from shapely.geometry import Point

# 读取社区SHP，自动加载几何列和坐标系
communities = gpd.read_file("你的社区文件.shp")
# 先确认坐标系，比如输出 EPSG:4326 或 EPSG:3857
print("社区坐标系:", communities.crs)

2. 处理命案CSV并转为GeoDataFrame #

CSV里的坐标列（比如lon、lat）必须是数值类型，还要转成Point对象，关键是坐标系要和社区SHP完全一致：

# 读取命案CSV
murders = pd.read_csv("命案数据.csv")

# 强制转坐标列为数值，转失败的会变成NaN
murders["lon"] = pd.to_numeric(murders["lon"], errors="coerce")
murders["lat"] = pd.to_numeric(murders["lat"], errors="coerce")

# 删掉坐标无效的行（避免后续报错）
murders = murders.dropna(subset=["lon", "lat"])

# 生成几何列，注意顺序是 (lon, lat)，对应GPS常用的WGS84坐标系
murders["geometry"] = murders.apply(lambda row: Point(row["lon"], row["lat"]), axis=1)

# 转为GeoDataFrame，先设置初始坐标系（假设你的CSV坐标是WGS84）
murders_gdf = gpd.GeoDataFrame(murders, crs="EPSG:4326")

# 核心操作：把命案数据的坐标系转成和社区SHP一模一样
murders_gdf = murders_gdf.to_crs(communities.crs)

3. 执行空间连接（点在多边形内） #

用predicate="within"明确判断逻辑，避免默认参数的歧义：

# 空间连接，保留命案所有列，添加社区的属性列（比如社区名称，假设列名是community_name）
result = gpd.sjoin(murders_gdf, communities, how="left", predicate="within")

如果还是有NaN，检查：

• 是不是把lat和lon写反了？比如应该是Point(lon, lat)你写成了Point(lat, lon)
• 社区SHP的多边形是否有效？可以用communities.is_valid检查，无效的话用communities = communities.buffer(0)修复

二、解决添加列后无法写入CSV的问题 #

写入CSV时要注意避开GeoDataFrame的几何列（CSV不支持几何对象），同时确保路径有权限：

# 先删掉几何列（不需要的话）
result_for_csv = result.drop(columns=["geometry"])

# 写入CSV，指定编码避免乱码，不要写索引列
result_for_csv.to_csv("带社区的命案数据.csv", index=False, encoding="utf-8-sig")

如果还是写不进去：

• 换个输出路径试试（比如不要写到C盘根目录，用用户文件夹）
• 可以加errors="replace"参数处理特殊字符：to_csv(..., errors="replace")

完整可运行示例代码 #

import geopandas as gpd
import pandas as pd
from shapely.geometry import Point

# 1. 读取社区边界SHP
communities = gpd.read_file("社区边界.shp")
print("社区坐标系:", communities.crs)

# 2. 处理命案坐标CSV
murders = pd.read_csv("命案坐标.csv")
# 转换坐标列为数值
murders["lon"] = pd.to_numeric(murders["lon"], errors="coerce")
murders["lat"] = pd.to_numeric(murders["lat"], errors="coerce")
# 过滤无效坐标行
murders = murders.dropna(subset=["lon", "lat"])
# 创建几何列
murders["geometry"] = murders.apply(lambda x: Point(x["lon"], x["lat"]), axis=1)
# 转为GeoDataFrame并匹配坐标系
murders_gdf = gpd.GeoDataFrame(murders, crs="EPSG:4326").to_crs(communities.crs)

# 3. 空间连接匹配社区
joined_data = gpd.sjoin(murders_gdf, communities, how="left", predicate="within")

# 4. 写入最终CSV
joined_data.drop(columns=["geometry"]).to_csv("命案带社区信息.csv", index=False, encoding="utf-8-sig")

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