54系到CGCS2000坐标转换：七参数、四参数与高斯投影-平芜编程栈

简介：面向GIS数据处理与测绘定位开发者的坐标转换工具集，聚焦大地坐标与平面坐标的互换问题：可在BJ54、Xian80与WGS84之间完成基准转换，并借助七参数、三参数、四参数及仿射变换等模型适配不同精度需求，同时涵盖高斯投影正反算等平面换算环节。压缩包共12个文件，以5个m脚本承载核心算法、4个txt存放待转换点与参数示例、3个asv为编辑备份，整体仅9KB，轻量易读，便于逐行对照公式理解实现过程。目前已有457人学习下载。对需要核验转换参数、比对不同模型精度的读者，可直接替换示例坐标文件复算；也可作为课程设计或项目中坐标模块的参考骨架，高斯正反算与四（六）参数求解等分支划分清晰，便于按需裁剪复用。

1. 一张老图斑压在新区底图上偏了一百多米：这项目到底在解决什么

把一份旧 54 系的存量管线图直接叠到 2000 系的影像底图上，图形整体往西北挪了一百多米，形状没变、边长没变，就是位置不对——这是大地坐标转换里最常见的一类现场问题。项目把这件事拆成三块：同一个椭球下大地坐标、空间直角坐标、高斯平面坐标之间的互换；不同基准之间靠七参数、四参数搭桥；再把旧坐标系的历史成果整批升级到 2000 系的换代流程。适合手里压着一堆老图、又必须跟新技术标准对齐的测绘和 GIS 开发，也适合只想弄明白「为什么我转完还是偏」的人。

2. 项目骨架：椭球、基准、投影三层拆开，五个模块各管一段

2.1 目录结构：谁负责算，谁负责编排

拿到源码先看目录，这个项目的分层基本一眼能看懂：

gtf/ ├── gtf/ │ ├── ellipsoid.py # 椭球定义与派生量：b、e²、e'² │ ├── coords.py # 大地坐标 ↔ 空间直角坐标 │ ├── datum.py # 七参数 / 四参数 / 三参数模型 │ ├── projection.py # 高斯-克吕格正反算、分带换算 │ ├── grid.py # 格网残差改正 │ ├── pipeline.py # 转换链路编排，对外只暴露 convert() │ └── cli.py # 批量转换命令行入口 ├── data/ │ ├── params.yaml # 坐标系定义 + 转换参数 │ └── sample.csv # 十来个点，够跑通链路 ├── tests/ │ └── test_roundtrip.py └── requirements.txt

ellipsoid.py只存常量和一个Ellipsoid数据类，不参与任何坐标计算，它的职责是把 a、1/f 这些「人填的」参数派生成 b、e²、e'² 这些「算的时候要用的」参数，避免每个模块各算一遍、各有各的精度。coords.py只认椭球，不认基准，给它 B、L、H 和椭球，它就返回 X、Y、Z，反过来也一样。这一层刻意做成无状态的纯函数，好处是单元测试特别好写，随便挑一组已知值就能对。

datum.py是唯一知道「参数」这个概念的地方，七参数、四参数、三参数的系数都从这里进。projection.py只管把大地坐标压到平面上，中央子午线、假东、带号前缀都在这里处理。pipeline.py是编排层，它按顺序调用前面三层，把「一串点从 A 系到 B 系」这件事封装成一个函数。cli.py只做参数解析和文件读写，一行坐标算法都不写——这条边界守住，后面换数据格式（CSV 换 Shapefile、换数据库游标）根本不用碰核心代码。

2.2 三层解耦的理由：换椭球不等于换基准

测绘里最容易混的三个词就是椭球、基准、投影。椭球是个纯几何体，两个数（长半轴、扁率）就定死了；基准是「椭球 + 定位 + 定向 + 尺度」，说白了是这个椭球摆在什么位置、朝哪个方向、多大；投影是把球面上的东西摊平到纸面上的数学规则。三件事互相独立，但代码里一旦写混，bug 会以最恶心的方式出现——结果看着「差不多对」，就差那么几十米。

