在通信工程和 GIS 应用中一个常见的需求是给定一批经纬度坐标点判断它们落在哪些区域多边形内或者距离某个地物有多远。这就是空间匹配要解决的问题。TableGIS 的匹配工具就是为此设计的。它支持将表格中的经纬度数据与 GeoJSON / Shapefile / KML 格式的图层进行空间匹配自动判断点与面的包含关系、计算点到地物的距离并将匹配结果导出为表格。本文将从原理层面介绍这个工具的核心算法并给出 Python 版本的实现代码。一、工具功能概览匹配工具的工作流程很简单输入一张带经纬度列的表格 一个参考图层文件GeoJSON / SHP / KML 输出原表格数据 图层属性字段 距离信息核心功能点功能说明多格式图层解析支持 GeoJSON、Shapefile、KML 三种主流格式点-面匹配判断经纬度点是否落在多边形内缓冲距离支持设置缓冲区半径米正值为向外扩展负值为向内收缩距离计算未匹配时自动计算到最近地物的距离三种匹配模式单个匹配 / 多值合并 / 多行展开三种匹配模式单个每个点只取第一个匹配结果输出一行多值合并一个点匹配到多个地物时属性值用逗号拼接仍为一行多行展开一个点匹配到 N 个地物时输出 N 行二、核心算法2.1 Haversine 公式球面距离计算这是最基础的地理距离算法。地球是近似球体不能直接用平面几何算距离。Haversine 公式通过两点的经纬度计算它们在球面上的大圆距离。原理a sin²(Δlat/2) cos(lat1) × cos(lat2) × sin²(Δlon/2) c 2 × atan2(√a, √(1−a)) d R × c其中 R 为地球半径约 6371000 米。Python 实现importmath EARTH_RADIUS6371000.0# 地球半径单位米defhaversine_distance(lat1,lon1,lat2,lon2):计算两个经纬度点之间的球面距离米lat1_rmath.radians(lat1)lat2_rmath.radians(lat2)d_latmath.radians(lat2-lat1)d_lonmath.radians(lon2-lon1)a(math.sin(d_lat/2)**2math.cos(lat1_r)*math.cos(lat2_r)*math.sin(d_lon/2)**2)amin(a,1.0)# 防止浮点误差导致 NaNc2*math.atan2(math.sqrt(a),math.sqrt(1-a))returnEARTH_RADIUS*c# 示例北京到上海大约 1068 公里print(haversine_distance(39.9042,116.4074,31.2304,121.4737))# 输出: ~1068000 米这个公式在匹配工具中用于计算点到点地物的距离、缓冲距离判断。2.2 射线法Ray Casting点在多边形内判断判断一个点是否在多边形内部是空间匹配的核心问题。TableGIS 使用的是经典的射线法也叫奇偶规则法。原理从待测点向任意方向通常取水平向右发射一条射线计算这条射线与多边形边界的交点数量交点数为奇数→ 点在多边形内部交点数为偶数→ 点在多边形外部┌──────────┐ │ │ P ───●───● │ 射线向右与边界交 2 次 → 偶数 → 外部 │ │ │ │ └──────┘ ┌──────────┐ │ │ │ P ────●──● 射线向右与边界交 1 次 → 奇数 → 内部 │ │ └──────────┘Python 实现defpoint_in_polygon(px,py,polygon): 射线法判断点 (px, py) 是否在多边形内 polygon: [(x1,y1), (x2,y2), ...] 顶点列表 nlen(polygon)insideFalsejn-1foriinrange(n):xi,yipolygon[i]xj,yjpolygon[j]# 判断射线是否穿过这条边if((yipy)!(yjpy))and\(px(xj-xi)*(py-yi)/(yj-yi)xi):insidenotinside jireturninside# 示例判断点是否在一个矩形内polygon[(0,0),(10,0),(10,10),(0,10)]print(point_in_polygon(5,5,polygon))# True - 在内部print(point_in_polygon(15,5,polygon))# False - 在外部这个算法的时间复杂度是 O(n)n 为多边形顶点数效率很高。2.3 点到多边形的距离当点不在多边形内部时需要计算它到多边形边界的最短距离。思路是遍历多边形的每条边计算点到每条线段注意是线段不是直线的距离取最小值。点到线段的距离defpoint_to_segment_distance(px,py,ax,ay,bx,by):点 (px,py) 到线段 (ax,ay)-(bx,by) 的距离度为单位dxbx-ax dyby-ayifdx0anddy0:# 线段退化为点returnmath.sqrt((px-ax)**2(py-ay)**2)# 计算投影参数 t限制在 [0, 1] 范围内t((px-ax)*dx(py-ay)*dy)/(dx*dxdy*dy)tmax(0.0,min(1.0,t))# 最近点坐标closest_xaxt*dx closest_yayt*dyreturnmath.sqrt((px-closest_x)**2(py-closest_y)**2)defpoint_to_polygon_distance(lon,lat,polygon):计算点到多边形边界的最短距离近似米min_distfloat(inf)nlen(polygon)foriinrange(n):j(i1)%n dpoint_to_segment_distance(lon,lat,polygon[i][0],polygon[i][1],polygon[j][0],polygon[j][1])ifdmin_dist:min_distd# 将度距离近似转换为米lat_radmath.radians(lat)deg_to_meterEARTH_RADIUS*math.pi/180.0*math.cos(lat_rad)returnmin_dist*deg_to_meter2.4 缓冲区匹配匹配工具支持设置缓冲距离这在通信场景中很实用——比如判断一个测点是否在某基站覆盖范围的 500 米以内。缓冲区的实现逻辑正缓冲点在多边形内部或者点到多边形的距离 ≤ 缓冲距离 → 匹配成功负缓冲点不仅要在多边形内部而且到边界的距离必须 ≥ |缓冲距离| → 相当于把多边形向内收缩defmatch_with_buffer(lon,lat,polygon,buffer_distance): 带缓冲区的匹配 buffer_distance 0: 向外扩展 buffer_distance 0: 向内收缩 buffer_distance 0: 严格在多边形内 is_insidepoint_in_polygon(lon,lat,polygon)ifbuffer_distance0:# 正缓冲在内部 或 距离边界在缓冲范围内ifis_inside:returnTrue,0.0distpoint_to_polygon_distance(lon,lat,polygon)ifdistbuffer_distance:returnTrue,distreturnFalse,distelse:# 负缓冲必须在内部且距边界至少 |buffer| 米ifnotis_inside:returnFalse,-1edge_distpoint_to_polygon_distance(lon,lat,polygon)ifedge_distabs(buffer_distance):returnTrue,0.0returnFalse,edge_dist三、图层文件解析匹配工具支持三种图层格式下面简要说明各自的解析思路。3.1 GeoJSONGeoJSON 是基于 JSON 的地理数据格式结构清晰{type:FeatureCollection,features:[{type:Feature,geometry:{type:Polygon,coordinates:[[[116.0,39.0],[117.0,39.0],[117.0,40.0],[116.0,40.0],[116.0,39.0]]]},properties:{name:区域A,type:城区}}]}Python 解析非常简单importjsondefparse_geojson(file_path):withopen(file_path,r,encodingutf-8)asf:datajson.load(f)features[]forfeatindata.get(features,[]):geom_typefeat[geometry][type]propsfeat.get(properties,{})coordsfeat[geometry][coordinates]ifgeom_typePoint:features.append({type:Point,lon:coords[0],lat:coords[1],properties:props})elifgeom_typePolygon:polygon[(c[0],c[1])forcincoords[0]]features.append({type:Polygon,polygon:polygon,properties:props})returnfeatures3.2 KMLKML 是 Google Earth 使用的 XML 格式核心结构是Placemark标签importxml.etree.ElementTreeasETdefparse_kml(file_path):treeET.parse(file_path)roottree.getroot()# KML 有命名空间需要处理ns{kml:http://www.opengis.net/kml/2.2}features[]forpminroot.findall(.//kml:Placemark,ns):props{}namepm.find(kml:name,ns)ifnameisnotNone:props[name]name.text# 解析 Pointpointpm.find(.//kml:Point/kml:coordinates,ns)ifpointisnotNone:partspoint.text.strip().split(,)features.append({type:Point,lon:float(parts[0]),lat:float(parts[1]),properties:props})continue# 解析 Polygonpolygonpm.find(.//kml:Polygon//kml:coordinates,ns)ifpolygonisnotNone:coords[]forcoordinpolygon.text.strip().split():partscoord.split(,)coords.append((float(parts[0]),float(parts[1])))features.append({type:Polygon,polygon:coords,properties:props})returnfeatures3.3 ShapefileShapefile 是 ESRI 的二进制格式由.shp几何、.dbf属性、.shx索引三个文件组成。解析比较复杂通常使用pyshp库importshapefile# pip install pyshpdefparse_shapefile(file_path):sfshapefile.Reader(file_path,encodinggbk)features[]field_names[field[0]forfieldinsf.fields[1:]]# 跳过 DeletionFlagforsrinsf.shapeRecords():propsdict(zip(field_names,sr.record))shapesr.shapeifshape.shapeTypeshapefile.POINT:features.append({type:Point,lon:shape.points[0][0],lat:shape.points[0][1],properties:props})elifshape.shapeTypein(shapefile.POLYGON,shapefile.POLYLINE):# 处理多部件partslist(shape.parts)[len(shape.points)]foriinrange(len(shape.parts)):polygon[shape.points[j]forjinrange(parts[i],parts[i1])]features.append({type:Polygon,polygon:polygon,properties:props})returnfeatures四、完整匹配流程把上面的算法组合起来就是完整的匹配流程defspatial_match(point_data,layer_features,buffer_distance0,match_fieldNone,methodsingle): 空间匹配主流程 point_data: [{lon: 116.4, lat: 39.9, ...}, ...] layer_features: 图层要素列表 buffer_distance: 缓冲距离米 match_field: 要匹配的字段名None 表示全部 method: single / merge / expand results[]forpointinpoint_data:lon,latpoint[lon],point[lat]matches[]forfeatinlayer_features:matchedFalsedist-1iffeat[type]Polygon:ifpoint_in_polygon(lon,lat,feat[polygon]):ifbuffer_distance0:matched,distTrue,0.0else:edge_distpoint_to_polygon_distance(lon,lat,feat[polygon])ifedge_distabs(buffer_distance):matched,distTrue,0.0ifnotmatched:iffeat[type]Point:disthaversine_distance(lat,lon,feat[lat],feat[lon])eliffeat[type]Polygon:distpoint_to_polygon_distance(lon,lat,feat[polygon])if0distbuffer_distance:matchedTrueifmatched:matches.append({properties:feat[properties],distance:dist})# 根据模式组装输出ifmethodsingle:rowdict(point)ifmatches:row.update(matches[0][properties])row[距离(米)]matches[0][distance]else:row[距离(米)]-1results.append(row)elifmethodmerge:rowdict(point)merged{}min_distfloat(inf)forminmatches:ifm[distance]min_dist:min_distm[distance]fork,vinm[properties].items():ifknotinmerged:merged[k]velifv!merged[k]:merged[k],v row.update(merged)row[距离(米)]min_distifmin_dist1e17else-1results.append(row)elifmethodexpand:ifnotmatches:rowdict(point)row[距离(米)]-1results.append(row)else:forminmatches:rowdict(point)row.update(m[properties])row[距离(米)]m[distance]results.append(row)returnresults五、实际应用场景场景 1基站覆盖分析将路测数据经纬度与基站扇区图层匹配判断每个测点落在哪个扇区的覆盖范围内输入路测点表经度、纬度、RSRP 扇区 KML含扇区名、方位角、覆盖距离 匹配结果每个测点 → 所属扇区名 到扇区中心的距离场景 2行政区划归属将一批 GPS 坐标点匹配到行政区划面图层自动标注所属省/市/区输入GPS 轨迹表 行政区划 GeoJSON 匹配结果每个坐标点 → 所属行政区名称场景 3缓冲区搜索查找某坐标点周围 500 米内的 POI兴趣点输入目标点坐标 POI 点图层 缓冲距离500 米 匹配结果500 米内的所有 POI按距离排序六、性能与优化匹配工具面对的数据量通常是几千到几万行图层要素从几十到几千个。朴素的遍历匹配每个点遍历所有要素在万级数据量下可能需要几秒。如果需要优化常见的空间索引方案有方案适用场景说明R-Tree多边形/线要素用 bounding box 构建索引树快速排除不可能相交的要素KD-Tree点要素适合最近邻搜索网格索引均匀分布的要素将空间划分为网格每个格子记录包含的要素当前版本的匹配工具使用朴素遍历对于日常使用已经足够。如果后续需要处理更大数据量可以引入空间索引来加速。七、总结TableGIS 匹配工具的核心就是三个算法的组合Haversine 公式→ 球面上两点之间的距离射线法→ 点是否在多边形内点到线段距离→ 点到多边形边界的最近距离配合缓冲区机制和三种匹配模式可以覆盖大多数点-面空间匹配的场景。工具支持 GeoJSON / Shapefile / KML 三种主流图层格式在 Excel 和 WPS 两个版本中都可以使用。如果你也有类似的空间匹配需求希望本文的算法讲解和 Python 代码能给你一些参考。