项目RFC 6：基于三角剖分的变换

作者

甚至鲁奥

联系

even.rouault@spatialys.com

状态

采用

实施目标

项目7.2

最后更新

2020-09-02

总结

这个RFC增加了一个新的转换方法， tinshift （TIN表示不规则三角网）

这项工作的动机是能够处理芬兰国家土地调查局创造的官方转变，用于：

  • KKJ和ETRS89水平基准之间的水平转换

  • N43和N60高度以及N60和N2000高度之间的垂直转换。

这种变换在某种程度上与传统的基于网格的变换相关，只是校正值由三角剖分的顶点保持，而不是位于网格的节点上。

在许多情况下，三角剖分是测地平差的初始积，而栅格是导出的积。例如，瑞士网格有原始三角剖分的衍生产品。

基于网格的转换仍然非常方便使用，因为访问校正值非常简单和有效，因此基于三角剖分的转换并不意味着取代它们，而是作为一种补充，为了能够将经过官方审查的转换结果复制到毫米或更高的精度（这里说的是数值结果的再现性，而不是转换的物理精度，可能是厘米），这有时是必要的。通常可以用网格逼近三角剖分的结果，但这可能需要采用较小的网格步长，从而生成比原始三角剖分大得多的网格。

细节

变换

一种新的变换方法， tinshift ，已添加。这需要一个强制性的论点， file ，它指向一个JSON文件，其中包含三角剖分和关联的元数据。输入和输出坐标必须是地理坐标或投影坐标。根据JSON文件的内容，可以转换坐标的水平、垂直或两个分量。

转换的用法如下：

$ echo 3210000.0000 6700000.0000 0 2020 | cct +proj=tinshift +file=./triangulation_kkj.json

209948.3217     6697187.0009    0.0000     2020

该变换是可逆的，其计算复杂度与前向变换相同。

算法

内部， tinshift 将整个文件放入内存。人们认为三角测量应该足够小。上面提到的KKJ到ETRS89三角剖分适合65kb的JSON，用于1449个三角形和767个顶点。

当一个点被变换时，必须找到它所落入的三角形。我们不需要迭代所有的三角形，而是在内存中建立一个四叉树来加速候选三角形的识别。在上面提到的KKJ->ETRS89三角剖分上，这会将整个转换速度提高10倍。四叉树结构很好地兼顾了它带来的性能提升和实现的简单性（我们移植了来自GDAL的实现，继承了shapefile.spx空间索引的实现）。

为了确定一个点是否落入一个三角形，我们计算它的3 barycentric coordinates 从它的投影坐标来看， \(\lambda_i\) 对于 \(i=1,2,3\) . 它们是真实的价值（在 [0,1] 三角形内点的范围），给出三角形3个顶点中每个顶点的权重。