我见过最常见的错法，是把 WGS84 到 2000 系也套一组七参数。这两个坐标系都是地心系、椭球参数差不到 1 毫米，实际差异在厘米级，套参数属于「用药过量」，反而把本来干净的成果搅浑。项目里对这种情况的处理很干脆：在params.yaml里给同一椭球、同一基准的两套坐标定义一个identity转换器，链路照样走，只是中间那步什么都不做。这样上层调用方不用改代码，也不会有人在里面偷偷塞参数。

三层在代码里的体现就是下面这几行引用，很朴素，但把依赖方向钉死了：

# gtf/pipeline.py 顶部：三个层各自独立，只有编排层同时引用它们 from .ellipsoid import ELLIPSOIDS from .datum import DatumShift from .projection import gauss_forward, gauss_inverse CRS = { # 大地坐标形式：只要椭球 "old54_3d_39": dict(ell=ELLIPSOIDS["krasovsky"], kind="geodetic"), "cgcs2000_3d_39": dict(ell=ELLIPSOIDS["cgcs2000"], kind="geodetic"), # 平面坐标形式：椭球 + 投影参数 "cgcs2000_gk39": dict(ell=ELLIPSOIDS["cgcs2000"], kind="projected", lon0=117.0, fe=500000.0, zone_prefix=39), }

kind字段决定了链路里要不要插一步投影正算或者反算，lon0决定中央子午线，zone_prefix决定东坐标前面要不要叠带号。把这三样写成配置而不是写死进函数，是为了后面加坐标系时只改 YAML——这一点在第 6 章还会再讲。

2.3 环境与依赖：numpy 够用，pyproj 只当裁判

跑起来需要的依赖很少：Python 3.9 以上、numpy、pyyaml，测试里额外装一个 pyproj。核心算法全部手写，不依赖任何投影库，原因有三个。第一是批量场景要可控，几百万个点跑一晚上，中间出问题得知道是哪一步、哪个点；外部库一旦抛异常，堆栈里全是它的内部细节，排查成本很高。第二是很多转换要拿中间量，比如先看空间直角坐标的增量是不是合理，外部库只给你最终结果，中间过程是个黑匣子。第三是区域格网改正这类东西，公共库不一定支持你要的那套参数。

pyproj 只出现在测试里，作用是当裁判：同一组点在两边各转一遍，差值超过阈值就报错。这种做法我很推荐，自己写的投影正反算最怕的就是「公式抄错一个系数」，肉眼看不出来，但跟权威实现一比，几厘米的差立刻现形。

python -m venv .venv && source .venv/bin/activate pip install -r requirements.txt # numpy、pyyaml pip install pyproj # 只在跑测试时需要 pytest tests/ -q

requirements.txt里把 numpy 的版本卡在大版本号上就行，这个项目对 numpy 的用法很浅，只是拿它做批量数组运算，不涉及什么冷门 API，卡太死反而在别人机器上装不上。用python -m gtf.cli而不是直接python cli.py启动，是为了让包内相对导入正常工作，这点新手特别容易翻车，直接跑文件会报ImportError。

3. 核心链路：一次大地坐标转换从角度到平面坐标的完整走法

3.1 椭球参数与角度单位：第一道关

所有转换的起点都是一个Ellipsoid。这个类写得很短，但每个属性都是后面反复要用的：

# gtf/ellipsoid.py from dataclasses import dataclass @dataclass(frozen=True) class Ellipsoid: name: str a: float # 长半轴，米 inv_f: float # 扁率倒数 1/f @property def f(self) -> float: return 1.0 / self.inv_f @property def b(self) -> float: # 短半轴 return self.a * (1.0 - self.f) @property def e2(self) -> float: # 第一偏心率平方：算卯酉圈曲率半径 N 时用 return 2.0 * self.f - self.f ** 2 @property def ep2(self) -> float: # 第二偏心率平方 e'²：高斯投影级数里出现得比 e² 还频繁 return self.e2 / (1.0 - self.e2) ELLIPSOIDS = { "krasovsky": Ellipsoid("旧 54 系椭球", 6378245.0, 298.3), "iag75": Ellipsoid("旧 80 系椭球", 6378140.0, 298.257), "cgcs2000": Ellipsoid("2000 大地坐标系椭球", 6378137.0, 298.257222101), "wgs84": Ellipsoid("WGS84", 6378137.0, 298.257223563), }

