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

本文介绍了在使用GDAL库进行经纬度到高斯坐标的转换时遇到的'inf'结果问题,该问题在GDAL V3.0以后出现。作者通过研究官方GitHub上的issue找到了解决方案,即设置转换策略。并提供了使用osgeo和pyproj库的封装代码,以解决坐标转换问题。代码经过测试,可正常工作。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

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

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

OSGeoEarth是基于Google Earth API的一个开源库,它允许用户在地图上显示地理信息。经纬度坐标(经度longitude纬度latitude)通常表示地球表面某一点在球面上的位置,而笛卡尔坐标(Cartesian coordinates)则是在二维平面上的直角坐标系统。 将经纬度转换为笛卡尔坐标的过程通常是通过三维空间中的大地投影(Geographic Projection)完成的,这个过程涉及到数学计算。常见的大地投影如Mercator、UTM等。以下是一个简单的例子,展示了如何使用EPSG4326(WGS84,通常用于Web地图)到Web Mercator(也叫Spherical Mercator)的转换: ```java import org.osgeo.proj4j.CoordinateReferenceSystem; import org.osgeo.proj4j.CoordinateTransform; import org.osgeo.proj4j.ProjCoordinate; import org.osgeo.proj4j.Transformer; public class LonLatToCartesian { private static final String FROM_EPSG = "EPSG:4326"; // WGS84 private static final String TO_WEB_MERCATOR = "EPSG:3857"; // Web Mercator public static void main(String[] args) throws Exception { Transformer transformer = TransformerFactory.createTransformer(FROM_EPSG, TO_WEB_MERCATOR); double longitude = 120.123; // 经度 double latitude = 31.234; // 纬度 ProjCoordinate fromCoord = new ProjCoordinate(longitude, latitude); ProjCoordinate toCoord = new ProjCoordinate(); transformer.transform(fromCoord, toCoord); double x = toCoord.x; // 笛卡尔坐标中的x轴值 double y = toCoord.y; // 笛卡尔坐标中的y轴值 System.out.println("经纬度:" + longitude + ", " + latitude + " 转换为笛卡尔坐标:" + x + ", " + y); } } ``` 在这个示例中,我们创建了一个`Transformer`对象,然后将经纬度传递给`transform`方法,得到的就是在Web Mercator下的笛卡尔坐标。注意,这只是一个基本示例,实际应用中可能需要处理精度问题其他细节。
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值