1. 为什么处理空间数据要先玩明白shapely
1.1 空间数据这件事,绕不开的几何基础
做地理信息相关项目的人,迟早会遇到一套感觉有点像CAD、又有点像数据库的东西:点、线、面、距离、相交、缓冲区。不管你是做POI选址、路网分析、地块规划,还是单纯想判断两个小区围栏有没有重叠,本质上都是在问同一个问题——几何对象之间到底长什么样、位置关系如何。
Python生态里处理空间数据,最底层的那块砖就是shapely。它不做地图渲染,不下载卫星影像,也不管坐标系转换,它只做一件事:在给定坐标系里表示几何对象,并且对这些对象做精确的布尔运算和关系判断。换句话说,把空间分析拆开看,其他库负责“拿到数据”和“画出来”,shapely负责“算出来”。
很多初学者有个误区:跑空间分析第一件事就去学渲染或者调接口。实际上,如果几何层的计算逻辑都不稳,后面画图再好看也是白搭。比如你要计算“某个点离某条道路的最短距离”,shapely直接一行.distance()给结果;你要判断“两个规划地块是否重叠”,一行.intersects()返回布尔值。这些最基础的操作,恰恰是很多上层GIS软件的底层内核。
1.2 shapely到底解决什么问题
一句话版本:shapely是Python中用于空间几何对象建模、操作和关系分析的库。它基于GEOS(一个C++写的成熟几何引擎),所以计算精度和稳定性比不少纯Python自写的几何函数靠谱得多。
它解决的问题分成两类:
第一类是“单对象操作”:比如一个多边形的面积、周长、中心点、外接矩形;一条线的长度、起终点;一个点的坐标;一个面被切掉一块后的剩余部分。
第二类是“多对象关系”:两个图形是否相交、是否包含、是否相邻;它们公共部分长什么样;合并之后的轮廓是什么;从一个点往外画一个半径1000米的圆,覆盖了多少个现有设施点。
这两类操作基本覆盖了90%的日常空间分析需求。再往上,例如最近邻分析、空间插值、路径规划,shapely可能不直接提供,但它常常是这些算法处理几何输入时真正干活的工具。所以,不管最后用PostGIS也好、GeoPandas也好,底层处理几何的引擎都是同一个思路。搞清楚shapely,等于搞清楚了空间计算的底层逻辑。
2. 几何对象的从0到1:创建与基础操作
2.1 默认的WKT/WKB坐标系,先别急着说坐标转换
shapely本身不区分“经纬度”还是“米”单位,它只把坐标当成数值来看。所以第一步必须明确:坐标的单位其实是由你自己“赋予”的。如果你输入的是lon/lat,那么distance算出来就是“度”,而不是“米”。这是一个极其常见的坑,后面详聊。
shapely里面最常接触的两种几何文本表示是WKT和WKB:
- WKT是文本,例如
POINT (116.3 39.9),人类可读。 - WKB是二进制,适合存储在数据库里,PostGIS就常用WKB传输几何对象。
创建几何对象最常用的方式是直接从WKT文本解析,或者直接从坐标元组构造。
2.2 点、线、面、集合:怎么构造和验证
安装就不用多说了,pip install shapely。要注意,shapely 2.x在API上有一些变化,但核心构造函数保持稳定。下面直接看例子:
from shapely.geometry import Point, LineString, Polygon from shapely import wkt # 点 p1 = Point(0, 0) p2 = Point(3, 4) # 线 line = LineString([(0, 0), (1, 1), (2, 0)]) # 面 polygon = Polygon([(0, 0), (4, 0), (4, 4), (0, 4)])这里有个需要注意的地方:面对象的构造要求顶点按顺序围成一个闭合区域,但不一定需要把第一个点重复一遍,shapely会自动认为首尾相连。如果顶点顺序形成“自相交”的图形,就会报错,比如:
bad_polygon = Polygon([(0, 0), (4, 4), (4, 0), (0, 4)])这种交叉“蝴蝶形”会被判定为无效几何,很多后续计算会失败。所以创建完毕一定要做验证:
print(polygon.is_valid) # True print(bad_polygon.is_valid) # False在我实际处理地块数据的时候,经常遇到从第三方数据源读出来的多边形是无效的,常见原因是重复点、小缺口或者顶点顺序错乱。一定要养成先is_valid再往下算的习惯,不然某个多边形面积算出来是负数,你不会第一时间想到是几何对象坏了。
2.3 基础属性:面积、长度、边界、中心点
拿到了几何对象,最基础的一批属性在shapely里都很直白:
line.length # 线长度 polygon.area # 面面积 polygon.length # 面边界长度,即周长 polygon.bounds # 外接矩形范围 (minx, miny, maxx, maxy) polygon.centroid # 几何中心点 polygon.representative_point() # 多边形内部一个保证落在内部的点这里必须重点想清楚一件事:area、length这类数值的单位和你输入坐标的单位一模一样。如果你的坐标是经纬度,那么polygon.area算出来是“平方度”,这个数值和实际平方公里之间没有简单的固定换算关系。真要算真实距离和面积,请务必把数据投影到合适的平面坐标系下,比如用 UTM 区域投影。
另外一个实用场景:centroid是几何重心,但某些凹多边形的重心可能落在多边形外部。这时候如果需要“一个肯定在多边形内部的点”,比如给某地块圆形范围内选投放点,representative_point()更可靠,它保证返回点一定落在几何内部。
3. 空间关系判断:两个几何对象之间到底什么关系
3.1 相交、包含、接触:空间谓词怎么用
在空间分析里,最常问的问题就是两个对象有没有交集。比如新建的规划地块和生态红线是否冲突,巡检路线是否穿越了某个禁入区。shapely给了一个个直观的谓词方法:
from shapely.geometry import Point, Polygon circle = Point(0, 0).buffer(5) inner = Point(1, 1).buffer(1) outer = Point(10, 10).buffer(1) print(circle.intersects(inner)) # True print(circle.intersects(outer)) # False print(circle.contains(inner)) # True print(circle.touches(outer)) # False(因为没有任何接触)contains的意思是一个几何对象完全包含另一个,边界上的点可以包含在内;covers类似但限制更宽松一点。touches判断两个对象是否只是在边界上接触而没有内部重叠。这几个词含义微妙,做合规性检查时一定要选准。
我在一个模拟项目里处理过商业选址的避让规则:要求距离小学边界200米范围内不能开设某些类型的店铺。这时候不是简单判断一个点是否在缓冲区内,还可能涉及“部分重叠”的情况。如果用intersects,只要店铺面和学校缓冲区沾边就算命中;如果用within,则要求店铺面完全落在缓冲区内。需求不同,选的方法完全不同,选错了会直接得出相反的结论。
3.2 从DE-9IM看空间关系的本质
有时候你会发现,光靠几个谓词不够用。比如,两个面共享一条边,却又不是完全包含,要如何精确描述?这个问题的本质,在计算几何里有一套标准模型:DE-9IM(Dimensionally Extended Nine-Intersection Model)。
它把每个几何对象分成内部、边界、外部三部分,然后比较两个对象这三部分之间的交集情况,得到一张3×3的矩阵。shapely里可以直接输出这个矩阵:
from shapely.geometry import LineString, Point line = LineString([(0, 0), (2, 2)]) pt = Point(1, 1) print(line.relate(pt)) # 输出一个9位字符串,如 "0F0FFF000"这类字符串也就是所谓的关系模式。比如某个字符串表示“点在线内部”,另一个表示“点在线的端点”。看不懂没关系,shapely已经封装好常用的谓词方法,但理解DE-9IM能帮你在遇到“没有现成谓词”的情况时,自己设计规则。尤其是遇到复杂集合对象、点集与面集之间做拓扑校验时,DE-9IM是最精确的通用语言。
3.3 实际业务里的关系判断场景
给一个真实的简化例子:某项目需要把一堆圆形候选点落在行政区内,并且要求这些点不能落在水域面里。候选点可能几万个,水域面几千个,如果两两判断intersects,复杂度接近几十亿次,直接跑死。
常规做法是先看外接矩形:
points_bounds_mask = [point.within(polygon.buffer(0.001)) for polygon in water_polys]但更合理的思路是用空间索引,比如STRtree。shapely自带一个STRtree,可以快速检索出“可能相交”的几何候选,大幅减少精确判断的次数。这也是后面性能优化里非常关键的一环,建议早些习惯。
4. 空间操作实战:缓冲区、联合、交集、差异
4.1 缓冲区(buffer)为什么是万金油
空间分析里最常用、最容易理解的操作就是缓冲区。给一个点一个半径,得到一个圆;给一条线一个宽度,得到一条长条面;给一个面一个扩展宽度,得到一个更大的面。shapely的buffer方法实现得很直接:
from shapely.geometry import Point, LineString pt = Point(0, 0) buf = pt.buffer(2.0) # buf是一个圆形面,半径约2 line = LineString([(0, 0), (10, 0), (20, 2)]) line_buf = line.buffer(1.0) # 一条两倍宽的长带这里有个注意点:buffer的精度由分辨率QuadSegs控制。比如:
pt.buffer(2.0, quad_segs=4)quad_segs表示每四分之一圆周用多少个线段逼近。数值越大,生成的圆越光滑,但顶点数越多。如果你处理的是全国路网,给上万条线段做缓冲,每段都用64段圆弧逼近,数据量会爆炸。所以要根据用途合理选择:计算面积时要求轮廓尽量精确,画图时可以适当降低精度。
更实用的技巧:buffer也可以接受负值,相当于向内收缩一个多边形。比如计算一个地块的“建筑退线范围”,只要对地块边界向内缓冲-10米,再拿缓冲结果和实际建筑轮廓做交集,就能快速判断是否违建。这个思路在规划审查项目里能省很多手工画平行线的功夫。
4.2 union/intersection/difference:集合操作怎么组合
逻辑上,空间几何可以类比为一组色块,它们之间能做并集、交集、差集。在shapely里就是union、intersection、difference这几个函数。
from shapely.geometry import box a = box(0, 0, 4, 4) b = box(2, 2, 6, 6) print(a.intersection(b)) # 重叠的矩形区域 print(a.union(b)) # L型组合区域 print(a.difference(b)) # a去掉和b重叠的部分注意,2.x版本里这些方法依然可用,但推荐直接使用函数式接口:
from shapely import intersection, union, difference result = intersection(a, b)如果对大量几何对象做批量“溶解”,比如把重叠地块合并成一个整体,更好的做法是用unary_union:
from shapely.ops import unary_union merged = unary_union([a, b, c])unary_union一次性把一堆几何对象融合成一个或几个多边形,自动处理互相重叠、嵌套的情况。这个过程实际中很常见:一个区域被多条道路切碎,你想知道“去掉道路带后的可用地块”,就是先合并道路带,再整体取差。
4.3 简化几何与平滑:处理数据噪声
从测绘或遥感数据里读出来的几何对象往往过于精细,顶点多到令人崩溃,很多顶点实际上只是噪声。比如一段本来应该是直线的小路,由于采集误差出现几十个抖动点,看着就烦。shapely提供一个simplify方法,它基于Douglas-Peucker算法,用容差参数控制简化程度:
from shapely.geometry import LineString line = LineString([(0, 0), (0.1, 0.01), (0.2, -0.01), (0.3, 0), (1, 0)]) simplified = line.simplify(tolerance=0.05)容差是“允许的最大形变距离”,设得越大,结果越粗糙,顶点越少。这个操作在做空间分析之前特别有用,因为很多算法复杂度跟顶点数量直接挂钩。比如给一个复杂多边形合并邻域时,先把顶点数降一个量级,性能能提升好几倍。
但注意:simplify主要用于“允许轻微改变形状”的场景。如果做合规审查,必须保留原始边界精度,就不能随意简化。另外,简化之后可能会产生一些极小缝隙或重叠,建议后续再做一次缓冲区“清洗”或者make_valid。
5. 从单对象到几何集合:MultiPolygon与GeometryCollection
5.1 为什么要有多几何对象
想象一个城市的地块集合,它可能包含很多不相连的多边形小岛,比如一个区被河流分隔成两块。在shapely里,这种“多个多边形组成一个逻辑整体”的对象,就是MultiPolygon。
from shapely.geometry import MultiPolygon poly1 = box(0, 0, 1, 1) poly2 = box(5, 5, 7, 7) mp = MultiPolygon([poly1, poly2]) print(mp.area) # 2.0,面积是两个子的和MultiPolygon特别适合表示“一个行政区划的多块飞地”。同样还有MultiLineString、MultiPoint,以及更通用的GeometryCollection,它是可以包含任意类型混合体的集合。
有了MultiPolygon,很多分析可以直接进行。比如判断某个点是否落在这个逻辑区域内,直接mp.contains(point),会自动对每个子面判断。当然,如果你希望对每个子面单独处理,可以遍历其geoms属性:
for g in mp.geoms: print(g)5.2 拆分与合并:在shapely里怎么玩
实际项目中,经常需要把一个MultiPolygon拆成多个单一多边形。很多人会想到遍历geoms,这在大多数场景下已经够用。但如果一个Polygon自带内环(岛洞),又需要把外部边界和洞拆开,那就稍微复杂一点:
mypoly = Polygon( shell=[(0, 0), (10, 0), (10, 10), (0, 10)], holes=[[(2, 2), (3, 2), (3, 3), (2, 3)]] )外部边界用mypoly.exterior,内部洞用mypoly.interiors。如果需要把带洞的多边形拆成“外圈 + 若干洞”的独立面,可以分别构造Polygon(shell)和Polygon(hole_coords)。
反过来,把零散几何合并成MultiPolygon或者进行整体融合,用unary_union更灵活。合并后结果可能是Polygon也可能是MultiPolygon,取决于各部件是否相连。
要注意的是:不要天真地对两个重叠的普通Polygon直接用MultiPolygon([p1, p2]),因为MultiPolygon的语义是“多个不相重叠的成员”。如果输入重叠,结果会产生无效对象。正确做法是先做union,或者用unary_union去融合。
6. 典型项目场景:从数据到答案的完整链路
6.1 场景一:选址分析中的“周边范围圈”
一个常见的需求是:给定一批候选位置点,要找出它们周围500米内有多少个便利店。初学者可能想循环所有候选点和所有便利店点,然后计算距离。这个思路没错,但数据量一旦上万,直接两层循环就成了悲剧。
正确做法是先给候选点做缓冲区,然后利用空间关系判断或空间索引查询。如果只是“统计数量”,更精确的做法是直接用距离阈值判断,但shapely里更自然的写法是:
from shapely.geometry import Point from shapely.strtree import STRtree stores = [Point(1, 1), Point(5, 5), Point(9, 9)] candidates = [Point(1.2, 1.1), Point(8, 8), Point(0, 0)] tree = STRtree(stores) for c in candidates: buffer_geom = c.buffer(1.0) possible_matches = tree.query(buffer_geom) # 此时 possible_matches 是可能位于缓冲区内的点 # 再做精确过滤 count = sum(1 for store in possible_matches if buffer_geom.contains(store))STRtree是shapely自带的R树空间索引,它其实不直接返回“精确命中”,而是返回“可能命中”的候选,所以记得再过滤一次。这里的关键思路是“先通过索引缩小范围,再精确判断”,几乎适用于任何空间关联计算。
6.2 场景二:路网简化和缓冲区叠加
再举个模拟例子:某项目要给一批高速公路线生成道路两侧各500米的影响范围,再和生态保护区叠加,找出重叠区域。整个过程可以很顺畅地写成:
from shapely.geometry import LineString from shapely.ops import unary_union from shapely import intersection road = LineString([(0, 0), (100, 50), (200, 0)]) eco_zone = unary_union([some_zone1, some_zone2]) road_buf = road.buffer(500) overlap = intersection(road_buf, eco_zone)这里有个性能优化点:如果eco_zone是一个巨大复杂多边形,建议先用simplify降维,或者先按要素的外接矩形做快速过滤,再做精确交集。真实项目里我不止一次遇到“就两三个多边形求相交,但跑了十分钟”的情况,原因就是用了一个几千星点的多边形去和一个几百万点的整体做交集。先条碎、再切割、最后合并,效率能快几个数量级。
6.3 性能问题:循环慢、坐标精度、大数据量怎么优化
空间数据量一大,性能就会变成最头疼的问题。三个常见优化点:
第一,尽量避免在Python层循环里逐个调用shapely方法,尤其是遍历上千个对象时。如果只是“把每个点做缓冲区”,可以用列表推导,但最好考虑用shapely的向量化接口。shapely 2.x中,shapely.buffer可以直接接受数组形式:
import shapely points_arr = shapely.points([(0, 0), (1, 1), (2, 2)]) buffers = shapely.buffer(points_arr, 1.0)底层做了并行化和内存优化,比循环快不少。
第二,统一精度。shapely默认浮点双精度,这对很多业务已经够了。如果做大数据量分析,可以把坐标通过round或四舍五入降到合适的小数位,减少后续计算复杂度,但这种方法需谨慎使用,因为它可能破坏拓扑关系。
第三,使用STRtree做空间连接。很多“统计一个点附近有多少个面”的问题,本质是空间连接。shapely的STRtree.query是单侧查询,如果要做批量连接,可以考虑配合GeoPandas的sjoin,但底层仍然是在用相同的R树逻辑。
7. 常见问题与避坑指南
7.1 明明相交为什么返回False?再看一眼坐标系
这两个几何对象一个用经纬度坐标,另一个用投影米制坐标,互相之间差了十万八千里,计算时intersects自然为False。这属于数据类型混乱,不是shapely的bug。
还有一种情况:两个对象物理上靠得很近,但由于浮点误差没有实际相交。比如一条线和一个面之间的距离是0.000001米,intersects返回False。这时可以用buffer(0.01)给其中一个对象加极小的缓冲,再判断关系。
7.2 面对象的顶点顺序:面积出现负数?道理在哪
shapely构造的普通Polygon默认顶点顺序应该按外边界逆时针,内环顺时针。如果你不按这个顺序给坐标,shapely不会直接报错,但部分操作可能怪异。实际上面积永远是正数,因为area取了绝对值。真正会让你困惑的是WKT输出的时候方向不同,以及某些标准格式对方向有严格校验。所以建议统一规范:外部边界逆时针,内部洞顺时针。
如果需要修正,可以使用shapely.geometry.polygon.orient:
from shapely.geometry import Polygon from shapely.geometry.polygon import orient p = Polygon([(0, 0), (4, 0), (4, 4), (0, 4)]) oriented_p = orient(p, sign=1.0)signed_area在简单多边形里可以判断方向。尤其处理从第三方接口拉来的GeoJSON时,很多人会忽略这个,等到做数据入库时被数据库告知“坐标顺序非法”才回头改,印象会很深。
7.3 数据精度与舍入问题:能直接用浮点数比较吗
shapely基于GEOS,用的是双精度浮点。实际中坐标小数位过少会导致两个本该重叠的面边界差了一点点,进而difference之后留下一条极细的残片。这类残片非常典型,比如两个相邻地块各自由不同数据源绘制,边界线在数值上差了0.0000000001,导致合并后出现难以察觉的缝隙。
解决方案有两个方向:一是用snap方法把两个几何对齐,即给定一个容差,将其中一个几何的顶点吸附到另一个几何的边上:
from shapely.ops import snap snapped = snap(polygon_a, polygon_b, tolerance=0.001)二是对差集结果做面积过滤,把微小的残差多边形直接丢弃。实际项目中我经常先用buffer(0)来“清洗”无效几何,这个操作可以解决很多自相交、极小缝隙问题,算是一个不起眼但极为实用的技巧。
7.4 小技巧:用buffer(0)修复无效几何
buffer(0)是一个神奇操作。当一个多边形因自相交或顶点顺序混乱被判定为无效时,buffer(0)常常能修复它:
faulty_polygon = Polygon([(0, 0), (4, 4), (4, 0), (0, 4)]) fixed = faulty_polygon.buffer(0) print(fixed.is_valid)原理是给几何加一个宽度为零的缓冲,等于对所有边界做了一次拓扑重建,自动处理重叠和交叉。但它也可能产生意料之外的效果,比如把一个有效多边形拆分成多个部分,建议使用后检查一下几何类型和面积是否符合预期。
8. 写在最后:我的实际操作感受
用shapely这些年,我最大的体会是:它的API看起来简单,但真正难的是“把你的业务问题翻译成几何关系”。比如“一个POI是否在营业范围内”,翻译成contains还是covers,结果可能有微妙差异;“周边500米内”翻译成缓冲区还是直接用距离阈值,也会影响速度和精度。
如果非要说一条最想分享的经验,就是别急着用高级框架,先把单个几何对象的创建、验证、属性计算练熟。空间数据的世界里,90%的奇怪结果都源于几何对象本身无效、坐标系不统一、单位含义不明确这三类问题。把这三关守住了,后面无论用shapely直接计算,还是接PostGIS、GeoPandas,都会顺很多。
最后一个实际建议:写空间分析代码的时候,多打印几条中间结果进行可视化和人工抽查。空间分析和普通数值计算不同,一个错误的几何往往不会抛出异常,而是静悄悄地给出一个看似合理、实则完全错误的结果。每处理完一个环节,都顺手plot一下或者打印一下wkt,看清楚再继续走,能帮你省下大把排查时间。