GeoPandas教程：空间连接sjoin怎么用？（附：空间索引优化技巧）

wowwwai GIS研习社 · 工具流程与项目排障

《GeoPandas教程：空间连接sjoin怎么用？（附：空间索引优化技巧）》这篇文章解决一个很常见的 GIS 数据处理问题：手里有点、线、面两类空间数据，想判断“点落在哪个面里”“道路穿过哪些行政区”“地块与规划区是否相交”，应该如何用 GeoPandas 的 sjoin 做空间连接，并尽量避免运行很慢。

引言：GeoPandas空间连接sjoin适合解决什么问题

在 GIS 工作中，普通表连接通常依赖一个共同字段，例如行政区代码、地块编号或站点 ID。空间连接不同，它依赖的是几何之间的空间关系，例如包含、相交、邻近、覆盖等。

GeoPandas 的 sjoin 就是 Python GIS 中最常用的空间连接工具之一。它可以把一个 GeoDataFrame 的属性，根据空间关系连接到另一个 GeoDataFrame 上。

典型场景包括：

• 把 POI 点位匹配到所在街道、区县或网格。
• 统计每个行政区内有哪些监测站点。
• 判断道路、河流、管线与哪些规划范围相交。
• 给采样点追加土地利用类型、保护区名称或人口栅格矢量化后的分区属性。

本文会用一个“点落入行政区面”的例子讲清楚 GeoPandas sjoin怎么用，同时解释 GeoPandas空间索引 为什么能提升速度，以及哪些情况下 sjoin 的结果容易出错。

背景：为什么不能直接用普通表连接

假设你有两个数据：

• points：一批学校点位，字段包括 school_id、name、geometry。
• districts：行政区面，字段包括 district_id、district_name、geometry。

如果学校点表里没有行政区编码，就无法直接用 merge 按字段连接。此时真正的判断条件是：学校点的坐标是否位于某个行政区面内部。

这类问题就应该使用 GeoPandas空间连接：

import geopandas as gpd

points = gpd.read_file("schools.gpkg")
districts = gpd.read_file("districts.gpkg")

result = gpd.sjoin(
    points,
    districts,
    how="left",
    predicate="within"
)

上面的代码表示：以 points 为主表，查找每个学校点位于哪个行政区面内，并把行政区属性追加到点表结果中。

原理：sjoin的核心是空间关系判断

sjoin 的本质是“空间关系连接”。它不是看两个表有没有相同字段，而是判断左表几何与右表几何之间是否满足某种空间谓词。空间谓词可以理解为几何关系条件。

常用的 predicate 参数包括：

predicate	含义	常见用途
within	左侧几何完全位于右侧几何内部	点匹配所属行政区、采样点匹配地类面
contains	左侧几何包含右侧几何	行政区面查找内部点位时使用
intersects	两个几何有任意交集	道路穿过行政区、面与面叠加筛选
touches	两个几何只在边界接触	相邻地块、相邻行政区边界检查
overlaps	两个同维度几何部分重叠但互不完全包含	面与面重叠检查

对 GIS 初学者来说，最容易混淆的是 within 和 contains。判断点落入面时，如果左表是点、右表是面，通常用 predicate="within"。如果左表是面、右表是点，则可以考虑 predicate="contains"。

简单记忆：sjoin(left, right) 中的 predicate 是从左表几何看右表几何。左边的点在右边的面内，就是 within。

步骤：GeoPandas空间连接sjoin怎么用

步骤1：读取空间数据并检查几何字段

首先读取两个空间数据。GeoPandas 支持 Shapefile、GeoPackage、GeoJSON、FileGDB 中的部分数据源，具体取决于本地 GDAL/Fiona 或 Pyogrio 环境。

import geopandas as gpd

points = gpd.read_file("data/schools.gpkg", layer="schools")
districts = gpd.read_file("data/districts.gpkg", layer="districts")

print(points.head())
print(districts.head())

print(points.geometry.name)
print(districts.geometry.name)

需要确认两点：

• 两个对象都是 GeoDataFrame，而不是普通 DataFrame。
• 两者都有有效的 geometry 字段。

步骤2：检查并统一坐标系

GeoPandas sjoin坐标系 是必须检查的重点。空间连接的几何判断是在当前坐标系下完成的。如果两个图层 CRS 不一致，即使图形在地图上看起来应该重叠，计算结果也可能为空或明显错误。

print(points.crs)
print(districts.crs)

