【541】shapely 相关功能

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

参考：Shapely用户手册

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

---

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

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

?

1 2 3 4 5 6 7 8 9 10 11 12 13

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. 多边形显示

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

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21

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')

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17

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')

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28

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 的具体信息。

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19

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 的思路

?

1 2 3 4 5 6 7 8 9 10

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

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14

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)