互联网常见坐标之间的转换（Python）

最新推荐文章于 2024-08-08 07:08:17 发布

置顶遥哥_GB_Fore

最新推荐文章于 2024-08-08 07:08:17 发布

阅读量1.1k

点赞数

分类专栏：地理信息

本文链接：https://blog.csdn.net/mymhj/article/details/117027453

版权

地理信息专栏收录该内容

9 篇文章 5 订阅

订阅专栏

1.场景描述

随着互联网的兴趣，在企业应用系统的开发过程中，几乎每一个APP（CS、BS、移动端）都会搜集用户或设备的位置数据，然后在相关地图（来自互联网电子地图或自建地图服务）上进行标注展示，如果对空间坐标系不是很了解的话，不管你Coding能力有多强都会被各种位置匹配问题（错误、偏离、飞出等）折磨的晕头转向，不知所措。首先针对几种常见的互联网坐标我们对其做一简单场景梳理

常见的互联网坐标：

地球坐标：WGS84，国际标准、从GPS设备获取的坐标的参考坐标系、国际地图提供商使用的坐标系。OSM、谷歌外国、ArcGISonline

火星坐标：GCJ-02，中国标准、也叫国测局坐标，从国行移动设备中定位获取的坐标的参考坐标系。国内出版的各种地图系统，必须至少采用GCJ-02对地理位置进行首次加密。高德地图、搜搜地图、阿里云地图、腾讯地图。

百度坐标：DB09，百度标准。百度地图、百度空间数据的参考坐标系都使用此坐标系。是在GCJ-02坐标系的基础上做了二次加密

备注：天地图使用的是CGC2000，参考椭球非常接近。当前阶段（星历变化可能导致差异变大），变率差异引的椭球面上的维度和高程变化量不大0.1mm,两者相容误差在cm级别，在互联网应用中可以将CGC2000就当作WGS84来用。

针对上述互联网坐标系，我将转换场景梳理为以下六种：

WGS84转GCJ02：WGS84_To_GCJ02

GCJ02转WGS84：GCJ02_To_WGS84

GCJ02转BD09：GCJ02_To_BD09

BD09转GCJ02：BD09_To_GCJ02

WGS84转BD09：WGS84_To_BD09

BD09转WGS84：BD09_To_WGS84

2.功能实现

参考资料

代码实现

开发工具：VScode 使用语言：Python

只实现经纬度坐标之间的转换，支持点、线、面。Web墨卡托的转换目前没设计到没有实现

# -*- coding: utf-8 -*-
import arcpy
import math


# define ellipsoid
x_pi = 3.14159265358979324 * 3000.0 / 180.0
pi = 3.1415926535897932384626  # π
a = 6378245.0  # 长半轴
ee = 0.00669342162296594323  # 偏心率平方


def _transformlat(lng, lat):
    ret = -100.0 + 2.0 * lng + 3.0 * lat + 0.2 * lat * lat + \
          0.1 * lng * lat + 0.2 * math.sqrt(math.fabs(lng))
    ret += (20.0 * math.sin(6.0 * lng * pi) + 20.0 *
            math.sin(2.0 * lng * pi)) * 2.0 / 3.0
    ret += (20.0 * math.sin(lat * pi) + 40.0 *
            math.sin(lat / 3.0 * pi)) * 2.0 / 3.0
    ret += (160.0 * math.sin(lat / 12.0 * pi) + 320 *
            math.sin(lat * pi / 30.0)) * 2.0 / 3.0
    return ret


def _transformlng(lng, lat):
    ret = 300.0 + lng + 2.0 * lat + 0.1 * lng * lng + \
          0.1 * lng * lat + 0.1 * math.sqrt(math.fabs(lng))
    ret += (20.0 * math.sin(6.0 * lng * pi) + 20.0 *
            math.sin(2.0 * lng * pi)) * 2.0 / 3.0
    ret += (20.0 * math.sin(lng * pi) + 40.0 *
            math.sin(lng / 3.0 * pi)) * 2.0 / 3.0
    ret += (150.0 * math.sin(lng / 12.0 * pi) + 300.0 *
            math.sin(lng / 30.0 * pi)) * 2.0 / 3.0
    return ret


def out_of_china(lng, lat):
    """
    判断是否在国内，不在国内不做偏移
    :param lng:
    :param lat:
    :return:
    """
    return not (lng > 73.66 and lng < 135.05 and lat > 3.86 and lat < 53.55)


