在生产中,会遇到将我们既有的坐标系转换到高德或者百度坐标系下,从而可以使用百度高德的基础地图。亦或者将获得的百度高德数据转换到普通的84坐标系下。本文基于Arcpy脚本实现了对shp或者gdb数据的坐标系之间的相互转换。
1、从常用坐标系到百度高德坐标系的转换
思路如下:
1)读取原始数据的坐标系;
2)利用Arcpy将原始坐标转换成WGS84经纬度;
3)对于线和面,利用Arcpy对数据进行修复;
4)逐个取点,将WGS84经纬度转换成高德或者百度坐标系。
# -*- coding:UTF-8 -*-
import arcpy, os
import ConvertLag_Lat as CVLL
import sys
reload(sys)
sys.setdefaultencoding('utf8')
#输入工作空间
in_features = arcpy.GetParameterAsText(0)
#输出的空间坐标系
outCSStr=arcpy.GetParameterAsText(1)
#输出工作空间
output_folder = arcpy.GetParameterAsText(2)
outCS=arcpy.SpatialReference(4326)
arcpy.BatchProject_management(in_features, output_folder, outCS)
arcpy.env.workspace = output_folder
transType=0
if(outCSStr=='wgs84'):
arcpy.AddMessage("transform succeed")
elif(outCSStr=='gaode'):
arcpy.AddMessage("target coords gaode")
transType='3'
elif(outCSStr=='baidu'):
transType='6'
arcpy.AddMessage("target coords baidu")
arcpy.AddMessage("start transforming......")
featureclasses = arcpy.ListFeatureClasses()
for fc in featureclasses:
cur= arcpy.UpdateCursor(fc)
SR= arcpy.Describe(fc).spatialReference
des = arcpy.Describe(fc)
arcpy.RepairGeometry_management(fc)
fc_count=int(arcpy.GetCount_management(fc).getOutput(0))
arcpy.AddMessage(fc_count)
arcpy.SetProgressor("step", "transform "+fc,0, fc_count, 1)
if des.shapeType.upper()=='POINT':
for r in cur:
geom = r.getValue("SHAPE")
r.setValue("SHAPE",CVLL.offsetPoint(geom,transType))
cur.updateRow(r)
arcpy.SetProgressorPosition()
del r,cur
elif des.shapeType.upper()=='POLYLINE':
for r in cur:
geom = r.getValue("SHAPE")
r.setValue("SHAPE",CVLL.offsetLine(geom,SR,transType))
cur.updateRow(r)
arcpy.SetProgressorPosition()
del r,cur
elif des.shapeType.upper()=='POLYGON':
for r in cur:
geom = r.getValue("SHAPE")
r.setValue("SHAPE",CVLL.offsetPolygon(geom,SR,transType))
cur.updateRow(r)
arcpy.SetProgressorPosition()
del r,cur
arcpy.AddMessage(fc+" transform finished")
2、百度高德与WGS84互转
import arcpy, os
import ConvertLag_Lat as CVLL
import sys
reload(sys)
sys.setdefaultencoding('utf8')
#input path of folder where all shapefiles reside in
in_fold = arcpy.GetParameterAsText(0)
output_folder = arcpy.GetParameterAsText(1)
outCSStr= arcpy.GetParameterAsText(2)
outCS=arcpy.SpatialReference(4326)
arcpy.env.workspace = in_fold
#copy data
featureclasses = arcpy.ListFeatureClasses()
for fc in featureclasses:
arcpy.Copy_management(fc,output_folder+"\\"+fc)
arcpy.env.workspace=output_folder
transType='0'
if(outCSStr=='mars2baidu'):
transType='1'
elif(outCSStr=='baidu2mars'):
transType='2'
elif(outCSStr=='mars2baidu'):
transType='3'
elif(outCSStr=='mars2wgs84'):
transType='4'
elif(outCSStr=='baidu2wgs'):
transType='5'
elif(outCSStr=='wgs2baidu'):
transType='6'
arcpy.AddMessage("transfor type is "+str(transType))
featureclasses = arcpy.ListFeatureClasses()
for fc in featureclasses:
cur= arcpy.UpdateCursor(fc)
SR= arcpy.SpatialReference(4326)
des = arcpy.Describe(fc)
arcpy.RepairGeometry_management(fc)
fc_count=int(arcpy.GetCount_management(fc).getOutput(0))
arcpy.SetProgressor("step", "transform "+fc,0, fc_count, 1)
if des.shapeType.upper()=='POINT':
for r in cur:
geom = r.getValue("SHAPE")
r.setValue("SHAPE",CVLL.offsetPoint(geom,transType))
cur.updateRow(r)
arcpy.SetProgressorPosition()
del r,cur
elif des.shapeType.upper()=='POLYLINE':
for r in cur:
geom = r.getValue("SHAPE")
r.setValue("SHAPE",CVLL.offsetLine(geom,SR,transType))
cur.updateRow(r)
arcpy.SetProgressorPosition()
del r,cur
elif des.shapeType.upper()=='POLYGON':
for r in cur:
geom = r.getValue("SHAPE")
r.setValue("SHAPE",CVLL.offsetPolygon(geom,SR,transType))
cur.updateRow(r)
arcpy.SetProgressorPosition()
del r,cur
arcpy.AddMessage(fc+" transform finished")
不足
这里有几个方面的不足
1)我使用的是git上给出的坐标转换方式,经测试转换后有米级的误差;
2)这里用的是取点在重组的方式,因此对于弧,多部件要素的支持不行,只能支持对简单的点线面进行转换。
这里给出源码,希望有需要的人可以继续优化完善。