简介：面向GIS坐标转换学习与小型工程实践的MATLAB实现包，重点处理1954北京坐标系、1980西安坐标系与WGS84之间的相互转换，并涉及平面坐标系投影变换，强调从坐标输入、参数选择到结果输出的完整链路，对理解坐标系换代与投影差异亦有帮助。内容覆盖七参数、三参数、四参数和仿射变换等常见模型，同时包含高斯正算与反算思路，适合测绘、地理信息相关课程实验、算法对照及坐标换代场景。压缩包共12个文件，其中m脚本与函数5个、txt坐标数据及说明4个、asv自动备份3个，整体约9KB；m文件承担转换函数与主流程求解，txt用于存放待转换点坐标及参数说明，asv便于回溯修改过程。已有457人学习下载，可作为参数求解、转换流程组织和结果核对的参考样例，帮助读者理解不同方法的适用条件、输入格式及精度控制要点。

1. 大地坐标转换的第一次翻车：平面坐标差 80 米，问题不在投影

一批历史测量点从旧参心系迁到 CGCS2000，平面坐标最大差 80 米。投影带、中央经线、单位都查过，最后发现只换了椭球参数，没有做基准转换。大地坐标转换在这里不是“经纬度转平面”这么简单，它至少分成基准变换、投影变换、轴序与单位处理三层；不同坐标系转换及换代，真正容易出问题的是基准和参数方向。这份大地坐标转换实现包把三层拆开：参数文件管基准，转换入口管投影与轴序，校验脚本管残差。适合手里有旧坐标要迁移的测绘、GIS、无人机航迹和工程测量从业者，也适合想自己写转换模块但不想从公式推起的人。

2. 坐标转换链路拆解：基准、投影、轴序三层各自管什么

2.1 基准变换和投影变换为什么要拆开

很多人拿到一组旧平面坐标，第一反应是改中央经线，或者把椭球换掉再算一遍。这样做的结果常常是“看起来变了，但怎么都对不上控制点”。原因在于，旧参心系到 CGCS2000 的差异首先是基准差异：椭球参数不同、地心位置不同、定向不同。基准没换，投影参数再调也只是在同一个错误基准上做平面挪动。

这份实现包把转换链路拆成三段。第一段是基准变换，负责把源坐标系的地心坐标或参心坐标归算到目标基准；第二段是投影变换，负责地理坐标和平面坐标之间的正反算；第三段是轴序与单位处理，负责经纬度顺序、度分秒转十进制度、弧度与度、米与毫米。三段分开后，排查时就能定位到底是哪一段出的问题，而不是把所有参数搅在一起试。

选型上，跨基准的转换优先用七参数布尔莎模型或格网方法；同一基准、小范围的平面坐标纠正才用四参数。四参数不能替代基准变换，这是血泪经验。资源包里用datum、projection、io三个层次组织代码，参数文件里也把基准参数和投影参数分开放，目的就是避免“改一个参数影响整条链路”。

2.2 七参数布尔莎模型：单位、旋转符号与方向

七参数布尔莎模型的核心是三个平移、三个旋转、一个尺度。平移单位通常是米，旋转单位通常是角秒，尺度单位通常是 ppm。写进代码前必须统一单位：角秒转弧度要乘π / (180 × 3600)，ppm 转比例因子要乘1e-6。这两个单位如果漏掉，结果会差得离谱。

旋转符号是第二个容易翻车的地方。位置矢量法和坐标框架法的旋转矩阵符号相反，不同软件、不同文档的约定不一样。资源包的参数文件里要求显式写rotation_convention: position_vector或frame，不允许靠默认值猜。另一个必须显式写的是方向：source_to_target还是target_to_source。正算和反算不是简单把七个参数取负，严格做法是用逆矩阵或逆管线。

注意：下面示例里的七参数只用于演示字段格式，数值不能直接用于生产，必须用你的控制点重新解算。

2.3 资源包里的调用链：参数、CRS、批量入口、校验

我拿到这份实现后先看三处：参数文件、转换入口、校验脚本。参数文件决定基准和投影，转换入口决定轴序和单位，校验脚本决定结果能不能验收。调用链一般是：读取参数文件，构建源 CRS 和目标 CRS，创建转换器，单点或批量转换，最后用检查点算残差。资源包没有把逻辑绑在某个 GIS 平台上，核心转换可以单独跑，这一点对批量任务比较友好。

