简介面向测绘、GIS及国土规划从业者的一份技术文献针对1980西安坐标系与2000国家大地坐标系的转换难题展开研究。文中系统梳理了两套坐标系的参考椭球参数、原点定义及差异重点解析七参数转换模型旋转、缩放、平移的构建思路并结合ArcGIS软件给出从数据导入、参数设置到坐标转换的完整操作办法可用于处理历史测绘成果向CGCS2000的基准统一。我国自2008年启用2000国家大地坐标系过渡期内大量基于西安80的存量数据面临转换需求该文对数据管理、地图制图与工程规划人员有直接借鉴意义。资源包内为1个PDF文档共514KB篇幅虽短但理论结合实践层次清晰。已有5413人学习/下载适合从事地理信息生产与坐标转换工作的技术人员参考有助于提升坐标转换工作的效率与准确性。1. 80西安坐标系转2000国家坐标系为什么是块硬骨头参心与地心之间隔着一整套基准我接过一个市政数据整合的活儿甲方送来的管线和地籍资料全是80西安坐标验收却要求提交2000国家坐标系成果。当时团队里有人说“在ArcGIS里把坐标系改成CGCS2000不就行了”结果第一批数据转完点位整体向东南方向漂了几十米局部地区甚至上百米。问题不在ArcGIS操作而在于80西安坐标系是参心坐标系2000国家坐标系是地心坐标系两套基准的椭球参数、坐标原点、定向方式都不一样不能靠改标注解决。ArcGIS里真正要做的是建立地理变换Geographic Transformation把一套基准下的坐标换算到另一套基准下。这篇文章就按我自己做这类项目的顺序来写先讲清两套坐标系的差异和转换模型怎么选再给出一条能在ArcGIS里完整落地的作业流程包括控制点处理、七参数解算、批量投影和最后的精度验证最后把几处最容易翻车的坑单独拿出来说。2. 坐标系差异与转换模型先搞清楚你在跟什么打交道2.1 椭球参数和坐标原点的差异一张表说清两套基准80西安坐标系ArcGIS里叫Xian 1980GCS_Xian_1980采用的是1975年国际大地测量与地球物理联合会推荐的椭球长半轴6378140米扁率分母298.2572000国家大地坐标系CGCS2000采用CGCS2000椭球长半轴6378137米扁率分母298.257222101。从椭球参数看两者差值很小但坐标原点差异是根本性的。80西安坐标系是参心坐标系坐标原点位于我国中部某个大地原点参考椭球面与我国大地水准面整体最佳拟合CGCS2000是地心坐标系原点在地球质心与ITRF框架对齐。原点和定向不同导致同一地面点在两套坐标系下的经纬度差异并不是固定常数而是随地区变化一般在我国境内表现为几十米量级的差异东西部还不一样。这个差异在ArcGIS里非常直观。你在ArcMap或ArcGIS Pro里新建一个空白地图先加载一个带80坐标的要素类再加载一个CGCS2000底图两套数据叠不到一起去偏移方向大多指向东南或西南视地区而定。真正做转换时椭球参数差异会被换算公式吸收一部分剩下由基准差异引起的部分则要靠转换参数来消掉。所以我一般建议先不要在ArcGIS里急着做Project而是把两套坐标系的差异和转换模型确定好再决定用几参数。对比项80西安坐标系2000国家大地坐标系坐标系性质参心坐标系地心坐标系椭球IAG-75椭球CGCS2000椭球长半轴 a6378140 米6378137 米扁率 1/f298.257298.257222101坐标原点参考椭球中心与我国大地水准面最佳拟合地球质心与ITRF框架对齐ArcGIS中地理坐标系名称GCS_Xian_1980GCS_China_Geodetic_Coordinate_System_20002.2 三参数、七参数、格网改正选型要看测区范围和精度指标转换模型的选择取决于测区大小和成果精度要求。最常见的是三参数和七参数。三参数只做三个坐标轴方向的平移相当于把两套地心直角坐标系的原点错位抹平适用于半径几十公里以内、精度要求不高的小范围工程比如单一地块的示意图或有备案需求的快速出图。七参数是在三个平移量之外再增加三个坐标轴旋转和一个尺度因子即布尔沙-沃尔夫模型能够吸收两套基准之间的旋转差和尺度差是80转2000最常用的模型。ArcGIS的自定义地理变换中提供多种方法做80到2000时一般选BURSA_WOLF。需要注意的是ArcGIS里还有COORDINATE_FRAME方法它在数学上和BURSA_WOLF几乎一样但旋转参数的符号约定相反。国内很多坐标转换软件输出的七参数默认是Bursa-Wolf约定若在ArcGIS里选成COORDINATE_FRAME旋转方向会反过来点位误差会显著变大这一点在后面避坑部分再展开。当测区范围大且精度要求高比如一个地级市以上的基础测绘成果整合七参数残差往往不够。更稳妥的做法是格网改正法先按七参数做整体转换再用剩余残差生成格网改正模型。ArcGIS本身不直接提供国家发布的格网改正文件但可以通过对控制点残差做空间插值来实现类似效果我把它放在最后一章专门说。另外还存在平面四参数加高程拟合的路线平面用四参数两个平移、一个旋转、一个尺度高程单独用一个高程拟合面处理但ArcGIS的自定义地理变换不支持四参数这条路通常要借助第三方程序完成平面转换再把高程改正叠加上去。3. 动手前先理顺数据ArcGIS里的坐标系识别与投影定义检查3.1 先确认数据是不是真的“西安80”别急着转我见过不少翻车现场问题不是转换参数而是源数据根本就不是西安80。有些采集设备默认输出WGS84经纬度现场记录时被标成了80坐标有些数据在流转过程中被反复Define Projection坐标系信息已经变成Unknown。所以第一步不是转坐标而是确认数据当前带的是什么坐标系。在ArcGIS里最简单的办法是看一眼要素类的空间参考属性。也可以直接用arcpy写两行代码批量检查import arcpy fc rD:\work\X80\land.shp sr arcpy.Describe(fc).spatialReference print(sr.name) print(sr.factoryCode)这段代码输出的是要素类的空间参考名称和工厂代码。如果输出类似“Xian_1980_3_Degree_GK_CM_117E”说明数据已经带了80的投影坐标系如果输出“Unknown”或“Unnamed”说明空间参考信息缺失需要先用Define Projection手工定义。注意factoryCode是一个整数编号ArcGIS里每个坐标系都有对应的WKID80的地理坐标系是4610CGCS2000的地理坐标系是4490看到这两个数字可以快速判断数据基准。还有一种常见情况是数据本身是80地理坐标经纬度只是被放到了一个通常为投影坐标系的项目里加载。这时要看sr.linearUnit是否为米如果线性单位是度DecimalDegrees说明数据还是地理坐标系得先投影到高斯-克吕格投影再做后续转换。我一般会把检查结果导成一个Excel清单逐层核对每个图层的空间参考确认无误后再统一进入转换流程。3.2 投影带号和中央经线Y坐标前两位就藏在带号里80西安坐标和CGCS2000在高斯-克吕格投影下的分带规则是相同的三度带和六度带的带号算法一致。但ArcGIS里带号选错造成的偏移比基准差异还夸张。一个常见场景是数据坐标Y值为8位前两位是20即六度带20带中央经线117度如果误定义成三度带117E或者定义成六度带19带点位东西方向能差出几十公里。判断带号不靠猜靠看Y坐标的位数和前缀。高斯-克吕格投影的Y坐标加带号时Y值会有7位或8位前一位或前两位就是带号。具体判断方法可以这样如果Y坐标约为8位且前两位在13到23之间多半是六度带带号中央经线按带号和经度的关系计算六度带中央经线等于带号乘6再减3如果Y坐标前三位在表示三度带带号三度带中央经线等于带号乘3实际作业中用带号推算中央经线最不容易出错。拿到数据后把源数据的Y坐标列一遍看是否带带号、带号是多少再和ArcGIS坐标系名称里的CM值比对。比如坐标系名Xian_1980_3_Degree_GK_CM_117E表示三度带中央经线117度如果数据实际是六度带20带两者就不能混用。import arcpy fc rD:\work\X80\land.shp sr arcpy.Describe(fc).spatialReference if sr.type Projected: print(中央经线:, sr.centralMeridian) print(投影名称:, sr.name) else: print(当前是地理坐标系先投影再转)这段代码用来判断投影坐标系的中央经线避免带号问题带到转换环节。如果打印出来的中央经线与数据实际带号对不上不要用Project而应该用Define Projection重新定义正确的坐标系定义正确后再做转换。这一步是我每次做批量转换前的固定动作数据量大时尤其值得花时间因为一旦带号错后面的七参数再怎么调都白搭。4. 从控制点到七参数ArcGIS完成80到2000转换的完整作业流程4.1 控制点整理与地心直角坐标换算七参数解算需要控制点在两套坐标系下同时具有坐标。控制点来源一般是测区内的已知控制点资料要求点位分布均匀至少3个实际作业我建议不少于6个再多留几个点做检查。这里先给一个最容易忽略的前提解算七参数用的是地心直角坐标X、Y、Z不是平面坐标也不是经纬度。所以拿到控制点的80平面坐标后要先把它转成80经纬度再连同椭球高一起换算成空间直角坐标。控制点的80平面坐标是高斯投影坐标要先用ArcGIS转成地理坐标。把控制点做成一个点要素类坐标系定义为80的投影坐标系然后用Project工具输出到GCS_Xian_1980得到经纬度。注意这个过程中不能用任何地理变换因为是同一个基准内部做投影反算ArcGIS会直接处理。转出经纬度后再用下面的公式把经纬度转成地心直角坐标import numpy as np def geo2cart(lon_deg, lat_deg, h, a, f): 经纬度和椭球高转地心直角坐标 lon_deg/lat_deg: 十进制度 h: 椭球高单位米 a: 椭球长半轴 f: 扁率 lon np.radians(lon_deg) lat np.radians(lat_deg) e2 2 * f - f * f N a / np.sqrt(1 - e2 * np.sin(lat) ** 2) X (N h) * np.cos(lat) * np.cos(lon) Y (N h) * np.cos(lat) * np.sin(lon) Z (N * (1 - e2) h) * np.sin(lat) return X, Y, Z这里的h必须是椭球高不是正常高。西安80坐标系下的控制点高程资料通常指正常高椭球高等于正常高加高程异常。如果手头没有高程异常模型一个替代做法是优先选择已知椭球高或GPS观测高程的控制点参与解算或者暂时忽略高程分量的影响只让控制点提供平面约束。需要说明的是忽略高程异常会给Z分量带来误差继而影响七参数中的旋转参数和尺度参数这是加密控制点时最容易出问题的地方也是后面避坑章要展开的一条。4.2 最小二乘解算七参数用Python脚本而不是手工算有了控制点在两套坐标系下的地心直角坐标下一步就是解算布尔沙-沃尔夫七参数。手工按公式计算几乎不现实常见做法是用坐标转换软件来算但如果想保证过程透明、参数可复现我更推荐直接用Python写一个最小二乘解算脚本。布尔沙-沃尔夫模型在小旋转角假设下可以线性化每个控制点能列三个误差方程方程组的系数阵由源坐标和待求参数组成用numpy的最小二乘求解即可。import numpy as np def bursa_wolf_lstsq(src_pts, dst_pts): 用最小二乘解算布尔沙-沃尔夫七参数 src_pts: 源坐标系控制点列表每个点为 (X, Y, Z) dst_pts: 目标坐标系控制点列表每个点为 (X, Y, Z) 返回平移(m)、旋转(角秒)、尺度(ppm) n len(src_pts) A np.zeros((3 * n, 7)) L np.zeros(3 * n) for i, (src, dst) in enumerate(zip(src_pts, dst_pts)): X, Y, Z src A[3 * i] [1, 0, 0, 0, -Z, Y, X] A[3 * i 1] [0, 1, 0, Z, 0, -X, Y] A[3 * i 2] [0, 0, 1, -Y, X, 0, Z] L[3 * i:3 * i 3] np.array(dst) - np.array(src) x, _, _, _ np.linalg.lstsq(A, L, rcondNone) tx, ty, tz, rx_rad, ry_rad, rz_rad, scale x rx_sec np.degrees(rx_rad) * 3600 ry_sec np.degrees(ry_rad) * 3600 rz_sec np.degrees(rz_rad) * 3600 return {TX: tx, TY: ty, TZ: tz, RX: rx_sec, RY: ry_sec, RZ: rz_sec, Scale: scale * 1e6}这段代码的逻辑是基于线性化后的误差方程把所有控制点代入后组成法方程用最小二乘得到七参数。旋转角输出单位换算成角秒尺度输出单位换算成ppm这两个单位也是ArcGIS自定义地理变换里要求的单位。脚本里有一个隐含假设是旋转角足够小线性化误差可忽略如果解出来的旋转角超过几十角秒就要考虑控制点本身有问题或者基准差异超出线性化范围不建议直接使用这组参数。参数解算完不能直接用要先检查各控制点残差。把每一对控制点代入解算结果算一下转换后坐标与已知目标坐标的差值残差大的点往往坐标本身有问题需要剔除后重新解算。我一般会把残差控制在3厘米以内平面检查点在测区内控制在5厘米以内如果达不到优先检查控制点坐标的可靠性。4.3 在ArcGIS里创建自定义地理变换并批量投影七参数解算完成后要把这组参数填进ArcGIS。ArcMap 10.x里通过Data Management Tools下的Create Custom Geographic Transformation创建方法选择BURSA_WOLF把平移、旋转、尺度依次填入。ArcGIS Pro中则在Project工具里点新建地理变换界面类似。这里有一个版本差异值得注意ArcGIS Pro较新版本里Project对话框会提供Calculate Transformation按钮可以直接选择源坐标系和目标坐标系再指定控制点要素类让软件自动解算七参数并生成变换省去外部脚本这一步。我习惯先用检查点反代一遍再创建变换。具体方式是拿一个控制点的80坐标用解出的七参数手工算一遍目标坐标和该点已知的2000坐标对比确认方向和数值无误后再往ArcGIS里填。因为哪怕参数只差一个符号转换结果也是几米甚至十几米的偏差这个反代验证能提前拦截大部分低级错误。批量转换用arcpy脚本是最省事的。假设已经创建好了名为X80_to_CGCS2000的自定义变换对所有要素类执行Projectimport arcpy import os arcpy.env.workspace rD:\work\X80 out_ws rD:\work\CGCS2000\result.gdb transformation X80_to_CGCS2000 out_sr arcpy.SpatialReference(CGCS2000_3_Degree_GK_CM_117E) if not arcpy.Exists(out_ws): arcpy.CreateFileGDB_management(os.path.dirname(out_ws), os.path.basename(out_ws)) for fc in arcpy.ListFeatureClasses(): name arcpy.Describe(fc).basename out_fc os.path.join(out_ws, name _2000) arcpy.Project_management(fc, out_fc, out_sr, transformation) print(完成:, name)输出坐标系要和变换配套。如果源数据是80三度带117E输出也应该是CGCS2000三度带117E不要跨带转换。Project工具在遇到源和目标坐标系带不同中央经线时会先转到地理坐标再转到目标投影过程中如果地理变换设置正确换带和换基准会一次完成。脚本里arcpy.Project_management的第四个参数就是前面创建的地理变换名称如果空着不填ArcGIS会用内置的默认变换但80到2000通常没有内置可用变换所以这个参数必须明确指定。5. 常见问题避坑指南转换结果偏了几十米大概率踩了这四个坑5.1 三个点为什么不够控制点数与分布问题有朋友问我“规范上说三个公共点就能解七参数我找了三个点怎么算出来的结果残差这么大”。现象是解算结果在控制点上符合得很好但压到其他检查点上误差达到几米。原因是三个点只能满足唯一解控制点本身的采集误差会直接进入参数没有多余观测来平差更严重的是如果三个点分布近似共线或都集中在测区一角法方程矩阵接近病态旋转参数和尺度参数解算不稳定外部检查点一验证就露馅。解决方法是控制点数量至少6个分布上要包住整个测区四个角点加中间区域都要有点。解算后逐个检查控制点残差把残差明显偏大的点剔除再解算一次。如果测区跨越大比如一个县或一个市建议分区解算七参数不要用一套参数硬套整个区域。我实践中会把控制点分成两组一组参与解算另一组只做检查确保检查点不参与参数计算这样才能反映转换的真实精度。5.2 旋转参数的单位和符号跨软件填数最容易翻车现象是用某个国产坐标转换软件算出的七参数填进ArcGIS后所有点位向同一方向弧线漂移从几十厘米到几米不等。原因有两类一是旋转参数单位没换算某些软件输出的是弧度ArcGIS要求角秒差一个3600倍的系数二是符号约定相反ArcGIS里BURSA_WOLF和COORDINATE_FRAME两种方法的旋转符号本就相反如果从别的软件复制参数时选的坐标系和ArcGIS的方法不一致旋转方向就反了。解决方法是统一单位旋转参数必须换算成角秒再填入尺度参数用ppm。填完后找一个已知控制点做反代验证把该点的80坐标输入用ArcGIS的地理变换转一遍看输出是否和已知2000坐标对上。这个验证步骤只需要几分钟能省掉后面整批数据返工的麻烦。反代验证时注意源坐标和目标坐标要用同一套椭球模型不要把80坐标直接当成经纬度喂给WGS84椭球。5.3 带号/中央经线填错根本不是坐标转换问题是投影定义问题现象是数据转换后整体向东西方向偏移几十公里远超基准差异的几十米量级。这类问题出在输入数据的坐标系定义阶段而不是转换参数。常见场景是源数据Y坐标带有带号如六度带20带但ArcGIS里坐标系被定义成三度带117E或者反过来。坐标转换工具会先按错误带号反算经纬度再用这个错误的经纬度去算目标坐标结果自然整体偏移。解决办法是在转换前核对Y坐标。如果Y坐标前两位是20对应六度带20带中央经线117度如果是三度带40带Y坐标前三位是40中央经线120度数据文件名或台账里一般会注明。发现定义错误时用Define Projection重新定义坐标系不要用Project去“纠正”。Project是在两套坐标系之间做变换定义错了要先纠正投影定义本身这个逻辑顺序不能反。批量处理前建议先抽样检查两个图层的坐标系定义别让带号问题混进大批量作业。5.4 高程异常没处理平面残差合格了高程还是对不上现象是平面转换后XY误差在2到3厘米看起来挺漂亮但检查点的高程差了30到40厘米或者用正常高直接参与七参数解算时平面残差出现系统性波动。原因在于七参数解算用的是空间三维坐标Z分量必须是椭球高而通常控制点提供的是正常高。正常高和椭球高之间隔着高程异常不同地区差异不同山区和平原差异也不同。把这个没有处理的高差带进解算旋转参数会被污染表面上控制点拟合得不错控制点以外的区域就会出现平面上的扭曲。解决方法是解算前先把控制点高程统一换算成椭球高。如果有省级似大地水准面精化模型或CORS站提供的高程异常优先使用没有模型时尽量选择已知椭球高或做过GPS静态观测的点做控制点并检查解算后残差中是否带有与高程相关的趋势。如果项目成果只需要平面坐标可以改用二维转换思路把平面转换和高程转换分开处理平面用四参数或七参数用椭球高尽量准确高程单独建立拟合模型不混在一套解算里。6. 精度验证与再提升残差统计、检查点评估与格网改正的进阶用法6.1 用检查点算残差与RMSE别拿转换参数自证转换完成后一定要独立验证。验证用的检查点不能是参与七参数解算的那些点否则只是自证。正确做法是把检查点的已知2000坐标和转换后的坐标放到一起逐个算平面点位误差再统计最大误差和均方根误差。ArcGIS里操作路径是先用Add XY Coordinates把转换后点要素的坐标写入属性表再用Join Field把检查点已知坐标合并进来最后用字段计算器算误差。import arcpy arcpy.env.workspace rD:\work\CGCS2000 pts check_pts_after arcpy.AddXY_management(pts) arcpy.JoinField_management(pts, PtID, check_pts_known, PtID, [Known_X, Known_Y]) arcpy.AddField_management(pts, dN, DOUBLE) arcpy.AddField_management(pts, dE, DOUBLE) arcpy.AddField_management(pts, Dist, DOUBLE) arcpy.CalculateField_management(pts, dN, !POINT_Y! - !Known_Y!, PYTHON3) arcpy.CalculateField_management(pts, dE, !POINT_X! - !Known_X!, PYTHON3) arcpy.CalculateField_management(pts, Dist, (!dN! ** 2 !dE! ** 2) ** 0.5, PYTHON3)这段代码先把转换后的检查点坐标提取出来再连接已知坐标最后计算南北向误差dN、东西向误差dE和平面点位误差Dist。注意AddXY用的POINT_X和POINT_Y是地理坐标还是投影坐标取决于要素类的坐标系如果要素类是投影坐标系输出的X和Y就是平面坐标可以直接用于误差统计。统计完看RMSE是否满足项目验收指标。我一般要求检查点平面中误差在5厘米以内如果达不到就回到参数解算环节查原因。6.2 残差插值改正七参数残差大时的补救手段七参数模型是整体拟合测区范围大或局部地区基准变化复杂时即使控制点残差很小中间区域的检查点仍可能出现厘米级甚至分米级偏差。这时可以做一个格网改正把所有检查点残差生成点要素用IDW或克里金插值成改正栅格再对转换后的坐标做二次改正。具体做法是先在属性表里把Dist字段的符号按实际偏差方向整理成dX和dY两个字段生成点要素集然后分别在dX和dY上执行插值得到两个改正栅格最后用栅格计算器把改正值加到转换后点坐标上。这个方法的适用范围是测区内部插值栅格不能外推到测区以外否则改正值不可信。如果残差呈现出明显的分区特征说明一套七参数已经撑不起整个测区更干脆的办法是把测区拆分成几个片区各自解算七参数后再接边。我习惯在每个转换项目里都留一套未参与解算的检查点数量不少于控制点总数的三分之一并且在成果提交时把检查点误差统计表一并交付。这样验收时对方可以直接看数据而不是听我口头保证。做这类转换项目多了以后最大的体会是坐标转换的难点从来不是ArcGIS工具本身而是输入数据是否干净、控制点是否可靠、参数单位是否一致。希望帮到你。本文还有配套的精品资源点击获取