Shapely

Geometric objects, predicates, and operations

Geometry

Geometry — Shapely 2.0.0 documentation

Geometry types

Shapely_-_SpatialDataModel.png

Point(*args) - A geometry type that represents a single coordinate with x,y and possibly z values.
LineString([coordinates]) - A geometry type composed of one or more line segments.
LinearRing([coordinates]) - A geometry type composed of one or more line segments that forms a closed loop.
Polygon([shell, holes]) - A geometry type representing an area that is enclosed by a linear ring.
MultiPoint([points]) - A collection of one or more Points.
MultiLineString([lines]) - A collection of one or more LineStrings.
MultiPolygon([polygons]) - A collection of one or more Polygons.
GeometryCollection([geoms]) - A collection of one or more geometries that may contain more than one type of geometry.

APIs

영역 계산: polygon.area
길이 계산: polygon.length
Bounding box: polygon.bounds
Polygon coords to list: list(polygon.exterior.coords)

difference

Shapely_-_difference.png

a = Point(1, 1).buffer(1.5)
b = Point(2, 1).buffer(1.5)
a.difference(b)

intersection vs symmetric_difference

Shapely_-_intersection-sym-difference.png

a = Point(1, 1).buffer(1.5)
b = Point(2, 1).buffer(1.5)
a.intersection(b)

a = Point(1, 1).buffer(1.5)
b = Point(2, 1).buffer(1.5)
a.symmetric_difference(b)

union

Shapely_-_union.png

a = Point(1, 1).buffer(1.5)
b = Point(2, 1).buffer(1.5)
a.union(b)

a.union(b).boundary
# ...
a.boundary.union(b.boundary)
# ...

intersection 결과로 GeometryCollection 를 받았을 경우

shapely.GeometryCollection — Shapely 2.0.0 documentation

동일한 결과는 아니지만, convex_hull을 사용하면 Polygon 으로 받을 수 있다.

from shapely.geometry import Polygon

def intersection_polygons(*polygons: Polygon) -> Polygon:
    pivot = polygons[0]
    for polygon in polygons[1:]:
        pivot = pivot.intersection(polygon)

    if pivot.geom_type == "GeometryCollection":
        return pivot.convex_hull
    elif pivot.geom_type == "MultiPolygon":
        return pivot.convex_hull
    elif pivot.geom_type == "Polygon":
        return pivot
    else:
        raise NotImplementedError(f"Unsupported geom type: {pivot.geom_type}")

convex_hull에 대한 자세한 내용은 Convex hull항목 참조.

Access raw datas

Examples

Merge Polygons

from shapely.geometry import Polygon
from shapely.ops import cascaded_union

polygon1 = Polygon([(0, 0), (5, 3), (5, 0)])
polygon2 = Polygon([(0, 0), (3, 10), (3, 0)])

polygons = [polygon1, polygon2]

u = cascaded_union(polygons)

도형과 직선의 교점 구하기

(Python) Shapely를 이용하여 도형과 직선의 교점 구하기

Shapely_-_circle_and_line_segments.png

from shaply.geometry import Point, LineString

line = LineString([(-1., -1.),(1., 1.)])
# (-1,-1)과 (1,1)을 잇는 직선을 뜻하는 line이라는 객체를 생성할 수 있습니다.

circle = Point(0., 0.).buffer(1.).boundary
# Point 함수로 원의 중점을 지정해주고 buffer라는 operation을 이용하여 반지름을 지정해주면 됩니다.
# 저는 속이 꽉 찬원이 아닌, 원의 경계선을 만들어 주고 싶었기 때문에 .boundary 옵션을 붙여주었습니다.

segment = line.difference(circle)[1]
# 위와 같이 .difference 함수를 사용하면 차집합의 원소가 구해지고, segment들이 직선이니까 양 끝점의 좌표를 통해 표시됩니다.
# 인덱스를 통해 각 원소에 접근 가능한데, 0번 인덱스는 위 사진에서 segment1, 1번 인덱스는 segment2, 2번 인덱스는 segment3을 뜻합니다.

print(segment.coords[0], segment.coords[1])
# .coords 옵션을 통해 1번 인덱스에 해당하는 좌표를 뽑아보면
# `(-0.7071067811865475, -0.7071067811865475) (0.7071067811865474, 0.7071067811865474)`
# 결과는 우리가 설정해준대로 unit circle과 y=x 직선의 교점이 출력됨을 알 수 있습니다.

Geometry coordinates

