INS-D 惯性导航 坐标系转换，由地理坐标系转换到投影坐标系 Pyproj

INS-D 组合导航系统输出的地理信息：经纬度（WGS84），ins-d 默认输出的是椭球高度，如果要输出海拔高的话，需要在配置软件里面，在 ins gui 里面 device option-GNSS receiver，将该选项改为MSL。

UTM(UNIVERSAL TRANSVERSE MERCARTOR GRID SYSTEM,通用横墨卡托格网系统)坐标是一种平面直角坐标,在UTM系统中，北纬84度和南纬80度之间的地球表面积按经度6度划分为南北纵带(投影带)。从180度经线开始向东将这些投影带编号，从1编至60(北京处于第50带)。每个带再划分为纬差8度的四边形。四边形的横行从南纬80度开始。用字母C至X(不含I和O)依次标记(第X行包括北半球从北纬72度至84度全部陆地面积，共12度)每个四边形用数字和字母组合标记。参考格网向右向上读取。

1 安装python包

在进行坐标转换之前，先简单的介绍下几个专业概念：
地理坐标系：地理坐标是用纬度、经度表示地面点位置的球面坐标
投影坐标系：简单点来讲就是通过投影的方式把地球地理坐标系这种三维空间坐标展开成平面坐标。

wgs84 坐标常用于我们身边的各种gps模块。从前面的文章（点这里查看文章）我们可以从gps模块中得到wgs坐标下的经纬度值，然而经纬度值有时候并不是以距离为单位的坐标，有时候更多的是用到以距离为单位的二维坐标，所以就需要从wgs84（三维坐标）投影转换为二维坐标。
python有一个模块叫做pyproj可以帮助我们完成坐标的转换
安装

pip install pyproj

Ubuntu或者Raspbian 也可以

sudo apt-get install python3-pyproj

2 找到EPSG码

https://epsg.io/?q=UTM+china


注意：在一个经度带内，以中央经线经度为变换结果的x轴的0坐标，然后为了使这个带的结果都为大于0的数，转换结果又加了500km。（如果不想要这个500km，减去就是了）
赤道纵坐标的起点。
单位m

其中要找的主要UTM投影带号
北半球地区，选择最后字母为“N”的带，
带数=（经度整数位/6）的整数部分+31
如：上海市经度范围121左右，带数=121/6+31=51，选51N（N表示北半球），即WGS 1984 UTM ZONE 51N，对应的EPSG代码 为 32651

3 进行转换

3.1 初始化坐标参考系统（CRS）

pyproj库提供了四种初始化参考坐标的方式，括号里面的数字是地图坐标系统在epsg中的代号，每一种地图坐标系统在epsg中都有一个唯一的数值与之对应，例如，4326对应的就是WGS84地理坐标系统，对应代码坐标系统的详细信息可以在https://epsg.io/查询。

from  pyproj  import  CRS
crs=CRS.from_epsg(4326)
crs = CRS.from_string("epsg:4326")
crs = CRS.from_proj4("+proj=latlon")
crs = CRS.from_user_input(4326)

3.2 检查坐标参考系统信息

直接输入刚才构造的crs对象，即可查看当前坐标系统的详细信息，在进行坐标转换之前，首先要检查坐标转换的范围和轴向信息，比如这里的wgs84地理坐标系的轴向信息是（纬度，经度），在后面输入经纬度坐标的时候切记是先纬度，后经度的顺序，不然整个坐标系统就完全乱了，转换出来的坐标可能会出现无穷大的错误值。

3.3 转换坐标

注意：在一个经度带内，以中央经线经度为变换结果的y轴的0坐标，
在转换坐标前，首先构造一个转换对象，确定是从哪个地理坐标向哪个投影坐标进行转换，
切记经纬度的输入顺序，输出的即是UTM坐标系统的xy系统坐标

from pyproj import CRS
from pyproj import Transformer
crs = CRS.from_epsg(4326)
crs_sh = CRS.from_epsg(32651) 
transformer_WGS84toUTM = Transformer.from_crs(crs, crs_sh)   
transformer_WGS84toUTM.transform(31.030848,121.443328)  
# crs是wgs84坐标参考系统
# crs_sh是上海市所在地区UTM投影坐标系（这里代号为51N带）

输出： (x: 351429.8875922263, y: 3434061.2439915854)
其中，x在地图上看水平向右，实际方向为正东方向；
y在地图上看垂直向上，实际方向为正北方向。
在地图上给x加上20,定位点向右（正东）移动了20m，在地图上给y加上20,定位点向上（正北）移动了20m

和在线计算得到的结果保持一致
https://epsg.io/transform#s_srs=4326&t_srs=32651&x=121.4433281&y=31.0308481