def WGS84_To_GCJ02(lng, lat):
    """
    WGS84转GCJ02(火星坐标系)
    :param lng:WGS84坐标系的经度
    :param lat:WGS84坐标系的纬度
    :return:
    """
    if out_of_china(lng, lat):  # 判断是否在国内
        return [lng, lat]
    dlat = _transformlat(lng - 105.0, lat - 35.0)
    dlng = _transformlng(lng - 105.0, lat - 35.0)
    radlat = lat / 180.0 * pi
    magic = math.sin(radlat)
    magic = 1 - ee * magic * magic
    sqrtmagic = math.sqrt(magic)
    dlat = (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * pi)
    dlng = (dlng * 180.0) / (a / sqrtmagic * math.cos(radlat) * pi)
    mglat = lat + dlat
    mglng = lng + dlng
    return arcpy.Point(mglng, mglat)


def GCJ02_To_WGS84(lng, lat):
    """
    GCJ02(火星坐标系)转GPS84
    :param lng:火星坐标系的经度
    :param lat:火星坐标系纬度
    :return:
    """
    if out_of_china(lng, lat):
        return [lng, lat]
    dlat = _transformlat(lng - 105.0, lat - 35.0)
    dlng = _transformlng(lng - 105.0, lat - 35.0)
    radlat = lat / 180.0 * pi
    magic = math.sin(radlat)
    magic = 1 - ee * magic * magic
    sqrtmagic = math.sqrt(magic)
    dlat = (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * pi)
    dlng = (dlng * 180.0) / (a / sqrtmagic * math.cos(radlat) * pi)
    mglat = lat + dlat
    mglng = lng + dlng
    return arcpy.Point(lng * 2 - mglng,lat * 2 - mglat)


def GCJ02_to_BD09(lng, lat):
    """
    火星坐标系(GCJ-02)转百度坐标系(BD-09)
    谷歌、高德——>百度
    :param lng:火星坐标经度
    :param lat:火星坐标纬度
    :return:
    """
    z = math.sqrt(lng * lng + lat * lat) + 0.00002 * math.sin(lat * x_pi)
    theta = math.atan2(lat, lng) + 0.000003 * math.cos(lng * x_pi)
    bd_lng = z * math.cos(theta) + 0.0065
    bd_lat = z * math.sin(theta) + 0.006
    return arcpy.Point(bd_lng, bd_lat)


def BD09_to_GCJ02(bd_lon, bd_lat):
    """
    百度坐标系(BD-09)转火星坐标系(GCJ-02)
    百度——>谷歌、高德
    :param bd_lat:百度坐标纬度
    :param bd_lon:百度坐标经度
    :return:转换后的坐标列表形式
    """
    x = bd_lon - 0.0065
    y = bd_lat - 0.006
    z = math.sqrt(x * x + y * y) - 0.00002 * math.sin(y * x_pi)
    theta = math.atan2(y, x) - 0.000003 * math.cos(x * x_pi)
    gg_lng = z * math.cos(theta)
    gg_lat = z * math.sin(theta)
    return arcpy.Point(gg_lng, gg_lat)


def bd09_to_wgs84(bd_lon, bd_lat):
    newPoint = arcpy.Point()
    newPoint = BD09_to_GCJ02(bd_lon, bd_lat)
    return GCJ02_To_WGS84(newPoint.X, newPoint.Y)


def WGS84_To_BD09(lon, lat):
    newPoint = arcpy.Point()
    newPoint = WGS84_To_GCJ02(lon, lat)
    return GCJ02_to_BD09(newPoint.X, newPoint.Y)


def TransPoint(geom, transType):
    newPoint = arcpy.Point()
    oldPoint = geom.getPart(0)
    if transType == 'WGS84_To_GCJ02':
        newPoint = WGS84_To_GCJ02(oldPoint.X, oldPoint.Y)
    elif transType == 'GCJ02_To_WGS84':
        newPoint = GCJ02_To_WGS84(oldPoint.X, oldPoint.Y)
    elif transType == 'GCJ02_to_BD09':
        newPoint = GCJ02_to_BD09(oldPoint.X, oldPoint.Y)
    elif transType == 'BD09_to_GCJ02':
        newPoint = BD09_to_GCJ02(oldPoint.X, oldPoint.Y)
    elif transType == 'WGS84_To_BD09':
        newPoint = WGS84_To_BD09(oldPoint.X, oldPoint.Y)
    elif transType == 'BD09_To_WGS84':
        newPoint = bd09_to_wgs84(oldPoint.X, oldPoint.Y)
    return newPoint


