局部（x,y）坐标 转 WGS84经纬度坐标

局部坐标转WGS84坐标

背景：

最近在做一个小项目，其中需要把局部坐标转换成wgs84坐标。我的局部地图是用激光雷达扫出来的一张图，就是下面这张图

用激光雷达建好图以后，就知道图上每一个点的局部坐标了，局部坐标是（x,y,z）的型式，我只用了x和y，没有用z。我们的目标就是通过这张图上的每个（x,y）坐标得到它实际上在地球上的经纬度。我没有找到局部xy坐标直接转经纬度的方法，因为局部坐标是（x，y）型式的，WGS84是经纬度，直接转换有些困难，所以我们先将局部xy坐标转换成UTM坐标，因为UTM坐标也是xy型式的，这样一来就方便多了，只要计算出两种xy坐标之间的转换矩阵就行了，而且UTM转WGS84有很多种方法。

前期准备：

1.记录局部地图上的几个坐标点（x,y）
2.测量你记录坐标点的UTM坐标和经纬度
要想得到准确的转换公式或者是转换程序，需要先拿几个点做一下试验。首先找到图片上的几个点的（x,y）坐标，再用打点仪器（测经纬度的仪器，比如千寻的打点器）测一下图片上对应点的UTM坐标和经纬度，是可以直接测出来的，最好多测几组数据，得出来的转换矩阵会比较准确。

局部xy和UTM xy转换矩阵的计算

局部坐标和UTM坐标转换是需要一个转换矩阵的，有了一些局部坐标点，有了这些局部坐标点的UTM坐标点，直接开始计算。程序使用python写的

import cv2
import numpy as np

if __name__ == '__main__' :
    points1 = np.array([[588682.56, 4074258.76], [588680.27, 4074255.52], [588682.30, 4074249.22], [588685.34, 4074248.69],[588689.00, 4074248.18], [588693.20, 4074247.44],[588696.24,4074246.92],[588698.66,4074245.95],[588700.08,4074246.96],[588701.67,4074248.64],[588702.09,4074250.75],[588700.74,4074251.63],[588698.68,4074252.49],[588695.91,4074253.46],[588693.76,4074253.99],[588691.17,4074254.52],[588688.84,4074255.05],[588684.80,4074249.13]])

    points2 = np.array([[13.10, 4.71],[16.82,3.15],[22.90,6.21],[22.59,9.27],[22.16,12.91],[21.63,17.08],[21.27,20.19],[20.55,23.48],[19.18,24.58],[17.13,25.58],[15.31,25.37],[14.43,24.33],[14.05,22.44],[13.75,19.32],[14.06,16.83],[14.37,14.46],[14.60,12.00],[15.15,7.69]])

    h1, status1 = cv2.findHomography(points1,points2)
    h2, status2 = cv2.findHomography(points2, points1)

    print(h1)

    print(h2)

解释一下：points1是很多个点的UTM坐标，points2是很多个点的局部坐标。h1是UTM转局部坐标的矩阵，h2是局部转UTM坐标的矩阵，在这里我们要用h2
计算出来是这样的：
[3.70929738e+05 7.82297943e+04 4.07430397e+06]
[5.35959224e+04 1.13062436e+04 5.88658249e+05]
[9.10430553e-02 1.92011805e-02 1.00000000e+00]

局部转UTM加UTM转WGS84程序

from pyproj import Transformer
import numpy as np

if __name__=="__main__":
    while True:
        LX = input("请输入局部坐标系X坐标：")
        lx = float(LX)
        LY = input("请输入局部坐标系Y坐标：")
        ly = float(LY)
        H = np.array([[3.70929738e+05,7.82297943e+04,4.07430397e+06],[5.35959224e+04,1.13062436e+04 ,5.88658249e+05],[9.10430553e-02 ,1.92011805e-02 ,1.00000000e+00]])  #H矩阵

        m = np.array([lx,ly,1])
        M = m.T  #m矩阵转置
        F = np.dot(H,M)   #矩阵相乘
        E = F/(F[2])    #矩阵除以第三行元素，将第三行元素变为1
        x = E[0]    #UTM的x坐标
        y = E[1]    #UTM的y坐标
        print(x,y)  #局部坐标转UTM坐标程序结束
    #################################################################################################################
        transformer = Transformer.from_crs(32650, 4326)
        transformer = Transformer.from_crs("EPSG:32650", "EPSG:4326")   #32650为华北地区UTM代码，4326为WGS84坐标系
        north, east = transformer.transform(y, x)   #对UTM坐标进行转换
        North = north+0.000021 #小数部分为误差补偿，没有误差就去掉
        East = east+0.00008  #小数部分为误差补偿，没有误差就去掉
        print(North, East)  #UTM坐标转WGS坐标完成，程序结束

这样一来局部坐标就可以转成经纬度了。需要注意的是，UTM的代码，华北地区UTM代码是50。还有非常重要的一点，就是程序倒数第二和第三行，加误差补偿那里。我这里加误差补偿的原因是因为自己测的局部坐标不是很精确，转成经纬度总是有一些误差，所以干脆手动添加误差补偿。一般不会有问题。

误差补偿

转换过程可能会存在一些误差，我存在误差的原因是因为局部坐标不精确，导致h矩阵并不是十分精确。下面说一下如何判断用程序计算出来的经纬度和实际的经纬度有没有误差。
在前期准备阶段，我们已经记录了一些局部坐标点的精确的经纬度值，把这个值当作标准值。我们拿用程序算出来的经纬度值与标准值做对比，加或减一个误差得到标准值。将经纬度的误差分别写到程序的倒数二三行就可以。一般打点的时候规范操作就不用加误差补偿。

总结

打点的时候要规范，最好多打几个点。最后分享一个输入经纬度显示定位的网址：
https://www.latlong.net/lat-long-utm.html
各位大佬可以再研究一下加上z转高度的部分。
有问题可以留言或者私信我。