一、前言
搞地图和自动驾驶的都知道,坐标转换是非常频繁的事情,有时候需要在各种坐标之间来回的转换,最近使用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