本文共 8185 字,大约阅读时间需要 27 分钟。
随着互联网的兴趣,在企业应用系统的开发过程中,几乎每一个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
参考资料
代码实现
开发工具: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/