1. 空间参考这件事为什么值得单独写一篇文章1.1 一次坐标全漂了的返工经历做GIS开发的这几年坐标系这个看似基础的概念实际坑过我不少回。印象最深的一次是接手一个第三方提交的规划数据属性表整整齐齐地块边界画得也对但一叠加到项目底图上整个地块跑到了几公里外的农田里。排查了一圈才发现数据本身的坐标没错错在坐标系信息丢掉了——对方交付时只拷贝了.dbf和.shp偏偏漏了那个不起眼的.prj文件平台端按默认坐标系解析经纬度被当成平面坐标用不漂才怪。这个经历让我意识到所有GIS项目里最难修的往往不是业务逻辑而是空间参考这个地基。地基不正后面叠加分析、投影转换、出图定位全都是错的。而要在GIS开发中真正把坐标系这件事搞清楚绕不开两个核心概念WKTWell-Known Text和EPSGEuropean Petroleum Survey Group编码。这篇文章我会把这两个东西彻底拆开讲明白再结合GDAL给出可以直接复用的实现方法。如果你是个GIS开发新手或者做了几年开发但一直对坐标系知其然不知其所以然这篇文章应该能帮你省掉不少走弯路的时间。老手也可以重点看第4章和第5章里面有GDAL 3.x版本切换过程中容易踩的几个隐蔽坑。1.2 先分清楚地理坐标系与投影坐标系在讲WKT和EPSG之前必须先把坐标系的大盘子立起来。很多开发者的困惑根源在于分不清地理坐标系和投影坐标系。地理坐标系描述的是地球椭球表面上的位置用经纬度表达单位是度。它本质上是球面坐标。你可以把地球想象成一个不太规则的土豆我们要用一个数学上可计算的椭球来近似它于是有了WGS84椭球、CGCS2000椭球这些说法。数据记录的是在这个椭球上的经度和纬度。投影坐标系则是把椭球面展开到平面上的结果。因为球面无论如何不能无损展开成平面所以数学家发明了一堆投影方法高斯-克吕格、横轴墨卡托UTM、兰伯特等角圆锥、Web墨卡托……每种方法都在保角度、保面积、保距离之间做取舍。投影坐标系你看到的坐标就是X、Y米可以直接量距离和面积。日常开发中经常遇到的情况是数据库里存的是经纬度WGS84但前端地图用的是Web墨卡托3857业务分析又需要转成当地的高斯投影。这一通操作背后的执行标准就是WKT和EPSG。1.3 WKT与EPSG在这个体系里的分工一句话总结这两者的关系EPSG是字典编号WKT是完整描述。EPSG体系给每个坐标系、基准面、椭球体分配一个全球唯一且稳定的整数编号比如EPSG:4326就是WGS84地理坐标系EPSG:3857是Web墨卡托。它解决的是标识问题就像身份证号。WKT则是一段结构化的文本把坐标系的所有必要信息完整地写出来椭球参数、基准面、本初子午线、投影方法、投影参数、单位等等。它解决的是表达问题让不同软件之间能传递同一个坐标系定义。一个shapefile旁边那个.prj文件里面存的就是一段WKT字符串。很多GIS软件、数据库和SDK都支持这两种表达方式的互相转换你给我一个EPSG码我能查出对应的WKT你给我一段WKT我也能反查它是不是某EPSG码的标准定义。GDAL在这里扮演的角色就是翻译器转换器。2. 把WKT拆开看每个字段都不是多余的2.1 WKT的演变从WKT1到WKT2WKT并不是一个静止不变的格式。最早由OGC在1999年左右发布的规范就是大家最熟悉的WKT1特征是以GEOGCS、PROJCS、DATUM这些标签开头。很长一段时间里它几乎是行业默认的坐标系文本格式。ArcGIS的.prj文件、GDAL早期版本导出的WKT基本都是这个风格。后来ISO在WKT1基础上做了大幅升级形成ISO 19162标准也就是WKT2。它修掉了WKT1里很多表达不清的地方比如支持动态坐标系、坐标历元、坐标轴顺序、参考系ensemble一个参考系的多个实例等开头的标签也变成了GEOGCRS、PROJCRS、DATUM这类更规范的名字。这里有个实际影响GDAL 3.x版本默认导出的WKT已经是WKT2系格式不同小版本可能是2015版也可能是2019版。如果你拿着GDAL 3导出的WKT字符串硬塞给只支持WKT1的老系统解析器很可能直接报错。这个问题我会在第5章专门展开先记住有这个坑就好。2.2 地理坐标系WKTGEOGCS、DATUM、SPHEROID、PRIMEM、UNIT先看一个最经典的地理坐标系EPSG:4326WGS84的WKT1长什么样GEOGCS[WGS 84, DATUM[WGS_1984, SPHEROID[WGS 84,6378137,298.257223563, AUTHORITY[EPSG,7030]], AUTHORITY[EPSG,6326]], PRIMEM[Greenwich,0, AUTHORITY[EPSG,8901]], UNIT[degree,0.0174532925199433, AUTHORITY[EPSG,9122]], AUTHORITY[EPSG,4326]]这段字符串一层套一层我来逐块解读第一层GEOGCS表示这是一个地理坐标系它的名字叫WGS 84。下面包着四块核心信息基准面、本初子午线、单位、权威码。DATUM定义基准面。基准面描述的是椭球体相对于地球实际表面怎么放置的。它内部又嵌套了SPHEROID也就是椭球体的具体参数6378137是长半轴单位是米298.257223563是反扁率。扁率简单理解就是地球的扁程度两极稍扁、赤道略鼓反扁率越大椭球越接近正球。有了这两个数椭球的形状就完全确定了。PRIMEM是本初子午线通常是Greenwich偏移角度为0。如果某个国家和地区有自己的零子午线这里就会出现非0的值。虽然多数业务用不到但WKT里必须表达清楚。UNIT是单位。地理坐标系的角度单位默认是度但WKT里它不是简单写个degree就完事第二个参数0.0174532925199433表示这个单位相对于弧度的换算系数。因为国际标准中角度的基础单位是弧度π/180约等于0.0174532925199433。这个细节很多人不看但如果你在程序中直接拿角秒、弧度混用坐标结果就会偏得离谱。AUTHORITY[EPSG,4326]则是把这一整段描述和EPSG编号对应起来相当于给这段WKT盖了个章我就是EPSG:4326。2.3 投影坐标系WKTPROJCS里多了哪些关键参数地理坐标系只回答在地球哪里投影坐标系还要回答如何把地球展开成平面。看一个投影坐标系的WKT1——EPSG:32650WGS84 / UTM zone 50NPROJCS[WGS 84 / UTM zone 50N, GEOGCS[WGS 84, DATUM[WGS_1984, SPHEROID[WGS 84,6378137,298.257223563, AUTHORITY[EPSG,7030]], AUTHORITY[EPSG,6326]], PRIMEM[Greenwich,0, AUTHORITY[EPSG,8901]], UNIT[degree,0.0174532925199433, AUTHORITY[EPSG,9122]]], PROJECTION[Transverse_Mercator], PARAMETER[latitude_of_origin,0], PARAMETER[central_meridian,117], PARAMETER[scale_factor,0.9996], PARAMETER[false_easting,500000], PARAMETER[false_northing,0], UNIT[metre,1, AUTHORITY[EPSG,9001]], AXIS[Easting,EAST], AXIS[Northing,NORTH], AUTHORITY[EPSG,32650]]和GEOGCS相比PROJCS最核心的增量是PROJECTION和一组PARAMETER。PROJECTION是投影方法这里是横轴墨卡托投影Transverse_Mercator。UTM本质上就是一种带比例因子的横轴墨卡托投影。投影参数逐个解析latitude_of_origin投影原点纬度UTM固定为0即赤道。central_meridian中央经线UTM 50带的中央经线是117°E。因为UTM把全球从西经180°开始每6度一个投影带第50带的中央经线是-180 6*50 - 3 117。scale_factor比例因子UTM固定为0.9996。它的意思是中央经线上的长度会被故意缩短到原来的0.9996倍换取整个带内长度变形尽量均匀。没有这个参数投影边缘的变形会大得多。false_easting东偏移量UTM固定为500000米。因为中央经线以西的横坐标会是负数加上500000米可以让带内所有横坐标变成正数工程上处理起来方便。false_northing北偏移量北半球UTM为0南半球UTM为10000000米避免南半球纵坐标出现负数。你会注意到投影坐标系的单位通常是米而不是度。这正是投影坐标系的直接目的——让坐标可以像普通平面几何那样参与计算。理解了这些参数你就能看懂绝大多数投影坐标系的WKT了。2.4 一个容易被忽略的事不同软件生成的WKT长得不一样同样是EPSG:4326GDAL导出的WKT、ESRI写入.prj的WKT、PostGIS中使用的WKT细节上可能各有差异。比较典型的是ESRI风格它的命名习惯是GCS_WGS_1984椭球名也可能写成D_WGS_1984而且ESRI的.prj通常不包含AUTHORITY节点。这意味着你在做系统对接时永远不要用字符串比较两个WKT是否相等。正确做法是通过库函数去判断坐标系的几何语义是否相同GDAL里就是IsSame()方法。比较字符串是新手最爱犯的错因为两个WKT序列化出来差几个空格或名称大小写不同就误判成两个坐标系了。3. EPSG其实它是一个全球统一的坐标系字典3.1 EPSG的身世从石油勘探组织到全球标准EPSG这三个字母的本义是European Petroleum Survey Group欧洲石油勘探组织。上世纪80年代一群石油勘探领域的地球科学家发现各家公司在不同矿区用的坐标系五花八门数据整合起来极度痛苦于是决定建立一个统一的坐标参考系统登记库给每个常见的坐标系、基准面、椭球体分配固定编号。后来这个维护权逐步移交给了国际油气生产商协会IOGP下属的地球科学委员会。数据库对全球免费开放直到今天还一直在更新现在版本已经到v9.x以上里面不仅有全球通用的WGS84、UTM分带也有各国的高斯投影、地方坐标系、动态参考框架。得益于这个历史渊源EPSG几乎成了GIS世界最广泛使用的坐标系字典。查询编码可以直接访问epsg.io或者到epsg.org下载官方数据。3.2 最常用的几个EPSG编码以及它们的适用场景我用一个表格把这些年项目里最高频使用的编码列出来方便你直接对照EPSG编码坐标系名称类型适用场景4326WGS 84地理坐标系GPS原始数据、全球经纬度存储3857WGS 84 / Pseudo-Mercator投影坐标系Web地图底图、前端切片4490CGCS2000地理坐标系中国国家国土数据、测绘成果4547CGCS2000 / 3-degree Gauss-Kruger zone 39投影坐标系中国区域内中央经线117°E的三度带32650WGS 84 / UTM zone 50N投影坐标系北半球UTM 50带覆盖北京附近32750WGS 84 / UTM zone 50S投影坐标系南半球UTM 50带实际使用中互联网地图服务几乎统一用3857因为它是墨卡托投影能保持角度不变形且全球覆盖范围规则测绘成果在国内常要求用CGCS2000及其高斯投影而遥感影像处理则习惯用UTM分带。这里要说一个认知EPSG编码只是约定它不代表这个坐标系就是正确的。同一个地块用4326存储、3857出图、4547计算的投影坐标数值完全不同但它们描述的是同一个地理位置。项目里最怕的是有人把4326的经纬度直接当成3857的平面坐标去算长度面积结果自然一塌糊涂。3.3 EPSG与WKT的一对多陷阱很多人以为一个EPSG编码一定对应唯一的WKT字符串其实不完全对。同一个EPSG:4326官方EPSG数据库中的定义、GDAL导出的WKT2、ESRI写入prj文件的WKT1文本内容都不一样。EPSG编码是稳定的键WKT是某个软件序列化出来的值。不同的软件、不同的WKT版本生成的字符串有差异很正常。反过来也一样一段自定义的WKT很可能找不到对应的EPSG编码。比如某个城市的独立坐标系原点是城市广场中央经线经过实地勘测确定再叠加椭球偏移参数这种坐标系在EPSG数据库里没有编号。它只能靠WKT或者PROJ字符串在软件之间传递不能靠EPSG码。另外EPSG数据库本身也会修订。WGS84这个名称下的坐标参考框架在老版本里是静态的新版本则被标记为动态参考系因为地球板块运动导致全球地心坐标系坐标值随历元漂移。这就出现了同一个EPSG:4326用2010年的数据和用2024年的数据理论上在地球上的位置有细微差异的情况。绝大多数普通GIS应用感知不到这个差异但在高精度测绘、卫星轨道计算中必须明确使用的框架和历元。4. GDAL实战从看懂到会改把坐标系握在自己手里4.1 环境准备Python、Windows与Java侧最省事的接法GDAL的安装一直是劝退新手的门槛我按开发场景分别给结论。Python我日常最推荐conda install -c conda-forge gdalconda会自动把C库、PROJ依赖、Python绑定一起搞定。如果你用pippip install gdal也可以但切记GDAL的Python包版本要和系统里的GDAL库版本一致否则运行时会报找不到动态库的错误。Windows Visual Studio 2022优先用vcpkg安装或者直接下载预编译的GDAL安装包。如果项目要求必须源码编译常规流程是解压源码后编辑nmake.opt确认WIN64YES设置好PROJ依赖路径然后依次执行nmake /f makefile.vc release和nmake /f makefile.vc devinstall。编译完成后把bin目录加入PATH设置GDAL_DATA指向gdal-data目录。大部分编译失败都是debug/release库混用、x64/x86平台不匹配导致的这两点先自查。JavaGDAL官方没有把Java绑定发布到Maven Central通常是在官网下载预编译的gdal jar包和对应平台的native库dll/so然后System.loadLibrary(gdal)并把dll目录加入java.library.path。有个细节经常被忽略GDAL运行时还要找proj.db所以光加dll不够还需要设置PROJ_LIB环境变量指向proj数据目录。三种环境按项目复杂度选一种就行。核心API逻辑是一致的。4.2 用 gdalsrsinfo 把EPSG和WKT的关系摸清楚命令行工具gdalsrsinfo是排查坐标系问题最顺手的工具没有之一。比如我想看EPSG:4326的WKT1、WKT2和PROJ字符串gdalsrsinfo -o wkt1 EPSG:4326 gdalsrsinfo -o wkt2 EPSG:4326 gdalsrsinfo -o proj4 EPSG:4326-o后面可以跟wkt1、wkt2、proj4、xml等格式。对于一段陌生的WKT你也可以把WKT当作输入传进去让GDAL帮你反解出对应的EPSG码gdalsrsinfo -o epsg GEOGCS[\WGS 84\,DATUM[\WGS_1984\,...]]这种先看定义、再查编码、最后转换的操作习惯能帮你避免很多低级错误。4.3 核心API读取、定义、导出坐标系在Python GDAL里坐标系相关操作集中在osr.SpatialReference对象上。先看最基本的读写流程from osgeo import osr, gdal # 1. 从栅格文件读取坐标系tif的投影信息存在文件头里 ds gdal.Open(raster.tif) wkt ds.GetProjection() print(文件里存的WKT, wkt) # 2. 把WKT字符串变成SpatialReference对象 srs osr.SpatialReference() srs.ImportFromWkt(wkt) # 3. 直接用EPSG码生成空间参考 srs2 osr.SpatialReference() srs2.SetFromUserInput(EPSG:4326) # 也支持WKT字符串、proj4字符串、prj文件路径 # 4. 导出成各种文本 print(srs2.ExportToWkt()) # GDAL 3.x默认输出WKT2系 print(srs2.ExportToPrettyWkt()) # 多行格式化方便人看 print(srs2.ExportToProj4()) # PROJ库的短字符串在读取shapefile时坐标系则存在.prj文件里通过layer.GetSpatialRef()拿到的就是SpatialReference对象。这个对象可以和栅格、矢量数据互相赋值。比如给一個没有定义坐标系的数据补坐标系就是先构造SpatialReference再写入数据集的投影信息。这里有个导出格式的重要提醒GDAL 3.x默认导出的WKT是WKT2格式以GEOGCRS、PROJCRS开头不是老工具的GEOGCS、PROJCS。如果你需要和老系统对接WKT1可以用命令行gdalsrsinfo -o wkt1或者在代码中设置环境变量OSR_WKT_FORMATWKT1_GDAL再导出。否则你写进.prj的WKT某些古老系统可能直接识别不了。4.4 投影转换坐标点和整体数据一起转坐标系定义清楚了下一步就是转换。GDAL里osr.CoordinateTransformation负责单个坐标点的投影变换from osgeo import osr src osr.SpatialReference() src.ImportFromEPSG(4326) dst osr.SpatialReference() dst.ImportFromEPSG(32650) # 关键步骤GDAL 3默认使用EPSG定义的轴序(纬度/经度) # 传统GIS习惯是经纬度(经度/纬度)这里必须显式指定 src.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER) dst.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER) ct osr.CoordinateTransformation(src, dst) # 以北京附近一个点为例输入是经度116.40、纬度39.90 lon, lat 116.40, 39.90 x, y, z ct.TransformPoint(lon, lat) print(fUTM 50N坐标X{x:.2f} 米, Y{y:.2f} 米)运行后X大约在45万米量级因为在50带的中央经线117°E以西约1.5度减去500000偏移后横坐标小于50万所以X大约是44万多Y大约是441万多米。这个坐标拿去做距离量算得到的单位就是米直接可用。如果你要转换整个文件不需要自己循环每一个要素用命令行工具最快。矢量转坐标ogr2ogr -t_srs EPSG:32650 output.shp input.shp栅格转坐标顺便重采样gdalwarp -t_srs EPSG:32650 -r bilinear input.tif output.tifogr2ogr和gdalwarp的内部逻辑和我们上面写的代码完全一致只是封装成了文件级操作。4.5 RPC正射校正场景为什么输出投影常选UTM再结合一个现实中常见的场景精准测绘或遥感影像处理中经常用到RPC正射校正。卫星影像自带RPC有理多项式参数gdalwarp可以做地形校正把倾斜影像转成正射影像。典型命令长这样gdalwarp -rpc -t_srs EPSG:32650 -to RPC_DEMdem.tif image.tif ortho.tif-rpc开启RPC纠正-to RPC_DEMdem.tif传入数字高程模型-t_srs EPSG:32650把输出结果投影到UTM 50N。这里有个值得思考的问题为什么输出用UTM而不是直接用WGS84经纬度核心原因是UTM在局部区域的长度变形极小而且均匀影像纠正后可以直接按米量测、做地块核算后续套合地形图也方便。同时UTM的投影参数是全球标准化的不用自己临时定义中央经线避免了这个自定义投影参数写错了的隐患。这也是遥感处理行业多年形成的习惯——只要影像范围在一个6度带内UTM就是最稳妥的输出投影。5. 实战中我踩过的坑坐标顺序、WKT版本和自定义参考系5.1 EPSG:4326到底是经纬度还是纬经度这是GDAL 3.x时代最隐蔽、最容易翻车的一个坑。EPSG官方数据库中4326定义的轴顺序是纬度在前、经度在后参考的是ISO 19111的坐标轴顺序规范。但在传统GIS习惯里几乎所有数据都用经度在前、纬度在后比如GeoJSON的[longitude, latitude]shp文件读写也默认这个顺序。GDAL 3.0之后很多API的行为都开始遵循数据源自己声明的轴序当你手动调用CoordinateTransformation.TransformPoint()时如果源SRS是4326且没做任何设置GDAL可能把第一个参数当作纬度、第二个当作经度。你明明传入的是(116.40, 39.90)它却当成(纬度116.40, 经度39.90)来处理转换结果直接甩到北极圈附近。这就是我文章开头说的坐标全漂了的另一种常见原因。规避方法很简单就是我在4.4代码里写的那两行src.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER) dst.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)把轴序策略切换回传统GIS顺序经度/纬度。只要是手动传坐标、手动接收坐标的代码都要加上这个设置。文件级转换一般不用管因为驱动程序已经处理了轴序问题但裸坐标转换一定要自己控制。5.2 WKT1与WKT2老系统对接的兼容性前前后后遇到过几次GDAL导出坐标写回.prjArcGIS打不开的案例。原因就是GDAL 3导出的WKT2格式老版本ArcGIS不认识。如果你要对接老平台最简单的办法是让GDAL输出WKT1。命令行靠gdalsrsinfo -o wkt1代码里设置环境变量OSR_WKT_FORMATWKT1_GDAL后再ExportToWkt()。设置之后导出的文本就是老朋友GEOGCS开头的WKT1格式了兼容性立刻上来了。另外一个相关坑是ESRI的WKT风格。ESRI在很多软件里用的并不是标准ISO WKT而是自己的一套写法标签名、命名规则都有差异比如把坐标系写成GCS_WGS_1984而不是WGS 84。在做跨软件互换时不要想当然认为EPSG:4326在任何软件里导出的WKT都完全一致。判断两个坐标系是否等价应该使用库函数GDAL里就是SpatialReference.IsSame()。5.3 自定义坐标系没有EPSG编号怎么办地方规划、矿区测量里经常遇到一种坐标系椭球是参照国家椭球的但中央经线是本地某一条经线还叠加了独特的平移参数。这类坐标系在EPSG数据库里根本没有编号唯一的通行凭证就是WKT。处理办法是把坐标系定义写成标准WKT存入项目数据库或配置文件在需要的时候通过SetFromUserInput(wkt)加载。GDAL还接受PROJ字符串例如srs osr.SpatialReference() srs.SetFromUserInput(projtmerc lat_00 lon_0114 k1 x_0500000 y_00 ellpsGRS80 unitsm no_defs)但我个人建议跨系统传递时优先用WKT而不是PROJ字符串。PROJ字符串是PROJ库的串行化格式表达能力有限而WKT可以完整表达基准面、椭球、投影参数、轴序、动态框架等全部信息是行业通用的完整描述。另外提醒一件事如果自定义坐标系项目里很多人共同使用最好在协作规范里固定WKT的版本和生成工具比如统一用gdalsrsinfo -o wkt1生成并保存避免不同人导出不同WKT历史数据对不上。5.4 UTM分带计算与输出前的检查清单最后给一个能直接用的技能如何判断一个经纬度落在UTM哪个带。UTM全球从西经180°开始每6度一个带带号计算公式是zone floor((lon 180) / 6) 1以北京为例经度116.4°zone floor((116.4 180) / 6) 1 floor(49.4) 1 50所以北京落在UTM 50带。北半球用EPSG:32650南半球用EPSG:32750。如果是其他区域再用同样的公式和南北半球范围选编号。不过这里有个冷知识挪威和斯瓦尔巴群岛附近的UTM分带和标准公式不一致因为这两个区域做了带合并和偏移处理相关的EPSG编码并不等于纯公式计算结果。如果你处理的是这两个区域的遥感数据要留心。总结一下做坐标系输出前建议按这个清单自查输入数据的坐标系信息是否完整.prj/srs是否丢失输出的EPSG编码对应的带号、中央经线是否正确手动坐标转换有没有设置轴序策略导出WKT的版本是否兼容下游系统投影单位确认是米而不是度转换后的坐标数值是否在该区域合理范围比如北京UTM X应在40万到50万之间Y应在440万左右。我在实际项目里已经形成了一个习惯所有入库数据先过一遍GDAL的SpatialReference检查给每份数据生成一张WKT快照存档防止哪天有人把.prj删了导致排查半天。坐标系这种地基问题宁可前期多花五分钟也别等数据铺开了再返工。
