Python GNSS、地理坐标以及投影坐标之间的转换

最近捣鼓了一些室外机器人的路线，位置信息来源于GPS接收器拿到的信息，之前接收项目的时候他们采用了相对坐标散点图画路径，现在想基于经纬度利用现成的高德地图接口实现，在提取数据的过程中发现了一些问题，记录下来。

1. GNSS

卫星导航系统（Global Navigation Satellite System,GNSS）是覆盖全球的自主地利空间定位的卫星系统，允许小巧的电子接收器确定它的所在位置（经度、纬度和高度），并且经由卫星广播沿着视线方向传送的时间信号精确到10米的范围内。接收机计算的精确时间以及位置，可以作为科学实验的参考。

截至2020年6月，只有美国的全球定位系统 （GPS；共由24颗卫星组成）、俄罗斯的格洛纳斯系统（GLONASS）和中国的北斗卫星导航系统（BDS）[1]覆盖全球。欧洲联盟的伽利略定位系统则为在初期部署阶段的全球导航卫星系统，预定最早到2020年才能够充分的运作[2]。一些国家，包括法国、日本和印度[3]，都在发展区域导航系统。

每个覆盖全球的系统通常都是由20-30颗卫星组成的卫星集群，以中地球轨道分布在几个轨道平面上。实际的系统各自不同，但是使用的轨道倾斜都大于50°，和轨道周期大约都是12小时（高度大约20,000千米（12,000英里））。

不同的GNSS

GNSS的工作原理分很多部分，比如：GNSS接收器是如何工作的？；GPS 接收器如何知道我在哪里？；GPS信号；计算出到卫星的距离等等

这并不是本篇文章所要探讨的，具体可以看GNSS一系列相关问题

2. 经纬度

经纬度这个概念想必不用我说了吧。其实我还真忘记了，太难了~

小学是没地理这门课的，那肯定就是初中知识。

这边主要讲讲经纬度的三种方式：

经纬度格式分为三种：度、度-分、度-份-秒

1.） ddd.ddddd °【度 . 度 格式】的十进制小数部分（5位）

2.） ddd°mm.mmm’ 【度 . 分 . 分 格式】的十进制小数部分（3位）

3.) ddd°mm’ss’’ 【度 . 分 . 秒 格式】

Google 使用的是第三种格式 度°分’秒’’

度分转换：
将度分单位数据转换为度单位数据
度=度+分/60
例如：
经度 = 116°20.12’
纬度 = 39°12.34’
经度 = 116 + 20.12 / 60 = 116.33533°
纬度 = 39 + 12.34 / 60 = 39.20567°

度分秒转换：
将度分秒单位数据转换为度单位数据
度 = 度 + 分 / 60 + 秒 / 60 / 60
例如：
经度 = 116°20’43”
纬度 = 39°12’37”
经度 = 116 + 20 / 60 + 43 / 60 / 60 = 116.34528°
纬度 = 39 + 12 / 60 + 37 / 60 / 60 = 39.21028°

3. GPGGA信息转换为经纬度

NMEA-0183是美国国家海洋电子协会（National Marine Electronics Association）为海用电子设备制定的标准格式。目前业已成了GPS导航设备统一的RTCM（Radio Technical Commission for Maritime services ）标准协议。

NMEA0183的六种输出协议：GPGGA、GPGLL、GPGSA、GPGSV、GPRMC、GPVTG。

而我们这里要解析的是**$GPGGA(GPS定位信息)**

pynmea2 是一个用来处理 NMEA 0183 协议的第三方模块，pynmea2 模块兼容 Python 2 和 Python 3，能够解析 GSA、GGA、GSV、RMC、VTG、GLL 等 NMEA 0183 协议定义的各类数据，功能强大。该模块目前以 MIT 协议开源并托管在 Github 网站上。

注意，解析字符串中 NMEA 0183 协议的数据，可以使用 pynmea2.parse(data, check=False) 方法，其中的 check 参数指定是否对消息中的检校字段进行检查。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
>>> import pynmea2
>>> msg = pynmea2.parse("$GPGGA,184353.07,1929.045,S,02410.506,E,1,04,2.6,100.00,M,-33.9,M,,0000*6D")
>>> msg
<GGA(timestamp=datetime.time(18, 43, 53), lat='1929.045', lat_dir='S', lon='02410.506', lon_dir='E', gps_qual='1', num_sats='04', horizontal_dil='2.6', altitude=100.0, altitude_units='M', geo_sep='-33.9', geo_sep_units='M', age_gps_data='', ref_station_id='0000')>