这里有个取舍：自己实现布尔莎公式能看清每一步，但容易在轴序、单位、反算上出错；直接用 PROJ/pyproj 能借力成熟管线，但参数文件必须写清楚，否则默认值会把你带偏。我一般会保留资源包自己的七参数实现做对照，同时用 pyproj 构建一条独立管线做交叉校验。两条链路结果差得多，先查参数方向和旋转约定，不要先怀疑数据。

2.4 单点与批量最小示例

先跑单点，确认轴序和基准变换方向。下面这段用 pyproj 构建源和目标 CRS，always_xy=True强制输入顺序为经度、纬度。

from pyproj import CRS, Transformer # 源：旧参心系，towgs84 仅演示字段顺序，数值不可用于生产 src = CRS.from_proj4( "+proj=longlat +ellps=krass " "+towgs84=-10,20,30,1.2,-0.8,0.5,2.0 +no_defs" ) # 目标：CGCS2000 地理坐标 dst = CRS.from_epsg(4490) tr = Transformer.from_crs(src, dst, always_xy=True) lon, lat = 114.123456, 30.654321 x, y = tr.transform(lon, lat) print(x, y)

这段代码只做基准变换，输出仍是经纬度。towgs84的七个值依次是 dx、dy、dz、rx、ry、rz、s，旋转单位是角秒，尺度单位是 ppm。always_xy=True很关键，不写它时不同 CRS 的轴序可能按纬度、经度解释，批量数据会整批错位。单点跑通后，再换投影目标，比如 CGCS2000 三度带投影EPSG:4547，看平面坐标是否落在合理范围。

批量转换时不要逐行创建Transformer，复用同一个转换器，并且把异常坐标隔离出来。

import csv from pyproj import Transformer tr = Transformer.from_crs("EPSG:4490", "EPSG:4547", always_xy=True) with open("in.csv", newline="", encoding="utf-8") as f, \ open("out.csv", "w", newline="", encoding="utf-8") as g: reader = csv.DictReader(f) writer = csv.DictWriter(g, fieldnames=["id", "x", "y"]) writer.writeheader() for row in reader: # 输入列必须是经度、纬度，不要写成纬度、经度 x, y = tr.transform(float(row["lon"]), float(row["lat"])) writer.writerow({"id": row["id"], "x": f"{x:.3f}", "y": f"{y:.3f}"})

批量示例里EPSG:4490是 CGCS2000 地理坐标，EPSG:4547是三度带投影。输出保留三位小数，是为了和常见成果表对齐；如果后续还要参与计算，建议多保留几位，落表时再取位。生产环境还要加坐标范围校验、空值捕获和拒绝文件，这些放到第 4 章展开。

3. 参数配置实战：七参数、四参数与 EPSG 校验怎么落到文件

3.1 参数来源分级：哪些七参数不能直接抄

七参数不是通用常量。不同区域、不同控制点网、不同解算方案得到的参数不一样，甚至同一区域不同年份的成果也会有差异。我一般把参数来源分成四级：A 级是测区控制点联测后自己解算的；B 级是同区域权威成果里明确给出适用范围的；C 级是软件默认值或历史图纸里抄来的；D 级是网上随手找的。跨基准换代至少要用 A 级或 B 级，C 级只能做初步试试，D 级不要进生产。

控制点数量也有讲究。三个点能解七参数，但只能勉强闭合；五个点以上、分布能包围测区，残差才有参考价值。如果控制点全挤在一条线或一个小角上，解出来的旋转和尺度会很不稳定，表面残差小，外推出去却差得远。资源包的校验脚本会输出 RMS 和最大残差，我一般要求最大残差也落在可接受范围内，而不是只看平均值。

3.2 七参数配置文件的字段设计

资源包用 YAML 管参数，字段设计比代码更重要。下面这份配置只演示结构，数值仍是示例。

datum: name: legacy_to_cgcs2000 source_ellipsoid: krass target_ellipsoid: GRS80 transform: type: helmert_7 rotation_convention: position_vector direction: source_to_target dx: -10.0 dy: 20.0 dz: 30.0 rx_arcsec: 1.2 ry_arcsec: -0.8 rz_arcsec: 0.5 scale_ppm: 2.0 projection: type: gauss_kruger lon_0: 114.0 k0: 1.0 x_0: 500000.0 y_0: 0.0 zone_width: 3 validation: check_points: check_points.csv max_rms_m: 0.05

