Shapely简介

作者 : RY    标签: polygon shapely 几何

Shapely是一个很方便的几何图形处理扩展库,本文通过几个实例介绍它的一些基本功能。

Shapely是一个基于几何引擎GEOS的Python扩展库。使用它可以很方便的处理多线和多边形等几何形状。本文通过一些实例介绍Shapely的一些基本用法。

Linux下输入如下命令:

sudo apt-get install libgeos-dev
sudo pip install shapely

Windows可以从下面的页面下载:

http://www.lfd.uci.edu/~gohlke/pythonlibs/

In [1]:
%pylab inline
import pylab as pl
Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.kernel.zmq.pylab.backend_inline].
For more information, type 'help(pylab)'.

Shapely使用ctypes调用GEOS库,如果安装了libgeos-dev,则可以编译speedups模块,提高Python和GEOS的接口速度。下面的程序启动提速功能。

In [2]:
from shapely import speedups
speedups.enable()

Shapely只能处理由线段组成的几何图形,对于椭圆这样的曲线图形,需要用许多细分的线段近似表示。在Shapely中用LineString对象表示多条线段连接而成的多线图形。

下面的ellipse_polyline()使用椭圆的参数方程,计算出椭圆上n个点的坐标,将这些坐标点用线段连接起来,就可以近似表示椭圆:

In [3]:
import numpy as np
def ellipse_polyline(ellipse, n=100):
    t = np.linspace(0, 2*np.pi, n)
    st = np.sin(t)
    ct = np.cos(t)
    x0, y0, a, b, angle = ellipse
    angle = np.deg2rad(angle)
    sa = np.sin(angle)
    ca = np.cos(angle)
    p = np.empty((n, 2))
    p[:, 0] = x0 + a * ca * ct - b * sa * st
    p[:, 1] = y0 + a * sa * ct + b * ca * st
    return p

让我们创建两个椭圆,并将其绘制出来,可以看到这两个椭圆最多可以有4个交点:

In [4]:
a = ellipse_polyline((1, 1, 2, 1, 45))
b = ellipse_polyline((2, 0.5, 5, 1.5, -30))
fig, ax = subplots()
ax.set_aspect("equal")
ax.plot(a[:,0], a[:,1])
ax.plot(b[:,0], b[:,1]);

数组a和b的形状为(n, 2),这样的数组可以直接通过LineString()转换成LineString对象:

In [5]:
from shapely.geometry.polygon import LineString
line_a = LineString(a)
line_b = LineString(b)

LineString对象提供了数组接口,可通过numpy.asarray()将其转换回数组:

In [6]:
np.allclose(a, np.asarray(line_a))
Out[6]:
True

调用intersection()方法可以计算两个几何对象之间的交点,返回的是一个表示多点的MultiPoint对象,同样可以通过numpy.asarray()将其转换成数组:

In [7]:
points = np.asarray(line_a.intersection(line_b))
points
Out[7]:
array([[-0.53075832,  0.39651172],
       [ 0.65717929, -0.43986125],
       [ 1.43832647,  2.47828182],
       [ 2.57965213,  1.89525131]])
In [8]:
ax.plot(points[:, 0], points[:, 1], "o")
fig
Out[8]:

下面的程序创建随机的50个点,我们希望得到由这些线条包围的区域的面积。

In [9]:
np.random.seed(0)
points = np.random.normal(size=(50, 2))
points = np.vstack((points, points[:1,:]))
fig, ax = subplots()
ax.plot(points[:,0], points[:, 1]);

首先将这些点转换成LineString对象:

In [10]:
lines = LineString(points)

LineString对象中的线是没有宽度的。我们可以通过buffer()方法用一个带宽度笔沿着这个路径描绘,得到的是一个Polygon对象:

In [11]:
polygon = lines.buffer(0.05)

Shapely中的Polygon对象包含两个区域:

  • exterior: 一个LinearRing对象,LinearRing是一组组成环形的线段,线段是不能自相交的。这个exterior属性保存的是Polygon外部的边缘。

  • interiors: 一个LinearRingSequence对象,其中的每个LinearRing对象都从Polygon中挖出一块。

下面让我们用matplotlib将这个多边形画出来,其中:

  • 蓝色线条表示Polygon的外部边缘:exterior,红色线段表示内部边缘,绿色的点线表示原始的线段。
In [12]:
fig, ax = subplots()
p = pl.Polygon(np.asarray(polygon.exterior), facecolor="none", edgecolor="blue")
ax.add_patch(p)
for ring in polygon.interiors:
    p = pl.Polygon(np.asarray(ring), facecolor="none", edgecolor="red")
    ax.add_patch(p) 
ax.plot(points[:,0], points[:, 1], "g:")
ax.relim()
ax.autoscale()

我们可以将外部边缘转换成一个新的Polygon对象,然后再计算其面积:

In [13]:
from shapely.geometry.polygon import Polygon

Polygon(polygon.exterior).area
Out[13]:
7.8145147170093185

当然前面的代码为了演示效果,使用了较粗的笔触,因此得到的面积比实际值要大一些,为了得到更准确的面积,可以用非常细的笔触:

In [14]:
polygon = lines.buffer(1e-6)
area = Polygon(polygon.exterior).area
print area
6.20552431436
In [15]:
fig, ax = subplots()
p = pl.Polygon(np.asarray(polygon.exterior), facecolor="red", edgecolor="none", alpha=0.5)
ax.add_patch(p)
ax.plot(points[:,0], points[:, 1], "b");

Shapely提供了几何图形的集合运算,可以很方便地求多个多边形的并集、交集以及差集。下面我们对一系列旋转的椭圆求并集:

In [16]:
ellipses = [Polygon(ellipse_polyline((0, 0, 3, 1, angle)))
    for angle in np.linspace(0, 360, 7, endpoint=False)]

union()方法可以计算两个几何图形的并集:

In [17]:
p = np.asarray(ellipses[0].union(ellipses[1]).exterior)
subplot(111, aspect="equal")
plot(p[:,0], p[:,1], lw=2);

为了求多个多边形的并集,可以使用cascaded_union():

In [18]:
from shapely.ops import cascaded_union
union = cascaded_union(ellipses)
p = np.asarray(union.exterior)
subplot(111, aspect="equal")
plot(p[:,0], p[:,1], lw=2);

loading...