如果两个图层坐标系不同，需要把其中一个转换到另一个坐标系：

districts = districts.to_crs(points.crs)

注意，to_crs 是坐标转换，不是简单设置坐标系。如果数据本身缺少 CRS，但你知道它实际使用的是 WGS84，可以先设置：

points = points.set_crs("EPSG:4326")

不要把 set_crs 当成投影转换使用。错误设置 CRS 是 GeoPandas 空间连接结果异常的高频原因。

步骤3：选择正确的how参数

how 控制连接结果保留哪一侧的数据。常用值有 left、right、inner。

how	结果含义	适用场景
left	保留左表所有要素，匹配不到右表时右表字段为空	给点追加所在行政区，且不想丢失任何点
inner	只保留匹配成功的要素	只关心落在研究区范围内的数据
right	保留右表所有要素	较少使用，适合以右表为主的特殊整理需求

在“学校点匹配行政区”的例子中，推荐先使用 how="left"，这样可以检查哪些点没有匹配到行政区。

步骤4：执行sjoin空间连接

下面是一个完整的 GeoPandas空间连接 示例：

import geopandas as gpd

points = gpd.read_file("data/schools.gpkg", layer="schools")
districts = gpd.read_file("data/districts.gpkg", layer="districts")

if points.crs != districts.crs:
    districts = districts.to_crs(points.crs)

joined = gpd.sjoin(
    points,
    districts[["district_id", "district_name", "geometry"]],
    how="left",
    predicate="within"
)

print(joined[["school_id", "name", "district_id", "district_name"]].head())

这里特意只保留了右表中的 district_id、district_name 和 geometry 字段。这样做可以减少无关字段进入结果表，避免字段名冲突，也有助于控制内存占用。

步骤5：检查未匹配记录

空间连接完成后，不要马上导出结果。建议先检查未匹配数量：

unmatched = joined[joined["district_id"].isna()]
print("未匹配点数量：", len(unmatched))

如果未匹配点很多，通常需要检查：

• 点和面是否使用相同 CRS。
• 点是否确实落在行政区面范围内。
• 面数据是否存在几何破损或缝隙。
• 是否选错了 predicate，例如应该用 intersects 却用了 within。

步骤6：处理边界点问题

如果点刚好落在行政区边界线上，within 可能不会把它视为位于面内部。此时可以根据业务含义改用 intersects：

joined = gpd.sjoin(
    points,
    districts[["district_id", "district_name", "geometry"]],
    how="left",
    predicate="intersects"
)

但 intersects 也可能带来一个问题：边界点可能同时匹配两个相邻面，导致一条点记录变成多条结果。行政区归属、网格归属这类业务通常需要额外制定边界归属规则。

步骤7：导出结果

检查结果没有明显问题后，可以导出为空间文件或普通表。

joined.to_file("output/schools_with_district.gpkg", layer="schools_joined", driver="GPKG")

joined.drop(columns="geometry").to_csv(
    "output/schools_with_district.csv",
    index=False,
    encoding="utf-8-sig"
)

如果后续还要在 QGIS、ArcGIS Pro 或其他 GIS 软件中查看，推荐导出 GeoPackage。它比 Shapefile 更适合保存较长字段名和中文属性。

常见坑：sjoin结果为空、重复或很慢怎么办

坑1：两个图层CRS不一致

这是 GeoPandas sjoin结果为空 的首要原因。空间连接不会自动帮你判断两个坐标系是否应该转换。执行前必须检查：

print(points.crs)
print(districts.crs)
print(points.total_bounds)
print(districts.total_bounds)

total_bounds 可以快速查看数据范围。如果一个图层范围类似 120, 30，另一个类似 13000000, 3500000，基本可以判断它们不在同一坐标系统中。

坑2：predicate选错

点落入面，通常使用 within。线穿过面、面与面筛选，通常使用 intersects。如果把所有问题都写成 intersects，结果可能过宽；如果把所有问题都写成 within，结果可能漏掉边界和相交对象。

坑3：字段名冲突导致结果难读

如果左右表有相同字段名，GeoPandas 会给字段加后缀。为了减少混乱，建议在连接前先筛选右表字段：

right_cols = ["district_id", "district_name", "geometry"]

joined = gpd.sjoin(
    points,
    districts[right_cols],
    how="left",
    predicate="within"
)

坑4：一个要素匹配到多个对象