已知这些权重后，插值目标水平坐标就是将这些权重与三角形3个顶点处的目标水平坐标进行线性组合 (\(Xt_i\) 和 \(Yt_i\) ）：

\[X{target}=Xt\u 1 * \lambda_1 + Xt_2 * \lambda_2+Xt_3*\lambda_3 * \lambda_1 + Yt_2 * \λ2+钇3*\λ3\]

此插值在三角剖分的顶点处是精确的，并且在每个三角形内具有线性特性。它完全等价于三角插值的其他公式，例如

\[ \begin{align}\begin{aligned}X_{目标}=A+X_{源} * B + Y_{{source}} * C\\Y_{目标}=D+X_{源} * E + Y_{{source}} * F\end{aligned}\end{align} \]

式中，A、B、C、D、E、F常数（对于给定三角形）是通过求解2个3线性方程组，由三角形3个顶点的源坐标对和目标坐标对约束得到的：

\[Xt_i=A+Xs_i * B + Ys_i * C级 * E + Ys_i * F\]

备注

从实验中可以看出，当插值幅度为1e5/1e6量级的投影坐标时，使用重心坐标的插值在数值上更为稳健，因为计算涉及坐标的差异。而含有A、B、C、D、E、F的公式中，A和D常数的值往往较大，C和E的值子句为0，B和F的值子句接近1。然而，两种方法之间的差异在实际应用中可以忽略不计（低于微米精度）

与垂直坐标变换类似，其中 \(Zoff_i\) 是三角形每个顶点的垂直偏移：

\[Z_{目标}=Z_{源}+ZOFF_1 * \lambda_1 + Zoff_2 * \lambda_2+Zoff_3*\lambda_3\]

三角剖分的约束条件

没有检查三角测量的一致性。强烈建议三角形之间不要重叠（在考虑源坐标或正变换时，或反变换的目标坐标时），否则将选择哪个三角形是未指定的。除此之外，三角剖分不需要有特定的属性（比如Delaunay三角剖分）

文件格式

据我们所知，没有现成的文件格式将大地坐标变换表示为三角剖分。存储tin的潜在类似格式有 ITF 或 XMS . 它们都需要扩展，以便处理基准面偏移信息，因为它们都主要用于DEM。

因此，我们提出了一种基于文本的格式，使用JSON作为序列化。与二进制格式相比，使用基于文本的格式可能被认为是性能方面的限制，但是对于所考虑的三角剖分的大小（几千个三角形/顶点），没有问题。加载这样的文件大约需要20毫秒。作为参考，加载约115000个三角形和71000个顶点的三角剖分需要450ms。

使用JSON提供了通用的格式和解析规则，并方便从Python脚本创建它作为示例。这也可以很容易地由不支持JSON的编写者生成。

对于通用元数据，我们紧密地重用用于 Deformation model master file

下面是一个最小的示例，从KKJ到ETRS89转换，只有一个三角形：

{
    "file_type": "triangulation_file",
    "format_version": "1.0",
    "name": "Name",
    "version": "Version",
    "publication_date": "2018-07-01T00:00:00Z",
    "license": "Creative Commons Attribution 4.0 International",
    "description": "Test triangulation",
    "authority": {
        "name": "Authority name",
        "url": "http://example.com",
        "address": "Address",
        "email": "test@example.com"
    },
    "links": [
    {
        "href": "https://example.com/about.html",
        "rel": "about",
        "type": "text/html",
        "title": "About"
    },
    {
        "href": "https://example.com/download",
        "rel": "source",
        "type": "application/zip",
        "title": "Authoritative source"
    },
    {
        "href": "https://creativecommons.org/licenses/by/4.0/",
        "rel": "license",
        "type": "text/html",
        "title": "Creative Commons Attribution 4.0 International license"
    },
    {
        "href": "https://example.com/metadata.xml",
        "rel": "metadata",
        "type": "application/xml",
    "title": " ISO 19115 XML encoded metadata regarding the triangulation"
        }
    ],
    "input_crs": "EPSG:2393",
    "target_crs": "EPSG:3067",
    "transformed_components": [ "horizontal" ],
    "vertices_columns": [ "source_x", "source_y", "target_x", "target_y" ],
    "triangles_columns": [ "idx_vertex1", "idx_vertex2", "idx_vertex3" ],
    "vertices": [ [3244102.707, 6693710.937, 244037.137, 6690900.686],
                  [3205290.722, 6715311.822, 205240.895, 6712492.577],
                  [3218328.492, 6649538.429, 218273.648, 6646745.973] ],
    "triangles": [ [0, 1, 2] ]
}

因此，在通用元数据之后，我们定义了输入和输出CRS（仅供参考），并且转换会影响坐标的水平分量。我们命名 vertices 和 triangles 数组。我们定义了每个顶点的源坐标和目标坐标，并通过引用每个顶点的索引来定义一个三角形 vertices 数组。

对于文件，更具体的项目是：

input_crs

标识顶点中源坐标的CR的字符串。典型的 EPSG:XXXX . 如果转换是针对垂直分量，则这应该是复合CRS的代码（可以是EPSG:XXXX+YYYY其中XXXX是水平CRS的代码，YYYY是垂直CRS的代码）。例如，对于KKJ->ETRS89转换，这是EPSG:2393 (KKJ / Finland Uniform Coordinate System ). 假设输入坐标按“可视化标准化”/“地理信息系统友好”顺序传递，即地理坐标为经度、纬度，投影坐标为东距、北距。

output_crs

标识顶点中目标坐标的CR的字符串。典型的 EPSG:XXXX . 如果转换是针对垂直分量，则这应该是复合CRS的代码（可以是EPSG:XXXX+YYYY其中XXXX是水平CRS的代码，YYYY是垂直CRS的代码）。例如，对于KKJ->ETRS89转换，这是爱普生：3067（“ETRS89/TM35FIN（E，N）”）。输出坐标将以“可视化标准化”/“地理信息系统友好”的顺序返回，即地理坐标为经度、纬度，投影坐标为东距、北距。

transformed_components

可以包含一个或两个字符串的数组：坐标的水平分量变换时为“水平”和/或垂直分量变换时为“垂直”。

vertices_columns

指定中行的列的名称 vertices 数组。一定有很多元素 vertices_columns 就像在一排 vertices . 以下名称具有特殊含义： source_x ， source_y ， target_x ， target_y ， source_z ， target_z 和 offset_z . source_x 和 source_y 是强制性的。 source_x 表示源经度（度）或东距。 source_y 表示源纬度（度）或北距。 target_x 和 target_y 是强制性的 horizontal 在中指定 transformed_components . (source_z 和 target_z 或 offset_z 是强制性的 vertical 在中指定 transformed_components

triangles_columns

指定中行的列的名称 triangles 数组。一定有很多元素 triangles_columns 就像在一排 triangles . 以下名称具有特殊含义： idx_vertex1 ， idx_vertex2 ， idx_vertex3 . 它们是强制性的。

顶点

一种数组，其项本身就是数组，列数与中所述的相同 vertices_columns .

三角形

一种数组，其项本身就是数组，列数与中所述的相同 triangles_columns . 的值 idx_vertexN 列必须是索引（介于0和len之间） (vertices -1） 项目的 vertices 数组。

代码影响

在src/transformations中添加了以下新文件：

  • tinshift.cpp文件：用于定义新操作的项目特定代码。在需要时负责输入和输出坐标转换（在输入坐标系和三角测量坐标系、源坐标系和目标坐标系和输出坐标系之间）。

  • tinshift.hpp公司：基于标头的实现。此文件包含API。

  • tinshift公司_例外.hpp：文件分析期间可能引发的异常

  • tinshift公司_impl.hpp公司：实现文件加载、三角搜索和插值。

这是变形模型实现所遵循的方法，它使单元测试变得更容易。

src公司/四叉树.hpp包含四叉树实现。

性能指标

在Intel（R）Core（TM）i7-6700HQ CPU@2.60GHz上测试，可转换400万个点

对于KKJ到ETRS89的变换（1449个三角形和767个顶点），可以每秒变换440万个点。

相比之下，基于Helmert的KKJ到ETRS89转换的速度为160万点/秒。

一个包含约115000个三角形和71000个顶点的三角剖分以220万点/秒的速度运行（由于三角剖分的初始加载在这里是不可忽略的，所以在更多点上的吞吐量会更好）

向后兼容性

新功能完全向后兼容。

测试

PROJ测试套件将被增强以测试新的转换，用一个新的.GIE文件和一个C++单元测试来测试较低的级别。

文档

  • tinshift方法将记录在案。

  • JSON格式将记录在https://proj.org/specifications/

  • 还将提供JSON模式。

拟议实施

初始实现可在https://github.com/rouault/PROJ/tree/tinshift

工具书类

Finnish coordinate transformation (automated translation to English)

收养状况

RFC于2020-09-02通过，并由以下PSC成员提供+1

  • 克里斯蒂安·埃弗斯

  • 查尔斯·卡尼

  • 托马斯·克努森

  • 甚至鲁奥

基金

这项工作是由 National Land Survey of Finland