rotation_convention和direction必须显式写，不能靠默认值。lon_0是中央经线，x_0是东伪偏移，y_0是北伪偏移，zone_width表示三度带还是六度带。validation里的max_rms_m是验收阈值，不是装饰字段；批量任务跑完后要拿检查点回算，超过阈值就停下来查参数，不要继续出成果。

3.3 四参数与高程拟合的边界

四参数适合同一基准下的小范围平面坐标纠正，典型形式是平移、旋转、尺度。它解决不了基准差异，所以不能拿四参数把旧参心系直接拽到 CGCS2000。判断边界时看两点：源和目标是不是同一基准；测区范围是不是足够小。范围一大，四参数的旋转和尺度会外推失真，常见做法是分区给参数，或者退回七参数加投影的链路。

高程要单独说。七参数输出的是椭球高，不是正常高；平面坐标准，不代表高程能用。把椭球高当正常高直接交付，差十几米并不罕见。资源包里把height_mode单独拆出来，支持不做高程、加高程异常、或接高程拟合模型。高程拟合需要高程异常控制点，不能用平面控制点凑。二维七参数只解决平面基准，不要指望它顺带把高程也修好。

3.4 用 EPSG 管线做交叉校验

自有参数跑通后，用 EPSG 管线做结构校验。下面这段把经纬度转地心、做七参数、再转回地理坐标并投影，目的是验证管线顺序，不是验证示例数值。

from pyproj import CRS, Transformer own = Transformer.from_pipeline( "+proj=pipeline " "+step +proj=unitconvert +xy_in=deg +xy_out=rad " "+step +proj=cart +ellps=krass " "+step +proj=helmert +x=-10 +y=20 +z=30 " "+rx=1.2 +ry=-0.8 +rz=0.5 +s=2.0 +convention=position_vector " "+step +proj=cart +ellps=GRS80 +inv " "+step +proj=tmerc +lat_0=0 +lon_0=114 +k=1 " "+x_0=500000 +y_0=0 +ellps=GRS80" ) epsg = Transformer.from_crs(CRS.from_epsg(4490), CRS.from_epsg(4547), always_xy=True) lon, lat = 114.123456, 30.654321 print("own:", own.transform(lon, lat)) print("epsg:", epsg.transform(lon, lat))

unitconvert先把度转弧度，cart做地理坐标到地心坐标，helmert做七参数，cart +inv转回地理坐标，tmerc做横轴墨卡托投影。convention=position_vector必须和参数文件一致。两条管线结果不一致时，先核对旋转符号和参数方向，再核对中央经线和东伪偏移。这个脚本适合放进 CI，每次改参数都跑一遍，避免手工改配置把方向改反。

4. 批量转换工程化：CSV、PostGIS 与服务接口的三条落地路径

4.1 CSV 批量转换：异常隔离与追溯字段

单点跑通后，批量任务最怕两件事：一条坏数据打断整批，以及出了问题找不到用了哪套参数。我的做法是每条记录独立捕获异常，坏数据写拒绝文件，好数据继续跑；输出里带params_id、src_crs、dst_crs，方便回溯。

import csv import logging from pyproj import Transformer logging.basicConfig(level=logging.INFO, format="%(levelname)s %(message)s") tr = Transformer.from_crs("EPSG:4490", "EPSG:4547", always_xy=True) PARAMS_ID = "legacy_to_cgcs2000_2024A" def convert(row): lon, lat = float(row["lon"]), float(row["lat"]) if not (-180 <= lon <= 180 and -90 <= lat <= 90): raise ValueError(f"坐标越界: {lon},{lat}") x, y = tr.transform(lon, lat) return {"id": row["id"], "x": f"{x:.3f}", "y": f"{y:.3f}", "params_id": PARAMS_ID, "src_crs": "4490", "dst_crs": "4547"} with open("in.csv", newline="", encoding="utf-8") as f, \ open("out.csv", "w", newline="", encoding="utf-8") as g, \ open("rejected.csv", "w", newline="", encoding="utf-8") as r: reader = csv.DictReader(f) writer = csv.DictWriter(g, fieldnames=["id", "x", "y", "params_id", "src_crs", "dst_crs"]) reject = csv.DictWriter(r, fieldnames=["id", "lon", "lat", "reason"]) writer.writeheader() reject.writeheader() for row in reader: try: writer.writerow(convert(row)) except Exception as exc: reject.writerow({"id": row.get("id"), "lon": row.get("lon"), "lat": row.get("lat"), "reason": str(exc)}) logging.warning("reject %s: %s", row.get("id"), exc)

