本文介绍了Python库Shapely的基础几何对象操作，如点、线和多边形的创建、属性获取及距离计算，并展示了如何进行地理分析，如判断点是否在多边形内，结合GeoJson数据处理地理围栏。此外，还通过示例展示了在实际应用中如何判断地理位置是否在特定地理区域内。
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
1368
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