e2用2f - f²而不是(a²-b²)/a²来算，是因为扁率是精确已知量，直接推比先算 b 再回头算 e² 少一次浮点舍入。ep2看着像个可有可无的派生量，但高斯投影正算的级数里，η² = e'²cos²B 这一项决定了 6 度带边缘的精度，漏掉它边缘点能差十几厘米，属于典型的「公式抄一半」。

角度单位是这一层真正的坑。原始数据里经纬度经常是116°23'29.5"这种字符串，或者116.23295这种把度分秒硬拼成小数的伪十进制。项目里统一在入口处解析成十进制度，再在计算函数内部转弧度：

# gtf/coords.py import math def dms_to_deg(s: str) -> float: """'116°23\'29.5"' -> 116.391527... 负号只作用在「度」上，南纬西经要显式判断，不能整体取负""" s = s.strip() neg = s.startswith("-") t = (s.lstrip("-").replace("°", " ") .replace("'", " ") .replace('"', " ")) d, m, sec = (float(x) for x in t.split()) v = d + m / 60.0 + sec / 3600.0 return -v if neg else v def deg_to_rad(x: float) -> float: return x * math.pi / 180.0

dms_to_deg里那个「负号只作用在度上」是血泪经验。西经 75°30' 写成-75.5是对的，但如果解析时把整个字符串取负再按度分秒展开，-75°30'会被算成-(75 + 30/60) = -75.5，其实也对——真正出问题的是-0°30'这种，度数为 0，符号丢了，结果变成正东经半小时。所以判断符号必须看原始字符串的开头，不能等算完再补。

3.2 大地坐标与空间直角坐标的互化

这一对换算是整个链路的枢纽。基准转换（七参数）必须在空间直角坐标下做，投影正算必须在大地坐标下做，所以任何跨基准的转换都要先上到 XYZ 再下来：

# gtf/coords.py（续） def geodetic_to_ecef(B, L, H, ell): """B、L 单位度，H 单位米；返回 X、Y、Z（米）""" B = deg_to_rad(B) L = deg_to_rad(L) sinB, cosB = math.sin(B), math.cos(B) # 卯酉圈曲率半径 N = ell.a / math.sqrt(1.0 - ell.e2 * sinB * sinB) X = (N + H) * cosB * math.cos(L) Y = (N + H) * cosB * math.sin(L) Z = (N * (1.0 - ell.e2) + H) * sinB return X, Y, Z

正算没有迭代，一行一个问题。反算就有讲究了，因为有闭式解但精度不稳，工程上普遍用迭代：

def ecef_to_geodetic(X, Y, Z, ell, tol=1e-12, max_iter=10): """返回 (B, L, H)，B/L 单位度，H 单位米""" p = math.hypot(X, Y) L = math.atan2(Y, X) # 极点保护：p 趋近 0 时 cosB 也趋近 0，下面除 cosB 会炸 if p < 1e-9: b = ell.b B = math.pi / 2 if Z >= 0 else -math.pi / 2 return math.degrees(B), math.degrees(L), abs(Z) - b B = math.atan2(Z, p * (1.0 - ell.e2)) # 初值，误差在米级以内足够 for _ in range(max_iter): sinB = math.sin(B) N = ell.a / math.sqrt(1.0 - ell.e2 * sinB * sinB) H = p / math.cos(B) - N B_new = math.atan2(Z, p * (1.0 - ell.e2 * N / (N + H))) if abs(B_new - B) < tol: B = B_new break B = B_new return math.degrees(B), math.degrees(L), H

迭代 4 到 5 次就能收敛到 1e-12 弧度，换算到地面是微米级，max_iter=10是保险。初值那行atan2(Z, p*(1-e²))不要写成atan2(Z, p)，后者在高纬度收敛慢，虽然多迭代两次也能出结果，但在批量场景里会白白多花时间。