空间连接不是一对一连接。一个点可能落在多个重叠面中，一条道路可能穿过多个行政区，一个面也可能与多个面相交。因此 sjoin 后结果行数可能大于左表行数。

如果业务要求每个点只保留一个结果，需要明确规则，例如：

• 优先选择面积更大的面。
• 优先选择行政级别更高的面。
• 按某个业务字段排序后去重。
• 先清理重叠面，再做空间连接。

坑5：几何无效导致空间关系异常

面数据自相交、空几何、破碎几何都可能导致空间关系判断异常。可以先检查无效几何：

print("points空几何：", points.geometry.is_empty.sum())
print("districts空几何：", districts.geometry.is_empty.sum())
print("districts无效几何：", (~districts.geometry.is_valid).sum())

对于面数据，可以尝试修复：

districts["geometry"] = districts.geometry.make_valid()

修复后仍建议在 GIS 软件中抽查图形，因为自动修复可能改变复杂几何的结构。

方法比较：sjoin、overlay、clip和普通merge怎么选

GeoPandas 中有多个看起来相似的操作。选错方法，会导致结果字段、几何或性能不符合预期。

方法	主要依据	是否改变几何	适合问题
merge	字段值相等	不改变几何	按行政区代码、ID、名称连接属性表
sjoin	空间关系	通常保留一侧几何	点落面、线与面相交、按空间关系追加属性
overlay	空间叠加	会生成新的几何切割结果	面叠加分析、相交区域面积计算
clip	裁剪范围	会裁剪几何	按研究区范围裁剪点线面数据

如果你的目标是“把右表属性追加到左表”，优先考虑 sjoin。如果你的目标是“生成相交后的新面并计算面积”，应考虑 overlay。如果你的目标是“保留研究区内的数据范围”，可以考虑 clip。

步骤：空间索引优化技巧

为什么sjoin需要空间索引

如果有 10 万个点和 1 万个面，最笨的方法是每个点都和每个面做一次空间关系判断，这会产生非常多的几何计算。GeoPandas空间索引 的作用是先用几何外包矩形快速筛选候选对象，再对少量候选几何做精确判断。

可以把空间索引理解为 GIS 里的“快速目录”。它不能改变最终的空间关系规则，但可以减少不必要的计算。

技巧1：确认空间索引可用

GeoPandas 会在很多空间操作中使用空间索引。你可以通过 sindex 检查空间索引对象：

print(points.sindex)
print(districts.sindex)

如果环境缺少必要依赖，某些空间索引功能可能不可用。建议使用较新的 GeoPandas、Shapely 及其兼容环境，并尽量通过 conda 或稳定的 Python 环境统一安装。

技巧2：先裁剪研究区，减少参与连接的数据

空间索引能减少候选几何，但如果原始数据范围过大，仍然会浪费读取、索引和匹配时间。可以先用边界框筛选数据：

minx, miny, maxx, maxy = districts.total_bounds

points_sub = points.cx[minx:maxx, miny:maxy]

joined = gpd.sjoin(
    points_sub,
    districts[["district_id", "district_name", "geometry"]],
    how="left",
    predicate="within"
)

这种方法适合全国点数据匹配某个城市或某个项目区的面数据。先缩小范围，再做 GeoPandas sjoin空间索引优化，通常更稳妥。

技巧3：只保留必要字段

很多空间数据带有几十个甚至上百个字段。空间连接时，如果右表字段全部进入结果，会增加内存压力，也会让结果难以检查。

districts_small = districts[["district_id", "district_name", "geometry"]].copy()

joined = gpd.sjoin(
    points,
    districts_small,
    how="left",
    predicate="within"
)

这是简单但有效的优化。尤其在点数量较大时，字段越少，后续导出和检查越方便。

技巧4：避免在地理坐标系下做距离类判断

如果需要用 sjoin_nearest 做最近邻空间连接，或者要使用距离阈值，建议先投影到米制坐标系。经纬度坐标系的单位是度，不适合直接解释为米。

points_proj = points.to_crs("EPSG:3857")
facilities_proj = facilities.to_crs("EPSG:3857")

nearest = gpd.sjoin_nearest(
    points_proj,
    facilities_proj[["facility_id", "facility_name", "geometry"]],
    how="left",
    max_distance=1000,
    distance_col="dist_m"
)

如果项目对距离精度要求较高，应选择适合本区域的投影坐标系，而不是简单套用 Web Mercator。

技巧5：分块处理超大数据

