Pandas 如何使用Pandas和geopandas来判断点是否在多边形内

在本文中，我们将介绍如何使用Pandas和geopandas来判断点是否在多边形内。这是一个常见的GIS问题，例如我们需要判断一个地址是否在某个特定的区域内，或者我们需要判断一个犯罪案件是否发生在某个特定的行政区域内。

数据准备

首先，我们需要准备两个数据集：点数据和多边形数据。假设我们有一个包含10000个点的数据集，我们可以用Pandas来读取和处理它。

import pandas as pd

# 读取点数据
points = pd.read_csv('points.csv')

# 查看数据
points.head()

点数据应该至少包含两列：经度和纬度。如果不包含这些列可以通过下面的代码新增加这些列。

# 将经度和纬度列命名为 'x' 和 'y'
points = points.rename(columns={'longitude':'x', 'latitude':'y'})

# 查看新增列后的数据
points.head()

接下来，我们需要读取包含所有多边形的数据集。最常见的数据格式是ESRI Shapefile，但geopandas也支持其他格式，例如GeoJSON和KML。

import geopandas as gpd

# 读取多边形数据
polygons = gpd.read_file('polygons.shp')

# 查看多边形数据
polygons.head()

多边形数据应该包含至少两列：几何列和任何其他需要使用的列。在ESRI Shapefile中，几何列通常被命名为“geometry”，它包含多边形的几何信息。

Geopandas point in polygon

现在我们已经准备好了点和多边形数据，我们可以使用geopandas中的sjoin函数来判断每个点是否在多边形内。sjoin函数可以完美地处理空间查询为主的数据分析问题，它可以识别两个GeoDataFrame中基于点、线和面的空间关系，例如点是否在多边形内。

# 使用 sjoin 函数查询点是否在多边形内
result = gpd.sjoin(points, polygons, op='within')

# 查看结果
result.head()

此时，result数据集将包含两个数据集的所有列，加上一个名为index_right的新列，它标识了每个点所在的多边形。如果点不在任何多边形内，则值为NaN。

现在我们可以使用value_counts()方法统计每个多边形内的点数。

# 统计点落入每个多边形内的数量
count = result['index_right'].value_counts()

# 查看结果
print(count)

性能优化

如果我们有很多点和/或很多多边形，上面的代码可能会花费很长的时间来执行。幸运的是，geopandas提供了一种简单的方法来优化这个过程：将多边形数据集转换为一个凸多边形集合。

一组凸多边形具有以下优点：

减少计算量，特别是当数据集包含很多多边形时。
统计计算更容易，比如可以使用快速点位于多边形凸壳集合内的算法来快速识别匹配的多边形。

# 将多边形转换为凸多边形
convex_polygons = polygons.convex_hull

# 重新命名几何列为 'geometry'
convex_polygons = convex_polygons.rename(columns={'geometry':'geometry'})

接下来，我们用新的多边形数据集再次运行sjoin函数，并计算每个多边形内的点数：

# 使用 sjoin 函数查询点是否在多边形内result = gpd.sjoin(points, convex_polygons, op='within')

# 统计点落入每个多边形内的数量
count = result['index_right'].value_counts()

# 查看结果
print(count)

示例

我们可以使用一个具体的示例来演示如何使用Pandas和geopandas判断点是否在多边形内。这里，我们将使用一个包含美国各州的多边形数据集，然后判断50个随机经纬度点是否在美国某个州内。

import random

# 读取美国各州的多边形数据集
us_states = gpd.read_file('us_states.shp')

# 生成50个随机经纬度点
random.seed(123)
points = pd.DataFrame({
    'x': [random.uniform(-130, -65) for i in range(50)],
    'y': [random.uniform(25, 50) for i in range(50)]
})

# 将点转换为GeoDataFrame
geometry = gpd.points_from_xy(points['x'], points['y'])
points = gpd.GeoDataFrame(points, geometry=geometry)

# 判断每个点是否在美国某个州内
result = gpd.sjoin(points, us_states, op='within')

# 查看结果
print(result.head())