【转载】osgeo和pyproj：经纬度坐标和高斯坐标互相转换

一、前言
搞地图和自动驾驶的都知道，坐标转换是非常频繁的事情，有时候需要在各种坐标之间来回的转换，最近使用python代码处理地图数据，在使用osgeo库中的gdal时，发现了gdal v2和V3的一些不同之处，研究了一下，这里分享出来。

二、问题描述
经纬度转高斯的过程中，发现3.0一直出现的转换结果是’inf’，经过查看官方github上的issue，才知道，gdal V3.0以后，转换需要设置转换策略，具体看后面代码中的说明，现象截图如下：

三、解决后封装的代码
下面代码使用osgeo库和pyproj库分别实现，底层调用都是gdal，代码经过测试。

#!/usr/bin/python3

__author__ = 'ISmileLi'

from osgeo import gdal, ogr, osr
from pyproj import Transformer

'''
osgeo底层坐标转换使用的库还是proj,下面函数中的espg值需要根据自己的需求进行修改，
下文测试使用的是wgs84与中国区高斯-克吕格EPSG码为21460区的转换
'''

def lonLat_to_gauss(lon, lat, from_epsg=4326, to_epsg=21460):
    '''
    经纬度转高斯
    :param lon:
    :param lat:
    :param from_epsg:
    :param to_EPSG:
    :return:
    '''

    from_spa = osr.SpatialReference()
    '''
    gdal版本大于3.0以后必须设置转换策略才能正确显示结果，否则结果将会输出'inf'
    可以了解官方的这个issue说明：https://github.com/OSGeo/gdal/issues/1546
    '''
    if int(gdal.__version__[0]) >= 3:
        from_spa.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
    from_spa.ImportFromEPSG(from_epsg)
    to_spa = osr.SpatialReference()
    to_spa.ImportFromEPSG(to_epsg)
    coord_trans = osr.CoordinateTransformation(from_spa, to_spa)

    t = coord_trans.TransformPoint(lon, lat)
    return t[0], t[1]

def gauss_to_lonLat(x, y, from_epsg=21460, to_epsg=4326):
    '''
    高斯转经纬度
    :param x:
    :param y:
    :param from_epsg:
    :param to_EPSG:
    :return:
    '''

    from_spa = osr.SpatialReference()
    #if int(gdal.__version__[0]) >= 3:
        #from_spa.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
    from_spa.ImportFromEPSG(from_epsg)
    to_spa = osr.SpatialReference()
    to_spa.ImportFromEPSG(to_epsg)
    coord_trans = osr.CoordinateTransformation(from_spa, to_spa)

    t = coord_trans.TransformPoint(x, y)
    return t[0], t[1]


def lonLat_to_gauss_proj(lon, lat, from_epsg="EPSG:4326", to_epsg="EPSG:21460"):
    '''
    使用proj库经纬度转高斯
    :param lon:
    :param lat:
    :param from_epsg:
    :param to_epsg:
    :return:
    '''
    transfromer = Transformer.from_crs(from_epsg, to_epsg,always_xy=True)  # WGS-84对应码->EPSG:4326, 中国高斯对应码：EPSG:21460
    x, y = transfromer.transform(lon, lat)
    print('lonLat_to_gauss_proj x, y:',x, y)
    return x, y

def gauss_to_lonLat_proj(x, y, from_epsg="EPSG:21460", to_epsg="EPSG:4326"):
    '''
    使用proj库高斯转经纬度
    :param x:
    :param y:
    :param from_epsg:
    :param to_epsg:
    :return:
    '''
    transfromer = Transformer.from_crs(from_epsg, to_epsg, always_xy=True)  # WGS-84对应码->EPSG:4326, 中国高斯对应码：EPSG:21460
    lon, lat = transfromer.transform(x, y)
    print('lonLat_to_gauss_proj lon, lat:', lon, lat)
    return lon, lat

if __name__ == '__main__':
    lon = 116.2446370442708300
    lat = 40.0670713975694400
    x, y = lonLat_to_gauss(lon, lat)
    print('x, y: ', x, y)
    lat_t, lon_t = gauss_to_lonLat(x, y)
    print('lon_t, lat_t: ', lon_t, lat_t)

    '''
    这里要注意pyproj的转换会交换x/y返回，可以对比osgeo使用打印结果看出来，
    详细了解可以参考官网文档：https://pyproj4.github.io/pyproj/stable/api/transformer.html
    '''
    lon_t = 116.2446370442708300
    lat_t = 40.0670713975694400
    x_t, y_t = lonLat_to_gauss_proj(lon_t, lat_t)
    gauss_to_lonLat_proj(x_t, y_t)

运行截图：

 四、下载地址：
你也可以直接从github上下载文中代码：https://github.com/toby-King/MapFormat

————————————————
版权声明：本文为CSDN博主「ISmileLi」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/toby54king/article/details/117599464