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.功能实现
参考资料
- https://github.com/hujiulong/gcoord
- https://github.com/wandergis/coordtransform
代码实现
开发工具：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效果验证
本文章为转载内容，我们尊重原作者对文章享有的著作权。如有内容错误或侵权问题，欢迎原作者联系我们进行内容更正或删除文章。