这段代码把异常隔离成拒绝文件，params_id让结果能追溯到具体参数版本。批量转换不要用for循环里反复from_crs，转换器初始化一次就够。输出坐标保留三位小数是落表习惯，内部计算不要提前取整，否则多条链路叠加后误差会放大。

4.2 Shapefile/GeoJSON 与 PostGIS 的原地转换

文件型数据可以用 GDAL 的ogr2ogr快速跑，但跨基准时不要只写-s_srs和-t_srs两个 EPSG 代码，除非源 CRS 已经带上了towgs84或格网参数。更稳的做法是给完整 PROJ 管线。

ogr2ogr -f GeoJSON out.geojson in.shp \ -s_srs "+proj=longlat +ellps=krass +towgs84=-10,20,30,1.2,-0.8,0.5,2.0 +no_defs" \ -t_srs "EPSG:4547" \ -lco RFC7946=YES

-s_srs是源 CRS，-t_srs是目标 CRS，-lco RFC7946=YES让 GeoJSON 按经纬度输出。跨基准时源 CRS 里必须带基准转换参数，否则 GDAL 可能只做椭球替换。批量跑完先抽几个检查点对比，不要直接覆盖原文件。

数据库里我一般先加新列，不直接改原列。PostGIS 的自定义 SRID 可以挂 PROJ 字符串，再用ST_Transform原地转换。

INSERT INTO spatial_ref_sys (srid, auth_name, auth_srid, proj4text) VALUES (990001, 'LOCAL', 990001, '+proj=longlat +ellps=krass +towgs84=-10,20,30,1.2,-0.8,0.5,2.0 +no_defs') ON CONFLICT (srid) DO NOTHING; ALTER TABLE survey_points ADD COLUMN geom_cgcs2000 geometry(Point, 4490); UPDATE survey_points SET geom_cgcs2000 = ST_Transform(ST_SetSRID(geom, 990001), 4490);

自定义 SRID 用990001这类本地段，避免和正式 EPSG 代码冲突。ST_SetSRID只改元数据，不做转换；ST_Transform才走 PROJ。先加新列、双跑对比，确认无误后再建索引和切读。直接UPDATE原列会让回滚很麻烦，这是后悔药。

4.3 服务接口的幂等与精度回传

服务化时，接口不能只返回 x、y。至少要带request_id、params_id、目标 CRS 和精度状态。request_id用于幂等，同一批请求重复提交时直接返回缓存结果；params_id用于追溯参数版本；精度状态告诉调用方是“检查点通过”还是“仅几何转换”。

from fastapi import FastAPI, HTTPException from pydantic import BaseModel app = FastAPI() PARAMS_ID = "legacy_to_cgcs2000_2024A" class PointIn(BaseModel): id: str lon: float lat: float class BatchIn(BaseModel): request_id: str points: list[PointIn] @app.post("/transform") def transform(batch: BatchIn): if not batch.points: raise HTTPException(400, "points empty") out = [] for p in batch.points: x, y = tr.transform(p.lon, p.lat) out.append({"id": p.id, "x": x, "y": y, "params_id": PARAMS_ID}) return {"request_id": batch.request_id, "crs": "EPSG:4547", "accuracy": "checked", "data": out}

接口里tr是启动时初始化好的全局转换器，不要每次请求重建。accuracy字段不要写死checked，应该由检查点残差任务回填。批量大小要限流，超时要有拒绝策略，坏点要返回逐条错误而不是整批 500。服务化最容易忽略的是参数版本，结果表里没有params_id，半年后根本说不清这批坐标是用哪套参数算的。

5. 避坑与排查：转换结果差几十米时先查这五件事

5.1 现象：经纬度颠倒，平面坐标直接飞到海里

现象很直观：输入点明明在内陆，转换后的平面坐标却跑到海域或几百公里外。原因通常是轴序。CSV 列名写的是lat,lon，代码却按lon,lat读；或者 pyproj 没有设always_xy=True，某些 CRS 按纬度、经度解释。另一个变种是度分秒没有转十进制度，114°12′30″被当成114.1230。

