WKT 与 EPSG 如何表达空间参考坐标系？附 GDAL 实现

2 空间参考坐标系的表达通过第2章介绍的地理空间参考系统的知识我们可以知道一个空间参考坐标系的内容是非常丰富的需要描述参考椭球参数、基准面参数等如果是投影坐标系还需要描述投影参数并且每个投影算法的参数都不同。而在实际的GIS应用中对空间参考坐标系的表达是复用程度和通用程度都很高的功能那么有没有一种可能采用一种约定俗成的方式去描述各种空间参考坐标系呢。答案是肯定的那就是WKT字符串和EPSG编码。3.2.1 WKT字符串WKT(Well-known Text)是开放地理空间信息联盟OGCOpen Geospatial Consortium制定的一种文本标记语言用于表示矢量对象和地理空间参考系统。提示开放地理空间信息联盟OGC是一个非盈利的国际标准化组织致力于提供地理信息行业软件和数据及服务的标准化工作。不同的GIS厂商可以按照这些标准定义空间数据存储的编码和开放服务的接口从而可以保证空间数据的互操作。OGC的成员包括了像ESRI、Google、Oracle等业界强势企业可以认为是国际上一种比较“权威”的标准化机构。OGC目前提供的标准很多我们在后面也会接触到一些。根据最新的OGC标准使用WKT来表达的空间参考系统的标准叫做“WKT CRS”也就是坐标参考系统的WKT表示。在这个标准规范中详细定义了如何通过WKT来表示各种空间参考系统。例如我们最常用的WGS84坐标系其WKT字符串为GEOGCS[GCS_WGS_1984, DATUM[D_WGS84, SPHEROID[WGS84,6378137,298.257223563]], PRIMEM[Greenwich,0], UNIT[Degree,0.017453292519943295]]结合第2章的内容我们很容易识别其描述的具体含义GEOGCS表明其是一个地理坐标系DATUM表示大地基准面SPHEROID则定义了一个参考椭球体每一个参数后面的方括号内部则是其具体的参数。而一个投影坐标系则形如PROJCS[Xian_1980_3_degree_Gauss_Kruger_zone_39, GEOGCS[GCS_XIAN 1980, DATUM[D_XIAN 1980, SPHEROID[Xian_1980,6378140,298.2569978029111]], PRIMEM[Greenwich,0], UNIT[Degree,0.017453292519943295]], PROJECTION[Transverse_Mercator], PARAMETER[scale_factor,1], PARAMETER[central_meridian,117], PARAMETER[latitude_of_origin,0], PARAMETER[false_easting,39500000], PARAMETER[false_northing,0], UNIT[Meter,1]]PROJCS标识其是一个投影坐标系定义投影坐标系的首要内容当然是需要先定义一个地理坐标系这里的GEOGCS定义了一个XIAN80坐标系PROJECTION则标识其采用的横轴墨卡托投影算法PARAMETER则是其投影参数。由于参数false_easting在500000前面加了带号39中央经线central_meridian为117所以这是一个高斯-克吕格投影坐标系使用的地理坐标系为XIAN80坐标系采用3度带进行投影带号为39。使用WKT字符串的好处是易于识别可以让用户和开发人员清楚的知道我们使用的到底是哪一个空间参考系统。但是让我们创建一个这么长的字符串就比较麻烦了这需要非常了解WKT对坐标参考的具体规范和标准才行。简单总结来说就是易读不易写。3.2.2 EPSG编码与WKT字符串易读不易写的特性相反使用EPSG编码来表达空间坐标参考系统的特点则是易写不易读。使用EPSG编码的原理在于世界上空间坐标参考系统固然很多但是通用的经常使用的个数确是有限的我们可以将这些常用的空间坐标参考系统都分配一个编号这个编号就是EPSG编码。提示EPSG的是英文European Petroleum Survey Group的缩写中文名称为欧洲石油调查组织。该组织维护了一个包含所有空间坐标参考系统相关的实体及其EPSG编码的在线数据库https://epsg.io/ 。这个数据集在GIS业界内得到了广泛的使用。EPSG编码非常简单一般来说就是4位数字及其以上编码。在下表3.3中展示了一些经常使用的空间参考坐标系的EPSG编码。可以看到一些常用的分带投影坐标系都被赋予了EPSG编码表3.3 常用空间参考坐标系及其EPSG编码名称EPSG编码Beijing 19544214Xian 19804610WGS844326CGCS20004490CGCS2000 / Gauss-Kruger zone 154493WGS 84 / UTM zone 44N32644Web Mercator3857比较有意思的是Web墨卡托投影坐标系的EPSG编码。我们在第2.5.3节中介绍过Web墨卡托投影坐标系在进行投影计算时将地球的参考椭球体近似成正球体处理是近似等角而不是严格等角。所以Web墨卡托投影是一种伪墨卡托投影由Google最先发明造成Web墨卡托投影长期得不到GIS领域业内承认。不过随着Web墨卡托的广泛使用EPSG最终还是给予了其编码3857。在这中间的过程中不同的GIS厂商各自使用不同的编码来表达Web墨卡托投影坐标系900913、3587、54004、41001、102113、102100、3785。现在这些编码都已废弃统一使用3857。3.3 空间参考坐标系的代码实现虽然已经有WKT字符串和EPSG编码来表达空间坐标参考系统但是这并不意味着我们需要从头开发一个程序来实现空间参考坐标系的相互转换。我们可以借助第三方开源库来实现这个功能也就是本节要介绍的PROJ和GDAL。3.3.1 开源工具PROJ和GDALPROJ是一个通用的坐标转换开源软件最主要的功能就是支持坐标从一个空间坐标参考系转换到另一个坐标参考系。PROJ还公开了一个应用程序编程接口API以便于让用户在自己的程序中使用PROJ的功能而无需自己实现类似的功能。PROJ是一个C库优点是强大的坐标转换算法。不过对空间参考坐标系的表达不是很直观它使用了自己定义的字符串来表达空间参考坐标系这个字符串在别的地方用的不是很多。所以很多时候是将PROJ集成进入GDAL来使用。GDALGeospatial Data Abstraction Library是一个开源的GIS数据抽象库包括栅格数据、矢量数据、空间坐标参考等并且还提供了对这些GIS数据的读取、写出、处理和转换。基本上来说这个库是做GIS开发的必备库集成了很多GIS的基础功能功能强大的同时最难能可贵的是这些都是开源免费的。GDAL可以通过WKT字符串和EPSG编码实现对空间坐标系的表达在集成PROJ之后GDAL就能实现空间坐标转换的功能。3.3.2 具体实现在Visual Studio 2019中新建一个控制台程序我们按照第3.1节中论述的方法配置好开发环境引入GDAL第三方库写出我们第一个示例程序主要代码如下例3.1所示例3.1 空间参考坐标系的表达#include ogr_spatialref.h #include iostream using namespace std; void OutputGcs() { OGRSpatialReference spatialReference; spatialReference.importFromEPSG(4326); // WGS84 // spatialReference.importFromEPSG(4214);//BeiJing54 // spatialReference.importFromEPSG(4610);//XIAN80 // spatialReference.importFromEPSG(4490);//CGCS2000 char* pszWKT nullptr; spatialReference.exportToPrettyWkt(pszWKT); cout pszWKT endl; CPLFree(pszWKT); pszWKT nullptr; } void OutputPcs() { OGRSpatialReference spatialReference; spatialReference.importFromEPSG(4490); // CGCS2000 spatialReference.SetTM(0, 114, 1.0, 38500000, 0); char* pszWKT nullptr; spatialReference.exportToPrettyWkt(pszWKT); cout pszWKT endl; CPLFree(pszWKT); pszWKT nullptr; } int main() { //设置Proj数据 string projDataPath getenv(GISBasic); projDataPath /share/proj; CPLSetConfigOption(PROJ_LIB, projDataPath.c_str()); cout 地理坐标系WGS84坐标系 endl; OutputGcs(); cout 投影坐标系高斯克吕格投影坐标系 endl; OutputPcs(); }这段代码输出了两种有代表性的空间参考坐标系第一种是WGS84坐标系是我们最常用的地理坐标系很显然去传入各种地理坐标系的描述是非常麻烦了最简单的就是直接传入具体的EPSG代码。第二种是国内专业地图最常用的高斯克吕格投影坐标系不过很可惜GDAL没有封装专门的接口。但是正如第2章中所述高斯克吕格投影坐标系其实就是横轴墨卡托投影因此我们可以直接用横轴墨卡托投影的接口。如代码所述第二种投影坐标系描述的其实是高斯克吕格投影坐标系3度带第38度带中央经线为114。代码运行结果如下所示输出了两种空间参考坐标系的WKT字符串地理坐标系WGS84坐标系 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]], AXIS[Latitude,NORTH], AXIS[Longitude,EAST], AUTHORITY[EPSG,4326]] 投影坐标系高斯克吕格投影坐标系 PROJCS[unnamed, GEOGCS[China Geodetic Coordinate System 2000, DATUM[China_2000, SPHEROID[CGCS2000,6378137,298.257222101, AUTHORITY[EPSG,1024]], AUTHORITY[EPSG,1043]], PRIMEM[Greenwich,0, AUTHORITY[EPSG,8901]], UNIT[degree,0.0174532925199433, AUTHORITY[EPSG,9122]], AUTHORITY[EPSG,4490]], PROJECTION[Transverse_Mercator], PARAMETER[latitude_of_origin,0], PARAMETER[central_meridian,114], PARAMETER[scale_factor,1], PARAMETER[false_easting,38500000], PARAMETER[false_northing,0], UNIT[metre,1, AUTHORITY[EPSG,9001]], AXIS[Easting,EAST], AXIS[Northing,NORTH]]