几种坐标系转换简述

港控/mmm° 2022-10-23 08:29 224阅读 0赞

# 1.大地主题解算 #

## 1.1 大地主题正解 ##

由点A的经纬度,以及与B的方位角和相对距离,计算点B的经纬度,见demo.

## 1.2 大地主题反解 ##

由点A,点B的经纬度,计算其前向方向角,后向方向角以及相对距离,见demo.

## 1.3 大地正反解demo ##

import pyproj
    geod = pyproj.Geod(ellps="WGS84")
    # ======= 大地主题正解 =======
    # 由lon, lats, az, dist 计算另一个点的经纬度 (azimuth:方位角 0度正北,180度正南)
    lon, lat, az = geod.fwd(120, 0, 0, 100000)
    print(lon, lat, az)
    # 输出：(120.0, -0.9043687229127633, 0.0)
    
    # ======= 大地主题反解 =======
    # 由lon1,lats1,lon2,lats2 计算距离和方位
    forward_az, back_az, dist = geod.inv(120, 0, lon, lat)
    print(forward_az, back_az, dist)
    # 输出：(180.0, 0.0, 99999.99999999999)

# 2.UTM转WGS #

## 2.1 定义 ##

## 2.1.1 UTM ##

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

每一四边形划分为很多边长为1000 000米的小区，用字母组合系统标记。在每个投影带中，位于带中心的经线，赋予横坐标值为500 000米。对于北半球赤道的标记坐标值为0，对于南半球为10000000米，往南递减。大比例尺地图UTM方格主线间距离一般为1KM,因此UTM系统有时候也被称作方里格。因为UTM系统采用的是[横墨卡托投影][Link 1]，沿每一条南北格网线(带中心的一条格网线为经线)比例系数为常数，在东西方向则为变数。沿每一UTM格网的中心格网线的比例系数应为0．99960(比例尺较小)，在南北纵行最宽部分(赤道)的边缘上，包括带的重叠部分，距离中心点大约363公里，比例系数为 1.00158。----(Baidu百科)

# 2.1.2 WGS #

WGS-84坐标系（World Geodetic System一1984 Coordinate System）

一种国际上采用的[地心坐标系][Link 2]。坐标原点为[地球质心][Link 3]，其地心[空间直角坐标系][Link 4]的Z轴指向BIH （国际时间服务机构）1984.0定义的协议地球极（CTP）方向，X轴指向BIH 1984.0的零[子午面][Link 5]和CTP赤道的交点，Y轴与Z轴、X轴垂直构成[右手坐标系][Link 6]，称为1984年世界[大地坐标系统][Link 7]。----(Baidu百科)

## 2.2 demo ##

注:UTM坐标系,又常称为XY坐标系单位为m,容易引起误导.

以demo为例,其中x+1000m,但经纬度均发生了改变,故与大地主题正解并不是同一回事.

import math
    from pyproj import Proj
    
    
    def transform_utm_into_lat_lon(x, y, zone, hemisphere):
        # verify the hemisphere
        h_north = False
        h_south = False
        if (hemisphere == 'N'):
            h_north = True
        elif (hemisphere == 'S'):
            h_south = True
        else:
            print("Unknown hemisphere: " + hemisphere)
        proj_in = Proj(proj='utm', zone=zone, ellps='WGS84', south=h_south, north=h_north, errcheck=True)
        lon, lat = proj_in(x, y, inverse=True)
        # just printing the floating point number with 9 decimal points will round it
        lon = math.floor(lon * 1000000) / 1000000
        lat = math.floor(lat * 1000000) / 1000000
        lon = "%.6f" % lon
        lat = "%.6f" % lat
        return lon, lat
    
    
    if __name__ == '__main__':
        # 定义初始化UTM,取上海的UTM坐标点
        x = 326474.4650588095
        y = 3392563.1296714176
        lon, lat = transform_utm_into_lat_lon(x, y, 51, 'N')
        print(lon, lat)
        # 输出:(121.189000 30.653162)
        lon, lat = transform_utm_into_lat_lon(x+1000, y, 51, 'N')
        print(lon, lat)
        # 输出:(121.199432 30.653307)

[Link 1]: https://baike.baidu.com/item/%E6%A8%AA%E5%A2%A8%E5%8D%A1%E6%89%98%E6%8A%95%E5%BD%B1/991817
[Link 2]: https://baike.baidu.com/item/%E5%9C%B0%E5%BF%83%E5%9D%90%E6%A0%87%E7%B3%BB/5155013
[Link 3]: https://baike.baidu.com/item/%E5%9C%B0%E7%90%83%E8%B4%A8%E5%BF%83/7013007
[Link 4]: https://baike.baidu.com/item/%E7%A9%BA%E9%97%B4%E7%9B%B4%E8%A7%92%E5%9D%90%E6%A0%87%E7%B3%BB/10743618
[Link 5]: https://baike.baidu.com/item/%E5%AD%90%E5%8D%88%E9%9D%A2/9688202
[Link 6]: https://baike.baidu.com/item/%E5%8F%B3%E6%89%8B%E5%9D%90%E6%A0%87%E7%B3%BB/2562697
[Link 7]: https://baike.baidu.com/item/%E5%A4%A7%E5%9C%B0%E5%9D%90%E6%A0%87%E7%B3%BB%E7%BB%9F/6131526