简介面向测绘、地理信息与导航定位从业者及测绘专业学习者的坐标转换资料包围绕国家大地2000CGCS2000与西安80之间的高精度转换需求汇集转换工具、参数配置与示例截图可用于历史数据兼容处理、GIS数据处理和地图制作等场景。压缩包为zip格式共97个文件、约4.26MB以49个cod坐标转换定义文件为核心辅以csv参数表、txt说明、xml配置及jpg/png/gif操作截图另有zgf、ggf、cfg等辅助数据便于按坐标系统与转换模型查阅。资料中配套COORD GM2.0相关程序与配置可支持七参数、三参数等转换思路的对照学习。已有3176人学习下载适合需要理解CGCS2000与西安80转换流程、核对参数文件组织方式并借助示例截图排查转换问题的读者参考。1. 一批西安80 的老成果要并入 CGCS2000 库卡住的从来不是代码手上接过一批测绘老成果要入库控制点还是西安80 的成果表里 X、Y 是 6 位加 8 位的一串数甲方要的是 CGCS2000 下的经纬度和平面坐标。两套坐标系同一个物理点平面位置差几十米到一百多米直接叠加就对不上图纸和实地也就对不上。这套坐标转换器解决的就是这条链路把西安80 的平面成果还原成大地坐标走七参数转进 CGCS2000再按目标带号做高斯投影落回平面。适合测绘、GIS 数据治理、不动产与管线入库这类必须动真参数的场景也适合只会调库、没自己解过七参数的人补一次原理课。2. 两套椭球的硬参数与转换路径为什么不能直接套一个 towgs84选型这件事上翻过车的人不少拿一段网上抄来的七参数扔进towgs84就跑跑完看着差不多就交了回头发现测区边缘差了三十多米。根子在于西安80 和 CGCS2000 不是同一个椭球两套坐标系的原点、定向、尺度也不一致转换参数跟测区强相关不存在一套全国通用的数。2.1 椭球参数差在哪几个数上先把账算清楚。CGCS2000 采用的是 GRS80 那一套椭球参数西安80 采用的是 IAG-75 椭球两者长半轴差 3 米扁率倒数差约 0.0002看着极小落到大地坐标上却是可观的偏差。项目CGCS2000西安80差值长半轴 am6378137.06378140.03.0扁率倒数 1/f298.257222101298.257约 0.0002221第一偏心率平方 e²0.006694380022900.00669438499959约 4.98e-9常用 EPSG地理坐标44904610—常用 EPSG地心直角4479——短半轴差 3 米带来的直接后果是同一个空间直角坐标在两套椭球下反算出来的经纬度不一样同一个经纬度在两套椭球下正算出来的平面坐标也不一样。工程上这个差异的量级通常在几十米再叠加两套坐标系的定向与尺度差异最终点位差落在几十米到一百多米之间具体取决于测区位置。所以任何零参数直接转的做法只能用在精度要求极低的示意性场景。真正需要盯住的是三个量椭球参数、七参数三个平移、三个旋转、一个尺度、投影带号与中央子午线。三者错一个结果就废。2.2 三参数、四参数、七参数各自的边界转换模型的选型不是越复杂越好是看你的公共点条件和测区尺度。模型未知数最少公共点适用场景主要局限三参数平移31精度要求低、小范围示意忽略旋转与尺度边缘误差大四参数平面42小范围平面成果跨度 30 km 以内只处理平面高程另算七参数布尔莎-沃尔夫73跨椭球、跨基准的三维转换主流做法公共点分布差时参数外推不稳定七参数莫洛琴斯基73与布尔莎等价公式形式不同同样依赖公共点质量我一般的判断顺序是跨椭球一律上七参数除非测区小到几公里且只做平面示意四参数留着做同一椭球下不同投影带之间的平面挪移三参数只当调试用的临时手段不进正式成果。七参数至少 3 个公共点才够解每个点给 3 个方程7 个未知数需要 9 个方程但工程上少于 4 个点我基本不敢用因为 3 个点解出来的参数没有多余观测残差恒为零你根本不知道自己解对没解对。公共点还得在测区里铺开四个角加中心最好全挤在一条线上必翻车。2.3 把转换器的运行环境搭起来这套东西依赖很轻Python 3.8 以上加两个库就够。先建个干净环境别在系统 Python 里折腾坐标转换这种活儿最怕版本玄学。python -m venv venv source venv/bin/activate # Windows 下用 venv\Scripts\activate pip install pyproj3.4 numpy1.23 python -c import pyproj, numpy; print(pyproj.__version__, numpy.__version__)pyproj是 PROJ 的 Python 绑定底层负责椭球定义、投影正反算、基准转换管线numpy负责七参数的最小二乘求解。版本上给个下限就行pyproj3.4 之后Transformer.from_pipeline的行为比较稳定再老的版本在inv步骤上有过不一致。工程目录我习惯这么切coord-convert/ ├── params/ # 七参数、带号配置按测区一份 ├── data/ # 输入的西安80 成果、公共点对 ├── core/ # 椭球、投影、七参数、批量转换 └── check/ # 检核脚本与残差报告冒烟测试别省先拿一个已知点跑通全链路再上批量。转换器本身不复杂复杂的是参数和边界所以入口脚本要能把参数显式打出来别藏在代码里。3. 七参数求解布尔莎模型的矩阵形式与 numpy 实现七参数是整个转换链路里唯一需要解的部分其余都是确定性的公式。这一章把求解过程拆开代码可以直接抄。3.1 大地坐标转空间直角入口必须对齐七参数是在空间直角坐标系地心坐标系里定义的不是在大地坐标里。所以第一步是把两套坐标系下的公共点各自按自己的椭球参数转成空间直角坐标。这一步错了后面全错。import numpy as np # 两套椭球的硬参数手抄容易错集中放一处 ELLIPSOIDS { CGCS2000: {a: 6378137.0, f: 1 / 298.257222101}, XIAN80: {a: 6378140.0, f: 1 / 298.257}, } def geodetic_to_ecef(lon_deg, lat_deg, h, ellps): 大地坐标 (B, L, H) - 空间直角坐标 (X, Y, Z)角度单位为度高程单位米 a ellps[a] f ellps[f] e2 f * (2 - f) # 第一偏心率平方 B np.radians(lat_deg) L np.radians(lon_deg) N a / np.sqrt(1 - e2 * np.sin(B) ** 2) # 卯酉圈曲率半径 X (N h) * np.cos(B) * np.cos(L) Y (N h) * np.cos(B) * np.sin(L) Z (N * (1 - e2) h) * np.sin(B) return np.array([X, Y, Z]) def ecef_to_geodetic(xyz, ellps, tol1e-10, max_iter100): 空间直角坐标 - 大地坐标用迭代法反算收敛判据是纬度变化量 a ellps[a] f ellps[f] e2 f * (2 - f) X, Y, Z xyz lon np.arctan2(Y, X) p np.hypot(X, Y) lat np.arctan2(Z, p * (1 - e2)) # 初值给球面近似 for _ in range(max_iter): N a / np.sqrt(1 - e2 * np.sin(lat) ** 2) h p / np.cos(lat) - N lat_new np.arctan2(Z, p * (1 - e2 * N / (N h))) if abs(lat_new - lat) tol: lat lat_new break lat lat_new N a / np.sqrt(1 - e2 * np.sin(lat) ** 2) h p / np.cos(lat) - N return np.degrees(lon), np.degrees(lat), h几个参数得说清楚。a是长半轴f是扁率e2 f * (2 - f)是第一偏心率平方这是由扁率推偏心率的恒等式别写成2f - f²之外的花样。N是卯酉圈曲率半径代表该纬度处椭球面法线方向的曲率。高程h用的是大地高不是正常高如果输入是 1985 国家高程基准的正常高严格来说要先做高程异常改正但在几十公里的测区里高程异常对平面转换的影响很小做平面成果时可以忽略做三维成果时必须补。反算用迭代初值取球面近似一般 3 到 5 次就收敛到 1e-10 弧度。别用那个闭式解公式边界条件容易出错。3.2 布尔莎-沃尔夫模型的矩阵形式布尔莎模型说的是目标坐标系下的空间直角坐标等于源坐标系坐标经过平移、旋转、尺度缩放之后的结果。写成小角度近似下的线性形式X2 dx (1 m)·X1 - rz·Y1 ry·Z1 Y2 dy rz·X1 (1 m)·Y1 - rx·Z1 Z2 dz - ry·X1 rx·Y1 (1 m)·Z1其中dx, dy, dz是三个平移量单位米rx, ry, rz是三个旋转角单位弧度m是尺度因子量级在 1e-6 左右习惯上按 ppm 说。把每个公共点写成三个误差方程就得到V B·x - l的形式x是七个未知数B是系数阵l是坐标差观测向量。系数阵的每一行是B_X [1, 0, 0, 0, Z1, -Y1, X1] B_Y [0, 1, 0, -Z1, 0, X1, Y1] B_Z [0, 0, 1, Y1, -X1, 0, Z1]l [X2-X1, Y2-Y1, Z2-Z1]单位米。注意旋转项的符号这套符号对应 PROJ 里的position_vector约定用coordinate_frame约定的话三个旋转角全体反号这是最常见的翻车点后面避坑章节会专门说。严格讲布尔莎模型中(1m)会乘到旋转项上这里按小角度近似省略了。省略带来的误差是 1e-6 乘以旋转量再乘以坐标量量级在 1e-11 米完全在测量噪声之下可以忽略。3.3 用 numpy 解七参数最小二乘与残差公共点超过 3 个以后方程数大于未知数用最小二乘解。def solve_helmert_7params(src_ecef, dst_ecef): src_ecef: (n, 3) 源坐标系如西安80空间直角坐标单位米 dst_ecef: (n, 3) 目标坐标系如 CGCS2000空间直角坐标单位米 返回: dict(平移, 旋转(角秒), 尺度(ppm)), 残差向量 n len(src_ecef) if n 3: raise ValueError(公共点少于 3 个七参数无唯一解) B np.zeros((3 * n, 7)) l np.zeros(3 * n) for i, ((x1, y1, z1), (x2, y2, z2)) in enumerate(zip(src_ecef, dst_ecef)): r 3 * i B[r] [1, 0, 0, 0, z1, -y1, x1] B[r 1] [0, 1, 0, -z1, 0, x1, y1] B[r 2] [0, 0, 1, y1, -x1, 0, z1] l[r], l[r 1], l[r 2] x2 - x1, y2 - y1, z2 - z1 x, *_ np.linalg.lstsq(B, l, rcondNone) # 最小二乘 v B x - l # 残差单位米 return { dx: x[0], dy: x[1], dz: x[2], # 米 rx: np.degrees(x[3]) * 3600, # 弧度 - 角秒 ry: np.degrees(x[4]) * 3600, rz: np.degrees(x[5]) * 3600, s: x[6] * 1e6, # 无量纲 - ppm }, v def residual_report(v): 残差体检整体 RMS、各分量最大残差、是否有离群点 rms np.sqrt(np.mean(v ** 2)) mx np.max(np.abs(v)) per_axis np.max(np.abs(v.reshape(-1, 3)), axis0) return {rms_m: rms, max_m: mx, max_xyz: per_axis}np.linalg.lstsq用的是 SVD 分解求解比直接算(BᵀB)⁻¹Bᵀl稳公共点少或者分布差的时候差别很明显。解出来的x[6]是无量纲的尺度因子乘以 1e6 就是 ppmx[3..5]是弧度乘180/π再乘 3600 才是角秒。残差向量v的用途比参数本身更大。RMS 落在厘米级说明公共点质量好、模型合适RMS 到分米级就要回头看公共点坐标的精度如果某一个点的残差明显大于其他点大概率是那个点的坐标抄错了或者点之记对不上把它剔掉重解参数会明显改善。工程上我会做一次留一验证轮流剔掉一个公共点用剩下的解参数再把这个点转过去看偏差。七个点里如果有哪个被剔掉后参数大幅变化说明这个点权重过大或者数据有问题这对测区边缘的精度判断很关键。4. 高斯投影与带号把平面坐标还原成经纬度七参数管的是椭球之间的事投影管的是椭球面到平面的事。这两件事分开做中间用大地坐标衔接是最不容易出错的链路。4.1 带号与中央子午线先从 Y 值里读出来国内基本比例尺地形图用的是高斯-克吕格投影分 3 度带和 6 度带。带号和中央子午线的关系6 度带带号 N floor(L / 6) 1中央子午线 L0 6N - 3 3 度带带号 n round(L / 3) 中央子午线 L0 3nL是经度单位度。举个具体的某点经度 116.4°6 度带带号是floor(19.4)1 20中央子午线 117°3 度带带号是round(38.8) 39中央子午线也是 117°。同一个点两种分带的中央子午线可能一样但带号不一样写进成果表就错位了。从平面坐标反推带号靠 Y 值。加了带号的 Y 值一般是 8 位前两位就是带号比如39500000.000表示 39 带、假东 500 km 处。没加带号的 Y 是 6 位。这个判据简单但极容易搞错因为老成果表里加不加带号、加 3 度带号还是 6 度带号各家单位习惯不一样。def guess_zone(y): 从 Y 坐标粗略判断带号与分带方式返回 (带号, 分带, 是否含带号) s f{abs(y):.3f} int_part s.split(.)[0] if len(int_part) 8: zone int(int_part[:2]) # 3 度带带号范围 25~456 度带 13~23重叠区间要人工确认 if 25 zone 45: return zone, 3, True if 13 zone 23: return zone, 6, True return zone, None, True return None, None, False这个函数只能给参考重叠区间比如带号 25~45 里也有一部分是 6 度带号必须人工确认。我在实际项目里的做法是拿到成果先抽 3 个点用假定带号算一遍中央子午线反算经纬度看是否落在测区范围内落在范围内才认。4.2 pyproj 做投影正反算投影这块不用自己写公式pyproj足够稳但要把字符串拼对。from pyproj import CRS, Transformer def make_tm_crs(ellps_name, lon0, false_easting500000.0, k01.0): 构造高斯-克吕格横轴墨卡托投影 CRSlon0 为中央子午线度 proj4 ( fprojtmerc lat_00 lon_0{lon0} k{k0} fx_0{false_easting} y_00 ellps{ellps_name} unitsm no_defs ) return CRS.from_proj4(proj4) # 示例西安8039 带中央子午线 117° xian80_tm make_tm_crs(IAU76, 117.0) cgcs2000_tm make_tm_crs(GRS80, 117.0) # 大地坐标 - 平面坐标 tr_fwd Transformer.from_crs(CRS.from_proj4(projlonglat ellpsIAU76 no_defs), xian80_tm, always_xyTrue) x, y tr_fwd.transform(116.4, 39.9) # 注意 always_xy 决定了输入顺序是 (lon, lat) # 平面坐标 - 大地坐标 tr_inv Transformer.from_crs(xian80_tm, CRS.from_proj4(projlonglat ellpsIAU76 no_defs), always_xyTrue) lon, lat tr_inv.transform(x, y)ellps参数是关键西安80 在 PROJ 里叫IAU76长半轴 6378140扁率倒数 298.257CGCS2000 用GRS80长半轴 6378137扁率倒数 298.257222101。这两个名字写反了出来的坐标会差几十米而且不会有任何报错属于典型的静默错误。always_xyTrue这个参数必须加。不加的话pyproj会按 CRS 的轴序来地理坐标是纬度在前很多人不知道这件事结果经纬度整个对调反算出来的点跑到另一个半球去了。k1.0是中央子午线处的投影尺度因子高斯投影标准取值就是 1.0别改成 0.9996那是 UTM 的。4.3 批量转换的完整链路把前面几块串起来一条完整的转换链路是这样的def xian80_to_cgcs2000(lon, lat, h, params, zone_lon0): 西安80 经纬度 - CGCS2000 经纬度params 为上一章解出的七参数 # 1. 西安80 大地坐标 - 西安80 空间直角 xyz1 geodetic_to_ecef(lon, lat, h, ELLIPSOIDS[XIAN80]) # 2. 布尔莎七参数正算 dx, dy, dz params[dx], params[dy], params[dz] rx, ry, rz (np.radians(params[k] / 3600) for k in (rx, ry, rz)) m params[s] * 1e-6 X1, Y1, Z1 xyz1 X2 dx (1 m) * X1 - rz * Y1 ry * Z1 Y2 dy rz * X1 (1 m) * Y1 - rx * Z1 Z2 dz - ry * X1 rx * Y1 (1 m) * Z1 # 3. CGCS2000 空间直角 - CGCS2000 大地坐标 return ecef_to_geodetic(np.array([X2, Y2, Z2]), ELLIPSOIDS[CGCS2000])三个步骤各自独立可以单独验证。第一步对不上就用已知点的空间直角坐标直接比对第二步对不上就检查符号约定第三步对不上就检查椭球参数。这种拆法比一整块糊在一起好排查得多。批量的时候加一层循环和失败记录别让一个坏点把整个文件中断import csv def batch_convert(in_csv, out_csv, params): ok, fail 0, [] with open(in_csv, newline, encodingutf-8) as f, \ open(out_csv, w, newline, encodingutf-8) as g: reader csv.DictReader(f) writer csv.writer(g) writer.writerow([id, lon, lat, h, status]) for row in reader: try: lon, lat, h xian80_to_cgcs2000( float(row[lon]), float(row[lat]), float(row.get(h) or 0), params, None) writer.writerow([row[id], f{lon:.9f}, f{lat:.9f}, f{h:.3f}, OK]) ok 1 except Exception as e: writer.writerow([row[id], , , , fFAIL: {e}]) fail.append((row[id], str(e))) return ok, fail输出保留 9 位小数是因为 1e-9 度大约对应 0.1 毫米足够覆盖任何测量精度保留更多位没有意义反而让文件变大。失败点单独回写到status列跑完直接筛出来重跑比中断重来省事。5. 避坑与排查五条踩过的记录坐标转换这活儿公式都对、代码都跑得通结果还是错八成栽在下面这几处。每条我都按现象、原因、解决写遇到时照着对号入座。5.1 参数层单位与符号约定的三处陷阱第一条旋转角单位混用。现象是转出来的点整体偏移几十米而且偏移方向跟测区位置有关看着像参数不准。原因是七参数里的旋转角论文和成果表里普遍用角秒布尔莎公式里必须用弧度中间少了π/180/3600这一步。解决是把单位转换写成独立函数输入统一按角秒收进公式前立刻转弧度全程不再出现裸的旋转数值。我现在的习惯是参数字典里键名直接带单位rx_arcsec、m_ppm写错的时候一眼能看出来。第二条符号约定反了。现象是转换后残差 RMS 从厘米级突然跳到米级但平移量看着又合理。原因是布尔莎模型有两套等价的符号体系PROJ 里叫position_vector和coordinate_frame两者三个旋转角全体反号。你解参数时用的是哪套套公式时就必须用同一套。解决是全程锁死一套约定我在项目里固定用position_vector用pyproj的projhelmert时显式写上conventionposition_vector不写就默认默认值随版本变过写上最保险。第三条尺度因子小数点错位。现象是转换后的点位随距离中心越远偏差越大中心附近几乎没差。原因是尺度因子按 ppm 记公式里要乘 1e-6直接拿 ppm 的数值当无量纲用等于把尺度放大了 100 万倍。这个错误在小测区里可能只表现为边缘差几厘米容易被忽略等到跨几百公里就彻底暴露。解决是残差检验时同时看中心和边缘的点只盯着中心几个点看不出尺度问题。5.2 数据层带号与公共点分布第四条3 度带和 6 度带混用。现象是同一批成果里大部分点转换正常个别点偏了几百公里。原因是成果表拼凑自不同年代、不同单位有的加 3 度带号有的加 6 度带号还有的根本没加带号只给了个 6 位的 Y。解决是先把所有 Y 值按位数分组8 位的按前两位判带6 位的按测区经度反推判不出来的单独列出来人工确认。我一般会写一个带号体检脚本把每个点的 Y 值、推断带号、推断中央子午线、反算经度打成一张表人眼扫一遍比什么都快。第五条公共点分布不均匀导致参数外推发散。现象是公共点附近的检核点都对得上测区边缘的点差了半米以上。原因是七个公共点全挤在测区东侧西侧全靠外推七参数在小范围内是线性的外推出去误差被放大。解决是选公共点时按测区范围铺开四个角各一个、中心一个多出来的点均匀插在中间。如果实在没有边缘的公共点那就老实承认边缘精度不够在成果说明里写清楚适用范围别硬扛。5.3 流程层那些不会报错但会出事的细节还有一类问题不属于参数也不属于数据是流程上的。最常见的是拿同一批点既解参数又做检核。这批点解出来的参数回代残差当然是零看着完美实际上没有任何检核能力。正确做法是留出至少 1 到 2 个点不参与求解专门做独立检核。我在项目里会把公共点分成求解集和检核集两组比例大概 4:1检核集的中误差才是这份参数的真实精度。另一个是坐标顺序。pyproj的always_xy不写输入按纬度在前处理经纬度对调csv里 X、Y 列的顺序各家不同代码里写死row[lon]之前得确认这一列真的是经度。这类错误不会抛异常只会安静地产出一个错点靠人眼看是看不出来的只能靠已知点回代。6. 把转换固化成可复现的验证流程参数解出来、代码跑通只完成了一半。真正让这套转换器能交出去、能在半年后别人接手还能复现的是验证环节。我现在固定跑四道检核缺一道不交成果。第一道是参数自检七个参数打印出来跟历史参数比一眼平移量在几十到几百米、旋转角在角秒量级、尺度在 ppm 量级超出这个范围基本可以判定数据有问题。第二道是残差体检求解集的 RMS、最大残差、每个分量的最大残差RMS 超过 5 厘米就回头看公共点。第三道是独立检核检核集的点转过去之后跟已知值比对中误差和最大偏差都记下来这是唯一有说服力的精度指标。第四道是边界检核取测区四个角附近的点看转换后是否出现异常偏移专门用来抓外推发散。四道检核的结果我习惯打成一张表跟参数一起归档检核项指标合格线实测示例求解集残差 RMS米≤ 0.050.021求解集最大残差米≤ 0.100.043独立检核中误差米≤ 0.080.035独立检核最大偏差米≤ 0.150.067边界点最大偏移米≤ 0.200.112表里的合格线是我在几个常规测区里用出来的经验值跨基准转换能做到 5 厘米以内的中误差已经算不错硬要压到 2 厘米以下得先确认公共点本身的精度撑不撑得住。还有一个具体技巧值得单独说把整条链路写成pipeline字符串交给pyproj一次性跑比手写三步更不容易错也方便跟别人交换参数。from pyproj import Transformer def build_pipeline(p): 用 PROJ 管线把 西安80大地 - CGCS2000大地 一步串起来 return ( projpipeline step projcart ellpsIAU76 # 西安80 大地 - 空间直角 fstep projhelmert x{p[dx]} y{p[dy]} z{p[dz]} frx{p[rx]} ry{p[ry]} rz{p[rz]} s{p[s]} conventionposition_vector # 与求解约定保持一致 step inv projcart ellpsGRS80 # 空间直角 - CGCS2000 大地 ) tr Transformer.from_pipeline(build_pipeline(params)) lon2, lat2, h2 tr.transform(lon1, lat1, h1)管线里的convention必须跟求解时用的约定一致ellps必须一头IAU76一头GRS80顺序不能反。这条管线的好处是可移植把参数填进去就能交给别人复现不用带上整份 Python 代码。从那以后我每次交坐标转换成果都会强制走一遍参数打印、残差体检、独立检核、边界抽点这四步哪怕测区只有几个点。因为坐标转换最坑的地方不是算不出来是算错了它不报错。希望帮到你。本文还有配套的精品资源点击获取
