【541】shapely 相关功能

参考：Python | 任意多边形间的重叠面积计算

参考：Shapely用户手册

参考：Python地理数据处理库shapely支持函数总结

---

1. 构建集合图形以及获取集合图形点信息

　　得到点集合，适合后面的处理。

poly_1 = Polygon([(0, 0), (0, 2), (2, 2), (2, 0)])

# 获取多边形外边坐标信息，最后是闭合的
# 类似列表，不过里面是 tuple，按照 xy 顺序显示
poly_1.exterior.coords[:]

# 输出结果如下：
# [(0.0, 0.0), (0.0, 2.0), (2.0, 2.0), (2.0, 0.0), (0.0, 0.0)]

pts = poly_1.exterior.coords
type(pts) 
# 输出结果
# shapely.coords.CoordinateSequence

2. 多边形显示

　　最后一句保证显示按照比例

from matplotlib import pyplot
from shapely.geometry import Point
from descartes import PolygonPatch
import numpy as np 

fig = pyplot.figure(1, dpi=90)
ax = fig.add_subplot(121)

patch1 = PolygonPatch(poly_1, alpha=0.5, zorder=1)
ax.add_patch(patch1)

patch2 = PolygonPatch(poly_2, alpha=0.5, zorder=1)
ax.add_patch(patch2)

patchc = PolygonPatch(poly_1.intersection(poly_2) , alpha=0.5, zorder=2)
ax.add_patch(patchc)

pyplot.xlim((-1, 4))
pyplot.ylim((-1, 3))

ax.set_aspect('equal', adjustable='box') 

from matplotlib import pyplot
from shapely.geometry import Point
from descartes import PolygonPatch
import numpy as np 

fig = pyplot.figure(1, dpi=90)
ax = fig.add_subplot(121)

ax.add_patch(PolygonPatch(poly_1.buffer(0.5), alpha=0.5, zorder=1))

patch1 = PolygonPatch(poly_1, alpha=0.5, zorder=1)
ax.add_patch(patch1)

pyplot.xlim((-1, 3))
pyplot.ylim((-1, 3))

ax.set_aspect('equal', adjustable='box')

from matplotlib import pyplot
from shapely.geometry import Point
from descartes import PolygonPatch
import numpy as np 

fig = pyplot.figure(1, dpi=90)

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

# 1
ax = fig.add_subplot(121)

patch1 = PolygonPatch(a, alpha=0.5, zorder=1)
ax.add_patch(patch1)
patch2 = PolygonPatch(b, alpha=0.5, zorder=1)
ax.add_patch(patch2)
c = a.union(b)
patchc = PolygonPatch(c, alpha=0.5, zorder=2)
ax.add_patch(patchc)

#pyplot.xlim((-1, 4))
#pyplot.ylim((-1, 3))

pyplot.xticks(np.arange(-1, 5, 1))
pyplot.yticks(np.arange(-1, 4, 1)) 

ax.set_aspect('equal', adjustable='box')

3. 多边形分割

　　可以将多边形通过折线来分割

　　参考：Cut a polygon with two lines in Shapely

　　参考：shapely官方文档——Splitting

　　切割后得到一个多边形集合，通过遍历可以获取每一个 geometry 的具体信息。

from shapely.geometry import Polygon, Point, LineString 
poly_1 = Polygon([(0,0), (0,2), (2,2), (2,0)]) 
line_1 = LineString([(-1,0.5), (-1,1.5), (3,1.5), (3,0.5), (-1,0.5)])

import shapely 
# 通过线可以将多边形进行切割
result = shapely.ops.split(poly_1, line_1) 
result.wkt 

# output:
# 'GEOMETRYCOLLECTION (POLYGON ((0 0, 0 0.5, 2 0.5, 2 0, 0 0)), POLYGON ((0 0.5, 0 1.5, 2 1.5, 2 0.5, 0 0.5)), POLYGON ((0 1.5, 0 2, 2 2, 2 1.5, 0 1.5)))'

for g in result: 
    print(g.exterior.coords[:]) 

# output:
# [(0.0, 0.0), (0.0, 0.5), (2.0, 0.5), (2.0, 0.0), (0.0, 0.0)]
# [(0.0, 0.5), (0.0, 1.5), (2.0, 1.5), (2.0, 0.5), (0.0, 0.5)]
# [(0.0, 1.5), (0.0, 2.0), (2.0, 2.0), (2.0, 1.5), (0.0, 1.5)]

 

4. Polygon 被 MultiLineString 切割

　　典型的应用场景就是路网切割，路网是一段段的折线（LineString）组成的，因此需要切割指定的多边形，从而生成多个多边形，不能通过上面的方法实现。

　　参考：Divide a polygon into multiple small polygons using a MultiLineString

　　具体实现的思路是通过密集的 LineString 集合然后分裂转成 Polygon 的思路

from shapely.geometry import Polygon, LineString
from shapely.ops import linemerge, unary_union, polygonize

poly = Polygon([(0,0), (4,0), (4,4), (0,4)])
lines = [LineString([(0,1), (1,1), (2,1), (5,1)]), LineString([(1,1), (2,1), (3,3), (5,3)])]

lines.append(poly.boundary)
lines = unary_union(lines)
lines = linemerge(lines)
polygons = polygonize(lines)

　　路网数据  

　　

　　切割后的数据

　　

　　根据路网，切割望京区域（选取边界区域的线路，然后按照上面的操作获取多个区域，选取面积最大的）　　

　　

　　

 

　　然后根据包含关系，获取望京内部的AOI信息，需要做个buffer，不然有些边界的区域会出错

　　

5. merge 多个多边形

　　使用下面的方法

• 显示构建 list
• 然后再用 cascaded_union

def get_merge_poly(poly_dict):
    """
    获取给定的所有 polys 的合并结果，避免面积重复计算
    输入：字典，名称+数组型poly 
    输出：shapely poly 
    """
    from shapely.ops import cascaded_union
    # 先获取 shapely 格式 list
    polys = []
    for hull in poly_dict.values():
        poly = Polygon([(pt[1], pt[0]) for pt in hull])
        polys.append(poly) 
        
    return cascaded_union(polys)