>>> msg.timestamp
datetime.time(18, 43, 53)
>>> msg.lat
'1929.045'
>>> msg.lat_dir
'S'
>>> msg.lon
'02410.506'
>>> msg.lon_dir
'E'

这边的经纬度属性以python浮点数 (DD, “decimal degrees”)存在而不是DDDMM.MMMM (“Degrees, minutes, seconds”)，但latitude_minutes, latitude_seconds, longitude_minutes, 和 longitude_seconds这些属性能轻松创建不同格式的位置字符串 

1
2
3
4
5
6
7
8
>>> msg.latitude
-19.4840833333
>>> msg.longitude
24.1751
>>> '%02d°%07.4f′' % (msg.latitude, msg.latitude_minutes)
'-19°29.0450′'
>>> '%02d°%02d′%07.4f″' % (msg.latitude, msg.latitude_minutes, msg.latitude_seconds)
"-19°29′02.7000″"

4. 地理坐标系

地理坐标系又可分为 参心坐标系 和 地心坐标系，常见的参心坐标系北京54、西安80，常见的地心坐标系有WGS84、GCJ-02、BD-09、GCS2000，这篇文章就简单讲讲地心坐标系。

4.1 地心坐标系分类

通过上面的Python的GPGGA信息解析转换为经纬度之后，我想验证设备的精确度，当我把经纬度信息放在Google earth的时候，位置很准确，但当我把它放到Google map进行查询时，却发现位置有较大偏差。当时我就纳闷极了，同样是Google怎么位置还不一样，后来研究了一番，发现我进入的Google map是大陆的，使用的坐标系不一样，需要进行转换。

主流地图在各个地区使用的坐标系：

地图 大陆/港/澳 台湾省 海外
高德 GCJ-02 WGS84 WGS84
Google GCJ-02 WGS84 WGS84
百度 BD-09 / GCJ-02 BD-09 / GCJ-02 WGS84

而我们通过GPS拿到的GPGGA信息在解析到的经纬度是属于WGS84坐标系。WGS84是为 GPS 全球定位系统建立的坐标系统，是世界上第一个统一的地心坐标系，因此也被称为大地坐标系、原始坐标系。一般通过GPS记录仪记录下来的经纬度，就是基于WGS84坐标系的数据。

4.2 坐标系转换

国测局规定：互联网地图在国内必须至少使用 GCJ02 进行首次加密，不允许直接使用 WGS84 坐标下的地理数据，同时任何坐标系均不可转换为 WGS84 坐标。因此不存在将 GCJ-02 坐标转换为 WGS84 坐标的官方转换方法。

各大地图开发平台其实都有其官方的坐标系转换的接口，我这就说说Python的，直接使用coord-convert库即可

1
2
3
4
5
from coord_convert.transform import wgs2gcj, wgs2bd, gcj2wgs, gcj2bd, bd2wgs, bd2gcj 
lon, lat = 120, 40
gcj_lon, gcj_lat = wgs2gcj(lon, lat)
bd_lon, bd_lat = wgs2bd(lon, lat)
print(gcj_lon, gcj_lat) # the result should be: 120.00567568355486 40.0013047896019

相关函数说明如下：

1
2
3
4
5
6
wgs2gcj : convert WGS-84 to GCJ-02
wgs2bd  : convert WGS-84 to DB-09
gcj2wgs : convert GCJ-02 to WGS-84
gcj2bd  : convert GCJ-02 to BD-09
bd2wgs  : convert BD-09 to WGS-84
bd2gcj  : convert BD-09 to GCJ-02

5. 投影坐标系

投影坐标系在二维平面中进行定义。与地理坐标系不同，在二维空间范围内，投影坐标系的长度、角度和面积恒定。投影坐标系始终基于地理坐标系，而后者则是基于球体或旋转椭球体的。

