从 Voronoi 镶嵌到 Shapely 多边形答案

【问题标题】：From Voronoi tessellation to Shapely polygons从 Voronoi 镶嵌到 Shapely 多边形
【发布时间】：2015-02-17 08:17:39
【问题描述】：

根据一组点，我使用scipy 构建了 Voronoi 镶嵌：

from scipy.spatial import Voronoi
vor = Voronoi(points)

现在我想从 Voronoi 算法创建的区域构建一个Polygon in Shapely。问题是 Polygon 类需要一个逆时针顶点列表。虽然我知道如何order these vertices，但我无法解决问题，因为通常这是我的结果：

（重叠多边形）。这是代码（一个随机示例）：

def order_vertices(l):
    mlat = sum(x[0] for x in l) / len(l)
    mlng = sum(x[1] for x in l) / len(l)

    # https://stackoverflow.com/questions/1709283/how-can-i-sort-a-coordinate-list-for-a-rectangle-counterclockwise
    def algo(x):
        return (math.atan2(x[0] - mlat, x[1] - mlng) + 2 * math.pi) % 2*math.pi

    l.sort(key=algo)
    return l

a = np.asarray(order_vertices([(9.258054711746084, 45.486245994138976),
 (9.239284166975443, 45.46805963143515),
 (9.271640747003861, 45.48987234571072),
 (9.25828782103321, 45.44377372506324),
 (9.253993275176263, 45.44484395950612),
 (9.250114174032936, 45.48417979682819)]))
plt.plot(a[:,0], a[:,1])

我该如何解决这个问题？

【问题讨论】：

这看起来不像顶点的逆时针顺序。你确定你正确地实现了排序算法吗？
@Kevin 我认为我正确地实现了它......我用一个例子更新了这个问题
是生成图表的代码吗？或者它是一个单独的例子？
@Kevin 该代码生成了图表 :)

标签： python gis voronoi shapely

【解决方案1】：

如果您只是在收集多边形，则无需预先订购点来构建它们。

scipy.spatial.Voronoi 对象有一个 ridge_vertices 属性，其中包含构成 Voronoi 脊线的顶点索引。如果索引是-1，那么山脊会趋于无穷大。

首先从一些随机点开始构建 Voronoi 对象。

import numpy as np
from scipy.spatial import Voronoi, voronoi_plot_2d
import shapely.geometry
import shapely.ops

points = np.random.random((10, 2))
vor = Voronoi(points)
voronoi_plot_2d(vor)

您可以使用它来构建 Shapely LineString 对象的集合。

lines = [
    shapely.geometry.LineString(vor.vertices[line])
    for line in vor.ridge_vertices
    if -1 not in line
]

shapely.ops 模块有一个 polygonize，它返回一个用于 Shapely Polygon 对象的生成器。

for poly in shapely.ops.polygonize(lines):
    #do something with each polygon

或者，如果您想要由 Voronoi 镶嵌所包围的区域形成单个多边形，您可以使用 Shapely unary_union 方法：

shapely.ops.unary_union(list(shapely.ops.polygonize(lines)))

【讨论】：

如何保存带有多边形的 shapefile？
@ValerioD.Ciotti 在多边形上使用shapely.geometry.mapping 将它们转换为geojson 对象，然后使用fiona 之类的库将geojson 写入新的shapefile。如果您遇到问题，请在 http://gis.stackexchange.com/ 上提问。
@ValerioD.Ciotti 或者您可以从多边形制作geopandas.GeoSeries 并将其直接保存到 shapefile 中。
你也可以直接使用GeoPandas已经依赖的fiona库。
查看this answer 了解处理无限 Voronoi 区域的方法。

【解决方案2】：

正如其他人所说，这是因为您必须根据索引正确地从结果点重建多边形。尽管您有解决方案，但我想我应该提到还有另一个 pypi 支持的镶嵌包，称为 Pytess（免责声明：我是包维护者），其中 voronoi 函数返回完全为您构建的 voronoi 多边形。

【讨论】：

【解决方案3】：

该库可以生成坐标的有序列表，您只需要使用提供的索引列表：

import numpy as np
from scipy.spatial import Voronoi

...

ids = np.array(my_points_list)
vor = Voronoi(points)
polygons = {}
for id, region_index in enumerate(vor.point_region):
    points = []
    for vertex_index in vor.regions[region_index]:
        if vertex_index != -1:  # the library uses this for infinity
            points.append(list(vor.vertices[vertex_index]))
    points.append(points[0])
    polygons[id]=points

polygons 字典中的每个多边形都可以导出到 geojson 或带入 shapely，我能够在 QGIS 中正确渲染它们

【讨论】：

【解决方案4】：

您实现的函数（order_vertices()）在您的情况下无法使用，因为它只需要一个已经排序的坐标序列，该序列构建一个矩形，并反转多边形的方向（并且可能仅适用于矩形。 ..)。但是你有一个无序的坐标序列

一般来说，你不能从任意序列的无序顶点构建多边形，因为凹多边形没有唯一的解决方案，如下例所示：https://stackoverflow.com/a/7408711/4313133

但是，如果您确定您的多边形始终是凸的，您可以使用以下代码构建一个凸包：https://stackoverflow.com/a/15945375/4313133（现已测试，对我有用）

也许您也可以使用 scipy 构建凸包，但我尝试对其进行测试：scipy.spatial.ConvexHull

【讨论】：