互联网常见坐标之间的转换（Python）
本文共 8185 字，大约阅读时间需要 27 分钟。
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/wandergis/coordtransform
代码实现
开发工具：VScode 使用语言：Python
只实现经纬度坐标之间的转换，支持点、线、面。Web墨卡托的转换目前没设计到没有实现
# -*- coding: utf-8 -*-import arcpyimport math# define ellipsoidx_pi = 3.14159265358979324 * 3000.0 / 180.0pi = 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 retdef _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 retdef 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 newPointdef 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 newLinedef 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.spatialReferenceif des.shapeType.upper() == 'POINT': for r in cur: geom = r.getValue("SHAPE") r.setValue("SHAPE", TransPoint(geom, transType)) cur.updateRow(r) del r, curelif 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, curelif des.shapeType.upper() == 'POLYGON': for r in cur: geom = r.getValue("SHAPE") r.setValue("SHAPE", TransPolygon(geom, transType)) cur.updateRow(r) del r, cur 效果验证
转载地址：http://zotqj.baihongyu.com/
你可能感兴趣的文章
Cannot find class [org.apache.commons.dbcp.BasicDataSource（亲测有效）
查看>>
小白学习[leetcode]之402移掉k位数字
查看>>
DNN和RNN和DNN之间的区别
查看>>
小白学习[leetcode]之215数组中的第K个最大元素
查看>>
Java 将PDF 转为Word、图片、SVG、XPS、Html、PDF/A(亲测有效）
查看>>
小白学习之程序员代码面试指南之单调栈结构
查看>>
小白之java虚拟机工作学习面试准备（持续更新）
查看>>
小白学习[leetcode]之455分发饼干 （贪心算法）
查看>>
小白学习[leetcode]之135分发糖果 （贪心算法）
查看>>
小白读论文之. Show and Tell Lessons learned from the 2015 MSCOCO Image Captioning Challenge
查看>>
小白学习[leetcode]之435无重叠区间（贪心算法）
查看>>
java二维数组使用lambda表达式进行排序
查看>>
java图解从编译到运行（详细）
查看>>
为什么java的局部变量要初始化而全局变量不用
查看>>
校园一卡通的实现机制(图解）
查看>>
依然是springMVC+mysql的用户登陆小项目（全注释图解超详细）
查看>>
小白学习[leetcode]之605种花问题（贪心算法）
查看>>
解决Intellij IDEA 构建Maven没有无java及resources等文件
查看>>
Unknown system variable ‘query_cache_size‘
查看>>
小白学习[leetcode]之452用最少数量的箭引爆气球（贪心算法）
查看>>