- 1Visual Effect Graph——块和操作符_visual effects graph cellular noise
- 2书籍推荐：《Java数据结构与算法》
- 3怎样切换git地址_git切换地址
- 4【机器学习】使用R语言进行机器学习特征选择_r语言特征选择
- 5容易内存泄露的几种情况_c# 哪些代码容易引起内存泄漏
- 6程序员修炼之道（一）_卫语句 程序员修炼之道
- 7内网渗透信息收集总结_渗透dmz区域内网网段获取ip
- 8win10 + Qt5.14.0（mingw73_64） 下的 OpenCV4.5.1源码编译 注意事项_missing:python_executable interpreter
- 9c# try catch finaly用法_c#中try catch的用法
- 10python毕设选题 - 大数据B站数据分析与可视化 - python 数据分析 大数据
Python之几何集合操作模块Shapely的简介_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)
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
输出结果如下：
POINT (1 1) POINT (2 2)
MULTIPOINT (1 1, 2 2)
- 1
- 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)
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
输出结果如下：
0.0
0.0
[(1.0, 1.0)]
(1.0, 1.0, 1.0, 1.0)
- 1
- 2
- 3
- 4
计算两个点之间的欧式距离代码如下：
# -*- 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
- 2
- 3
- 4
- 5
- 6
输出结果为：
1.4142135623730951
- 1
线（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)
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
输出结果如下：
0.0
1.4142135623730951
[(1.0, 1.0), (2.0, 2.0)]
(1.0, 1.0, 2.0, 2.0)
- 1
- 2
- 3
- 4
多边形（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
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
输出结果如下：
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
- 2
- 3
- 4
- 5
可以看到，该三角形的面积为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)))
- 1
- 2
- 3
- 4
- 5
- 6
输出结果如下：
True
False
- 1
- 2
地理对象的实现与操作
Shapely模块还支持对地理信息的操作。在此之前，我们先需要了解下GeoJson。
GeoJson
GeoJson是用于编码地理数据结构的格式，是Json的子集，参考下面的例子：
{
"type": "Feature",
"geometry": {
"type": "Point",
"coordinates": [125.6, 10.1]
},
"properties": {
"name": "Dinagat Islands"
}
}
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
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'}
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
- 48
- 49
- 50
- 51
- 52
- 53
- 54
- 55
- 56
- 57
- 58
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))))
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
输出结果为：
True
False
- 1
- 2
总结
本文介绍了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