[https://geobgu.xyz/py/shapely.html#geometry-coordinates](https://geobgu.xyz/py/shapely.html#geometry-coordinates)

가장 일반적으로 사용되는 5가지 지오메트리 유형의 좌표에 액세스할 수 있는 방법:

Geometry type	Access method
Point coordinates	`.coords`
LineString coordinates	`.coords`
MultiLineString coordinates	`.geoms[0].coords`
Polygon coordinates	`.exterior.coords`, `.interiors[0].coords`
MultiPolygon coordinates	`.geoms[0].exterior.coords`, `.geoms[0].interiors[0].coords`

Raw array

순수 파이썬 자료구조로 변환하는 방법

# -*- coding: utf-8 -*-

from itertools import chain
from typing import List, Optional

from shapely import (
    GeometryCollection,
    LinearRing,
    LineString,
    MultiLineString,
    MultiPoint,
    MultiPolygon,
    Point,
    Polygon,
)
from shapely.geometry.base import BaseGeometry

NumberT = TypeVar("NumberT", int, float)

PointT = Tuple[NumberT, NumberT]


def raw_points(base: Optional[BaseGeometry]) -> List[PointT]:
    if base is None:
        return list()
    elif base.is_empty:
        return list()
    elif base.geom_type == "GeometryCollection":
        assert isinstance(base, GeometryCollection)
        return list(chain(*[raw_points(geom) for geom in base.geoms]))
    elif base.geom_type == "MultiPolygon":
        assert isinstance(base, MultiPolygon)
        raise NotImplementedError
    elif base.geom_type == "Polygon":
        assert isinstance(base, Polygon)
        raise NotImplementedError
    elif base.geom_type == "MultiLineString":
        assert isinstance(base, MultiLineString)
        return list(chain(*[raw_points(geom) for geom in base.geoms]))
    elif base.geom_type == "LinearRing":
        assert isinstance(base, LinearRing)
        raise NotImplementedError
    elif base.geom_type == "LineString":
        assert isinstance(base, LineString)
        return [(c[0], c[1]) for c in base.coords]
    elif base.geom_type == "MultiPoint":
        assert isinstance(base, MultiPoint)
        return list(chain(*[raw_points(geom) for geom in base.geoms]))
    elif base.geom_type == "Point":
        assert isinstance(base, Point)
        return [(base.x, base.y)]
    else:
        raise TypeError(f"Unsupported geom type: {base.geom_type}")


def filter_polygons(base: Optional[BaseGeometry]) -> List[Polygon]:
    if base is None:
        return list()
    elif base.is_empty:
        return list()
    elif base.geom_type == "GeometryCollection":
        assert isinstance(base, GeometryCollection)
        return list(chain(*[filter_polygons(geom) for geom in base.geoms]))
    elif base.geom_type == "MultiPolygon":
        assert isinstance(base, MultiPolygon)
        return list(chain(*[filter_polygons(geom) for geom in base.geoms]))
    elif base.geom_type == "Polygon":
        assert isinstance(base, Polygon)
        return [base]
    elif base.geom_type == "MultiLineString":
        assert isinstance(base, MultiLineString)
        return list()
    elif base.geom_type == "LinearRing":
        assert isinstance(base, LinearRing)
        return list()
    elif base.geom_type == "LineString":
        assert isinstance(base, LineString)
        return list()
    elif base.geom_type == "MultiPoint":
        assert isinstance(base, MultiPoint)
        return list()
    elif base.geom_type == "Point":
        assert isinstance(base, Point)
        return list()
    else:
        raise TypeError(f"Unsupported geom type: {base.geom_type}")

OpenCV Contour

OpenCV의 Contour 와 상호 변환 방법

from typing import List, Tuple, TypeVar

from numpy import array, int32
from numpy.typing import NDArray
from shapely import LineString, MultiPolygon, Polygon
from shapely.geometry.base import BaseGeometry

NumberT = TypeVar("NumberT", int, float)
PointT = Tuple[NumberT, NumberT]
LineT = Tuple[PointT, PointT]
RectT = Tuple[NumberT, NumberT, NumberT, NumberT]


def cvt_contour2polygon(contour: NDArray) -> Polygon:
    assert len(contour.shape) == 3
    assert contour.shape[0] >= 4
    assert contour.shape[1] == 1
    assert contour.shape[2] == 2
    return Polygon(contour[:, 0, :].tolist())


def cvt_contour2linestring(contour: NDArray) -> LineString:
    assert len(contour.shape) == 3
    assert contour.shape[0] >= 4
    assert contour.shape[1] == 1
    assert contour.shape[2] == 2
    return LineString(contour[:, 0, :].tolist())


def cvt_line2linestring(line: LineT) -> LineString:
    return LineString([*line])


def cvt_polygon2contour(polygon: Polygon, dtype=int32) -> NDArray:
    if polygon.is_empty:
        raise ValueError("Polygon must not be empty")

    coords = polygon.exterior.coords
    shape = len(coords), 1, 2
    return array(coords, dtype=dtype).reshape(shape)


def cvt_multipolygon2contours(multipolygon: MultiPolygon, dtype=int32) -> List[NDArray]:
    if multipolygon.is_empty:
        raise ValueError("Polygon must not be empty")

    return [cvt_polygon2contour(polygon, dtype) for polygon in multipolygon.geoms]


def cvt_geometry2contours(geom: BaseGeometry, dtype=int32) -> List[NDArray]:
    if geom.is_empty:
        raise ValueError("Geometry must not be empty")

    if geom.geom_type == "Polygon":
        assert isinstance(geom, Polygon)
        return [cvt_polygon2contour(geom, dtype)]
    elif geom.geom_type == "MultiPolygon":
        assert isinstance(geom, MultiPolygon)
        return cvt_multipolygon2contours(geom, dtype)
    else:
        raise TypeError(f"Unsupported geom type: {geom.geom_type}")


def cvt_roi2contour(roi: RectT, dtype=int32) -> NDArray:
    x1, y1, x2, y2 = roi
    shell = (x1, y1), (x2, y1), (x2, y2), (x1, y2), (x1, y1)
    raw = [p for xy in shell for p in xy]
    contour = array(raw, dtype=dtype).reshape((5, 1, 2))

    assert len(contour.shape) == 3
    assert contour.shape[0] == len(shell)
    assert contour.shape[1] == 1
    assert contour.shape[2] == 2

    return contour


def cvt_roi2polygon(roi: RectT) -> Polygon:
    x1, y1, x2, y2 = roi
    return Polygon([(x1, y1), (x2, y1), (x2, y2), (x1, y2)])


AreaType = Union[NDArray, RectT]


def cvt_area2polygon(area: AreaType) -> Polygon:
    if isinstance(area, ndarray):
        return cvt_contour2polygon(area)
    elif isinstance(area, (tuple, list)):
        if len(area) == 4:
            return cvt_roi2polygon(area)
    raise TypeError(f"Unsupported area type: {type(area).__name__}")