Python之几何集合操作模块Shapely的简介
Shapely是基于笛卡尔坐标的几何对象操作和分析的Python库，底层基于GEOS和JTS库。Shapely不关心数据格式或坐标系，但可以很容易地与这些文件包集成。
Shapely的pypi访问网址为:https://pypi.org/project/Shapely/，官方文档访问网址为：https://shapely.readthedocs.io/en/latest/。
本文将会介绍Shapely的一些基础几何对象的操作以及其在地理分析上的简单应用，其中基础几何对象包括点（Point）、线（LineString）、多边形（Polygon）。
基础几何对象的实现与操作
点（Point）
我们以下图中的点A(1,1)与点B(2,2)为例，演示Shapely中的点（Point）的实现与相关操作。
在Shapely模块中，使用
shapely.geometry.Point代表表，使用shapely.geometry.MultiPoint代表多个点。实现上述A、B两点的Python代码如下：# -*- coding: utf-8 -*-
from shapely.geometry import Point
from shapely.geometry import MultiPoint
point_a = Point(1, 1)
point_b = Point(2, 2)
points = MultiPoint([[1, 1], [2, 2]])
print(point_a, point_b)
print(points)
输出结果如下：
POINT (1 1) POINT (2 2)
MULTIPOINT (1 1, 2 2)
点（Point）也有面积（area）、长度（length）、坐标（coords）和边界（bounds）属性，其中面积和长度的值为0，如下所示：
# -*- coding: utf-8 -*-
from shapely.geometry import Point
point_a = Point(1, 1)
print(point_a.area)
print(point_a.length)
print(list(point_a.coords))
print(point_a.bounds)
输出结果如下：
0.0
0.0
[(1.0, 1.0)]
(1.0, 1.0, 1.0, 1.0)
计算两个点之间的欧式距离代码如下：
# -*- coding: utf-8 -*-
from shapely.geometry import Point
point_a = Point(1, 1)
point_b = Point(2, 2)
print(point_a.distance(point_b))
输出结果为：
1.4142135623730951
线（LineString）
我们以下图中的点A(1,1)与点B(2,2)组成的线段为例，演示Shapely中的线段（LineString）的实现与相关操作。
在Shapely模块中，使用
shapely.geometry.LineString代表线段，使用shapely.geometry.MultiLineString代表多个线段。实现上述AB线段及其相关属性的Python代码如下：# -*- coding: utf-8 -*-
from shapely.geometry import LineString
line_ab = LineString([[1, 1], [2, 2]])
print(line_ab.area)
print(line_ab.length)
print(list(line_ab.coords))
print(line_ab.bounds)
输出结果如下：
0.0
1.4142135623730951
[(1.0, 1.0), (2.0, 2.0)]
(1.0, 1.0, 2.0, 2.0)
多边形（Polygon）
我们以下图中的点A(1,1)、点B(2,2)、点C(1,3)组成的三角形为例，演示Shapely中的多边形（Polygon）的实现与相关操作。
在Shapely模块中，使用
shapely.geometry.Polygon代表多边形，使用shapely.geometry.MultiPolygon代表多个多边形。实现上述三角形ABC及其相关属性的Python代码如下：# -*- coding: utf-8 -*-
from shapely.geometry import Polygon
polygon_1 = Polygon(((1, 1), (2, 2), (3, 1)))
print(polygon_1.area)
print(polygon_1.length)
print(list(polygon_1.exterior.coords))
print(list(polygon_1.interiors))
print(list(polygon_1.bounds))
输出结果如下：
1.0
4.82842712474619
[(1.0, 1.0), (2.0, 2.0), (3.0, 1.0), (1.0, 1.0)]
[]
[1.0, 1.0, 3.0, 2.0]
可以看到，该三角形的面积为1，边长为4.828，外部的坐标点为[(1.0, 1.0), (2.0, 2.0), (3.0, 1.0), (1.0, 1.0)]，无内部左边点，边界（(minx, miny, maxx, maxy)）为[1.0, 1.0, 3.0, 2.0]。
另外，多边形几何对象还可以判断某个点是否在其内部，示例Python代码如下：
# -*- coding: utf-8 -*-
from shapely.geometry import Polygon, Point
polygon_1 = Polygon(((1, 1), (2, 2), (3, 1)))
print(polygon_1.contains(Point(2, 1.5)))
print(polygon_1.contains(Point(1.5, 2)))
输出结果如下：
True
False
地理对象的实现与操作
Shapely模块还支持对地理信息的操作。在此之前，我们先需要了解下GeoJson。
GeoJson
GeoJson是用于编码地理数据结构的格式，是Json的子集，参考下面的例子：
{
"type": "Feature",
"geometry": {
"type": "Point",
"coordinates": [125.6, 10.1]
},
"properties": {
"name": "Dinagat Islands"
}
}
GeoJson支持的几何对象为：Point, LineString, Polygon, MultiPoint, MultiLineString和 MultiPolygon。几何对象加上其属性为Feature对象。Feature对象的集合为FeatureCollection对象。
地理围栏
地理围栏指的是某个地理区域的边界信息，比如澳门圣方济各堂区（数据来源于：http://datav.aliyun.com/portal/school/atlas/area_selector）的地理围栏数据为：
{'features': [{'geometry': {'coordinates': [[[[113.554447, 22.107312],
[113.559643, 22.106216],
[113.561046, 22.106214],
[113.576328, 22.108147],
[113.578953, 22.108929],
[113.580194, 22.109694],
[113.597564, 22.125115],
[113.603681, 22.132371],
[113.603873, 22.132865],
[113.601946, 22.138113],
[113.589178, 22.144281],
[113.584062, 22.137923],
[113.583056, 22.136473],
[113.582311, 22.135821],
[113.581374, 22.135379],
[113.580612, 22.135238],
[113.579329, 22.135134],
[113.57802, 22.135011],
[113.577257, 22.134855],
[113.576591, 22.134553],
[113.575827, 22.134051],
[113.574848, 22.133293],
[113.574068, 22.132806],
[113.573184, 22.13248],
[113.572257, 22.132332],
[113.571408, 22.132305],
[113.568618, 22.132386],
[113.568324, 22.132349],
[113.568107, 22.132224],
[113.568003, 22.131992],
[113.567889, 22.128884],
[113.567856, 22.128712],
[113.567778, 22.128615],
[113.56763, 22.128485],
[113.567474, 22.128416],
[113.567266, 22.128386],
[113.564798, 22.128485],
[113.564651, 22.128477],
[113.564529, 22.128445],
[113.56439, 22.128297],
[113.56407, 22.127443],
[113.564018, 22.127342],
[113.563949, 22.127289],
[113.553916, 22.125216],
[113.553597, 22.117626],
[113.553389, 22.112663],
[113.554447, 22.107312]]]],
'type': 'MultiPolygon'},
'properties': {'acroutes': [100000, 820000],
'adcode': 820008,
'center': [113.5599542, 22.12348639],
'centroid': [113.577403, 22.123389],
'childrenNum': 0,
'level': 'district',
'name': '圣方济各堂区',
'parent': {'adcode': 820000}},
'type': 'Feature'}],
'type': 'FeatureCollection'}
http://geojson.io/网站是GeoJson可视化网址，上述围栏数据的可视化结果为：
判断地理位置是否在围栏中
我们以黑沙海滩和大东海广场两个地点为例，判断他们是否在澳门圣方济各堂区。注意，黑沙海滩的经纬度为(113.57, 22.12)，位于澳门圣方济各堂区，而大东海广场的经纬度为(109.52, 18.22)，位于三亚。查询经纬度的网站为：http://www.bigemap.net/city-208645.html。
使用Shapely模块可以帮助我们实现判断，示例代码如下：
# -*- coding: utf-8 -*-
import json
from shapely.geometry import Polygon, Point
with open('圣方济各堂区.json', 'r', encoding='utf-8') as f:
content = json.loads(f.read())
# 读取圣方济各堂区地理围栏数据，并封装成Polygon
macau_coords = content["features"][0]['geometry']['coordinates'][0][0]
polygon = Polygon(macau_coords)
# 判断黑沙海滩是否位于圣方济各堂区
print(polygon.contains(Point((113.57, 22.12))))
# 判断大东海广场是否位于圣方济各堂区
print(polygon.contains(Point((109.52, 18.22))))
输出结果为：
True
False
总结
本文介绍了Shapely的一些基础几何对象的操作以及其在地理分析上的简单应用，其中我们不难找到一些我们感兴趣的应用点，这也是Shapely模块的价值所在。
感谢大家阅读~
2021.12.12于上海长宁区
参考网址
Shapely官网文档地址: https://shapely.readthedocs.io/en/latest/manual.html
Shapely的Pypi地址：https://pypi.org/project/Shapely/
GeoJSON介绍网址：https://geojson.org/
GeoJSON可视化网址: http://geojson.io/
阿里云地理围栏数据网址: http://datav.aliyun.com/portal/school/atlas/area_selector
在线地图查询经纬度网站：http://www.bigemap.net/city-208645.html