RFC 73:WKT2的PROJ6集成、后期绑定功能、时间支持和统一CRS数据库

作者:

甚至鲁奥

联系方式:

even.rouault@spatialys.com网站

起动:

2019年1月8日

上次修改时间:

2019年5月2日

状态:

在GDAL 3.0中实现

总结

该文件描述了PROJ 6与GDAL集成的相关工作,GDAL增加了不同的功能:对CRS WKT 2版本的支持,CRS之间坐标转换的“后期绑定”功能,坐标操作的时间维度支持,以及统一CRS数据库的使用。

动机

动机是那些暴露在 https://gdalbarn.com/#why ,此处复制

GDAL、PROJ和libgeotiff中的坐标系缺少现代功能,需要彻底重构:

  • PROJúu LIB和GDALúu数据中可怕的特殊CSV数据库让用户感到沮丧,给开发人员带来挑战,并阻碍定义的互操作性。

  • GDAL和PROJ缺少OGC WKT2支持。

  • PROJ 5.0+不再需要通过WGS84进行数据转换枢轴,WGS84会引入高达2米的误差,但其他工具不会利用这一点。

CSV数据库

在EPSG和其他定义中使用基于SQLite的数据库将允许项目添加更多功能(区域感知验证),将项目的自定义特殊数据结构转换为更通用的可消费结构,并促进许多坐标系处理软件工具之间的定义互操作性。

工作周2

OGC WKT2 修复了长期存在的互操作性坐标系定义差异。WKT2包含用于描述时间相关坐标系的工具。PROJ 5+现在能够进行依赖时间的转换,但是GDAL和其他工具还不支持它们。

一些国家正在更新其大地测量基础设施,以纳入与时间有关的坐标系。例如,澳大利亚和美国分别在2020年和2022年调整了与时间相关的坐标系。北美熟悉的NAD83和NAVD88被NATF2022和NAPGD2022所取代,行业迟早要适应这些挑战。

WGS84枢轴

项目之前需要通过7参数转换通过WGS84旋转的数据转换。这个轴心是一个实用的解决方案,但它会引入大约两米的误差,许多传统的基准不能用WGS84定义。PROJ 5+现在提供了支持后期绑定的工具 transformation pipeline framework ,但是GDAL和其他工具还不能使用它。更高精度的转换避免了步进WGS84,并消除了来自当地大地测量机构的边车数据的额外转换步骤。

提案

第三方库要求

GDAL master(future 3.0)需要PROJ master(future PROJ 6.0)和libgeotiff master(future libgeotiff 1.5或2.0)来构建和执行。

关于PROJ,GDAL master中不会嵌入PROJ的内部副本。支持旧版本的PROJ是不可行的,因为OGRSpatialReference类已经被大量重写,以利用从GDAL完全移到PROJ的功能:PROJ字符串导入和导出、WKT字符串导入和导出、EPSG数据库利用。为了能够在复杂的设置中更容易地使用GDAL master和PROJ master,其中一些GDAL依赖项使用系统提供的libproj,如果将天真的PROJ master和这个旧的libproj混合在一起会导致运行时崩溃,可以使用CFLAGS/CXXFLAGS=-DPROJ_RENAME_符号来构建PROJ master,并将GDAL能够使用此自定义生成。请注意,这并不打算长期使用,因为正确的打包解决方案最终将使用PROJ 6来重建其所有的反向依赖关系。还应注意,PROJ在configure/nmake时是必需的,即通过dlopen()/LoadLibrary()在运行时动态加载不再可用。

关于libgeotiff,frmts/gtiff/libgeotiff中的内部副本已经用上游libgeotiff master的内容刷新。

所有持续集成系统(Travis CI和AppVeyor)都已更新,作为GDAL构建的一部分来构建PROJ master。

OGRSpatialReference重写

OGRSpatialReference类是GDAL/OGR中所有坐标参考系(CRS)操作的中心。在GDAL 2.4之前,这个类主要包含WKT 1表示的OGR-SRSNode根节点,所有getter和setter都操作这个树表示。作为这项工作的一部分,OGRSpatialReference内部包含的主对象现在是一个PROJ PJ对象,并且方法在此PJ对象上调用PROJ C API getter和setter。这使得,主要是(*),表示独立。

