如何使用Shapely获取MultiPolygon内每个点的列表

10

我有以下的 MultiPolygon

MULTIPOLYGON (
(
(10.8849956 49.8901705, 10.8849507 49.8902499, 10.884969 49.8902588, 10.8851033 49.8903298, 10.8851183 49.8903132, 10.88512882654868 49.8903054, 10.8851246 49.8903054, 10.8851246 49.8902754, 10.8851546 49.8902754, 10.8851546 49.89028643275958, 10.8853289 49.8901612, 10.885421 49.8901035, 10.8854414638889 49.8900896, 10.8854205 49.8900896, 10.8854205 49.8900596, 10.8854505 49.8900596, 10.8854505 49.89008346226415, 10.885527 49.8900315, 10.885519 49.8899952, 10.8854851 49.8899903, 10.8853164 49.8899957, 10.8852419 49.8899981, 10.8851711 49.8899919, 10.8851165 49.8899814, 10.8850728 49.8899652, 10.8850692 49.8899713, 10.8849925 49.8900275, 10.8850251 49.890083, 10.8850275 49.8901159, 10.8850185 49.8901733, 10.8849956 49.8901705),
(10.8852028 49.8901715, 10.8852328 49.8901715, 10.8852328 49.8902015, 10.8852028 49.8902015, 10.8852028 49.8901715),
(10.8852889 49.8900884, 10.8853146 49.8900884, 10.8853146 49.8901184, 10.8853078 49.8901184, 10.8853078 49.8901337, 10.8852808 49.8901337, 10.8852808 49.8901463, 10.8852508 49.8901463, 10.8852508 49.8901346, 10.8852239 49.8901346, 10.8852239 49.8901046, 10.8852305 49.8901046, 10.8852305 49.8900815, 10.8852589 49.8900815, 10.8852589 49.8900812, 10.8852889 49.8900812, 10.8852889 49.8900884),
(10.8851133 49.890201, 10.8851433 49.890201, 10.8851433 49.890231, 10.8851133 49.890231, 10.8851133 49.890201),
(10.8849849 49.8902202, 10.8850149 49.8902202, 10.8850149 49.8902502, 10.8849849 49.8902502, 10.8849849 49.8902202)
),

(
(10.8852605 49.8901112, 10.8852605 49.8901115, 10.8852539 49.8901115, 10.8852539 49.8901163, 10.8852778 49.8901163, 10.8852778 49.8901112, 10.8852605 49.8901112)
)
)

我该如何获得一个包含所有唯一点的平面列表?列表项不必是Shapely Points,而可以是元组。我不太理解如何迭代这个结构。


2个回答

18

遗憾的是,Shapely并不提供从MultiPolygon对象立即提取所有点功能。相反,你需要首先迭代MultiPolygon中的每个多边形,然后提取每个的各个点

有不同的方法来解决这个问题。例如,如果你知道你的多边形没有孔,你可以简单地执行以下操作:

points = []
for polygon in multipolygon:
    points.extend(polygon.exterior.coords[:-1])

请注意[:-1],它可以避免复制第一个顶点。如果您希望拥有更干净的语法并且不介意每个多边形有一个重复点,则可以将其删除。
这也可以使用具有两个循环的列表推导式在一行中编写:

points = [point for polygon in multipolygon for point in polygon.exterior.coords[:-1]]

或者使用itertools.chain.from_iterable帮助:

from itertools import chain

points = list(chain.from_iterable(polygon.exterior.coords[:-1] for polygon in multipolygon))

通常情况下,如果多边形内部包含空洞,我们可以编写以下函数从内部环中提取坐标:

def to_coords(multipolygon):
    for polygon in multipolygon:
        yield from polygon.exterior.coords[:-1]
        yield from chain.from_iterable(interior.coords[:-1] for interior in polygon.interiors)

使用示例:

mp = MultiPolygon([Polygon([(0, 0), (1, 0), (1, 1), (0, 1)]),
                   Polygon([(2, 0), (3, 0), (3, 1), (2, 1)], 
                           holes=[[(2.25, 0.25), (2.75, 0.25), (2.75, 0.75), (2.25, 0.75)]])])

在此输入图像描述

points = list(to_coords(mp))
# [(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0),
#  (2.0, 0.0), (3.0, 0.0), (3.0, 1.0), (2.0, 1.0), 
#  (2.25, 0.25), (2.75, 0.25), (2.75, 0.75), (2.25, 0.75)]

甚至可以进一步将其推广到任何输入几何体(Python ≥3.7):

from functools import singledispatch
from itertools import chain
from typing import (List, 
                    Tuple,
                    TypeVar)

from shapely.geometry import (GeometryCollection,
                              LinearRing,
                              LineString,
                              Point,
                              Polygon)
