QGIS WKT自定义坐标系研究及其与ArcGIS的对比
我国全国版图常用的投影为双标准纬线圆锥投影,如等积Albers投影或者等角Lambert投影,所使用的地理坐标参照系则依据数据生产时间而有所不同,如早期的北京54坐标系和西安80坐标系,以及近年来所普遍采用的国家2000大地坐标系。本文以获取自中国科学院资源环境科学与数据中心的2015年省级行政区划数据为示范,该数据为SHP格式,采用的就是基于北京54坐标系的Albers投影。
这种双标准纬线、根据制图区域选定中央经线、地理坐标系、投影方式的组合,在大多数GIS软件的预定义坐标系中普遍不存在,需要自行定义坐标系,否则该图层将被标记为未识别坐标系的数据,在使用过程中经常出现计算错误或者无法与其他来源的图层叠加显示等问题。
QGIS的自定义投影功能位于菜单【设置】->【自定义投影】下,点击后弹出自定义坐标参照系对话框:
点击右上角的
QGIS自定义坐标系支持WKT和Proj两种描述格式,由于Proj过于简洁,是一种有损描述,因此QGIS官方建议使用WKT语言来描述坐标系定义。
熟知文本(Well-known text ,简称WKT) 是由开放地理空间联盟(OGC)制订的一种文本标记语言,用便于人们阅读的简洁语法表达地理空间中的几何对象。这种以文字表达几何图形的方式因其方便人们理解、具有很高的可读性、易于存储和易于共享等优势而深受欢迎,进而广泛应用到坐标系定义、坐标变换参数描述等方面。
本文所使用的术语“WKT坐标系定义”,特指以WKT标记语言描述的坐标系定义字符串,若无特殊说明,本文中的“WKT坐标系定义”一词都是此含义。
但是,由上图可以看到,QGIS的自定义坐标系界面做得比较简陋,操作体验不太友好。
自定义坐标系本身就是一个复杂的过程,需要设置的参数包括椭球体、基准面、投影方式、标准纬线、中央经线等。商用GIS软件通常将这个过程分解为易于使用的几个操作步骤向导,用户只需要点击和输入对应的参数即可。例如,ArcMap 10.2自定义坐标系过程如下:
对普通用户来说,相比手动输入代码的方式,操作向导界面能极大提高用户体验。因为无论使用WKT格式还是Proj格式,构造一个新的坐标系定义,都需要了解其中的语法结构、关键字和取值。这些对于不了解相关语言规范又想用QGIS快速作图的人来说几乎无从下手。
问题最后集中在如何获取WKT坐标系定义描述字符串,除了希望随着QGIS版本更新,开发团队提供更加友好的操作方式进行自定义坐标系之外,本文通过对比QGIS与ArcGIS WKT坐标系定义的异同,得出通过导入ArcGIS PRJ文件中的WKT坐标系定义字符串,解决当前状况下QGIS自定义坐标系操作中WKT坐标系定义字符串难以构造的问题。
将ArcGIS WKT坐标系定义导入QGIS
如果一个图层的坐标系在ArcGIS中进行了正确定义,其SHP文件同名的.PRJ文件即为符合WKT语言规范的坐标系定义文本文件。可以用文本编辑器打开并查看其内容,例如2015年行政区划数据,其.PRJ文件的内容为:
这个WKT坐标系定义是否与QGIS的WKT坐标系定义兼容?如果兼容,则可用它作为QGIS自定义投影的坐标系WKT参数。为了节约大家阅读时间,我先把结论放在这里:QGIS兼容ArcGIS的WKT坐标系定义字符串,可以直接复制.PRJ文件的内容作为QGIS新坐标系的WKT参数。
具体操作步骤如下:
点击菜单【设置】->【自定义投影】,打开自定义坐标系对话框,导入过程如下图:
返回地图主窗口,单击右下角的坐标系区域,弹出坐标系选择对话框,在用户自定义区域找到刚刚新建的坐标系,点击【OK】,QGIS就可以识别该坐标系了。
接下来,我们看看为什么QGIS能够兼容ArcGIS的WKT坐标系定义字符串。要回答这个问题,先从二者的不同之处说起。
QGIS与ArcGIS WKT坐标系定义的差异
从QGIS的坐标系选择器打开一个投影坐标系,观察两个软件所生成的WKT坐标系定义的结构和取值(此处选择示范数据所带PRJ导入后生成的自定义坐标系进行比较):
对比两个字符串,可以看到,虽然都是WKT坐标系定义,二者却有很多不同之处,表现在:
1、关键字不同
在ArcGIS的PRJ文件中,WKT坐标系定义以“PROJCS”开头,表示该坐标系为投影坐标系。与之相对应的QGIS投影坐标系使用关键字“PROJCRS”,后者多了一个字母“R”。每个投影坐标系均以地理坐标系为基础,在ArcGIS中,地理坐标系使用关键字“GEOGCS”,QGIS则使用关键字“BASEGEOGCRS”,详细对比如下表:
2、结构不同
ArcGIS的WKT坐标系定义中,投影坐标系“PROJCS”的属性包括名称、地理坐标系(基准面、椭球体)、中央子午线、投影名称和投影参数、坐标系单位。
QGIS的WKT坐标系定义中,投影坐标系“PROJCRS”的属性包括名称、地理坐标系(基准面、椭球体)、中央子午线、投影方式名称、投影方法、投影参数、平面坐标系、适用范围说明,而且在每个数值型参数后面均增加了所使用的单位说明,例如椭球体的半径使用长度单位米,在字符串表现为“LENGTHUNIT["metre",1]”。
3、名称不同
WKT坐标系定义中很多元素带有名称属性,如投影坐标系名称、地理坐标系名称、基准面名称,各个GIS厂商对名称的处理方式不同。早期WKT1标准中对名称的取值没有明确定义,导致各个厂商的实现方式出现不一致,影响了数据交换,新版标准WKT2则要求名称取值以EPSG为基准。
总的来说,QGIS与ArcGIS的WKT坐标系定义不同,主要是由于依据的WKT坐标系定义标准不同引起。
目前GIS软件处理WKT坐标系定义基于如下标准:
【1】ISO 19162:2019,Geographic information — Well-known text representation of coordinate reference systems,WKT坐标系定义的现行标准,对ISO 19111:2019中WKT 坐标系定义部分做了合并,简称WKT2。
【2】ISO 19111:2019,Geographic information — Referencing by coordinates,Proj坐标系定义与转换开源项目所依据标准,关于WKT坐标系定义部分已合并到ISO 19162:2019。
【3】OGC 01-009/ISO 19125-1:2004,Coordinate Transformation Services,Revision 1.00,2001年发布,是OGC 99-049的升级版本,简称WKT1。
【4】OGC 99-049,Simple Features Specification For SQL,Revision 1.1,WKT坐标系定义标准第一版。
根据查阅相关资料得知,ArcGIS主要依据WKT1坐标系定义标准,即OGC 01-009/ISO 19125-1:2004;QGIS则符合WKT2坐标系定义标准,即ISO 19162:2019。
WKT2标准在向后兼容方面做出明确规定:要求符合WKT2标准的GIS软件能够读取和导入WKT1坐标系定义,但不要求输出WKT1坐标系定义,这也就解释了为什么QGIS可以读取和导入ArcGIS版本的WKT坐标系定义。
常见GIS软件WKT坐标系定义所遵循的标准
QGIS/GDAL/OGR,以WKT语言为内置坐标系定义描述,导入力图兼容新旧标准以及ESRI格式,依据WKT2(ISO 19162:2019)标准转换和和导出。
ESRI,符合WKT1(OGC 01-009)标准。
PostGIS,WKT坐标系定义用OGR生成并存储于spatial_ref_sys表中,因此其遵循的标准与OGR相同。
FME,基于OGR实现WKT坐标系定义的读写操作,与QGIS/GDAL/OGR一样,力图兼容新老标准和ESRI格式的WKT坐标系定义。
超图SuperMap桌面版,坐标系导出格式为自定义的XML文档,尚未支持WKT坐标系定义导出。(也可能是我没找到,如果此处与实际不符,请留言指正。)
QGIS和ArcGIS对WKT坐标系定义的处理机制分析
QGIS对WKT坐标系定义的处理机制
QGIS对WKT坐标系定义的处理分成导入(Import)和导出(Export)两个部分,导入接收WKT格式字符串,经过拆分和识别,将其转换为QGIS的坐标参照系对象(QgsCoordinateReferenceSystem);导出则是将QGIS的坐标参照系对象转换为符合WTK2标准的字符串返回给用户。
A、导入处理过程
从WKT坐标系定义字符串到QGIS坐标参照系对象的转换主要是由QgsCoordinateReferenceSystem 类的createFromWktInternal()实现,需要传入WKT坐标系定义字符串和该字符串对应的描述信息,即我们看到的坐标系名称,流程如下:
1.判断WKT坐标系定义字符串是否为空,如果为空,则导入失败。
2.在缓存查找WKT坐标系定义字符串对象。
3.如果缓存中找到匹配,将对应的坐标参照系对象赋值给当前坐标参照系实例。
4.如果缓存中没有匹配,则判断WKT坐标系定义的来源厂商,目前兼容的厂商为:epsg|esri|osgeo|ignf|zangi|iau2000|postgis|internal|user,根据来源厂商分别用对应的标准解析传入的WKT坐标系定义字符串,生成坐标参照系对象。
5.如果无法匹配来源厂商,则按照一般WKT坐标系定义字符串处理。判断当前使用的Proj版本,主版本大于等于6,则使用Proj4库查找与传入WKT坐标系定义匹配的对象。
6.如果Proj版本小于6,则用ORG库查找传入的WKT坐标系定义。
7.如果以上过程没有找到与传入WKT坐标系定义匹配的预定义坐标系,则创建新的自定义坐标参照系,存入本地自定义坐标系数据库。
仍以2015年省级行政区划数据的坐标系定义导入过程为例,对应的PRJ文件WKT字符串如下:
PROJCS["Krasovsky_1940_Albers",
GEOGCS["GCS_Krasovsky_1940",
DATUM["D_Krasovsky_1940",
SPHEROID["Krasovsky_1940",6378245.0,298.3]],
PRIMEM["Greenwich",0.0],
UNIT["Degree",0.0174532925199433]],
PROJECTION["Albers"],
PARAMETER["False_Easting",0.0],
PARAMETER["False_Northing",0.0],
PARAMETER["Central_Meridian",105.0],
PARAMETER["Standard_Parallel_1",25.0],
PARAMETER["Standard_Parallel_2",47.0],
PARAMETER["Latitude_Of_Origin",0.0],
UNIT["Meter",1.0]]
以该WKT坐标系定义作为参数,QGIS先对该参数做基本的判别(是否为空),然后开始导入过程:
第一步,在缓存中查找该字符串,此处缓存匹配失败。
第二步,判断来源厂商。QGIS直接用字符串查找来匹配来源厂商,例如,在传入的WKT坐标系定义字符串中如果找到字符子串“ESRI”,则判断该定义来自ArcGIS,按照ArcGIS的标准解析,得到坐标参照系对象。此处没有出现所兼容的厂商标识字符串,所以按照一般WKT坐标系定义进行解析。
第三步,如果当前使用的Proj版本,主版本大于等于6,则用Proj库将WKT坐标系定义转换为Proj定义格式,生成坐标参照系对象。当前Proj库的最新版本为7.2,只要经常更新QGIS,一般Proj库的版本均大于等于6,也就是说,没有标明来源厂商的WKT坐标系定义一般由Proj库来处理。当Proj版本低于6时,则用GDAL/OGR库来解析传入的WKT坐标系定义参数。
无论是Proj库还是GDAL/OGR库,当前遵循的WKT标准均为最新的WKT2(ISO 19162:2019)。
第四步,如果以上过程没有找到与传入WKT坐标系定义匹配的预定义坐标系,则创建新的自定义坐标参照系,存入本地自定义坐标系数据库。
第五步,如果上述过程均不能解析传入的WKT坐标系定义,则导入失败,并在消息窗口打印出详细的错误信息。
导入成功后,在坐标系选择对话框的自定义坐标系分组中将看到解析成功后按照WTK2标准输出的WKT坐标系定义:
PROJCRS["Krasovsky_1940_Albers",
BASEGEOGCRS["Unknown datum based upon the Krassowsky 1940 ellipsoid",
DATUM["Not specified (based on Krassowsky 1940 ellipsoid)",
ELLIPSOID["Krassowsky 1940",6378245,298.3,
LENGTHUNIT["metre",1]],
ID["EPSG",6024]],
PRIMEM["Greenwich",0,
ANGLEUNIT["Degree",0.0174532925199433]]],
CONVERSION["unnamed",
METHOD["Albers Equal Area",
ID["EPSG",9822]],
PARAMETER["Latitude of false origin",0,
ANGLEUNIT["Degree",0.0174532925199433],
ID["EPSG",8821]],
PARAMETER["Longitude of false origin",105,
ANGLEUNIT["Degree",0.0174532925199433],
ID["EPSG",8822]],
PARAMETER["Latitude of 1st standard parallel",25,
ANGLEUNIT["Degree",0.0174532925199433],
ID["EPSG",8823]],
PARAMETER["Latitude of 2nd standard parallel",47,
ANGLEUNIT["Degree",0.0174532925199433],
ID["EPSG",8824]],
PARAMETER["Easting at false origin",0,
LENGTHUNIT["metre",1],
ID["EPSG",8826]],
PARAMETER["Northing at false origin",0,
LENGTHUNIT["metre",1],
ID["EPSG",8827]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]]]
由于示范数据所导入的WKT坐标系定义无法匹配原有预定义坐标系,因此产生了新的自定义坐标系,所创建自定义坐标系在本地的数据库中也能找到对应的记录。
本地数据库的位置一般存储在C:\Users\用户\AppData\Roaming\QGIS\QGIS3\profiles\default\qgis.db,这是一个SQLite数据库,打开可以看到一共有四个表:tbl_bookmarks(存储空间书签)、tbl_ellipsoid(预定义椭球体)、tbl_projection(存储投影方法名称、别名等信息)、tbl_srs(存储用户自定义坐标系),查询tbl_srs表可以看到我们自定义的坐标系:
B、导出处理过程
导出的流程相对简单,根据需要输出的WKT标准版本、是否有换行和缩进空格数三个参数,将坐标参照系的各个属性用关键字和取值的方式组合成符合WKT标准的字符串返回给用户。
支持导出的WKT标准
出于兼容性考虑,QGIS将WKT字符串分为11种类型:
ArcGIS对WKT坐标系定义的处理
上文已经提过,ArcGIS的WKT输出格式为WKT1,即符合OGC 01-009标准。该标准于1999年发布,2001年修改,存在定义不明确的地方,导致各个GIS在实施该标准中出现偏差。
例如:
ESRI WKT中,节点名称的取值一般用下划线连接,而GDAL/OGR的WKT1版本用空格分隔。
基准面(DATUM)的名称在EPSG名称的基础上增加“D_”开头。
PARAMETER 元素的名称不同,例如第一标准纬线在GDAL/OGR名称为“Latitude of 1st standard parallel”,而在ESRI中参数名称变为“standard_parallel_1”。另外,GDAL/OGR的参数名称首字母大写,ESRI全部用小写字母,实际上根据OGC标准,关键字不区分大小写,但是参数取值(引号中的内容)则区分大小写。
单位UNIT元素的取值在GDAL/OGR中为首字母大写,例如“Degree”,而ESRI则全部小写“degree”。
……
还有一些未提及的差异,例如对Lambert投影,ESRI将单标准纬线和双标准纬线视为一个投影类型,只是在参数上区分,而GDAL/OGR则区分为2个投影:LCC 1SP和LCC 2SP。
Melita Kennedy是负责开发ESRI地图投影和基准面模块的主要成员,也是CRS WKT 2.0标准起草委员会成员之一,根据她在多种场合对ESRI WKT处理方式的说明,可以大致得到ArcGIS对WKT处理:
WKT1标准的定义本身就存在许多不明确的地方,因此阻碍了WTK坐标系的共享;她希望能改进WKT的实现方式,也确实在ArcGIS 10.5.1和ArcGIS pro 2.0完成了对WKT2的支持,但是后面备注写道:仅在GeoPackage中实现。
ESRI WKT在对象名称和某些参数的取值虽然增加了下划线以方便阅读,但是在字符串比较中(如参数、名称字符串)则是忽略下划线和前缀字符串“GCSE_”和“D_”,且不区分大小写。
为了兼容其他厂商的WKT坐标系,ESRI维护了一个名称与别名的对照列表,但是没有公开该列表。
小结
就自定义坐标系而言,QGIS在兼容性和开放性上做得非常优秀,这也是后来者的优势:在成熟的标准上开发代码,自然会避免前辈们踩过的坑。但是,在进行自定义坐标系操作的用户体验方面,QGIS与商用GIS软件相比还有一定差距。
ArcGIS作为世界GIS行业的领先者,拥有庞大的用户群体,生产并管理着海量的GIS数据。很多时候,ArcGIS对一些问题的处理方式即为GIS行业的事实标准,其他GIS软件在兼容性方面普遍需要考虑对ArcGIS的兼容,这种状况应该还会持续一段时间。与此同时,随着OGC等机构主导的开放标准日趋完善,对最新开放标准的支持将成为所有GIS软件的必修课,这也意味着GIS行业必将走向更加开放和标准的未来。
版权声明
本文欢迎转载,转载时请注明出处。