解决分三步。第一，输入列名强制写lon、lat，不写x、y，避免歧义。第二，pyproj 的Transformer.from_crs加always_xy=True。第三，转换前做范围校验，经度不在 -180 到 180、纬度不在 -90 到 90 的直接拒绝。度分秒在入口统一转十进制度，不要留到转换函数里。这个坑我踩过不止一次，批量数据里只要有一条列顺序错，整批成果表都不可信。

5.2 现象：同一区域两批点相差几百公里

两批点来自同一区域，一批正常，另一批平面坐标差几百公里。原因多半是投影带选错。三度带和六度带的带号、中央经线不一样；跨带批量时，如果统一用一个中央经线，边缘点会严重变形。另一种情况是带号写进了 y 坐标前两位，代码却把带号当成了坐标值。

解决时先按经度算带号和中央经线，再决定是否分带。

def gk_zone(lon, width=3): if width == 3: zone = round(lon / 3) lon0 = zone * 3 elif width == 6: zone = int((lon + 1.5) / 6) + 1 lon0 = zone * 6 - 3 else: raise ValueError("width must be 3 or 6") return zone, lon0 print(gk_zone(114.123456, 3))

三度带用round(lon / 3)，六度带用int((lon + 1.5) / 6) + 1，中央经线分别是zone * 3和zone * 6 - 3。批量任务里如果测区跨带，不要硬用一个投影，按带分区输出，或者统一转地理坐标再分发。带号是否写入 y 坐标要在成果表里写清楚，不要靠猜。

5.3 现象：正反算不闭合，来回差几厘米到几米

正算再反算，回不到原点，差值几厘米到几米。原因通常是参数方向反了，或者旋转符号约定不一致。把源到目标的七参数直接取负当反算，在小角度下近似能看，但严格来说不闭合；旋转符号如果和参数文件约定相反，误差会随点位放大。

解决时把方向写进配置，反算用逆管线或逆矩阵，不要手写取负。下面用 PROJ 管线加+inv做反算，和正算共用同一套参数。

from pyproj import Transformer fwd = Transformer.from_pipeline( "+proj=pipeline " "+step +proj=unitconvert +xy_in=deg +xy_out=rad " "+step +proj=cart +ellps=krass " "+step +proj=helmert +x=-10 +y=20 +z=30 " "+rx=1.2 +ry=-0.8 +rz=0.5 +s=2.0 +convention=position_vector " "+step +proj=cart +ellps=GRS80 +inv" ) inv = Transformer.from_pipeline( "+proj=pipeline " "+step +proj=unitconvert +xy_in=deg +xy_out=rad " "+step +proj=cart +ellps=GRS80 " "+step +proj=helmert +x=-10 +y=20 +z=30 " "+rx=1.2 +ry=-0.8 +rz=0.5 +s=2.0 +convention=position_vector +inv " "+step +proj=cart +ellps=krass +inv" ) lon, lat = 114.123456, 30.654321 lon2, lat2 = inv.transform(*fwd.transform(lon, lat)) print(lon2 - lon, lat2 - lat)

正算和反算的helmert段一个正向、一个+inv，convention保持一致。闭合差应该在数值精度范围内，如果差到米级，先查旋转符号，再查尺度单位是不是漏了1e-6。

5.4 现象：平面很准，高程差十几米

平面检查点都对得上，高程却差十几米。原因通常是把七参数输出的椭球高当成了正常高，或者二维七参数根本没有处理高程异常。旧参心系和 CGCS2000 的高程基准不是简单平移，高程异常随位置变化，不能用平面参数顺带修。

解决时把高程模式单独配置：none表示只出平面，ellipsoid表示输出椭球高，geoid表示接高程异常模型，fit表示用高程控制点做拟合。交付成果表里要写清楚高程类型，不要把椭球高标成正常高。高程拟合控制点要覆盖测区，不能用平面控制点代替。平面精度和高程精度分开验收，这是基本习惯。

5.5 现象：批量跑到一半失败，重跑结果不一致

批量任务跑到某一行报错退出，重跑后输出行数不一致；或者并发跑两次，结果文件互相覆盖。原因可能是空值、

本文还有配套的精品资源，点击获取