WKT1,WKT2,ESRI WKT,PROJ字符串导入和导出现在委托给PROJ。从EPSG数据库导入CRS也一样,现在它依赖于proj.db SQLite数据库。因此所有的数据/ * 包含CRS相关信息的.csv文件已从GDAL中删除。需要注意的是,当导入WKT时,来自ESRI WKT的“变形”现在是自动完成的。

虽然IsSame()或FindMatches()等方法的一般语义保持不变,但底层实现却大不相同,这在某些情况下可能会导致与以前的GDAL版本不同的结果。在FindMatches()的情况下,由于增强了数据库中的查询功能,CRS到EPSG条目的标识通常会得到改进。

(*)这里有“大多数”的精度,因为在每个地方重写并不实际。因此,对于某些方法,仍然会执行内部WKT1导出。对于以SRS节点的路径(如“GEOGCS | UNIT”)作为参数的方法,或者一些期望OGC WKT1特定名称的方法(如SetProjection()、GetProjParm())。那些被认为主要是用来驾驶的。将它们改为EPSG名称会影响许多驱动程序,其中一些驱动程序很少测试SRS支持,而且它们大多只支持WKT1表示。

坐标变换

由于GDAL 2.3和初始PROJ 5支持,在两个CRS之间转换时,我们仍然依赖源CRS和目标CRS的PROJ.4字符串导出来创建坐标操作管道。因此,这仅限于“早期绑定”操作,即通过to WGS84或nadgrids PROJ关键字使用WGS84枢轴。现在,PROJ使用了在两个CRS之间找到合适的坐标操作的新功能,提供了“后期绑定”功能,以考虑WGS84或使用区域以外的其他轴心点。

OGRCreateCoordinateOperation()现在使用额外的可选参数来定义选项。

其中一个选项是定义在搜索候选操作时要考虑的感兴趣区域。如果多个操作匹配,将选择“最佳”(根据项目排序标准)。注意:即使以后调用Transform()时使用的坐标超出了最初关注的区域,也会系统地使用它。

另一个选项是指定要应用的坐标操作的能力,以便替代将自动计算的GDAL/PROJ,可以是PROJ字符串(通常是a+PROJ=pipeline),也可以是WKT坐标操作/串联操作。用户通常可以使用新的PROJ projinfo实用程序选择特定的坐标操作,该实用程序可以从源/目标元组返回候选操作。

当没有指定选项时,GDAL将使用PROJ列出所有候选坐标操作。对于每次调用Transform(),它将计算输入坐标的平均坐标,并使用它从候选坐标中确定最佳坐标操作。

Transform()方法现在需要一个额外的参数来包含与时间相关的坐标操作的坐标纪元(通常是十进制年份值)。与此相关的是,gdalwarp通常使用的GDALTransform机制的转换选项现在接受用于相同目的的坐标纪元。

OGRSpatialReference在GDAL中的应用

目前,GDAL数据集接受并返回WKT 1字符串来描述SRS。为了更独立于实际编码,并且例如允许GeoPackage栅格数据集能够使用WKT 2,最好能够附加不依赖于表示的SRS(WKT 1或WKT 2),因此使用OGRSpatialReference对象而不是const char * 字符串。