在投影坐标系中，通过格网上的 x,y 坐标来标识位置，其原点位于格网中心。每个位置均具有两个值，这两个值是相对于该中心位置的坐标。一个指定其水平位置，另一个指定其垂直位置。这两个值称为 x 坐标和 y 坐标。采用此标记法，原点坐标是 x = 0 和 y = 0。

5.1 地图投影

无论将地球视为球体还是旋转椭球体，都必须变换其三维曲面以创建平面地图图幅。此数学变换通常称作地图投影。理解地图投影如何改变空间属性的一种简便方法就是观察光穿过地球投射到表面（称为投影曲面）上。想像一下，地球表面是透明的，其上绘有经纬网。用一张纸包裹地球。位于地心处的光会将经纬网投影到一张纸上。现在，可以展开这张纸并将其铺平。纸张上的经纬网形状与地球上的形状不同。地图投影使经纬网发生了变形。

地图投影使用数学公式将地球上的球面坐标与平面坐标关联起来。

5.2 墨卡托投影

墨卡托投影（Mercator Projection），又称麦卡托投影、正轴等角圆柱投影，是一种等角的圆柱形地图投影法。本投影法得名于法兰德斯出身的地理学家杰拉杜斯·墨卡托，他于1569年发表长202公分、宽124公分以此方式绘制的世界地图。在以此投影法绘制的地图上，经纬线于任何位置皆垂直相交，使世界地图可以绘制在一个长方形上。由于可显示任两点间的正确方位，航海用途的海图、航路图大都以此方式绘制。在该投影中线型比例尺在图中任意一点周围都保持不变，从而可以保持大陆轮廓投影后的角度和形状不变（即等角）；但墨卡托投影会使面积产生变形，极点的比例甚至达到了无穷大。

5.3 墨卡托投影范围

由于墨卡托投影在两极附近是趋于无限值得，因此它并没完整展现了整个世界，地图上最高纬度是85.05度。为了简化计算，采用球形映射，而不是椭球体形状。虽然采用墨卡托投影只是为了方便展示地图，需要知道的是，这种映射会给Y轴方向带来0.33%的误差。

墨卡托投影范围

X轴：由于赤道半径为6378137米，则赤道周长为2PIr = 20037508.3427892，因此X轴的取值范围：[-20037508.3427892,20037508.3427892]。

Y轴：当纬度φ接近两极，即90°时，Y值趋向于无穷。因此通常把Y轴的取值范围也限定在[-20037508.3427892,20037508.3427892]之间。

范围：因此在墨卡托投影坐标系（米）下的坐标范围是：最小为(-20037508.3427892, -20037508.3427892 )到最大坐标为(20037508.3427892, 20037508.3427892)。

墨卡托投影网格

5.4 经纬度与墨卡托投影转换

在Python中我们可以使用pyproj库来实现转换。pyproj是一个地图投影和坐标变换库，它可以将经度、纬度转换为原生地图投影的x,y坐标，反之亦然。

网上找的好多都是已经被废弃的函数，这边说下最新的，基于3.0版本的,从GPGGA解析到的经纬度转换为墨卡托平面坐标：

1
2
3
4
5
6
7
8
9
10
11
from pyproj import Transformer
#如果希望轴的顺序总是按 x、y 或 lon、lat 的顺序排列，你可以在创建转换器时使用 always_xy 选项。
# 转换为伪墨卡托坐标，，也叫web墨卡托坐标，用于谷歌地图等
transformer1 = Transformer.from_crs("EPSG:4326", "EPSG:3857", always_xy=True)
# 转换成中国地区的通用横轴墨卡托坐标(UTM投影)，用于计算等
transformer2 = Transformer.from_crs(4326, 32650, always_xy=True)

transformer1.transform(120.7, 31.3)
(13436262.538748119, 3671771.4487971356)
transformer1.transform(120.7, 31.3)
(852227.2443889615, 3468763.8018454057)

在国际上，每个坐标系统都会被分配一个 EPSG 代码，EPSG:4326 就是 WGS84 的代码。GPS是基于WGS84的

相关资源：

文章分类
Python GPS

本博客所有文章除特别声明外，均采用 CC BY-SA 3.0协议 。转载请注明出处！