简介:针对经纬度与大地坐标互相转换的实际需求,这份VC++工程源码为GIS、测绘及导航方向开发者提供了可直接运行的示例程序。资源内含完整对话框应用,涵盖经纬度坐标系与大地投影坐标系的转换算法实现,包括椭球参数设置、投影正反算等核心逻辑,适合需要理解坐标换算原理或快速集成到项目中的初中级开发者参考。压缩包共18个文件,以4个h头文件、3个cpp源文件为主,另含工程配置、资源脚本及ReadMe说明等配套文件,整体仅29KB,结构轻量清晰。已有3396人学习下载。通过阅读源码,可掌握从经纬度到平面直角坐标的正算过程以及由大地坐标反解经纬度的逆算过程,理解UTM、北京54、西安80等坐标系背后的数学变换,并能直接复用工程模板进行二次开发。
1. 经纬度和大地坐标互相转换:先搞清这套计算里的四个隐藏前提
把手机 GPS 读出来的经纬度换算成施工图纸上的「大地坐标」,是测绘、GIS、无人机航测里几乎天天都要做的操作。但互转不是套公式那么简单:椭球参数、投影带、中央子午线、500 公里偏移量,任何一个选错,结果偏出几百米,代码却不报错。这套互转经常被人当黑匣子用,正算反算各贴一段代码,跑通就不管了。
这里把「大地坐标」约定为高斯-克吕格投影后的平面坐标(北向 x、东向 y),也就是图纸上常见的那串米制数字。经纬度是球面角,平面坐标是投影后的距离,两者互相转换的核心就是投影正算和反算。
接下来从投影原理讲到可运行的 Python 正反算代码,再讲跨基准转换怎么选型、五个高频翻车点,以及一套能反复使用的精度验证方法。新手可以直接抄代码,熟手可以重点看第 4、5 章的边界和参数,踩坑大多集中在这些看起来不起眼的环节。
2. 经纬度转大地坐标:高斯正算公式与最小可跑的 Python 代码
2.1 为什么经纬度不能当米用:投影才是核心
经纬度本质是球面坐标,纬度是角、经度也是角。纬度一度对应约 111 公里,但经度一度对应的地面距离会随纬度变化:赤道约 111 公里,到纬度 60 度只剩约 55 公里。如果直接把经纬度当平面坐标用,距离、面积、方位角全都会失真,比例尺稍微大一点就完全没法施工。
所以工程上要把球面拉平,这就是地图投影。中国和很多国家采用高斯-克吕格投影,它属于等角横切椭圆柱投影,基本思路是:用一个椭圆柱横着套在地球上,让椭圆柱面与某一条经线相切,这条经线就是中央子午线,把中央子午线附近的经纬度一一映射到平面上。等角投影的好处是形状不扭曲、角度关系不变,适合大比例尺测图和工程放样。
投影必然带来长度变形,离中央子午线越远,变形越大。为了把变形控制在毫米级或厘米级,就用分带投影:每隔 3 度或 6 度换一条中央子午线。这就是为什么工程坐标里总带着「带号」这个信息——不告诉你带号,光给一堆坐标数字,外人根本不知道这个点在地球上哪一带。
2.2 高斯正算分步公式:从 B、L 到 x、y
高斯正算的任务是:已知大地纬度 B、大地经度 L、中央子午线经度 L0,求平面坐标 x(北向)、y(东向)。注意在测量惯例里 x 是北向、y 是东向,和数学坐标系相反,代码里我会保持这个约定,避免和图纸对不上。
先算经差 l = L - L0,单位必须是弧度。随后需要几个辅助量:
- 第一偏心率平方 e² = f(2 - f),其中 f 是椭球扁率
- 第二偏心率平方 η² = e'² cos²B,其中 e'² = e² / (1 - e²)
- 卯酉圈曲率半径 N = a / √(1 - e² sin²B)
- 子午线弧长 X,即从赤道沿子午线到纬度 B 的椭球弧长
X 没有初等闭式,工程上常用幂级数展开:
X = a(1 - e²)(A0·B + A2·sin2B + A4·sin4B + A6·sin6B + A8·sin8B)
其中 A0、A2、A4、A6、A8 是只与 e² 有关的系数,代码里动态算,换成任意椭球都能直接用。
有了这些辅助量,正算主项是:
x = X + N/2·sinB·cosB·l² + N/24·sinB·cos³B·(5 - t² + 9η² + 4η⁴)·l⁴ + N/720·sinB·cos⁵B·(61 - 58t² + t⁴ + 270η² - 330η²t²)·l⁶
y = N·cosB·l + N/6·cos³B·(1 - t² + η²)·l³ + N/120·cos⁵B·(5 - 18t² + t⁴ + 14η² - 58η²t²)·l⁵
这里的 t = tanB。公式前面还省略了 l⁸ 级的高阶项,对正常经差(3 度带内 l 不超过 1.5 度)来说,保留到 l⁶ 已经能把误差压到毫米级以下。
这段公式初看眼花,但核心逻辑不复杂:x 的主体是子午线弧长,再加上经差 l 的偶次项修正;y 的主体是 N·cosB·l,再加上奇次项修正。
2.3 Python 正算实现:函数、单位与返回值
下面这份代码我平时写项目也这样写,唯一要提醒的是单位:函数外面用度,函数内部全部转弧度,返回值是米。
import math def meridian_arc_length(B, a, e2): """子午线弧长:从赤道到纬度 B 的椭球弧长,B 为弧度""" e4 = e2 * e2 e6 = e4 * e2 e8 = e6 * e2 A0 = 1 + 3/4*e2 + 45/64*e4 + 175/256*e6 + 11025/16384*e8 A2 = -3/4*e2 - 15/16*e4 - 525/512*e6 - 2205/2048*e8 A4 = 15/64*e4 + 105/256*e6 + 2205/4096*e8 A6 = -35/512*e6 - 315/2048*e8 A8 = 315/16384*e8 return a * (1 - e2) * (A0*B + A2*math.sin(2*B) + A4*math.sin(4*B) + A6*math.sin(6*B) + A8*math.sin(8*B)) def geo_to_gauss(B_deg, L_deg, L0_deg, a=6378137.0, f=1/298.257223563): """经纬度 -> 高斯平面坐标自然值 B_deg: 纬度(度),北纬为正 L_deg: 经度(度),东经为正 L0_deg: 中央子午线经度(度) 返回: (x, y) = (北向, 东向),单位米 """ B = math.radians(B_deg) L = math.radians(L_deg) L0 = math.radians(L0_deg) l = L - L0 e2 = f * (2 - f) # 第一偏心率平方 e2_2 = e2 / (1 - e2) # 第二偏心率平方 sinB = math.sin(B) cosB = math.cos(B) tanB = math.tan(B) eta2 = e2_2 * cosB * cosB N = a / math.sqrt(1 - e2 * sinB * sinB) # 卯酉圈曲率半径 X = meridian_arc_length(B, a, e2) # 子午线弧长 # 高斯正算主项,x 为北向,y 为东向 x = (X + N/2 * sinB * cosB * l**2 + N/24 * sinB * cosB**3 * (5 - tanB**2 + 9*eta2 + 4*eta2**2) * l**4 + N/720 * sinB * cosB**5 * (61 - 58*tanB**2 + tanB**4 + 270*eta2 - 330*eta2*tanB**2) * l**6) y = (N * cosB * l + N/6 * cosB**3 * (1 - tanB**2 + eta2) * l**3 + N/120 * cosB**5 * (5 - 18*tanB**2 + tanB**4 + 14*eta2 - 58*eta2*tanB**2) * l**5) return x, y说明一下这段代码里的关键点。子午线弧长用了 8 次幂前的全部主项,A2 到 A8 按正负交替写进一个表达式里,A0 是常数项,这样在 e² 接近 0 时 X 会退化成 a·B,符合几何直觉。第二偏心率单独算,是为了下面的反算复用。
参数部分:a 默认取 WGS84 的长半轴 6378137.0,f 默认取 1/298.257223563。如果你做的是 2000 国家大地坐标系,只需要把 f 改成 1/298.257222101;做西安 80 就改成 a=6378140.0、f=1/298.257。函数返回的是自然值 y,还没有加 500 公里偏移,也没有带带号前缀,这是刻意的:内部计算一律用自然值,显示时才加工。
调用示例:
x, y = geo_to_gauss(30.5, 114.3333, 114.0) print(x, y)这个例子对应 3 度带 38 带、中央子午线 114 度附近的一个点。输出的是北向和东向自然值,如果要图纸上的常用坐标,还要做下一节里的两步加工。
2.4 加 500 公里偏移与带号前缀:让坐标符合图纸习惯
高斯平面坐标的自然值 y 在中央子午线以西会是负数,历史上为了避免负号,所有通用高斯坐标的东向自然值都统一加 500000 米,叫假东偏移(False Easting)。同时为了让人一眼看出这个点在哪个投影带,还把带号写在 y 值最前面,比如 38527665.32,表示 3 度带第 38 带,去掉带号得到 527665.32,再减 500000 得到自然值 27665.32。
具体加工规则是这样:6 度带带号 N = floor(L / 6) + 1,中央子午线 L0 = 6N - 3;3 度带带号 N = floor((L + 1.5) / 3),中央子午线 L0 = 3N。以 L = 120 度为例,6 度带 N = 20、L0 = 117 度;3 度带 N = 40、L0 = 120 度。同一组经纬度,落在 3 度带和 6 度带里的平面坐标差异可能上百公里,因为它们根本不是同一套网格。
所以正算之后要不要加偏移、加哪个带号,取决于接你数据的下游单位用什么规范。我的习惯是:函数只返回自然值,加工逻辑放到独立的展示层,避免在核心算法里塞业务规则,这样换项目时只需要改展示层配置。
3. 大地坐标转经纬度:反算迭代与收敛判据
3.1 反算难点:底点纬度的牛顿迭代
从 x、y 反推 B、L,不能直接反解正算的级数式。观察正算公式:x 的第一大项是子午线弧长 X,它只和 B 有关,其他都是经差 l 的小修正。反算的常规思路是先不管 l,用 x 解出一个中间纬度,叫底点纬度 Bf,再用 y 修正经差、用底点纬度附近的曲率修正纬度。
Bf 的定义是:如果从赤道向正北量一段子午线弧长恰好等于 x,那么这段弧长对应的纬度就是 Bf。问题是子午线弧长公式是 B 的复杂三角级数,没法直接反解,只能用数值迭代。最常见的牛顿迭代形式是:
Bf_new = Bf_old - (X(Bf_old) - x) / (dX/dB)
其中 dX/dB 是子午线弧长对纬度的导数,可以由弧长系数直接求导得到。初值取 x 除以赤道处的弧长导数 a(1-e²)·A0,通常两三次迭代就能收敛到 1e-12 弧度,对应精度在微米级。这个迭代很稳,因为 dX/dB 在 0 到 90 度范围内恒大于零,不会出现梯度为零的情况。这一步别偷懒换成固定迭代次数,收敛判据必须写,否则换椭球或换投影带时精度会莫名下降。
3.2 Python 反算实现:从 x、y 回到 B、L
反算代码继续复用前面的子午线弧长函数,同时需要新增导数函数。
def meridian_arc_derivative(B, a, e2): """子午线弧长对纬度的导数 dX/dB,B 为弧度""" e4 = e2 * e2 e6 = e4 * e2 e8 = e6 * e2 A0 = 1 + 3/4*e2 + 45/64*e4 + 175/256*e6 + 11025/16384*e8 A2 = -3/4*e2 - 15/16*e4 - 525/512*e6 - 2205/2048*e8 A4 = 15/64*e4 + 105/256*e6 + 2205/4096*e8 A6 = -35/512*e6 - 315/2048*e8 A8 = 315/16384*e8 return a * (1 - e2) * (A0 + 2*A2*math.cos(2*B) + 4*A4*math.cos(4*B) + 6*A6*math.cos(6*B) + 8*A8*math.cos(8*B)) def gauss_to_geo(x, y, L0_deg, a=6378137.0, f=1/298.257223563): """高斯平面坐标自然值 -> 经纬度 x: 北向(米) y: 东向自然值(米),注意:如果输入是带偏移和带号的坐标,先还原 L0_deg: 中央子午线经度(度) 返回: (B, L),单位度 """ e2 = f * (2 - f) e2_2 = e2 / (1 - e2) # 底点纬度初值:x 除以赤道处 dX/dB,再牛顿迭代 A0 = 1 + 3/4*e2 + 45/64*e2*e2 + 175/256*e2**3 + 11025/16384*e2**4 Bf = x / (a * (1 - e2) * A0) for _ in range(20): dX = meridian_arc_length(Bf, a, e2) - x ddB = meridian_arc_derivative(Bf, a, e2) Bf_new = Bf - dX / ddB if abs(Bf_new - Bf) < 1e-12: Bf = Bf_new break Bf = Bf_new sinBf = math.sin(Bf) cosBf = math.cos(Bf) tanBf = math.tan(Bf) etaf2 = e2_2 * cosBf * cosBf Nf = a / math.sqrt(1 - e2 * sinBf * sinBf) Mf = a * (1 - e2) / (1 - e2 * sinBf * sinBf) ** 1.5 y0 = y # 输入必须是自然值,偏移在函数外还原 B = (Bf - tanBf / (2 * Mf * Nf) * y0**2 + tanBf / (24 * Mf * Nf**3) * (5 + 3*tanBf**2 + etaf2 - 9*etaf2*tanBf**2) * y0**4 - tanBf / (720 * Mf * Nf**5) * (61 + 90*tanBf**2 + 45*tanBf**4) * y0**6) l = (y0 / (Nf * cosBf) - (1 + 2*tanBf**2 + etaf2) * y0**3 / (6 * Nf**3 * cosBf) + (5 + 28*tanBf**2 + 24*tanBf**4 + 6*etaf2 + 8*etaf2*tanBf**2) * y0**5 / (120 * Nf**5 * cosBf)) L = math.radians(L0_deg) + l return math.degrees(B), math.degrees(L)反算的核心是底点纬度迭代。初值我用赤道处的弧长导数直接除,比随便给一个中间纬度收敛快;迭代里如果前后两次的差小于 1e-12 弧度就退出,正常两三步就能满足,20 次上限只是保险。
后面的纬度修正公式里,Mf 是子午圈曲率半径,Nf 是卯酉圈曲率半径,都取底点纬度处的值。要注意 y0 必须是自然值:如果输入的是 38527665.32,你得先还原成 27665.32 再传进来,否则反算出的经度会偏到离谱。这个前置处理放到下一节专门讲。
代码里 A0 的计算重复了一遍,如果你不想写两遍,可以抽一个公共函数返回 A0 到 A8 系数,正算和反算共用。我的习惯是两个函数都能独立运行,避免改系数时相互影响。
3.3 反算前置处理:带号、偏移和中央子午线还原
反算最常出问题的不是公式,而是喂进来的坐标带了一堆附加信息没人拆。我见过最典型的错误:输入一个 38527665.32 的 y 值,中央子午线填 114,反算出来经度差了将近 6 度,因为代码把 38527665 当作自然值参与运算,相当于凭空加了 38527 公里。
正确的还原顺序是:先把 y 最前面的带号去掉,如果是 3 度带第 38 带,去掉 38 后还剩 527665.32;再减 500000,得到自然值 27665.32。然后根据带号确定中央子午线:3 度带 L0 = 3N,6 度带 L0 = 6N - 3。如果数据没带带号,就必须由业务方明确告知中央子午线,否则无解。
还有一个细节:北向 x 也可能有加假北(False Northing)的情况。南半球项目常见 y 从 10000 公里起算,这时要先把 x 减去 10000000 再进入反算函数。国内项目一般不加假北,x 就是从赤道量起的真实北向,但如果对接海外项目,这两条都要提前问清楚。
4. 坐标系基准与椭球参数:同样的经纬度,换套参数就偏差几百米
4.1 常见椭球参数对照:不同基准的经纬度本身就不一样
同一个物理点,在不同参考椭球下测量出来的经度、纬度是不一样的,因为椭球大小、定位和定向都不同。很多初学者以为「经纬度是绝对的」,其实经纬度只是相对于某个椭球的描述。下面这几个是国内和全球项目里最常见的椭球。
| 坐标系 | 椭球 | 长半轴 a (m) | 扁率分母 1/f | 备注 |
|---|---|---|---|---|
| WGS84 | WGS84 | 6378137.0 | 298.257223563 | GPS 原始输出,全球通用 |
| 2000国家大地坐标系 | CGCS2000 | 6378137.0 | 298.257222101 | 现行国家坐标框架 |
| 西安80 | IAG 75 | 6378140.0 | 298.257 | 老规范常用,存量数据多 |
| 北京54 | 克拉索夫斯基 | 6378245.0 | 298.3 | 历史遗留,已逐步淘汰 |
长半轴彼此只差几米到几十米,但如果把一套经纬度用错椭球的参数去投影,平面坐标能偏出去几百米。更隐蔽的是,即使 a 和 f 都用对了,如果输入的经纬度本来就在另一个椭球基准下,偏几公里都不奇怪。所以第一步永远是确认你的经纬度属于哪个基准。
2000 国家大地坐标系和 WGS84 在厘米级可以认为是兼容的,很多手持设备直接按 WGS84 处理问题不大。但西安80 和北京54 与 WGS84 之间没有简单换算关系,必须靠区域转换参数。存量图纸里这两种坐标非常常见,接手时先看数据字典里「坐标系」字段怎么写,别信文件名里的约定。
4.2 基准转换:四参数、七参数和「先转经纬度」的边界
如果你的经纬度是 WGS84 基准,而目标图纸是西安80 或 CGCS2000,直接投影是不够的,中间要插一次基准转换。常见有三种做法。
四参数转换适合小范围平面:拿到若干重合点,解出 2 个平移、1 个旋转、1 个尺度,公式简单,但只适用于平面坐标到平面坐标,范围一大就失真。它的前提是两个平面坐标系之间的投影参数几乎一致,否则转完也不准。
七参数转换(布尔沙模型)适合三维空间:把经纬度先转成空间直角坐标 X、Y、Z,再用 3 个平移、3 个旋转、1 个尺度,转到目标椭球的空间直角坐标,最后再转回目标椭球的经纬度。这才是跨大区域、跨基准转换的正路。
常见工作流是:WGS84 经纬度 -> WGS84 空间直角 -> 七参数 -> CGCS2000 空间直角 -> CGCS2000 经纬度 -> 高斯投影。这里的七参数通常是用测区里 3 个以上控制点解算出来的,不同地区的七参数不通用,不能拿 A 城市的参数套到 B 城市。解算后至少要检查控制点残差,残差超过厘米级就说明控制点本身有矛盾,参数不可信。
「先转经纬度再转目标平面」这句话只对了一半:如果你的经纬度和目标平面坐标的基准一致,比如都是 CGCS2000,那直接投影就行;如果基准不一致,就必须先做空间直角坐标系下的基准转换,而不是在经纬度上加减某个常数。经纬度是球面角,直接做平移毫无几何意义,这步做错就是第 5 章第一个坑的来源。
4.3 中央子午线和带号速查:动手前先算一遍
3 度带和 6 度带的带号规则看起来简单,实际项目里到处是坑。我一般会在需求评审阶段用一张小表把经度范围和投影带列出来,避免开发到一半才发现带号理解不一致。
| 经度范围 | 6度带带号 | 6度带中央子午线 | 3度带带号 | 3度带中央子午线 |
|---|---|---|---|---|
| 111°00' - 114°00' | 19 | 111° | 37 / 38 | 111° / 114° |
| 114°00' - 117°00' | 20 | 117° | 38 / 39 | 114° / 117° |
| 117°00' - 120°00' | 21 | 123° | 39 / 40 | 117° / 120° |
注意 6 度带跨了 3 度带的边界:同一组经纬度在 114°00' 附近,3 度带会跨两个带,应选择离中央子午线更近的那个 3 度带,而不是机械地按整除取。项目里如果图纸上标注了「38 带」但没有写是 3 度带还是 6 度带,要去查图廓或控制点资料,默认 3 度带是很危险的习惯。
中央子午线不能偷懒用「当地经度四舍五入到最近的 3 的倍数」来确定。比如经度 114.1 度,四舍五入到 3 度带可能是 114,但如果这个点归 38 带还是 39 带,取决于图幅编号和项目规范,而不是单纯数学取整。严谨做法永远是先确认带号,再由带号反推中央子午线。
5. 经纬度大地坐标互转避坑:5 个翻车现场与排查
下面这 5 个翻车场景,是我在实际项目里遇到最多、也最容易浪费一整天的问题。很多看起来像玄学,定位到最后都是参数问题。
5.1 点位整体偏出几百米:基准和椭球混用
现象:某个点用 WGS84 经纬度加 WGS84 投影参数算出来的平面坐标,和项目部给的标准坐标差 300 到 500 米,各个方向都有,不是单纯某个轴差固定值。
原因:标准坐标是 CGCS2000 或西安80 基准下的经纬度投影来的,输入的是 WGS84 经纬度,基准没有先转换。两个基准在同一个物理点上的经纬度差可以达到几十米弧长,投影后自然偏几百米。
解决:查清源数据基准和目标基准。如果基准不同,先做七参数空间基准转换,转成目标椭球的经纬度,再用目标椭球的参数投影。不要看到一个方向一致的偏差就调尺度或加常数,那只会让结果更乱。
提示:这类问题不要靠调整投影公式的系数来「校正」,根源在基准,基准转对了,偏差自然消失。
5.2 反算经度固定差 3 度:中央子午线填错
现象:正算反算互逆性测试全部通过,但和标准成果比经度差 3 度或 6 度,纬度几乎完全正确。
原因:中央子午线取错。6 度带 L0=117 的点,填成了 3 度带 L0=114 或 L0=120;或 y 坐标带号是 6 度带的 20,反算时却按 3 度带 L0=111 处理。
解决:先由带号反推 L0,再决定投影带。6 度带号 N 对应 L0 = 6N - 3,3 度带号 N 对应 L0 = 3N。最稳的做法是把带号解析出来,作为坐标属性随数据一起存,而不是每次人工输入。
5.3 y 值少 500000 或加了两次:偏移量重复处理
现象:批量转完坐标后,所有的东向坐标和标准值都差 500000 或 1000000,x 正常。
原因:正算结果没有加偏移,或反算时把已经带偏移的值又减了一次 500000。更隐蔽的是从数据库取数时,有的表存自然值、有的表存通用值,两种混在一起。
解决:内部一律使用自然值;输出平面坐标时统一加 500000 和带号;输入反算前先拆带号、再减偏移。建议在函数注释里把「输入/输出是否含偏移」写得明明白白,这个参数比公式本身更容易出错。
5.4 南半球或西经区符号处理失败
现象:海外项目布点,纬度输入 -33.9,反算出来的经度差 6 度或者纬度的绝对值不对。
原因:南半球高斯投影也可以直接用,但要注意北坐标可能加了假北;另外很多国内工具库默认北纬为正,把负纬度当无效值处理。
解决:先确认项目坐标系统的 False Northing 和 False Easting 约定。南半球常见 y 从 10000000 米起算,先还原 x 再反算。代码里不要把纬度取绝对值,应该保留符号并让函数按物理约定处理,比如gauss_to_geo(x - 10000000, y, L0),而不是在函数内部偷偷取正。
5.5 精度从厘米掉到分米:浮点类型踩坑
现象:同一份数据在本地跑互逆误差 1e-8,部署到另一台环境后误差变成零点几米,且不是系统性偏移。
原因:中间有人用了单精度浮点数存坐标。6000000 这个量级的浮点数,float32 的步长约 0.5 米,累加几次误差就到米级。
解决:坐标存储和计算统一用双精度 float64;接口传输时不要强转单精度;数据库表字段用 double precision,不要用 numeric(10,2) 存平面坐标,那会把毫米级精度直接砍掉。
6. 精度验证技巧:用已知控制点给你的互转代码校准
6.1 造一份最小已知点表,批量算残差
验证互转函数最靠谱的方法是找 3 到 5 个已知点,经纬度和平面坐标都是权威成果。把经纬度用你的正算函数转成平面坐标,和权威平面坐标求差,统计最大误差和均方根误差。误差在厘米级说明函数可用;如果误差大且有规律,直接对照第 5 章排查。
我常用的脚本框架:
points = [ # (纬度, 经度, 标准x, 标准y自然值) (30.5000, 114.3333, 3376000.0, 27665.0), (31.0000, 114.6667, 3432000.0, 56780.0), (29.8000, 113.9000, 3301000.0, -12345.0), ] for lat, lon, x_std, y_std in points: x, y = geo_to_gauss(lat, lon, 114.0) print(f"{lat:.4f} {lon:.4f} dx={x-x_std:10.4f} dy={y-y_std:10.4f}")这里的 x_std、y_std 必须统一到自然值再比较,如果标准成果本身带 500000 偏移,先减掉。标准成果一般来自控制点资料,点位分布要覆盖测区四角和中心,不要五个点全挤在一个角落。
如果 dx、dy 都小于 0.01 米,说明正算函数和参数选得对,可以放心用。如果某个方向差固定值,优先检查偏移和带号;如果差值和经差相关,检查中央子午线;如果全部点在同一个方向上差了几百米,回第 4 章检查基准。反算的验证同样简单:用正算结果喂给反算函数,看能否回到原始经纬度,互逆误差一般在 1e-8 度以内。
这套验证流程我每次接坐标转换需求都会先跑一遍,跑了这么多年,公式本身很少改,改的最多的反而是输入数据到底用的哪套参数。说白了,互转代码写法千奇百怪,但坑几乎都集中在参数边界上。希望帮到你。
本文还有配套的精品资源,点击获取