手把手教你用Python+Shapely解决实际问题：从判断快递配送范围到计算地块重叠面积

Python空间分析实战：用Shapely解决配送范围与地块重叠问题

空间数据分析正在成为各行各业的基础技能。无论是外卖平台的配送范围划定，还是房地产开发商的地块规划，亦或是农业种植区的权属划分，都离不开对地理空间关系的精确计算。本文将带你用Python的Shapely库解决两个典型业务场景：判断地址是否在配送范围内，以及计算相邻地块的重叠面积。

1. 环境准备与基础概念

在开始实战之前，我们需要先搭建好开发环境并理解几个核心概念。Shapely是一个专注于平面几何对象操作的Python库，它基于GEOS（Geometry Engine - Open Source）库构建，提供了点、线、面等几何对象的创建、操作和分析功能。

安装Shapely非常简单，只需执行以下命令：

pip install shapely

对于更复杂的地理空间分析，你可能还需要安装GeoPandas：

pip install geopandas

Shapely中的几个核心几何对象：

• Point：表示二维平面上的一个点，由(x,y)坐标定义
• LineString：由一系列点连接而成的线
• Polygon：由至少三个点定义的闭合多边形区域

这些几何对象都有一系列属性和方法，比如计算面积、判断包含关系、求交集等。下面是一个简单的示例：

from shapely.geometry import Point, Polygon # 创建一个点和一个多边形 point = Point(0.5, 0.5) polygon = Polygon([(0, 0), (0, 1), (1, 1), (1, 0)]) # 判断点是否在多边形内 print(point.within(polygon)) # 输出: True

注意：Shapely处理的是平面几何问题，不考虑地球曲率。对于大范围地理空间分析，需要考虑投影转换。

2. 判断地址是否在配送范围内

外卖和快递服务都需要精确判断客户地址是否在配送范围内。这是一个典型的"点是否在多边形内"的空间关系问题。让我们通过一个完整的例子来实现这个功能。

2.1 准备配送范围数据

假设某外卖平台的配送范围由以下坐标点定义的多边形表示：