from shapely.geometry.base import (BaseGeometry,
                                   BaseMultipartGeometry)

Geometry = TypeVar('Geometry', bound=BaseGeometry)


@singledispatch
def to_coords(geometry: Geometry) -> List[Tuple[float, float]]:
    """Returns a list of unique vertices of a given geometry object."""
    raise NotImplementedError(f"Unsupported Geometry {type(geometry)}")


@to_coords.register
def _(geometry: Point):
    return [(geometry.x, geometry.y)]


@to_coords.register
def _(geometry: LineString):
    return list(geometry.coords)


@to_coords.register
def _(geometry: LinearRing):
    return list(geometry.coords[:-1])


@to_coords.register
def _(geometry: BaseMultipartGeometry):
    return list(set(chain.from_iterable(map(to_coords, geometry))))


@to_coords.register
def _(geometry: Polygon):
    return to_coords(GeometryCollection([geometry.exterior, *geometry.interiors]))

使用示例:

from shapely.geometry import (MultiLineString,
                              MultiPoint,
                              MultiPolygon)

geometry_objects = [Point(0, 0),
                    LineString([(0, 0), (1, 1)]),
                    LinearRing([(0, 0), (1, 0), (1, 1)]),
                    Polygon([(0, 0), (1, 0), (1, 1), (0, 1)], 
                            holes=[[(0.25, 0.25), (0.75, 0.25), (0.75, 0.75), (0.25, 0.75)]]),
                    MultiPoint([(0, 0), (1, 1)]),
                    MultiLineString([LineString([(0, 0), (1, 1)]), LineString([(2, 0), (3, 1)])]),
                    MultiPolygon([Polygon([(0, 0), (1, 0), (1, 1), (0, 1)]),
                                  Polygon([(2, 0), (3, 0), (3, 1), (2, 1)], 
                                          holes=[[(2.25, 0.25), (2.75, 0.25), (2.75, 0.75), (2.25, 0.75)]])]),
                    GeometryCollection([Point(0, 0), LineString([(0, 0), (1, 1)])])]

for geometry in geometry_objects:
    print(f"For {geometry.wkt}\nwe got:\n"
          f"{to_coords(geometry)}\n")

输出:

For POINT (0 0)
we got:
[(0.0, 0.0)]

For LINESTRING (0 0, 1 1)
we got:
[(0.0, 0.0), (1.0, 1.0)]

For LINEARRING (0 0, 1 0, 1 1, 0 0)
we got:
[(0.0, 0.0), (1.0, 0.0), (1.0, 1.0)]

For POLYGON ((0 0, 1 0, 1 1, 0 1, 0 0), (0.25 0.25, 0.75 0.25, 0.75 0.75, 0.25 0.75, 0.25 0.25))
we got:
[(0.0, 1.0), (0.0, 0.0), (0.25, 0.25), (0.75, 0.25), (0.75, 0.75), (0.25, 0.75), (1.0, 0.0), (1.0, 1.0)]

For MULTIPOINT (0 0, 1 1)
we got:
[(0.0, 0.0), (1.0, 1.0)]

For MULTILINESTRING ((0 0, 1 1), (2 0, 3 1))
we got:
[(2.0, 0.0), (0.0, 0.0), (3.0, 1.0), (1.0, 1.0)]

For MULTIPOLYGON (((0 0, 1 0, 1 1, 0 1, 0 0)), ((2 0, 3 0, 3 1, 2 1, 2 0), (2.25 0.25, 2.75 0.25, 2.75 0.75, 2.25 0.75, 2.25 0.25)))
we got:
[(0.0, 1.0), (0.0, 0.0), (3.0, 0.0), (3.0, 1.0), (2.0, 1.0), (2.0, 0.0), (2.25, 0.25), (2.75, 0.25), (2.75, 0.75), (2.25, 0.75), (1.0, 0.0), (1.0, 1.0)]

For GEOMETRYCOLLECTION (POINT (0 0), LINESTRING (0 0, 1 1))
we got:
[(0.0, 0.0), (1.0, 1.0)]

1

根据从Shapely中提取多边形的点坐标,我自己解决了这个问题。我不确定它是否适用于更深层次的MultiPolygons,但它适用于我的使用情况。

import shapely as sh

def get_coords_from_polygon(shape):
    coords = set()

    if isinstance(shape, sh.geometry.Polygon):
        coords.update(shape.exterior.coords[:-1])
        for linearring in shape.interiors:
            coords.update(linearring.coords[:-1])
    elif isinstance(shape, sh.geometry.MultiPolygon):
        for polygon in shape:
            coords.update(get_coords_from_polygon(polygon))

    return coords

不知道为什么它会交换坐标的顺序,你有任何想法吗? - not2qubit

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接