一、前言
搞地圖和自動駕駛的都知道，坐標(biāo)轉(zhuǎn)換是非常頻繁的事情，有時候需要在各種坐標(biāo)之間來回的轉(zhuǎn)換，最近使用python代碼處理地圖數(shù)據(jù)，在使用osgeo庫中的gdal時，發(fā)現(xiàn)了gdal v2和V3的一些不同之處，研究了一下，這里分享出來。

二、問題描述
經(jīng)緯度轉(zhuǎn)高斯的過程中，發(fā)現(xiàn)3.0一直出現(xiàn)的轉(zhuǎn)換結(jié)果是’inf’，經(jīng)過查看官方github上的issue，才知道，gdal V3.0以后，轉(zhuǎn)換需要設(shè)置轉(zhuǎn)換策略，具體看后面代碼中的說明，現(xiàn)象截圖如下：

三、解決后封裝的代碼
下面代碼使用osgeo庫和pyproj庫分別實現(xiàn)，底層調(diào)用都是gdal，代碼經(jīng)過測試。

#!/usr/bin/python3__author__ = 'ISmileLi'from osgeo import gdal, ogr, osr from pyproj import Transformer''' osgeo底層坐標(biāo)轉(zhuǎn)換使用的庫還是proj,下面函數(shù)中的espg值需要根據(jù)自己的需求進行修改， 下文測試使用的是wgs84與中國區(qū)高斯-克呂格EPSG碼為21460區(qū)的轉(zhuǎn)換 '''def lonLat_to_gauss(lon, lat, from_epsg=4326, to_epsg=21460):'''經(jīng)緯度轉(zhuǎn)高斯:param lon::param lat::param from_epsg::param to_EPSG::return:'''from_spa = osr.SpatialReference()'''gdal版本大于3.0以后必須設(shè)置轉(zhuǎn)換策略才能正確顯示結(jié)果，否則結(jié)果將會輸出'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):'''高斯轉(zhuǎn)經(jīng)緯度: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庫經(jīng)緯度轉(zhuǎn)高斯:param lon::param lat::param from_epsg::param to_epsg::return:'''transfromer = Transformer.from_crs(from_epsg, to_epsg,always_xy=True) # WGS-84對應(yīng)碼->EPSG:4326, 中國高斯對應(yīng)碼：EPSG:21460x, y = transfromer.transform(lon, lat)print('lonLat_to_gauss_proj x, y:',x, y)return x, ydef gauss_to_lonLat_proj(x, y, from_epsg="EPSG:21460", to_epsg="EPSG:4326"):'''使用proj庫高斯轉(zhuǎn)經(jīng)緯度:param x::param y::param from_epsg::param to_epsg::return:'''transfromer = Transformer.from_crs(from_epsg, to_epsg, always_xy=True) # WGS-84對應(yīng)碼->EPSG:4326, 中國高斯對應(yīng)碼：EPSG:21460lon, lat = transfromer.transform(x, y)print('lonLat_to_gauss_proj lon, lat:', lon, lat)return lon, latif __name__ == '__main__':lon = 116.2446370442708300lat = 40.0670713975694400x, 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的轉(zhuǎn)換會交換x/y返回，可以對比osgeo使用打印結(jié)果看出來，詳細(xì)了解可以參考官網(wǎng)文檔：https://pyproj4.github.io/pyproj/stable/api/transformer.html'''lon_t = 116.2446370442708300lat_t = 40.0670713975694400x_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

————————————————
版權(quán)聲明：本文為CSDN博主「ISmileLi」的原創(chuàng)文章，遵循CC 4.0 BY-SA版權(quán)協(xié)議，轉(zhuǎn)載請附上原文出處鏈接及本聲明。
原文鏈接：https://blog.csdn.net/toby54king/article/details/117599464

總結(jié)

以上是生活随笔為你收集整理的【转载】osgeo和pyproj：经纬度坐标和高斯坐标互相转换的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。