GIS算法：7_python拓扑套件shapely

可用于GIS数据处理与空间计算的Python包有GDAL和shapely。

相对来说，GDAL功能更全面，还可以用于数据转换、瓦片地图发布、栅格影像分析等，shapely只能进行基本的数据处理和空间计算。

但大多数情况下shapely就完全够用了。

GDAL的官网：https://gdal.org/

Shapely的中文文档：https://www.osgeo.cn/shapely/

 

本例只讲shapely的用法。

可以通过cmd安装shapely，打开cmd，输入python，回车，查看python的版本号。

去下载shapely安装程序，https://www.lfd.uci.edu/~gohlke/pythonlibs/#shapely

根据python版本和电脑位数，选了Shapely-1.6.4.post2-cp37-cp37m-win32.whl下载。

安装的时候，打开cmd，运行语句：pip install D:\gismap\软件\Shapely-1.6.4.post2-cp37-cp37m-win32.whl，输入，回车。其中D:\gismap\软件\是放安装包的路径。

如果提示，

提示要升级pip，写啥就敲啥。python -m pip install --upgrade pip

升级完了，接着运行上一个语句，把shapely安装了。

 

 

构建RTree：生成50×50个点，根据这些点生成缓冲1.0的圆，用这些圆构建RTree，在RTree中查点[1,1]落在哪些圆里。
 

from shapely.geometry import Point
from shapely.strtree import STRtree
polygons=[]
for i in range(0,50):
    for j in range(0,50):
        polygon=Point([i,j]).buffer(1.0)
        polygon.__setattr__('id',str(i)+str(j))
        polygons.append(polygon)
tree=STRtree(polygons)
result=tree.query(Point(1,1))
for r in result:
    print(r.__getattribute__('id'))

 

 

 

线裁切面：如图所示，用红色线去裁切蓝色面，结果会生成3个面。

 

 

from shapely import wkt
from shapely.ops import polygonize
polygon=wkt.loads('POLYGON ((220 350, 400 440, 635 249, 380 80, 174 164, 179 265, 220 350))')
polyline=wkt.loads('LINESTRING (570 400, 392 315, 299 215, 430 140, 530 240, 450 360, 460 480)')
boundary=polygon.boundary
polyline=polyline.union(boundary)
lineList = []
for i in range(0,len(polyline.geoms)):
    lineList.append(polyline.geoms[i])
clipPolygon=polygonize(lineList)
buffer=polygon.buffer(1)
for c in clipPolygon:
    if buffer.contains(c):
        print(c)

 

 

将裁切线换成LINESTRING (320 330, 500 280, 400 150, 290 200, 400 360)，结果如下：

 

 

 

结果可以在JTS TestBuilder中查验。

在控制台input的文本框中输入wkt文本，点击load geometrys，就可以在面板中展示图形。

 

Shapely常用函数列表：

Geometric Objects
object.area	Returns the area (float) of the object.
object.bounds	Returns a (minx, miny, maxx, maxy) tuple (float values) that bounds the object.
object.length	Returns the length (float) of the object.
object.geom_type	Returns a string specifying the Geometry Type	类型
object.distance(other)	Returns the minimum distance (float) to the othergeometric object.
object.representative_point()	Returns a cheaply computed point that is guaranteed to be within the geometric object.	返回几何图形内的标志点，不是图心或质心，因为凹多边形的图心或质心可能在几何图形外
Unary Predicates		一元谓词关系判断</