在GDALDataset中添加了以下新方法:

  • 虚拟const OGRSpatialReference * GetSpatialRef()常量;

  • 虚拟cplersetspatialref(const OGRSpatialReference*);

  • 虚拟const OGRSpatialReference * GetGCPSpatialRef()常量;

  • 虚拟CPLErr SetGCPs(int nGCPCount,const GDAL_GCP pasGCPList,const OGRSpatialReference公司

为了简化转换,在GDALDataset中添加了以下非虚拟方法:

  • const OGRSpatialReference公司 * GetSpatialRefFromOldGetProjectionRef()常量;

  • CPLErr OldSetProjectionFromSetSpatialRef(const OGRSpatialReference * 波斯斯);

  • const OGRSpatialReference公司 * getgcpspatialrefromoldgetgcpproject()常量;

  • CPLErr OldSetGCPsFromNew(int nGCPCount,const GDAL_GCP * pasGCPList,const OGRSpatialReference公司 * poGCP_SRS);

前面的GetProjectionRef()、SetProjection()、GetGCPProjection()和SetGCPs()可用作投影虚拟方法,前缀为下划线

通过这种方式转换现有的驱动程序,需要将其GetProjectionRef()方法重命名为\u GetProjectionRef(),并添加:

const OGRSpatialReference* GetSpatialRef() const override {
    return GetSpatialRefFromOldGetProjectionRef();
}

默认WKT版本

不带选项的OGRSpatialReference::exportToWkt()将报告WKT 1(带显式轴节点)。请参见下面的“Axis order issues”段落)了解此表示的CRS兼容项,否则请使用WKT2:2018(通常用于地理3D CRS)。

exportToWkt()的增强版本接受用于指定所使用的确切WKT版本的选项(如果必须使用多行或单行输出等)。

或者,OSR_WKT_FORMAT configuration选项可用于修改exportToWk()使用的WKT版本(当exportToWkt()选项中未传递显式版本时)

gdalinfo、ogrinfo和gdalsrsinfo实用程序将默认输出WKT2:2018

轴订单问题

这是一个反复出现的痛点。这个RFC提出了一个新的方法(不假装完全解决它)来完成最初的 RFC 20: OGRSpatialReference Axis Support . 问题是,CRS官方定义使用的轴顺序与GIS应用程序中栅格或矢量数据的传统编码方式不一致。典型的例子是来自EPSG的地理“WGS 84”定义,EPSG:4326使用纬度作为第一轴,使用经度作为第二轴。RFC 20决定,默认情况下,当管理局的AXIS命令与地理信息系统友好命令不匹配时,AXIS定义将从WKT中删除(并使用自定义EPSGA管理局使WKT具有官方AXIS元素)

这在技术上是可能的,因为WKT 1语法定义了AXIS元素。然而,删除轴的定义是一个潜在的混淆源,因为不清楚实际使用的是哪个轴顺序。此外,在WKT2中,轴元素是必需的,内部PROJ表示还需要定义坐标系。所以有两个不令人满意的选择:

  • 返回带有地理信息系统友好顺序的官方定义修补版本,同时仍使用官方权限代码。因为我们保留了与源代码的链接,所以很实用,但是因为我们修改了它,所以是个谎言。用户不知道他们是否必须信任编码的订单,或者来自权威机构的正式订单。

  • 以地理信息系统友好顺序返回官方定义的修补版本,但不返回官方授权代码。这是符合的,但我们会失去与权威代码的链接。

本RFC提出的解决方案是添加一个“数据轴到SRS轴映射”的概念,这有点类似于在WCS DescribeCoverage response中所做的工作,以解释SRS轴如何映射到覆盖范围的网格轴

提取自 https://docs.geoserver.org/stable/en/user/extensions/wcs20eo/index.html 对于使用EPSG:4326的覆盖范围

<gml:coverageFunction>
  <gml:GridFunction>
    <gml:sequenceRule axisOrder="+2 +1">Linear</gml:sequenceRule>
    <gml:startPoint>0 0</gml:startPoint>
  </gml:GridFunction>
</gml:coverageFunction>

添加了一个类似的映射来定义地理变换矩阵或几何测量中的“x”和“y”分量如何映射到CRS定义所定义的轴。

这种映射是通过OGRSpatialReference中的一种新方法给出的

const std::vector<int>& GetDataAxisToSRSAxisMapping() const

为了解释它的语义,假设它返回2,-1,3。这被解释为:

  • 2: CRS的第一轴映射到数据的第二轴

  • -1: CRS的第二个轴映射到数据的第一个轴,值取反

  • 3: CRS的第三轴映射到数据的第三轴

