本文记录使用 Python 库 pyproj 实现地理坐标转换的流程。
简介
pyproj是一个Python库，用于执行坐标转换和投影变换。它基于Proj库，后者是一个C库，用于处理地图投影和坐标变换。pyproj提供了Python语言的接口，使得用户可以方便地使用这些功能。
坐标转换在地理信息系统（GIS）、地图制作、卫星导航、地震勘探以及许多其他科学和工程领域中都是基本需求。地球上的点的位置通常用一系列的数对来表示，这些数对称为坐标，可以是笛卡尔坐标、经纬度或其他形式。由于各种原因（如地图制作、技术限制等），这些坐标可能需要从一个系统转换到另一个系统。pyproj库正是用来执行这些转换的。
以下是pyproj的一些关键特点：
- 坐标系统支持广泛：
pyproj支持大量的坐标系统，包括各种国际和区域标准，如EPSG、ESRI、OGC等。
- 灵活性：它允许用户指定任意的源和目标坐标系统，以及相关的参数。
- 高性能：由于底层的Proj库是用C语言编写的，
pyproj在执行坐标变换时提供了较好的性能。
- 易用性：
pyproj的API设计简单直观，使得用户可以轻松地进行坐标变换和投影操作。
- 社区支持：
pyproj有一个活跃的社区，不断更新和改进库的功能。
- 跨平台：
pyproj可以在各种操作系统上运行，包括Windows、macOS和Linux。
Git 仓库：https://github.com/pyproj4/pyproj?tab=readme-ov-file
官方文档：https://pyproj4.github.io/pyproj/stable/
安装方法
地理坐标转换
pyproj 开发了很多 API ，这里简单介绍常用的地理坐标转换使用方法。
示例数据
地表上中国境内一点：
- 4538 坐标系下： (X, Y) = (548888, 4940000)
- 4326 坐标系下：(lat, lon) = (44.59390182452887, 87.6157036862893)
Pyproj1
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
| from pyproj import Proj, transformif __name__ == '__main__':
# 创建一个Transformer对象，用于坐标转换
from_coor = "epsg:4538" # 输入坐标系
to_coor = "epsg:4326" # 输出坐标系 transformer = Transformer.from_crs(from_coor, to_coor) # 输入坐标
x = 548888
y = 4940000 # 转换坐标
#pyproj1
input_proj = Proj(init='epsg:4538')
output_proj = Proj(init='epsg:4326') # WGS84坐标系 lon1, lat1 = transform(input_proj, output_proj, x, y)
print(f"pyproj1 转换后的坐标：({lon1}, {lat1})") -->
pyproj1 转换后的坐标：(87.6157036862893, 44.59390182452887)
|
Pyproj2
上代码
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
| from pyproj import Transformerif __name__ == '__main__':
# 创建一个Transformer对象，用于坐标转换
from_coor = "epsg:4538" # 输入坐标系
to_coor = "epsg:4326" # 输出坐标系 transformer = Transformer.from_crs(from_coor, to_coor) # 输入坐标
x = 548888
y = 4940000 # 转换坐标
# pyproj2
lat2, lon2 = transformer.transform(y, x)
print(f"pyproj2 转换后的坐标：({lon2}, {lat2})") -->
pyproj2 转换后的坐标：(87.6157036862893, 44.59390182452887)
|
两种实现比较
当前网上资料大多数以 pyoroj1 的实现方法为主，但是事实上 pyproj1 相比 pyproj2 速度慢很多
耗时对比：
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
| from pyproj import Proj, transform
from pyproj import Transformer
import timeif __name__ == '__main__':
# 创建一个Transformer对象，用于坐标转换
from_coor = "epsg:4538" # 输入坐标系
to_coor = "epsg:4326" # 输出坐标系 transformer = Transformer.from_crs(from_coor, to_coor) # 输入坐标
x = 548888
y = 4940000 # 转换坐标 run_times = 100 #pyproj1
input_proj = Proj(init='epsg:4538')
output_proj = Proj(init='epsg:4326') # WGS84坐标系 start = time.time()
for _ in range(run_times):
lon1, lat1 = transform(input_proj, output_proj, x, y)
end = time.time()
print(f"pyproj1 平均转换耗时：{(end - start)/run_times}s")
print(f"pyproj1 转换后的坐标：({lon1}, {lat1})") # pyproj2
start = time.time()
for _ in range(run_times):
lat2, lon2 = transformer.transform(y, x)
end = time.time()
print(f"pyproj2 平均转换耗时：{(end - start)/run_times}s")
print(f"pyproj2 转换后的坐标：({lon2}, {lat2})")
|
输出结果：
1
2
3
4
| pyproj1 平均转换耗时：0.09519411087036132s
pyproj1 转换后的坐标：(87.6157036862893, 44.59390182452887)
pyproj2 平均转换耗时：1.6641616821289064e-06s
pyproj2 转换后的坐标：(87.6157036862893, 44.59390182452887)
|
差了不是一点半点 …
官方建议
官方链接：https://pyproj4.github.io/pyproj/stable/gotchas.html#upgrading-to-pyproj-2-from-pyproj-1
Upgrading to pyproj 2 from pyproj 1
We recommended using the pyproj.transformer.Transformer and pyproj.crs.CRS in place of the pyproj.Proj and pyproj.transformer.transform().
pyproj 1 style:
1
2
3
4
5
6
| >>> from functools import partial
>>> from pyproj import Proj, transform
>>> proj_4326 = Proj(init="epsg:4326")
>>> proj_3857 = Proj(init="epsg:3857")
>>> transformer = partial(transform, proj_4326, proj_3857)
>>> transformer(12, 12)
|
pyproj 2 style:
1
2
3
| >>> from pyproj import Transformer
>>> transformer = Transformer.from_crs("EPSG:4326", "EPSG:3857")
>>> transformer.transform(12, 12)
|
结论
无脑用 pyproj2 就完了。
参考资料
文章链接：
https://www.zywvvd.com/notes/coding/python/pyproj/pyproj/