极点保护那段看着多余，实际上非常必要。有一批高纬度的数据，H 填的是 0 或者缺测值 -9999，算出来的 Z 很大、p 很小，H = p/cos(B) - N里 cos(B) 接近 0，H 直接变成 1e15 这种数，然后这个污染值继续往下传，整条链路的结果全废，日志里还看不出来是哪一步出的问题。

提示：写往返测试时不要只测中纬度。挑一个纬度 0°、一个 89.9°、一个经度跨 ±180° 的点，这三个位置是几乎所有实现的分支边界。

3.3 七参数与四参数：基准转换的两种口径

七参数是跨基准的主力模型，本质是把一个空间直角坐标做一次「平移 + 旋转 + 缩放」。项目里用的是布尔莎模型，写成矩阵是 X₂ = ΔX + k·R·X₁，代码里展开成标量，比调矩阵库快得多，也方便逐个参数打日志：

# gtf/datum.py import math def apply_7param(xyz, p): """xyz: (X, Y, Z) 米；p: 参数字典 旋转角单位是角秒，尺度单位是 ppm，这是行业里最常见的写法""" X, Y, Z = xyz k = 1.0 + p["scale_ppm"] * 1e-6 rx = math.radians(p["rx"] / 3600.0) ry = math.radians(p["ry"] / 3600.0) rz = math.radians(p["rz"] / 3600.0) # 小角度近似：cos≈1，sin≈角本身，此项误差在 1e-11 量级，可忽略 X2 = p["dx"] + k * (X - rz * Y + ry * Z) Y2 = p["dy"] + k * (rz * X + Y - rx * Z) Z2 = p["dz"] + k * (-ry * X + rx * Y + Z) return X2, Y2, Z2

参数单位这里必须写死并且写进配置文件的注释。旋转角用角秒、尺度用 ppm 是国内测绘参数的通行写法，但有些来源会直接给弧度和小数（0.0000012 而不是 1.2），两者差 2 万倍。我一般会在加载参数时做一次量级校验：如果abs(rx) > 60，几乎肯定是单位写错了，直接报错而不是硬算。

方向问题同样要命。同一组参数，从旧系到新系用正号，从新系回旧系就得取逆。严格逆要解矩阵，工程上因为旋转角和尺度都极小，直接取反号做近似逆就够了，残差在亚毫米级：

def invert_7param(p): """近似逆：平移取负并反向旋转，旋转、尺度取反号""" k = 1.0 + p["scale_ppm"] * 1e-6 return { "dx": -p["dx"] / k, "dy": -p["dy"] / k, "dz": -p["dz"] / k, "rx": -p["rx"], "ry": -p["ry"], "rz": -p["rz"], "scale_ppm": -p["scale_ppm"] / (1.0 + p["scale_ppm"] * 1e-6), }

四参数是平面的简化版，只有两个平移、一个旋转、一个尺度，用在高斯平面上做小范围（几十公里内）的坐标纠正。它的好处是不需要高程，坏处是一旦范围大了，地球曲率带来的误差就盖不住了，几十公里外能到分米级。所以params.yaml里给每个转换都标了model: bursa7 | helmert4和extent_km，超过范围加载时就打警告。这条规则救过我一次：有人把一组四参数用在跨两个地级市范围的数据上，跑完发现边缘偏了半米多，翻配置才发现extent_km明明写着 30。

3.4 高斯-克吕格投影正反算与分带

投影这块是整个项目里公式最长、最容易抄错的地方。核心是先把子午线弧长算出来，再往上叠加经差的高阶级数：