delivery_area = Polygon([ (116.300, 39.900), # 北京西单 (116.350, 39.900), # 北京王府井 (116.350, 39.950), # 北京朝阳门 (116.300, 39.950) # 北京金融街 ])

2.2 创建地址点对象

客户地址可以转换为经纬度坐标点。例如，某客户位于(116.320, 39.920)：

customer_location = Point(116.320, 39.920)

2.3 判断包含关系

使用within()方法可以轻松判断点是否在多边形内：

is_in_delivery_area = customer_location.within(delivery_area) print(f"该地址{'在' if is_in_delivery_area else '不在'}配送范围内")

2.4 处理实际业务场景

在实际应用中，我们通常需要处理更复杂的情况：

1. 多个配送区域：可能需要检查点是否在任意一个配送多边形内
2. 不规则形状：配送范围可能不是简单的矩形
3. 排除区域：某些区域虽然在大配送范围内，但可能有特殊限制

下面是一个更完整的实现示例：

def is_in_delivery_zones(point, delivery_zones, exclude_zones=None): """ 判断点是否在配送区域内，同时不在排除区域内 :param point: Point对象，客户位置 :param delivery_zones: 配送区域列表，每个元素是一个Polygon :param exclude_zones: 排除区域列表，默认为None :return: bool """ # 检查是否在任一配送区域内 in_delivery = any(point.within(zone) for zone in delivery_zones) if not in_delivery: return False # 如果在配送区域内，检查是否在排除区域 if exclude_zones: in_exclude = any(point.within(zone) for zone in exclude_zones) return not in_exclude return True

3. 计算地块重叠面积

另一个常见场景是计算两个地块的重叠面积，这在土地规划、农业种植和房地产开发中非常有用。Shapely提供了强大的空间关系分析方法来处理这类问题。

3.1 创建地块多边形

假设我们有两块相邻的土地：

from shapely.geometry import Polygon # 地块A plot_a = Polygon([ (0, 0), (0, 5), (5, 5), (5, 0) ]) # 地块B plot_b = Polygon([ (4, 4), (4, 7), (7, 7), (7, 4) ])

3.2 计算重叠区域

使用intersection()方法可以获取两个多边形的交集：

overlap = plot_a.intersection(plot_b)

3.3 计算重叠面积

交集区域也是一个多边形对象，可以直接计算其面积：

overlap_area = overlap.area print(f"重叠面积为: {overlap_area} 平方单位")

3.4 完整解决方案封装

为了便于复用，我们可以将上述逻辑封装成函数：

def calculate_overlap(plot1, plot2): """ 计算两个地块的重叠面积 :param plot1: 第一个地块的Polygon对象 :param plot2: 第二个地块的Polygon对象 :return: 重叠面积(float)，重叠区域(Polygon) """ if not plot1.intersects(plot2): return 0.0, None overlap_region = plot1.intersection(plot2) return overlap_region.area, overlap_region

这个函数返回两个值：重叠面积和重叠区域对象。这样既可以得到数值结果，也可以获取具体的重叠形状用于可视化。

4. 进阶应用与性能优化

掌握了基础用法后，让我们探讨一些更高级的应用场景和性能优化技巧。

4.1 处理复杂多边形

现实中的地块和配送区域往往不是简单的矩形。Shapely可以处理任意复杂的多边形，包括有"洞"的多边形。例如：

# 带孔洞的多边形 complex_plot = Polygon( [(0, 0), (0, 10), (10, 10), (10, 0)], # 外部边界 holes=[[(2, 2), (2, 8), (8, 8), (8, 2)]] # 内部孔洞 )

4.2 批量处理空间关系

当需要处理大量几何对象时，直接使用循环可能效率不高。结合GeoPandas可以显著提高处理效率：

import geopandas as gpd from shapely.geometry import Point # 创建包含多个点的GeoDataFrame points = gpd.GeoDataFrame(geometry=[ Point(1, 1), Point(2, 3), Point(4, 4) ]) # 创建多边形 poly = Polygon([(0, 0), (0, 5), (5, 5), (5, 0)]) # 批量判断点是否在多边形内 points['within'] = points.within(poly)

4.3 空间索引加速查询

对于大规模数据集，使用空间索引可以极大提高查询效率。Shapely本身不提供空间索引，但可以通过GeoPandas或rtree库实现：

import rtree # 创建空间索引 idx = rtree.index.Index() # 将多边形边界框插入索引 for i, polygon in enumerate(polygons): idx.insert(i, polygon.bounds) # 查询可能与目标多边形相交的多边形 potential_matches = list(idx.intersection(target_polygon.bounds))

4.4 可视化分析结果

将分析结果可视化有助于理解和验证。可以使用matplotlib进行简单的绘图：

import matplotlib.pyplot as plt fig, ax = plt.subplots() # 绘制地块A x, y = plot_a.exterior.xy ax.fill(x, y, alpha=0.5, fc='blue', label='地块A') # 绘制地块B x, y = plot_b.exterior.xy ax.fill(x, y, alpha=0.5, fc='green', label='地块B') # 绘制重叠区域 if overlap: x, y = overlap.exterior.xy ax.fill(x, y, alpha=0.5, fc='red', label='重叠区域') ax.legend() plt.show()

5. 常见问题与解决方案

在实际应用中，你可能会遇到一些典型问题。以下是几个常见问题及其解决方案。

5.1 坐标系统不一致

不同数据源可能使用不同的坐标系统，这会导致分析结果错误。解决方案：

• 统一使用相同的坐标参考系统(CRS)
• 必要时进行坐标转换
• 使用GeoPandas的to_crs()方法进行投影转换

import geopandas as gpd gdf = gpd.read_file('data.shp') gdf = gdf.to_crs('EPSG:4326') # 转换为WGS84坐标系统

5.2 几何对象无效

有时几何对象可能因为各种原因无效（如自相交、不闭合等）。可以使用is_valid属性检查：

if not polygon.is_valid: # 尝试修复无效几何 polygon = polygon.buffer(0)

5.3 精度问题

浮点数计算可能引入精度问题。可以设置精度或使用buffer(0)技巧：

from shapely import set_precision polygon = set_precision(polygon, 0.001) # 设置精度为0.001

5.4 大规模数据处理

处理大量几何对象时，内存可能成为瓶颈。解决方案：

• 使用空间索引减少不必要的计算
• 分批处理数据
• 考虑使用Dask-geopandas等分布式计算工具

import dask_geopandas as dgpd ddf = dgpd.read_parquet('large_dataset.parquet') result = ddf[ddf.within(polygon)].compute()

6. 扩展应用场景

Shapely的应用远不止于配送范围和地块分析。以下是一些其他有用的应用场景。

6.1 路径规划与缓冲区分析

通过创建线对象和缓冲区，可以进行简单的路径分析：

from shapely.geometry import LineString # 创建路径 path = LineString([(0, 0), (1, 1), (2, 1), (3, 2)]) # 创建路径缓冲区 buffer_zone = path.buffer(0.2) # 200米缓冲区 # 判断点是否在路径缓冲区内 check_point = Point(1.1, 1.1) print(check_point.within(buffer_zone)) # 输出: True

6.2 服务范围分析

结合Voronoi图，可以分析服务设施的覆盖范围：

from shapely.ops import voronoi_diagram from shapely.geometry import MultiPoint # 服务设施点 facilities = MultiPoint([(0, 0), (2, 3), (4, 1)]) # 生成Voronoi图 voronoi = voronoi_diagram(facilities) # 获取每个服务设施的服务区域 service_areas = list(voronoi.geoms)

6.3 地理围栏监控

实时监控移动物体是否进入或离开特定区域：

class GeoFenceMonitor: def __init__(self, fence_polygon): self.fence = fence_polygon self.previous_status = None def check_position(self, point): current_status = point.within(self.fence) if self.previous_status is None: self.previous_status = current_status return None if current_status and not self.previous_status: self.previous_status = current_status return "进入围栏" elif not current_status and self.previous_status: self.previous_status = current_status return "离开围栏" return None

6.4 农业地块管理

计算不同作物种植区的边界和重叠：

def calculate_crop_overlap(field1, crop1, field2, crop2): """计算两种作物种植区的重叠情况""" overlap_area, overlap_region = calculate_overlap(field1, field2) if overlap_area > 0: return { 'crop1': crop1, 'crop2': crop2, 'overlap_area': overlap_area, 'overlap_percentage_field1': overlap_area / field1.area, 'overlap_percentage_field2': overlap_area / field2.area } return None