30分鐘學(xué)會(huì)shapely空間幾何分析

公眾號(hào)后臺(tái)回復(fù)關(guān)鍵字:"源碼"，獲取本文全部代碼。

實(shí)踐證明，它的以下一些功能特性非常常用：

• 幾何對(duì)象可以和numpy.array互相轉(zhuǎn)換。

• 可以輕松求線的長(zhǎng)度(length)，面的面積（area)，對(duì)象之間的距離(distance),最小最大距離(hausdorff_distance)。

• 可以輕松求幾何對(duì)象之間的關(guān)系：相交(intersect)，包含(contain)，求相交區(qū)域(intersection)等。

• 可以輕松對(duì)幾何對(duì)象求幾何中心(centroid),緩沖區(qū)(buffer),最小旋轉(zhuǎn)外接矩形(minimum_rotated_rectangle)等。

• 可以求線的插值點(diǎn)(interpolate),可以求點(diǎn)投影到線的距離(project),可以求幾何對(duì)象之間對(duì)應(yīng)的最近點(diǎn)(nearestPoint)

• 可以輕松對(duì)幾何對(duì)象進(jìn)行旋轉(zhuǎn)(rotate)和縮放(scale)。

#安裝shapely
!pip install shapely 

from shapely import geometry as geo
from shapely import wkt 
from shapely import ops
import numpy as np 

一，Point對(duì)象

# 創(chuàng)建Point對(duì)象
pt1 = geo.Point([0,0])
print(pt1) 
coord = np.array([0,1])
pt2 = geo.Point(coord)
print(pt2)
pt3 = wkt.loads("POINT(1 1)")
print(pt3)

#批量可視化
geo.GeometryCollection([pt1,pt2,pt3])

# 常用屬性

print(pt1.x) 
print(pt1.y)
print(list(pt1.coords))  
print(np.array(pt1)) #可以和np.array互轉(zhuǎn)

# 常用方法

d = pt2.distance(pt1)
print(d)

二, LineString對(duì)象

# 創(chuàng)建LineString對(duì)象
line1 = geo.LineString([(0,0),(1,-0.1),(2,0.1),(3,-0.1),(5,0.1),(7,0)])
line1 

arr = np.array([(2,2),(3,2),(4,3)])
line2 = geo.LineString(arr)
line2 

line3 = wkt.loads("LineString(-2 -2,4 4)")
line3 

# 常用屬性
print(line2.length) 
print(list(line2.coords)) 
print(np.array(line2))  #可以和np.array互轉(zhuǎn)
print(line2.bounds) #坐標(biāo)范圍

center = line2.centroid #幾何中心
geo.GeometryCollection([line2,center])

bbox = line2.envelope #最小外接矩形
geo.GeometryCollection([line2,bbox]) 

rect = line2.minimum_rotated_rectangle #最小旋轉(zhuǎn)外接矩形
geo.GeometryCollection([line2,rect])

# 常用方法

d1 = line1.distance(line2) #線線距離
print(d1)
d2 = line1.distance(geo.Point([-1,0])) #線點(diǎn)距離
print(d2)
d3 = line1.hausdorff_distance(line2) #最小最大距離
print(d3)

pt_half = line1.interpolate(0.5,normalized=True) #插值
geo.GeometryCollection([line1,pt_half])

ratio = line1.project(pt_half,normalized=True) #投影
print(ratio)

line1_simplify = line1.simplify(0.5)  #化簡(jiǎn) DouglasPucker算法
print(line1)
print(line1_simplify)
line1_simplify 

buffer_with_circle = line2.buffer(0.2)  #端點(diǎn)按照半圓擴(kuò)展
geo.GeometryCollection([line2,buffer_with_circle])

buffer_without_circle = line2.buffer(0.2,cap_style=2) #端點(diǎn)不擴(kuò)展
geo.GeometryCollection([line2,buffer_without_circle])

buffer_with_square = line2.buffer(0.2,cap_style=3) #端點(diǎn)按照方形擴(kuò)展
geo.GeometryCollection([line2,buffer_with_square])

buffer_round_join = line2.buffer(0.2,join_style=1) #圓弧連接
geo.GeometryCollection([line2,buffer_round_join])

buffer_angle_join = line2.buffer(0.2,join_style=2) #折角連接
geo.GeometryCollection([line2,buffer_angle_join])