# gtf/projection.py import math from .coords import deg_to_rad def meridian_arc(B, ell): """从赤道到纬度 B 的子午线弧长，B 为弧度""" e2 = ell.e2 A0 = 1 + 3/4*e2 + 45/64*e2**2 + 175/256*e2**3 + 11025/16384*e2**4 A2 = 3/4*e2 + 15/16*e2**2 + 525/512*e2**3 + 2205/2048*e2**4 A4 = 15/64*e2**2 + 105/256*e2**3 + 2205/4096*e2**4 A6 = 35/512*e2**3 + 315/2048*e2**4 A8 = 315/16384*e2**4 return ell.a * (1 - e2) * ( A0 * B - A2 / 2 * math.sin(2*B) + A4 / 4 * math.sin(4*B) - A6 / 6 * math.sin(6*B) + A8 / 8 * math.sin(8*B) ) def gauss_forward(B, L, ell, lon0, fe=500000.0): """B、L、lon0 单位度；返回测绘习惯的 (x北, y东)""" B = deg_to_rad(B) l = deg_to_rad(L - lon0) # 经差，必须归一化到 ±3.5° 内 sinB, cosB = math.sin(B), math.cos(B) t = math.tan(B) t2 = t * t eta2 = ell.ep2 * cosB * cosB N = ell.a / math.sqrt(1 - ell.e2 * sinB * sinB) X = meridian_arc(B, ell) x = (X + N/2 * sinB * cosB * l**2 + N/24 * sinB * cosB**3 * (5 - t2 + 9*eta2 + 4*eta2**2) * l**4 + N/720 * sinB * cosB**5 * (61 - 58*t2 + t2**2) * l**6) y = (fe + N * cosB * l + N/6 * cosB**3 * (1 - t2 + eta2) * l**3 + N/120 * cosB**5 * (5 - 18*t2 + t2**2 + 14*eta2 - 58*eta2*t2) * l**5) return x, y

这里有两个「不写注释就会忘」的点。第一，测绘习惯的 x 是北方向、y 是东方向，和 GIS 里 x=经度、y=纬度完全相反，函数返回值一定要在文档字符串里写清楚，不然调用方一个(x, y)直接塞进几何库，图形会转 90 度。第二，l = L - lon0必须归一化，如果输入经度是 118.5、中央子午线是 117，差 1.5 度没问题；但如果数据里有 359 度这种表达，差值变成 242 度，级数展开直接发散，结果是天文数字。归一化到(-180, 180]还不够，得再收敛到(-3.5, 3.5]才能保证 6 度带内的级数精度。

反算走的是「由 x 反求底点纬度，再用 y 展开」的路子，底点纬度用牛顿迭代，比查表法干净：

def gauss_inverse(x, y, ell, lon0, fe=500000.0): """输入测绘习惯的 (x北, y东)，返回 (B, L) 十进制度""" y = y - fe # 1) 牛顿迭代求底点纬度 Bf，使 meridian_arc(Bf) == x Bf = x / (ell.a * (1 - ell.e2)) for _ in range(8): M = ell.a * (1 - ell.e2) / (1 - ell.e2 * math.sin(Bf)**2) ** 1.5 d = (x - meridian_arc(Bf, ell)) / M Bf += d if abs(d) < 1e-12: break # 2) 由 Bf 和 y 展开修正 sinBf, cosBf, tf = math.sin(Bf), math.cos(Bf), math.tan(Bf) tf2 = tf * tf etaf2 = ell.ep2 * cosBf**2 Nf = ell.a / math.sqrt(1 - ell.e2 * sinBf**2) Mf = ell.a * (1 - ell.e2) / (1 - ell.e2 * sinBf**2) ** 1.5 B = (Bf - tf/(2*Mf*Nf) * y**2 + tf/(24*Mf*Nf**3) * (5 + 3*tf2 + etaf2 - 9*etaf2*tf2) * y**4 - tf/(720*Mf*Nf**5) * (61 + 90*tf2 + 45*tf2**2) * y**6) l = (y/(Nf*cosBf) - 1/(6*Nf**3*cosBf) * (1 + 2*tf2 + etaf2) * y**3 + 1/(120*Nf**5*cosBf) * (5 + 28*tf2 + 24*tf2**2 + 6*etaf2 + <p> <a href="https://download.csdn.net/download/feifei199112/8601111" style="color:#ec7500;font-size:14px;"> 本文还有配套的精品资源，点击获取 </a> <img alt="menu-r.4af5f7ec.gif" src="https://csdnimg.cn/release/wenkucmsfe/public/img/menu-r.4af5f7ec.gif" style="width:16px;margin-left:4px;vertical-align:text-bottom;cursor:text;"> </p>