当点数据达到数百万级，单次读入和连接可能超出内存。可以考虑按行政区、网格、文件批次分块处理，再合并结果。

outputs = []

for city_code, points_part in points.groupby("city_code"):
    districts_part = districts[districts["city_code"] == city_code]

    if districts_part.empty:
        continue

    part_joined = gpd.sjoin(
        points_part,
        districts_part[["district_id", "district_name", "geometry"]],
        how="left",
        predicate="within"
    )

    outputs.append(part_joined)

joined_all = gpd.pd.concat(outputs, ignore_index=True)

上面的思路要求点表和面表都有可用于分块的字段。如果没有字段，也可以先构建规则网格或按空间范围切分。

检查清单：运行sjoin前后要确认哪些内容

做 GeoPandas空间连接sjoin 时，建议按下面清单逐项检查：

• 两个数据是否都是 GeoDataFrame。
• 两个数据是否都有有效的 geometry 字段。
• points.crs 和 districts.crs 是否一致。
• 是否使用了正确的 predicate。
• 是否选择了合适的 how。
• 右表字段是否已精简，避免无关字段进入结果。
• 是否检查了空几何和无效几何。
• 连接后行数是否符合预期。
• 未匹配记录数量是否合理。
• 是否存在一个对象匹配多个对象的情况。
• 大数据量时是否先做范围筛选或分块处理。
• 导出前是否抽样查看地图结果。

如果只记住一句话：先统一坐标系，再选对空间谓词，最后检查未匹配和重复匹配。

FAQ：GeoPandas sjoin常见问题

1. GeoPandas sjoin和Spatial Join是同一个概念吗？

是同一类 GIS 操作。ArcGIS Pro 中叫 Spatial Join，QGIS 中也有按空间位置连接属性的工具，GeoPandas 中常用 gpd.sjoin 实现。不同软件参数名称不同，但核心思想都是按空间关系连接属性。

2. 点落在面边界上，within为什么匹配不到？

within 强调左侧几何在右侧几何内部。边界点不一定被视为内部点。如果业务允许边界点归属到相邻面，可以尝试 predicate="intersects"，但要注意它可能产生多个匹配结果。

3. GeoPandas sjoin结果为什么比原始点表行数多？

因为空间连接可能是一对多关系。一个点可能落在多个重叠面内，一条线可能与多个面相交。结果行数变多通常不是程序错误，而是空间关系本身产生了多个匹配。

4. sjoin结果为空应该先检查什么？

先检查 CRS 和数据范围。执行 print(points.crs)、print(districts.crs)、print(points.total_bounds)、print(districts.total_bounds)。如果坐标系或范围明显不一致，需要先正确设置或转换坐标系。

5. GeoPandas空间索引需要手动创建吗？

多数情况下不需要手动创建。GeoPandas 在空间连接等操作中会使用空间索引机制。你可以通过 gdf.sindex 检查索引对象。真正需要优化时，更重要的是减少数据范围、精简字段、修复几何和分块处理。

6. sjoin可以用来计算每个区内有多少个点吗？

可以。先用 sjoin 把点匹配到区，再按区字段分组统计：

joined = gpd.sjoin(
    points,
    districts[["district_id", "district_name", "geometry"]],
    how="inner",
    predicate="within"
)

count_by_district = joined.groupby("district_id").size().reset_index(name="point_count")

如果最终还要得到行政区面图层，可以再把统计表按 district_id 合并回 districts。

7. sjoin和sjoin_nearest有什么区别？

sjoin 判断的是明确空间关系，例如相交、包含、位于内部。sjoin_nearest 判断的是最近对象，适合查找最近医院、最近道路、最近监测站等。做最近邻时应特别注意坐标系单位，距离分析通常需要投影坐标系。

结论：用好sjoin的关键不是代码长，而是检查到位

GeoPandas 的 sjoin 是 Python GIS 中非常实用的空间连接工具。对于点落面、线面相交、面面筛选这类任务，它比手写循环判断更简洁，也更符合 GIS 工作流。

实际项目中，GeoPandas教程：空间连接sjoin怎么用 的重点不只是记住一行代码，而是理解三个核心问题：左表和右表谁在前、空间谓词是否符合业务含义、坐标系是否一致。

如果遇到 GeoPandas sjoin结果为空、重复记录过多或运行很慢，优先按本文的检查清单排查。多数问题都能通过统一 CRS、修复几何、精简字段、使用合适的 GeoPandas空间索引 优化思路来解决。