Add RFC-6 text

2 files changed, 369 insertions, 0 deletions

diff --git a/docs/source/community/rfc/index.rst b/docs/source/community/rfc/index.rst
index ced162e2..bb8b1455 100644
--- a/docs/source/community/rfc/index.rst
+++ b/docs/source/community/rfc/index.rst
@@ -16,3 +16,4 @@ the project.
    rfc-3
    rfc-4
    rfc-5
+   rfc-6
diff --git a/docs/source/community/rfc/rfc-6.rst b/docs/source/community/rfc/rfc-6.rst
new file mode 100644
index 00000000..a4437962
--- /dev/null
+++ b/docs/source/community/rfc/rfc-6.rst
@@ -0,0 +1,368 @@
+.. _rfc6:
+
+====================================================================
+PROJ RFC 6: Triangulation-based transformations
+====================================================================
+
+:Author: Even Rouault
+:Contact: even.rouault@spatialys.com
+:Status: Adopted
+:Implementation target: PROJ 7.2
+:Last Updated: 2020-09-02
+
+Summary
+-------------------------------------------------------------------------------
+
+This RFC adds a new transformation method, ``tinshift`` (TIN stands for
+Triangulated Irregular Network)
+
+The motivation for this work is to be able to handle the official transformations
+created by National Land Survey of Finland, for:
+
+- horizontal transformation between the KKJ and ETRS89 horizontal datums
+- vertical transformations between N43 and N60 heights, and N60 and N2000 heights.
+
+Such transformations are somehow related to traditional grid-based transformations,
+except that the correction values are hold by the vertices of the triangulation,
+instead of being at nodes of a grid.
+
+Triangulation are in a number of cases the initial product of a geodesic adjustment,
+with grids being a derived product. The Swiss grids have for example
+derived products of an original triangulation.
+
+Grid-based transformations remain very convenient to use because accessing
+correction values is really easy and efficient, so triangulation-based transformations
+are not meant as replacing them, but more about it being a complement, that is
+sometimes necessary to be able to replicate the results of a officially vetted
+transformation to a millimetric or better precision (speaking here about reproducibility
+of numeric results, rather than the physical accuracy of the transformation that
+might rather be centimetric). It is always possible to approach the result of
+the triangulation with a grid, but that may require to adopt a small grid step,
+and thus generate a grid that can be much larger than the original triangulation.
+
+Details
+-------------------------------------------------------------------------------
+
+Transformation
+++++++++++++++
+
+A new transformation method, ``tinshift``, is added. It takes one mandatory
+argument, ``file``, that points to a JSON file, which contains the triangulation
+and associated metadata. Input and output coordinates must be geographic or projected.
+Depending on the content
+of the JSON file, horizontal, vertical or both components of the coordinates may
+be transformed.
+
+The transformation is used like:
+
+::
+
+    $ echo 3210000.0000 6700000.0000 0 2020 | cct +proj=tinshift +file=./triangulation_kkj.json
+
+    209948.3217     6697187.0009    0.0000     2020
+
+The transformation is invertible, with the same computational complexity than
+the forward transformation.
+
+Algorithm
++++++++++
+
+Internally, ``tinshift`` ingest the whole file into memory. It is considered that
+triangulation should be small enough for that. The above mentioned KKJ to ETRS89
+triangulation fits into 65 KB of JSON, for 1449 triangles and 767 vertices.
+
+When a point is transformed, one must find the triangle into which it falls into.
+Instead of iterating over all triangles, we build a in-memory quadtree to speed-up
+the identification of candidates triangles. On the above mentioned KKJ -> ETRS89
+triangulation, this speeds up the whole transformation by a factor of 10. The
+quadtree structure is a very good compromise between the performance gain it brings
+and the simplicity of its implementation (we have ported the implementation coming
+from GDAL, inherit from the one used for shapefile .spx spatial indices).
+
+To determine if a point falls into a triangle, one computes its 3
+`barycentric coordinates <https://en.wikipedia.org/wiki/Barycentric_coordinate_system#Conversion_between_barycentric_and_Cartesian_coordinates>`_
+from its projected coordinates, :math:`\lambda_i` for :math:`i=1,2,3`.
+They are real values (in the [0,1] range for a point inside the triangle),
+giving the weight of each of the 3 vertices of the triangles.
+
+Once those weights are known, interpolating the target horizontal
+coordinate is a matter of doing the linear combination of those weights with
+the target horizontal coordinates at the 3 vertices of the triangle (:math:`Xt_i` and :math:`Yt_i`):
+
+.. math::
+
+    X_{target} = Xt_1 * \lambda_1 + Xt_2 * \lambda_2 + Xt_3 * \lambda_3
+
+    Y_{target} = Yt_1 * \lambda_1 + Yt_2 * \lambda_2 + Yt_3 * \lambda_3
+
+This interpolation is exact at the vertices of the triangulation, and has linear properties
+inside each triangle. It is completely equivalent to other formulations of
+triangular interpolation, such as
+
+.. math::
+
+    X_{target} = A + X_{source} * B  + Y_{source} * C
+
+    Y_{target} = D + X_{source} * E  + Y_{source} * F
+
+where the A, B, C, D, E, F constants (for a given triangle) are found by solving
+the 2 systems of 3 linear equations, constraint by the source and target coordinate pairs
+of the 3 vertices of the triangle:
+
+.. math::
+
+    Xt_i = A + Xs_i * B + Ys_i * C
+
+    Yt_i = D + Xs_i * E + Ys_i * F
+
+
+.. note:: 
+
+    From experiments, the interpolation using barycentric coordinates is slightly
+    more numerically robust when interpolating projected coordinates of amplitude of the
+    order of 1e5 / 1e6, due to computations involving differences of coordinates.
+    Whereas the formulation with the A, B, C, D, E, F tends to have big values for
+    the A and D constants, and values clause to 0 for C and E, and close to 1 for
+    B and F. However, the difference between the two approaches is negligible for
+    practical purposes (below micrometre precision)
+
+Similarly for a vertical coordinate transformation, where :math:`Zoff_i` is the vertical
+offset at each vertex of the triangle:
+
+.. math::
+
+    Z_{target} = Z_{source} + Zoff_1 * \lambda_1 + Zoff_2 * \lambda_2 + Zoff_3 * \lambda_3
+
+Constraints on the triangulation
+++++++++++++++++++++++++++++++++
+
+No check is done on the consistence of the triangulation. It is highly
+recommended that triangles do not overlap each other (when considering the
+source coordinates or the forward transformation, or the target coordinates for
+the inverse transformation), otherwise which triangle will be selected is
+unspecified. Besides that, the triangulation does not need to have particular
+properties (like being a Delaunay triangulation)
+
+File format
++++++++++++
+
+To the best of our knowledge, there are no established file formats to convey
+geodetic transformations as triangulations. Potential similar formats to store TINs
+are `ITF <http://vterrain.org/Implementation/Formats/ITF.html>`_ or
+`XMS <https://www.xmswiki.com/wiki/TIN_Files>`_.
+Both of them would need to be extended in order to handle datum shift information,
+since they are both intended for mostly DEM use.
+
+We thus propose a text-based format, using JSON as a serialization. Using a text-based
+format could potentially be thought as a limitation performance-wise compared to
+binary formats, but for the size of triangulations considered (a few thousands triangles / vertices),
+there is no issue. Loading such file is a matter of 20 milliseconds or so. For reference,
+loading a triangulation of about 115 000 triangles and 71 000 vertices takes 450 ms.
+
+Using JSON provides generic formatting and parsing rules, and convenience to
+create it from Python script for examples. This could also be easily generated "at hand"
+by non-JSON aware writers.
+
+For generic metadata, we reuse closely what has been used for the
+`Deformation model master file <https://github.com/linz/deformation-model-format>`_
+
+Below a minimal example, from the KKJ to ETRS89 transformation, with just a
+single triangle:
+
+.. code-block:: json
+
+    {
+        "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": "Adress",
+            "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] ]
+    }
+
+So after the generic metadata, we define the input and output CRS (informative
+only), and that the transformation affects horizontal components of
+coordinates. We name the columns of the ``vertices`` and ``triangles`` arrays.
+We defined the source and target coordinates of each vertex, and define a
+triangle by referring to the index of its vertices in the ``vertices`` array.
+
+More formally, the specific items for the triangulation file are:
+
+input_crs
+  String identifying the CRS of source coordinates
+  in the vertices. Typically ``EPSG:XXXX``. If the transformation is for vertical
+  component, this should be the code for a compound CRS (can be EPSG:XXXX+YYYY
+  where XXXX is the code of the horizontal CRS and YYYY the code of the vertical CRS).
+  For example, for the KKJ->ETRS89 transformation, this is EPSG:2393
+  (``KKJ / Finland Uniform Coordinate System``). The input coordinates are assumed
+  to be passed in the "normalized for visualisation" / "GIS friendly" order,
+  that is longitude, latitude for geographic coordinates and
+  easting, northing for projected coordinates.
+
+
+output_crs
+  String identifying the CRS of target coordinates in the vertices.
+  Typically ``EPSG:XXXX``. If the transformation is for vertical component,
+  this should be the code for a compound CRS (can be EPSG:XXXX+YYYY where
+  XXXX is the code of the horizontal CRS and YYYY the code of the vertical CRS).
+  For example, for the KKJ->ETRS89 transformation, this is EPSG:3067
+  (\"ETRS89 / TM35FIN(E,N)\"). The output coordinates will be returned in
+  the "normalized for visualisation" / "GIS friendly" order,
+  that is longitude, latitude for geographic coordinates and
+  easting, northing for projected coordinates.
+
+
+transformed_components
+  Array which may contain one or two strings: "horizontal" when horizontal
+  components of the coordinates are transformed and/or "vertical" when the
+  vertical component is transformed.
+
+
+vertices_columns
+  Specify the name of the columns of the rows in the ``vertices``
+  array. There must be exactly as many elements in ``vertices_columns`` as in a
+  row of ``vertices``. The following names have a special meaning: ``source_x``,
+  ``source_y``, ``target_x``, ``target_y``, ``source_z``, ``target_z`` and
+  ``offset_z``.  ``source_x`` and ``source_y`` are compulsory.
+  ``source_x`` is for the source longitude (in degree) or easting.
+  ``source_y`` is for the source latitude (in degree) or northing. 
+  ``target_x`` and ``target_y`` are compulsory when ``horizontal`` is specified
+  in ``transformed_components``. (``source_z`` and ``target_z``) or
+  ``offset_z`` are compulsory when ``vertical`` is specified in ``transformed_components``
+
+
+triangles_columns
+  Specify the name of the columns of the rows in the
+  ``triangles`` array. There must be exactly as many elements in ``triangles_columns``
+  as in a row of ``triangles``. The following names have a special meaning:
+  ``idx_vertex1``, ``idx_vertex2``, ``idx_vertex3``. They are compulsory.
+
+
+vertices
+  An array whose items are themselves arrays with as many columns as
+  described in ``vertices_columns``.
+
+
+triangles
+  An array whose items are themselves arrays with as many columns as
+  described in ``triangles_columns``.
+  The value of the ``idx_vertexN`` columns must be indices
+  (between 0 and len(``vertices``-1) of items of the ``vertices`` array.
+
+Code impacts
+++++++++++++
+
+The following new files are added in src/transformations:
+
+- tinshift.cpp: PROJ specific code for defining the new operation. Takes care
+  of the input and output coordinate conversions (between input_crs and triangulation_source_crs,
+  and triangulation_target_crs and output_crs), when needed.
+- tinshift.hpp: Header-based implementation. This file contains the API.
+- tinshift_exceptions.hpp: Exceptions that can be raised during file parsing
+- tinshift_impl.hpp: Implementation of file loading, triangle search and interpolation.
+
+This is the approach that has been followed for the deformation model implementation,
+and which makes it easier to do unit test.
+
+src/quadtree.hpp contains a quadtree implementation.
+
+Performance indications
++++++++++++++++++++++++
+
+Tested on Intel(R) Core(TM) i7-6700HQ CPU @ 2.60GHz, transforming 4 million points
+
+For the KKJ to ETRS89 transformation (1449 triangles and 767 vertices),
+4.4 million points / sec can be transformed.
+
+For comparison, the Helmert-based KKJ to ETRS89 transformation operates at
+1.6 million points / sec.
+
+A triangulation with about 115 000 triangles and 71 000 vertices
+operates at 2.2 million points / sec
+(throughput on more points would be better since the initial loading of the
+triangulation is non-negligible here)
+
+Backward compatibility
+-------------------------------------------------------------------------------
+
+New functionality fully backward compatible.
+
+Testing
+-------------------------------------------------------------------------------
+
+The PROJ test suite will be enhanced to test the new transformation, with a
+new .gie file, and a C++ unit test to test at a lower level.
+
+Documentation
+-------------------------------------------------------------------------------
+
+- The tinshift method will be documented.
+- The JSON format will be documented under https://proj.org/specifications/
+- A JSON schema will also be provided.
+
+Proposed implementation
+-------------------------------------------------------------------------------
+
+An initial implementation is available at https://github.com/rouault/PROJ/tree/tinshift
+
+References
+-------------------------------------------------------------------------------
+
+`Finnish coordinate transformation (automated translation to English) <https://translate.google.fr/translate?sl=auto&tl=en&u=https%3A%2F%2Fwww.maanmittauslaitos.fi%2Fkartat-ja-paikkatieto%2Fasiantuntevalle-kayttajalle%2Fkoordinaattimuunnokset>`_
+
+Adoption status
+-------------------------------------------------------------------------------
+
+The RFC was adopted on 2020-09-02 with +1's from the following PSC members
+
+* Kristian Evers
+* Charles Karney
+* Thomas Knudsen
+* Even Rouault
+
+Funding
+-------------------------------------------------------------------------------
+
+This work is funded by `National Land Survey of Finland <https://www.maanmittauslaitos.fi/en>`_