这类似于PROJ axisswap操作: https://proj4.org/operations/conversions/axisswap.html

默认情况下,在新创建的OGRSpatialReference对象上,GetDataAxisToSRSAxisMapping()返回标识1,2 [,3] ,即符合管理局定义的轴顺序。

由于所有GDAL和绝大多数OGR驱动程序都依赖于使用“GIS轴映射”,因此添加了SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER或OAMS_AUTHORITY_COMPLIANT或OAMS_CUSTOM)方法,以使指定轴映射的工作更容易;

OAMS_TRADITIONAL_GIS_ORDER是指:

  • 对于地理二维CRS,

    • 对于纬度北,经度东(例如EPSG:4326),GetDataAxisToSRSAxisMapping()返回{2,1},这意味着数据顺序是经度,纬度

    • 对于东经、北纬(例如OGC:CRS84),返回{1,2}

  • 对于预计的CRS,

    • 对于东、北(即大多数预计的CRS),返回{1,2}

    • 北,东,返回{2,1}

    • 对于具有东/南、北/南的北极CRS,例如EPSG:5041(“WGS 84/UPS North(E,N)”),将返回{1,2}

    • 对于具有北/南、东/南的北极CRS,例如EPSG:32661(“WGS 84/UPS North(N,E)”),将返回{2,1}

    • 同样适用于南极CRS

    • 对于所有其他情况,返回{1,2}

OGRCreateCoordinateTransformation()现在执行数据轴到srs轴的映射。

注意:与我在上一封电子邮件中指出的相反,gdaltransform行为没有改变,因为gdaltransform机制在内部强制执行GIS友好顺序。

栅格数据集被修改为在OGRSpatialReference上调用SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER) * 它们返回,并在SetSpatialRef()中假定它(现在假定且未选中)

矢量层通常都调用OGRSpatialReference上的SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER) * 由GetSpatialRef()返回。在GML驱动程序的情况下,如果用户定义了反转轴顺序如果长打开选项,轴交换不会完成(如前所述),并且使用符合权限策略。当接收到OGRSpatialReference时 * 可能决定(大多数人都会这么做)更改轴映射策略。也就是说:如果它接收到一个符合权威指令的OGRSpatialReference,它可能会决定切换到传统的ogrspatialref,并且getspaceref()::GetDataAxisToSRSAxisMapping()将反映这一点。在这种情况下,修改ogr2ogr以进行几何轴交换。

与此更改相关,WKT 1导出现在总是返回AXIS元素,因此EPSG:xxxx的行为与EPSGA:xxxx相同

因此,这种方法的总结是,在正式的SRS定义中,我们不再对轴顺序进行减损,而是添加了一个额外的接口来描述如何使我们的匹配与SRS定义匹配。

驱动程序更改

通过GetProjectionRef()、SetProjection()、getgcpproject()和SetGCPs()方法将SRS作为WKT字符串返回/接受的栅格驱动程序已升级为使用新的虚拟方法,在大多数情况下使用兼容层。

GDALPamDataset(PAM.aux.xml文件)和GDAL VRT驱动程序已经完全升级,以支持新接口,并序列化/反序列化数据轴到SRS轴的映射值。

GeoPackage驱动程序现在完全支持官方的“gpkg_crs_wkt”扩展,该扩展用于将wkt 2字符串定义存储在gpkg_space_ref_sys表中。当SRS可以编码为WKT1字符串时,驱动程序尝试不使用扩展名,并且如果插入了需要WKT2(通常是地理3D CRS)的SRS,将自动将“definition_12_063”列添加到现有的gpkg_spatial_6refu sys表中。

公用设施的变化

  • 每当报告CRS时,gdalinfo和ogrinfo都会报告数据轴到CRS轴的映射。默认情况下,它们还将输出WKT2-u2018,除非指定“-wkt-u format wkt1”。

Driver: GTiff/GeoTIFF
Files: out.tif
Size is 20, 20
Coordinate System is:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["unknown"],
        AREA["World"],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1 <-- here
