ARTICLE DETAIL

资讯详情

深耕网站建设与运营推广的一线实战洞察。

经纬度坐标转CGCS2000:基准转换与七参数完整指南

经纬度坐标转CGCS2000:基准转换与七参数完整指南 简介面向ArcGIS用户的WGS84经纬度坐标转CGCS2000投影坐标系操作说明文档适合测绘、GIS数据处理及国土相关技术人员参考。文档完整梳理了从奥维地图导出shp文件到ArcGIS中进行投影与变换再到通过ITRF2000过渡并最终转换为CGCS2000坐标系的全流程涵盖WGS84、ITRF2000、CGCS2000三种坐标系的椭球参数与转换逻辑以及带号选择、3度带/6度带中央经线计算、坐标字段添加与十进制单位切换等关键细节。包体仅含1个doc文件大小3.06MB为图文步骤说明可直接对照ArcToolbox操作。该文档已有646人学习下载适合需要处理不同坐标系数据转换、规避常见误差的GIS使用者是一份实用且可复用的操作笔记。对于将奥维兴趣点、测量成果等外部数据转入国家2000坐标系的应用场景按此文档操作可减少转换错误步骤按菜单层级逐步展开也便于零基础读者对照完成。1. 经纬度坐标系转CGCS2000差的不是度数是基准许多人以为经纬度是全球通用的拿到GPS坐标直接当成CGCS2000用。实际上普通GPS接收机默认输出WGS84经纬度而CGCS2000是一个独立的地心坐标参考框架二者椭球参数、参考历元和框架实现都有细微差异造成的平面偏移在多数地区有0.51.5米。这在小比例尺地图上不明显但施工放样、地籍测量、CAD与GIS拼图时足以让线位错位。这篇文章只解决一件事如何把经纬度坐标系的数据正确转成CGCS2000坐标系包括原理、七参数计算、复用工具和最终验证。适合测绘内业、GIS开发者、CAD转GIS的数据处理人员。2. 经纬度转CGCS2000前先把椭球、基准和投影分开算2.1 经纬度本质是大地坐标必须先锁定参考椭球经纬度(B,L,H)本身只是表达方式不是基准。同一坐标数值放在不同椭球上对应的空间位置不同。CGCS2000使用CGCS2000椭球长半轴a6378137m扁率f1/298.257222101而WGS84椭球扁率为1/298.257223563长半轴相同扁率差很小但存在。转换第一步不是套公式而是确认源和目标各处于哪个椭球/框架。一般从在线地图或手持GPS获取的经纬度可视为WGS84国土、规划数据中的经纬度很可能已经是CGCS2000。要把经纬度转成空间直角坐标常见做法是先用椭球参数计算卯酉圈曲率半径N再按下面关系式分解X (N H) * cos(B) * cos(L) Y (N H) * cos(B) * sin(L) Z (N * (1 - e^2) H) * sin(B)其中e^2 2f - f^2H是大地高不是海拔。这段公式可以直接用Python实现def geodetic_to_cartesian(lat_deg, lon_deg, h_m, a6378137.0, rf298.257222101): from math import radians, sin, cos, sqrt f 1.0 / rf e2 2 * f - f * f b radians(lat_deg) l radians(lon_deg) n a / sqrt(1 - e2 * sin(b) ** 2) x (n h_m) * cos(b) * cos(l) y (n h_m) * cos(b) * sin(l) z (n * (1 - e2) h_m) * sin(b) return x, y, z这里把纬度、经度、大地高传进去默认使用CGCS2000椭球。函数返回的是以地心为原点的XYZ直角坐标单位米。rf是扁率倒数e2是第一偏心率平方这两个是椭球转换里最容易抄错的地方。如果你手里的是海拔高程还要在转换前加大地水准面差距否则转出的XYZ会整体偏高或偏低影响后续七参数。2.2 CGCS2000和WGS84只差一个基准转换椭球参数只决定了数学形状不决定地球框架。CGCS2000采用ITRF97参考框架参考历元2000.0WGS84在G1762版本后与ITRF对齐到厘米级。实际工作中很多地方直接忽略转换导致局部有规律性偏移。如果项目精度要求优于1米就必须做基准转换。国内常用的是布尔沙七参数模型。七参数包括平移参数DX、DY、DZ米旋转参数RX、RY、RZ角秒或弧度以及尺度因子mppm。转换关系可以写成矩阵形式小角度情况下常见做法是展开为X2 DX (1 m) * X RZ * Y - RY * Z Y2 DY (1 m) * Y - RZ * X RX * Z Z2 DZ (1 m) * Z RY * X - RX * Y也就是说把源坐标X、Y、Z做一次线性变换得到目标框架下的坐标。注意旋转参数单位要换算成弧度而且公式的旋转符号约定在不同地区可能相反需要小范围试算确认。七参数的来源要非常谨慎。示例参数如下仅为写法示意不能用于生产DX 1.20 m DY -0.80 m DZ -1.50 m RX -0.02 RY -0.01 RZ 0.03 m 0.5 ppm提示生产项目的七参数必须向当地自然资源或测绘主管部门索取或使用GNSS连续运行参考站提供的本地参数不要从网上随便抄一组。参数对不上时结果会偏离数百米。2.3 七参数怎么套进经纬度里手工搞定的套路是源经纬度转XYZ - 应用七参数得到目标XYZ - 目标XYZ反算目标经纬度。如果目标要投影平面坐标再对目标经纬度做高斯-克吕格投影。整个过程只有三步看起来不复杂但每一步都依赖椭球参数和投影带号其中一项选错后面全错。反算XYZ到经纬度需要迭代因为大地纬度出现在N的计算式里。常见做法是用初始纬度迭代几次收敛到毫米级就停。2.4 转到平面时要区分3度带和6度带当目标CGCS2000坐标是平面坐标如CAD图纸里的x3456789.12, y20612345.67实际是高斯-克吕格投影平面坐标。Y坐标前两位数“20”是6度带带号实际横坐标是612345.67东偏移500km后。3度带或6度带的中央经线计算公式6度带L0 6 * n - 3例如带号20中央经线L0 117°E3度带L0 3 * n例如带号39中央经线L0 117°E下表是东经114°附近几个带号的中央经线对照带号类型带号中央经线适用经度范围6度带19111°E108°E ~ 114°E6度带20117°E114°E ~ 120°E3度带38114°E112.5°E ~ 115.5°E3度带39117°E115.5°E ~ 118.5°E选择带号时按目标点经度算不要照抄图纸里任意坐标的带号尤其跨带图幅要额外说明。很多CAD转GIS的“6位坐标”问题本质上就是Y坐标省略了带号后面第4章会单独讲。3. 用 pyproj 和手动七参数把经纬度转成 CGCS20003.1 首选 pyproj两行代码先解决 WGS84 到 CGCS2000 的经纬度pyproj是GDAL生态里最顺手的坐标转换库。如果源数据是WGS84经纬度目标只需要CGCS2000经纬度直接定义Transformer即可。CGCS2000的地理坐标EPSG代码是4490WGS84地理坐标是4326。from pyproj import Transformer t Transformer.from_crs(4326, 4490, always_xyTrue) # 输入经度、纬度输出CGCS2000经纬度 lon, lat t.transform(117.123456, 39.654321) print(lon, lat)always_xyTrue表示输入输出都是经度在前、纬度在后避免和传统纬度在前习惯混在一起。转换时pyproj会按两个坐标系对应的椭球和基准做换算默认没有使用七参数时它做的是“忽略基准面差异”的转换极端情况下会产生米级偏差。如果你的数据源本身是CGCS2000下的经纬度直接跳过这步即可。3.2 立即转平面坐标高斯投影的 PROJ 字符串很多项目最终要的是CGCS2000平面坐标。比如把WGS84经纬度转成6度带20带的高斯平面坐标可以自定义一个PROJ坐标系使用projtmerc中央经线设置为117°E。这样既能同时完成基准转换和投影from pyproj import CRS, Transformer source CRS.from_epsg(4326) target CRS.from_proj(projtmerc lat_00 lon_0117 k1 x_0500000 y_00 ellpsGRS80 unitsm no_defs) t Transformer.from_crs(source, target, always_xyTrue) x, y t.transform(117.123456, 39.654321) print(x, y)这里的x是高斯平面北坐标y是横坐标注意x_0500000表示横坐标加500公里偏移输出Y在500000左右。如果你需要带带号在Y前拼接带号20例如结果y206123456.78。CGCS2000椭球比GRS80扁率略大但在毫米级精度要求不高的场景下用ellpsGRS80通常不影响工程应用若需要严格匹配把ellpsGRS80替换为a6378137 rf298.257222101。3.3 手动七参数转换不依赖 PROJ 也能算在某些离线环境里没有pyproj或者甲方只给了一组七参数我一般会直接手写转换函数。完整过程是源BLH转XYZXYZ过七参数再反算目标BLH。下面的函数实现了三步import math def blh_to_cgcs2000(lat_deg, lon_deg, h_m, dx, dy, dz, rx, ry, rz, m_ppm, rf_src, rf_dst): # 1. 源椭球BLH - XYZ x, y, z geodetic_to_cartesian(lat_deg, lon_deg, h_m, rfrf_src) # 2. 七参数转换 rx math.radians(rx / 3600.0) # 角秒转弧度 ry math.radians(ry / 3600.0) rz math.radians(rz / 3600.0) s 1 m_ppm * 1e-6 x2 dx s * (x rz * y - ry * z) y2 dy s * (-rz * x y rx * z) z2 dz s * (ry * x - rx * y z) # 3. XYZ - BLH迭代求纬度 lon math.atan2(y2, x2) p math.sqrt(x2 * x2 y2 * y2) lat math.atan2(z2, p * (1 - 2.0 / rf_dst)) for _ in range(5): n 6378137.0 / math.sqrt(1 - (1 - 1.0 / rf_dst) ** 2 * math.sin(lat) ** 2) h p / math.cos(lat) - n lat math.atan2(z2, p * (1 - (1 - 1.0 / rf_dst) * n / (n h))) lon_deg math.degrees(lon) lat_deg math.degrees(lat) return lat_deg, lon_deg, h七参数中的旋转单位在这里按角秒处理角秒转弧度要除以3600再乘π/180rx参数传入时是角秒值。第三步反算纬度用了迭代近似5次循环已足够收敛到毫米级。这个函数把源椭球扁率rf_src和目标椭球扁率rf_dst都显式传进去避免把大地高当海拔。3.4 转换结果对拍表到底差多少米下面用一组示意数据演示转换前后差异。输入WGS84经纬度CGCS2000近似等于其目标值平面坐标只是示例说明量级。输入WGS84经度输入WGS84纬度输出CGCS2000经度输出CGCS2000纬度6度带20带平面坐标Y117.12345639.654321117.12345239.65431820612344.12表格里的差值是示意性的实际差异由七参数决定。判断结果是否合理关键是看Y坐标是否落在带号对应的中央经线附近以及同一坐标在不同方法下回算残差是否小于0.01米。4. 不写代码的路径ogr2ogr、GeoHey 在线转换与 CAD 的 6 位坐标4.1 用 ogr2ogr 批量转矢量文件如果你手头是Shapefile或GeoPackage不想写Python用GDAL自带的ogr2ogr最直接。把矢量数据从WGS84经纬度转成CGCS2000地理坐标命令ogr2ogr -s_srs EPSG:4326 -t_srs EPSG:4490 output.shp input.shp-s_srs指定源坐标系-t_srs指定目标坐标系。要转到高斯平面坐标需要写更完整的参数ogr2ogr -s_srs EPSG:4326 -t_srs projtmerc lat_00 lon_0117 k1 x_0500000 y_00 ellpsGRS80 unitsm no_defs output.shp input.shp批量转换时建议先转一份小数据验证字段和坐标再用循环处理整个目录。ogr2ogr默认会重写OGR字段不会改变属性表结构。4.2 GeoHey 在线坐标转换做单点抽查需要快速核对单个经纬度时在线坐标转换是很多人的习惯做法。GeoHey在线坐标转换工具里通常要先选源坐标系和目标坐标系。这里有个容易忽略的细节工具里如果没有“CGCS2000地理位置”选项可以选“CGCS2000 / 3-degree Gauss-Kruger CM 117E”这类投影坐标系结果会直接是平面坐标。单点抽查只用来排查方向性错误不建议大批量生产。在线工具大多默认不公开具体七参数所以要确认结果是否满足你项目的平面精度要求。转换方式适合场景精度控制上手难度pyproj批量处理/脚本集成可显式设置七参数中ogr2ogr矢量文件批量转依赖-t_srs定义低GeoHey在线转换单点抽查工具内部实现很低手动七参数离线/学习原理参数可控高4.3 CAD到GIS6位坐标其实是投影坐标没带带号CAD图纸经常出现“6位坐标转换”的疑问CAD里标出的点比如X3456789.12Y612345.67导入GIS后跑到别的位置。常见原因是这个Y是高斯平面坐标去掉了带号。比如6度带20带的实际横坐标应为20612345.67Y612345.67只是EASTING的十位到个位部分。对应处理方法是先补带号再转经纬度或者补带号后直接定义成CGCS2000平面坐标。我一般先看图纸说明如果坐标范围只有6位很可能是省略了中央经线前的带号如果Y值大于500000则横坐标正常。补带号后在ArcGIS中定义坐标时选择CGCS2000 6度带对应的投影坐标系例如20带。如果图纸是3度带就把带号换成39之类的3度带带号。这样处理后CAD到GIS的6位坐标问题就解决了一半剩下需要确认长度单位是米还是毫米。5. 三个必查参数和一组回算验证转换结果5.1 必查参数源椭球、目标框架和转换参数转换结果异常时先检查三个参数源坐标的椭球和基准、目标坐标系有没有带投影带号、七参数是否匹配区域。按顺序排查多数问题会浮出来。源坐标如果是RTK导出的经纬度通常已落入CGCS2000如果是手机GPS基本是WGS84。目标坐标如果要求带带号的Y值就不要选不带带号的定义否则绘图软件会把它当普通坐标。5.2 用已知控制点做残差回算最可靠的验证方法是回算已知点。选一个当地已知的CGCS2000控制点平面坐标或经纬度把它和待转换点用同一套流程处理。具体做法是先把你的参数和代码转出一个值再用反函数转回原坐标系比较初始值和回算值的差。也可以取两个已知点一个做参数拟合一个做独立验证。编写一个小脚本计算残差for pt in known_points: lat2, lon2 transform(pt.lat, pt.lon) dx_m (lon2 - pt.lon_cgcs) * 111000 * math.cos(math.radians(pt.lat)) dy_m (lat2 - pt.lat_cgcs) * 110946 print(pt.name, dx_m, dy_m)这里用近似公式把经度纬度差换算成平面米数适合快速检查。正常残差应在厘米到分米级如果出现米级以上偏差回到5.1重新核对参数。验证完成后保留脚本和参数表作为交付文档的一部分。本文还有配套的精品资源点击获取
返回列表