print(line2.intersects(line3)) #線線關(guān)系，是否相交
print(line2.intersection(line3)) #線線交點(diǎn)
print(line2.contains(geo.Point(2.5,2))) #點(diǎn)線關(guān)系

三，Polygon對(duì)象

# 創(chuàng)建Polygon對(duì)象
poly1 = geo.Polygon([(0,0),(1,0),(1,1),(0,1),(0,0)]) #起點(diǎn)和終點(diǎn)相同
poly1

coords = np.array([(0,0),(1,0.1),(2,0),(1,2),(0,0)])
poly2 = geo.Polygon(coords)
poly2 

#第一個(gè)括號(hào)是外部坐標(biāo)，后面的是內(nèi)部空洞坐標(biāo)
poly3 = wkt.loads("POLYGON((0 0,2 0,2 2,0 2,0 0),(0.5 0.5,1.5 0.5,1.5 1.5,0.5 1.5,0.5 0.5))")
poly3 

#創(chuàng)建bbox對(duì)象
poly4 = geo.Polygon.from_bounds(xmin=0,ymin=0,xmax=20,ymax=20)
poly4 

#常用屬性
print(poly1.area)  #面積
print(poly1.length) #周長(zhǎng)
print(np.array(poly1.exterior))  #外圍坐標(biāo)點(diǎn)
print(poly3.bounds) #坐標(biāo)范圍

center = poly3.centroid #幾何中心
geo.GeometryCollection([center,poly3]) 

poly3.boundary #邊緣

rect = poly2.minimum_rotated_rectangle #最小外接矩形
geo.GeometryCollection([rect,poly2])

# 常用方法
r1 = poly2.contains(geo.Point(0,0)) #面點(diǎn)關(guān)系
print(r1)

r2 = poly2.intersects(geo.LineString([(0,0),(5,5)])) #面線關(guān)系
print(r2)

r3 = poly2.intersects(poly3) #面面關(guān)系
print(r3)

geo.GeometryCollection([poly1,line3])

inter = poly1.intersection(line3) #面線交集
geo.GeometryCollection([poly1,inter])

geo.GeometryCollection([poly1,poly2])

poly1.intersection(poly2) #面面交集

poly1.union(poly2) #面面并集

poly2.difference(poly1) #面面補(bǔ)集

poly2.simplify(0.5) #簡(jiǎn)化

print(poly2.area)
poly2_bigger = poly2.buffer(0.2) #外擴(kuò)面積變大
print(poly2_bigger.area)
poly2_smaller = poly2.buffer(-0.2) #內(nèi)擴(kuò)面積變小
print(poly2_smaller.area)
poly2_smaller 

# MultiPoint 多點(diǎn)

x = np.linspace(0,2*np.pi,10)
y = np.sin(x)
points = [geo.Point(i,j) for i,j in zip(x,y)]
multipoints = geo.MultiPoint(points ) 
multipoints

hull = multipoints.convex_hull  #凸包
geo.GeometryCollection([hull,multipoints]) 

# MultiLineString 多線
multilines = geo.MultiLineString([line1,line2])
multilines 

# MultiPolygon 多面

multipolys = geo.MultiPolygon([poly1,poly2])
multipolys 

# GeometryCollection 對(duì)象集合

geoms = [pt1,pt2,pt3,line3,poly3]
geo.GeometryCollection(geoms) #方便在jupyter 中對(duì)多個(gè)幾何對(duì)象可視化

以下是一些非常有用但是不屬于某個(gè)類(lèi)的方法的函數(shù)。

• ops.nearest_points 求最近點(diǎn)

• ops.split 分割線

• ops.substring  求子串

• affinity.rotate 旋轉(zhuǎn)幾何體

• affinity.scale 縮放幾何體

• affinity.translate 平移幾何體

from shapely import ops,affinity 

poly1 = geo.Polygon([(0,0),(2,0),(1,1),(0,0)])
poly2 = geo.Polygon([(4,0),(6,0),(6,2),(4,2),(4,0)])

p1,p2 = ops.nearest_points(poly1,poly2)

geo.GeometryCollection([poly1,poly2,p1,p2]) 

poly1_rot30 = affinity.rotate(poly1,30,origin = "centroid")
geo.GeometryCollection([poly1,poly1_rot30])

poly1_scale = affinity.scale(poly1,xfact=2.0,yfact=2.0)

geo.GeometryCollection([poly1,poly1_scale])