def TransPolyLine(geom, transType):
    newPoint = arcpy.Point()
    newArray = arcpy.Array()
    array = geom.getPart(0)
    for i in range(0, array.count):
        oldPoint = array[i]
        if transType == 'WGS84_To_GCJ02':
            newPoint = WGS84_To_GCJ02(oldPoint.X, oldPoint.Y)
        elif transType == 'GCJ02_To_WGS84':
            newPoint = GCJ02_To_WGS84(oldPoint.X, oldPoint.Y)
        elif transType == 'GCJ02_to_BD09':
            newPoint = GCJ02_to_BD09(oldPoint.X, oldPoint.Y)
        elif transType == 'BD09_to_GCJ02':
            newPoint = BD09_to_GCJ02(oldPoint.X, oldPoint.Y)
        elif transType == 'WGS84_To_BD09':
            newPoint = WGS84_To_BD09(oldPoint.X, oldPoint.Y)
        elif transType == 'BD09_To_WGS84':
            newPoint = bd09_to_wgs84(oldPoint.X, oldPoint.Y)
        newArray.add(newPoint)
    newLine = arcpy.Polyline(newArray, SR)
    newArray.removeAll()
    return newLine


def TransPolygon(geom, transType):
    newPoint = arcpy.Point()
    newArray = arcpy.Array()
    array = geom.getPart(0)
    for i in range(0, array.count):
        oldPoint = array[i]
        if transType == 'WGS84_To_GCJ02':
            newPoint = WGS84_To_GCJ02(oldPoint.X, oldPoint.Y)
        elif transType == 'GCJ02_To_WGS84':
            newPoint = GCJ02_To_WGS84(oldPoint.X, oldPoint.Y)
        elif transType == 'GCJ02_to_BD09':
            newPoint = GCJ02_to_BD09(oldPoint.X, oldPoint.Y)
        elif transType == 'BD09_to_GCJ02':
            newPoint = BD09_to_GCJ02(oldPoint.X, oldPoint.Y)
        elif transType == 'WGS84_To_BD09':
            newPoint = wgs84_to_bd09(oldPoint.X, oldPoint.Y)
        elif transType == 'BD09_To_WGS84':
            newPoint = bd09_to_wgs84(oldPoint.X, oldPoint.Y)
        newArray.add(newPoint)
    newLine = arcpy.Polyline(newArray, SR)
    newArray.removeAll()
    return newLine


# 输入工作空间
# in_Features = arcpy.GetParameter(0)
# out_Feature_FullPath = arcpy.GetParameter(1)
# transType = arcpy.GetParameter(2)
# arcpy.AddMessage(convert_Type)
transType = "WGS84_To_GCJ02"
inShpFile = "E:\\OutputData\\WGS84_R.shp"
outShpFile = "E:\\OutputData\\WGS_GCJ02_R.shp"
if arcpy.Exists(outShpFile):
    arcpy.Delete_management(outShpFile)
arcpy.CopyFeatures_management(inShpFile, outShpFile)
arcpy.RepairGeometry_management(outShpFile)
# arcpy.AddMessage(out_ShpFile)
cur = arcpy.UpdateCursor(outShpFile)
des = arcpy.Describe(outShpFile)
SR = des.spatialReference
if des.shapeType.upper() == 'POINT':
    for r in cur:
        geom = r.getValue("SHAPE")
        r.setValue("SHAPE", TransPoint(geom, transType))
        cur.updateRow(r)
    del r, cur
elif des.shapeType.upper() == 'POLYLINE':
    for r in cur:
        geom = r.getValue("SHAPE")
        newgeom = TransPolyLine(geom, transType)
        r.setValue("SHAPE", newgeom)
        cur.updateRow(r)
    del r, cur
elif des.shapeType.upper() == 'POLYGON':
    for r in cur:
        geom = r.getValue("SHAPE")
        r.setValue("SHAPE", TransPolygon(geom, transType))
        cur.updateRow(r)
    del r, cur

效果验证