Origin = (2.000000000000000,49.000000000000000)
Pixel Size = (0.100000000000000,-0.100000000000000)
  • gdalwarp、ogr2ogr和gdaltransform获得了一个-ct开关,高级用户可以使用它来指定坐标操作,可以是PROJ字符串(通常是a+PROJ=pipeline),也可以是WKT坐标操作/串联操作,如上面“ogrcoordinatedransformchanges”一段中所述。注意:管道必须考虑CRS的轴顺序,即使底层栅格/矢量驱动程序使用“地理信息系统友好”顺序。例如,当从EPSG:4326转换到EPSG:32631时,“+proj=pipeline+step+proj=axisswap+order=2,1+step+proj=unitconvert+xy-in=deg+xy-out=rad+step+proj=utm+zone=31+ellps=WGS84”。

  • gdalsrsinfo得到了增强,能够指定两个新的受支持的WKT变体:WKT2u 2015和WKT2u 2018。它将默认为输出WKT2U 2018

SWIG绑定更改

增强的ExportToWkt()和OGRCoordinateTransformation方法可通过SWIG绑定获得。对于非Python语言,可能需要额外的类型映射(特别是对于4dx,Y,Z,时间坐标的支持)

向后兼容性

这项工作的目的是 主要地 向后兼容,但不可避免的差异将被发现。例如,WKT 1和PROJ字符串导出已在PROJ中完全重写,因此虽然有望与GDAL 2.4或更早版本生成的字符串等效,但这并不完全相同:有效位数、PROJ参数的顺序、舍入等。。。

已更新MIGRATION_GUIDE.TXT以反映一些差异:

  • osrimportfromepg()考虑了官方轴顺序。

  • 删除OPTGetProjectionMethods()、OPTGetParameterList()和OPTGetParameterInfo()没有等效项。

  • 删除OSRFixup()和OSRFixupOrdering():不再需要,因为构造的对象始终有效

  • 删除OSRStripCTParms()。使用OSRExportToWktEx()而不是FORMAT=SQSQL选项

  • exportToWkt()输出轴节点

  • OSRIsSame():现在考虑数据轴到CRS轴的映射,除非IGNORE_data_axis_to_SRS_axis_mapping=YES被设置为OSRIsSameEx()的选项

  • ogr_srs_api.h:srs_WKT_WGS84宏在默认情况下不再声明,因为没有轴的WKT太不明确。首选修复方法:使用SRS_WKT_WGS84_LAT_LONG。或者在包含ogr_SRS_api.h之前定义USE_DEPRECATED_SRS_WKT_WGS84

引入新的虚拟方法GetSpatialRef()、SetSpatialRef()、GetGCPSpatialRef()和SetGCPs(…,const OGRSpatialReference)将影响树外栅格驱动程序 * 以及使用WKT字符串而不是OGRSpatialReference删除旧的等价项 * 实例。

文档

新方法已形成文件,现有方法的文件已在开发过程中适当更改。也就是说,需要更彻底的通过。教程也必须更新。

测试

由于多种原因(在WKT中导出的轴节点、WKT中的差异和PROJ字符串生成)导致预期结果发生更改,因此自动测试套件已在许多地方进行了调整。已经为新功能添加了新的测试。

需要注意的是,autotest不一定检查所有内容,问题已经通过手动测试发现并解决。“数据轴到CRS轴映射”概念的引入也很容易出错,因为它需要在许多不同的地方设置OAMS传统的GIS顺序策略。

因此,一旦这项工作在master中实现,就请用户和开发人员彻底测试GDAL。

实施

即使是鲁奥, Spatialys . 根据提供 PR 1185 资金来源 gdalbarn 赞助。

虽然它是作为多个提交提供给“更容易”审阅的,但由于在开发过程中发生了项目符号重命名(这将破坏可分性),中间步骤并不都是可构建的,因此它可能会被压缩到一个提交中以包含在master中。

投票历史

由PSC成员HowardB、JukkaR、DanielM和Ever+1采用

修改

2019年5月2日:GDAL 2.5改为GDAL 3.0