在两个包含Shapely对象的numpy数组上应用成对Shapely函数

Question

我有两个长度不同的数组。一个数组里是形状多边形，另一个数组里是形状点。我想对这两个数组里的每一个元素组合，使用一个叫做 a_polygon.contains(a_point) 的函数。

我在看 这篇帖子，因为如果能先建立一个两列的矩阵，把所有可能的组合放在行里，那可能是个不错的中间步骤。但是，'cartersian(arrays)' 函数里的循环在处理大量数据时可能会影响性能。

我尝试过对其中一个数组进行广播，然后应用这个 shapely 函数：

Polygons_array[:,newaxis].contains(Points_array)

但当然，这样做是不行的。我知道最近发布了 geopandas 库，但在我的 Canopy 安装环境中无法使用。

Answer 1

下面的代码展示了如何对两个长度不同的数组中的几何对象应用一个函数。这种方法避免了使用循环。我们需要用到Pandas的apply方法和Numpy的.vectorize以及广播功能。

首先，我们需要导入一些库，并创建以下两个数组：

import numpy as np
import pandas as pd
from shapely.geometry import Polygon, Point

polygons = [[(1,1),(4,3),(4,1),(1,1)],[(2,4),(2,6),(4,6),(4,4),(2,4)],[(8,1),(5,1),(5,4),(8,1)]]
points = [(3,5),(7,3),(7,6),(3,2)]

可以通过以下方式获取一个包含多边形和点的几何对象数组：

geo_polygons = pd.DataFrame({'single_column':polygons}).single_column.apply(lambda x: Polygon(x)).values
geo_points = pd.DataFrame({'single_column':points}).single_column.apply(lambda x: Point(x[0], x[1])).values
# As you might noticed, the arrays have different length.

接下来，我们定义一个要应用于这两个数组的函数，并将其向量化：

def contains(a_polygon, a_point):
    return a_polygon.contains(a_point)
contains_vectorized = np.vectorize(contains)

这样，函数就准备好可以应用到向量中的每个元素上。通过广播点数组，可以处理成对的评估：

contains_vectorized(geo_polygons, geo_points[:,np.newaxis])

这将返回以下数组：

array([[False,  True, False],
   [False, False, False],
   [False, False, False],
   [ True, False, False]], dtype=bool)

这个数组的列对应多边形，行对应点。数组中的布尔值表示，例如，第一个点是否在第二个多边形内。这个结果是正确的。将多边形和点进行映射可以验证这一点：

from descartes import PolygonPatch
import matplotlib.pyplot as plt
fig = plt.figure(1, figsize = [10,10], dpi = 300)
ax = fig.add_subplot(111)
offset_x = lambda xy: (xy[0] + 0.1, xy[1])
offset_y = lambda xy: (xy[0], xy[1] - 0.5)
for i,j in enumerate(geo_polygons):
    ax.add_patch(PolygonPatch(j, alpha=0.5))
    plt.annotate('polygon {}'.format(i + 1), xy= offset_y(tuple(j.centroid.coords[0])))
for i,j in enumerate(geo_points):
    ax.add_patch(PolygonPatch(j.buffer(0.07),fc='orange',ec='black'))
    plt.annotate('point {}'.format(i + 1), xy= offset_x(tuple(j.coords[0])))
ax.set_xlim(0, 9)
ax.set_ylim(0, 7)
ax.set_aspect(1)
plt.show()