Merge pull request #1756 from rouault/document_create_operations

Doc: add a operations_computation.rst describe the magic behind createOperations()

2 files changed, 749 insertions, 1 deletions

diff --git a/docs/source/operations/index.rst b/docs/source/operations/index.rst
index 87bfdc40..d4582049 100644
--- a/docs/source/operations/index.rst
+++ b/docs/source/operations/index.rst
@@ -19,4 +19,4 @@ exert a change in reference frame are called transformations.
    conversions/index
    transformations/index
    pipeline
-
+   operations_computation
diff --git a/docs/source/operations/operations_computation.rst b/docs/source/operations/operations_computation.rst
new file mode 100644
index 00000000..930044e3
--- /dev/null
+++ b/docs/source/operations/operations_computation.rst
@@ -0,0 +1,748 @@
+.. _operations_computation:
+
+================================================================================
+Computation of coordinate operations between two CRS
+================================================================================
+
+:Author: Even Rouault
+:Last Updated: 2019-11-28
+
+Introduction
+------------
+
+When using :command:`projinfo -s {crs_def} -t {crs_def}`,
+:command:`cs2cs {crs_def} {crs_def}` or the underlying
+:c:func:`proj_create_crs_to_crs` or :cpp:func:`proj_create_operations` functions,
+PROJ applies an algorithm to compute one or several candidate coordinate operations,
+that can be expressed as a PROJ :ref:`pipeline <pipeline>` to transform between the source
+and the target CRS. This document is about the description of this algorithm that
+finds the actual operations to apply to be able later to perform transform coordinates.
+So this is mostly about metadata management around coordinate operation methods,
+and not about the actual mathematics used to implement those methods.
+As a matter of fact with PROJ 6, there are about 60 000
+lines of code dealing with "metadata" management (including conversions between PROJ
+strings, all CRS WKT variants), to be compared to 30 000 for the purely computation part.
+
+This document is meant as a plain text explanation of the code for developers,
+but also as a in-depth examination of what happens under the hood for curious PROJ
+users. It is important to keep in mind that it is not meant to be the ultimate
+source of truth of how coordinate operations should be computed. There are clearly
+implementation choices and compromises that can be questionned.
+
+Let us start with an example to research operations between the NAD27 and NAD83
+geographic CRS:
+
+.. code-block:: shell
+
+    $ projinfo -s NAD27 -t NAD83 --summary --spatial-test intersects --grid-check none
+
+    Candidate operations found: 10
+    DERIVED_FROM(EPSG):1312, NAD27 to NAD83 (3), 1.0 m, Canada
+    DERIVED_FROM(EPSG):1313, NAD27 to NAD83 (4), 1.5 m, Canada - NAD27, at least one grid missing
+    DERIVED_FROM(EPSG):1241, NAD27 to NAD83 (1), 0.15 m, USA - CONUS including EEZ
+    DERIVED_FROM(EPSG):1243, NAD27 to NAD83 (2), 0.5 m, USA - Alaska including EEZ
+    EPSG:1462, NAD27 to NAD83 (5), 1.0 m, Canada - Quebec, at least one grid missing
+    EPSG:1573, NAD27 to NAD83 (6), 1.5 m, Canada - Quebec, at least one grid missing
+    EPSG:9111, NAD27 to NAD83 (9), 1.5 m, Canada - Saskatchewan, at least one grid missing
+    unknown id, Ballpark geographic offset from NAD27 to NAD83, unknown accuracy, World, has ballpark transformation
+    EPSG:8555, NAD27 to NAD83 (7), 0.15 m, USA - CONUS and GoM, at least one grid missing
+    EPSG:8549, NAD27 to NAD83 (8), 0.5 m, USA - Alaska, at least one grid missing
+
+The algorithm involves many cases, so we will progress in the explanation from
+the most simple case to more complex ones. We document here the working of this
+algorithm as implemented in PROJ 6.3.0.
+The results of some examples might also be quite sensitive to the content of the
+PROJ database and the PROJ version used.
+
+From a code point of view, the entry point of the algorithm is the C++
+:cpp:func:`osgeo::proj::operation::CoordinateOperation::createOperations` method.
+
+It combines several strategies:
+
+- look up in the PROJ database for available operations
+- consider the pair (source CRS, target CRS) to synthetize operations depending
+  on the nature of the source and target CRS.
+
+Geographic CRS to Geographic CRS, with known identifiers
+--------------------------------------------------------
+
+With the above example of two geographic CRS, that have an identified identifier,
+(:program:`projinfo` internally resolves NAD27 to EPSG:4267 and NAD83 to EPSG:4269)
+the algorithm will first search
+in the coordinate operation related tables of the ``proj.db`` if there are records
+that list direct transformations between the source and the target CRS. The
+transformations typically involve :ref:`Helmert <helmert>`-style operations or datum shift based on
+grids (more esoteric operations are possible).
+
+A request similar to the following will be emitted:
+
+.. code-block:: shell
+
+    $ sqlite3 proj.db "SELECT auth_name, code, name, method_name, accuracy FROM \
+                       coordinate_operation_view WHERE \
+                       source_crs_auth_name = 'EPSG' AND \
+                       source_crs_code = '4267' AND \
+                       target_crs_auth_name = 'EPSG' AND \
+                       target_crs_code = '4269'"
+
+    EPSG|1241|NAD27 to NAD83 (1)|NADCON|0.15
+    EPSG|1243|NAD27 to NAD83 (2)|NADCON|0.5
+    EPSG|1312|NAD27 to NAD83 (3)|NTv1|1.0
+    EPSG|1313|NAD27 to NAD83 (4)|NTv2|1.5
+    EPSG|1462|NAD27 to NAD83 (5)|NTv1|1.0
+    EPSG|1573|NAD27 to NAD83 (6)|NTv2|1.5
+    EPSG|8549|NAD27 to NAD83 (8)|NADCON5 (2D)|0.5
+    EPSG|8555|NAD27 to NAD83 (7)|NADCON5 (2D)|0.15
+    EPSG|9111|NAD27 to NAD83 (9)|NTv2|1.5
+    ESRI|108003|NAD_1927_To_NAD_1983_PR_VI|NTv2|0.05
+
+As we have found direct transformations, we will not attempt any more complicated
+research.
+One can note in the above result set that a ESRI:108003 operation was found,
+but as the source and target CRS are in the EPSG registry, and there are
+operations between those CRS in the EPSG registry itself, transformations from
+other authorities will be ignored (except if they are in the PROJ authority,
+which can be used as an override).
+
+As those results all involve operations that does not have a perfect accuracy and that
+does not cover the area of use of the 2 CRSs, a
+'Ballpark geographic offset from NAD27 to NAD83' operation is synthetized by PROJ.
+This operation is a sort of dummy operation that only takes into account potential
+difference of axis orders (long-lat vs lat-long), units (degree vs grads) and
+prime meridian (Greewich vs Paris/Rome/other historic prime meridians). It does
+not attempt any datum shift, hence the "ballpark" qualifier in its name.
+
+Filtering and sorting of coordinate operations
+----------------------------------------------
+
+The last step is to filter and sort results in order of relevance.
+
+The filtering takes into account the following criteria to decide which operations
+must be retained or discarded:
+
+* a minimum accuracy that the user might have expressed,
+* an area of use on which the coordinate operation(s) must apply
+* if the absence of grids needed by an operation must result in discarding it.
+
+The sorting algorithm determines the order of relevance of the operations we got.
+A comparison function compares pair of operations to determine which of the
+two is the most releavant. This is implemented by the :cpp:func:`operator ()`
+method of the SortFunction structure.
+When comparing two operations, the following criteria are used. The tests are
+performed in the order they are listed below:
+
+1. consider as more relevant an operation that can be expressed as a PROJ operation string
+   (the database might list operations whose method is not (yet) implemented by PROJ)
+2. if both operations evaluate identically with respect to the above criterion,
+   consider as more relevant an operation that does not include a synthetic
+   ballpark vertical transformation (occurs when there is a geoid model).
+3. if both operations evaluate identically with respect to the above criterion,
+   consider as more relevant an operation that does not include a synthetic
+   ballpark horizontal tranformation.
+4. consider as more relevant an operation that refers to shift grids that are locally available.
+5. consider as more relevant an operation that refers to grids that are available
+   in one of the proj-datumgrid packages, but not necessarily locally available
+6. consider as more relevant an operation that has a known accuracy.
+7. if two operations have unknown accuracy, consider as more relevant an operation
+   that uses grid(s) if the other one does not (grid based operations are assumed
+   to be more precise than operations relying on a few parameters)
+8. consider as more relevant an operation whose area of use is larger
+   (note: the computation of the are of use is approximate, based on a bounding box)
+9. consider as more relevant an operation that has a better accuracy.
+10. in case of same accuracy, consider as more relevant an operation that does
+    not use grids (operations that use only parameters will be faster)
+11. consider as more relevant an operation that involves less transformation steps
+12. and for completness, if two operations are comparable given all the above criteria,
+    consider as more relevant the one which has the shorter name, and if they
+    have the same length, consider as more relevant the one whose name comes first in
+    lexicographic order (obviously completely arbitrary, but a sorting
+    algorithm must be able to compare all entries)
+
+Geodetic/geographic CRS to Geodetic/geographic CRS, without known identifiers
+-----------------------------------------------------------------------------
+
+In a number of situations, the source and/or target CRS do not have an identifier
+(WKT without identifier, PROJ string, ..)
+The first step is to try to find in the ``proj.db`` a CRS of the same nature of
+the CRS to identify and whose name exactly matches the one provided to the
+:c:func:`createOperations` method. If there is exactly one match and that the CRS are
+"computationnaly" equivalent, then use the code of the CRS for further computations.
+
+If this search did not succeed, or if the previous case with known CRS identifiers
+did not result in matches in the database, the search will be based on the
+datums. That is, a list of geographic CRS whose datum matches the datum of the
+source and target CRS is searched for in the database (by querying the `geodetic_crs`
+database table). If the datum has a known
+identifier, we will use it, otherwise we will look for an equivalent datum in the
+database based on the datum name.
+
+Let's consider the case where the datum of the source CRS is EPSG:6171 "Reseau
+Geodesique Francais 1993" and the datum of the target CRS is EPSG:6258 "European
+Terrestrial Reference System 1989".
+For EPSG:6171, there are 10 matching (non-deprecated) geodetic CRSs: 
+
+- EPSG:4171, RGF93, geographic 2D
+- EPSG:4964, RGF93, geocentric
+- EPSG:4965, RGF93, geographic 3D
+- EPSG:7042, RGF93 (lon-lat), geographic 3D
+- EPSG:7084, RGF93 (lon-lat), geographic 2D
+- IGNF:RGF93, RGF93 cartesiennes geocentriques, geocentric
+- IGNF:RGF93GDD, RGF93 geographiques (dd),geographic 2D
+- IGNF:RGF93GEODD, RGF93 geographiques (dd), geographic 3D
+- IGNF:RGF93G, RGF93 geographiques (dms), geographic 2D
+- IGNF:RGF93GEO, RGF93 geographiques (dms), geographic 3D
+
+The first three entries from the EPSG dataset are typical: for each datum,
+one can define a geographic 2D CRS (latitude, longitude), a geographic 3D crs
+(latitude, longitude, ellipsoidal height) and a geocentric one. For that particular
+case, the EPSG dataset has also included two extra definitions corresponding to a
+longitude, latitude, [ellipsoidal height] coordinate system, as found in the official
+French IGNF registry. This IGNF registry has also definitions for a geographic 2D
+CRS (with an extra subtelty with an entry using decimal degree as unit and another
+one degree-minute-second), geographic 3D and geocentric.
+
+For EPSG:6258, there are 7 matching (non-deprecated) geodetic CRSs: 
+
+- EPSG:4258, ETRS89, geographic 2D
+- EPSG:4936, ETRS89, geocentric
+- EPSG:4937, ETRS89, geographic 3D
+- IGNF:ETRS89, ETRS89 cartesiennes geocentriques, geocentric
+- IGNF:ETRS89G, ETRS89 geographiques (dms), geographic 2D
+- IGNF:ETRS89GEO, ETRS89 geographiques (dms), geographic 3D
+- ESRI:104129, GCS_EUREF_FIN, geographic 2D
+
+So the 3 typical EPSG entries, 3 equivalent (with long, lat ordering for the
+geographic CRS) and one from the ESRI registry;
+
+PROJ can now test 10 x 7 different combinations of source x target CRSs, using
+the database searching method explained in the previous section. As soon as
+one of this combination returns at least one non-ballpark combination, the result
+set coming from that combination is used. PROJ will then add before that
+transformation a conversion between the source CRS and the first intermediate CRS,
+and will add at the end a conversion between the second intermediate CRS and the
+target CRS. Those conversions are conversion between geographic 2D and geographic 3D
+CRS or geographic 2D/3D and geocentric CRS.
+
+This is done by the :c:func:`createOperationsWithDatumPivot()` method.
+
+So if transforming between EPSG:7042, RGF93 (lon-lat), geographic 3D and
+EPSG:4936, ETRS89, geocentric, one get the following concatenated operation,
+chaining an axis order change, the null geocentric translation between
+RGF93 and ETRS89 (EPSG:1591), and a conversion between geographic and geocentric
+coordinates. This concatenated operation is assumed to have a perfect accuracy
+as both the initial and final operations are conversions, and the middle transformation
+accounts for the fact that the RGF93 datum is one realization of ETRS89, so they
+are equivalent for most purposes.
+
+.. code-block:: shell
+
+    $ projinfo projinfo -s EPSG:7042 -t EPSG:4936
+
+    Candidate operations found: 1
+    -------------------------------------
+    Operation n°1:
+
+    unknown id, axis order change (geographic3D horizontal) + RGF93 to ETRS89 (1) + Conversion from ETRS89 (geog2D) to ETRS89 (geocentric), 0 m, France
+
+    PROJ string:
+    +proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=cart +ellps=GRS80
+
+    WKT2:2019 string:
+    CONCATENATEDOPERATION["axis order change (geographic3D horizontal) + RGF93 to ETRS89 (1) + Conversion from ETRS89 (geog2D) to ETRS89 (geocentric)",
+        SOURCECRS[
+            GEOGCRS["RGF93 (lon-lat)",
+                [...]
+                ID["EPSG",7042]]],
+        TARGETCRS[
+            GEODCRS["ETRS89",
+                [...]
+                ID["EPSG",4936]]],
+        STEP[
+            CONVERSION["axis order change (geographic3D horizontal)",
+                METHOD["Axis Order Reversal (Geographic3D horizontal)",
+                    ID["EPSG",9844]],
+                ID["EPSG",15499]]],
+        STEP[
+            COORDINATEOPERATION["RGF93 to ETRS89 (1)",
+                [...]
+                METHOD["Geocentric translations (geog2D domain)",
+                    ID["EPSG",9603]],
+                PARAMETER["X-axis translation",0,
+                    LENGTHUNIT["metre",1],
+                    ID["EPSG",8605]],
+                PARAMETER["Y-axis translation",0,
+                    LENGTHUNIT["metre",1],
+                    ID["EPSG",8606]],
+                PARAMETER["Z-axis translation",0,
+                    LENGTHUNIT["metre",1],
+                    ID["EPSG",8607]],
+                OPERATIONACCURACY[0.0],
+                ID["EPSG",1591],
+                REMARK["May be taken as approximate transformation RGF93 to WGS 84 - see code 1671."]]],
+        STEP[
+            CONVERSION["Conversion from ETRS89 (geog2D) to ETRS89 (geocentric)",
+                METHOD["Geographic/geocentric conversions",
+                    ID["EPSG",9602]]]],
+        USAGE[
+            SCOPE["unknown"],
+            AREA["France"],
+            BBOX[41.15,-9.86,51.56,10.38]]]
+
+Geodetic/geographic CRS to Geodetic/geographic CRS, without direct transformation
+---------------------------------------------------------------------------------
+
+Still considering transformations between geodetic/geographic CRS, but let's
+consider that the lookup in the database for a transformation between
+the source and target CRS (possibly going through the "equivalent" CRS based on
+the same datum as detailed in the previous section) leads to an empty set.
+
+Of course, as most operations are invertible, one first tries to do a lookup switching
+the source and target CRS, and inverting the resulting operation(s):
+
+.. code-block:: shell
+
+    $ projinfo -s NAD83 -t NAD27 --spatial-test intersects --summary
+
+    Candidate operations found: 10
+    INVERSE(DERIVED_FROM(EPSG)):1312, Inverse of NAD27 to NAD83 (3), 1.0 m, Canada
+    INVERSE(DERIVED_FROM(EPSG)):1241, Inverse of NAD27 to NAD83 (1), 0.15 m, USA - CONUS including EEZ
+    INVERSE(DERIVED_FROM(EPSG)):1243, Inverse of NAD27 to NAD83 (2), 0.5 m, USA - Alaska including EEZ
+    INVERSE(DERIVED_FROM(EPSG)):1313, Inverse of NAD27 to NAD83 (4), 1.5 m, Canada - NAD27, at least one grid missing
+    INVERSE(EPSG):1462, Inverse of NAD27 to NAD83 (5), 1.0 m, Canada - Quebec, at least one grid missing
+    INVERSE(EPSG):1573, Inverse of NAD27 to NAD83 (6), 1.5 m, Canada - Quebec, at least one grid missing
+    INVERSE(EPSG):9111, Inverse of NAD27 to NAD83 (9), 1.5 m, Canada - Saskatchewan, at least one grid missing
+    unknown id, Ballpark geographic offset from NAD83 to NAD27, unknown accuracy, World, has ballpark transformation
+    INVERSE(EPSG):8555, Inverse of NAD27 to NAD83 (7), 0.15 m, USA - CONUS and GoM, at least one grid missing
+    INVERSE(EPSG):8549, Inverse of NAD27 to NAD83 (8), 0.5 m, USA - Alaska, at least one grid missing
+
+That was an easy case. Now let's consider the transformation between the Australian
+CRS AGD84 and GDA2020. There is no direct transformation from AGD84 to GDA2020, or
+in the reverse direction, even when considering alternative geodetic CRS based on
+the underlying datums. PROJ will then do a cross-join in the coordinate_operation_view
+table to find the tuples (op1, op2) of coordinate operations such that:
+
+- SOURCE_CRS = op1.source_crs AND op1.target_crs = op2.source_crs AND op2.target_crs = TARGET_CRS or
+- SOURCE_CRS = op1.source_crs AND op1.target_crs = op2.target_crs AND op2.source_crs = TARGET_CRS or
+- SOURCE_CRS = op1.target_crs AND op1.source_crs = op2.source_crs AND op2.target_crs = TARGET_CRS or
+- SOURCE_CRS = op1.target_crs AND op1.source_crs = op2.target_crs AND op2.source_crs = TARGET_CRS
+
+Depending on which case is selected, op1 and op2 should be reversed, before
+being concatenated.
+
+This logic is implement by the ``findsOpsInRegistryWithIntermediate()`` method.
+
+Assuming that the proj-datumgrid-oceania package is installed, we get the
+following results for the AGD84 to GDA2020 coordinate operations lookup:
+
+.. code-block:: shell
+
+    $ projinfo -s AGD84 -t GDA2020 --spatial-test intersects -o PROJ
+
+    Candidate operations found: 4
+    -------------------------------------
+    Operation n°1:
+
+    unknown id, AGD84 to GDA94 (5) + GDA94 to GDA2020 (1), 0.11 m, Australia - AGD84
+
+    PROJ string:
+    +proj=pipeline +step +proj=axisswap +order=2,1 \
+                   +step +proj=unitconvert +xy_in=deg +xy_out=rad \
+                   +step +proj=hgridshift +grids=National_84_02_07_01.gsb \
+                   +step +proj=push +v_3 \
+                   +step +proj=cart +ellps=GRS80 \
+                   +step +proj=helmert +x=0.06155 +y=-0.01087 +z=-0.04019 \
+                                       +rx=-0.0394924 +ry=-0.0327221 +rz=-0.0328979 \
+                                       +s=-0.009994 +convention=coordinate_frame \
+                   +step +inv +proj=cart +ellps=GRS80 \
+                   +step +proj=pop +v_3 \
+                   +step +proj=unitconvert +xy_in=rad +xy_out=deg \
+                   +step +proj=axisswap +order=2,1
+
+    -------------------------------------
+    Operation n°2:
+
+    unknown id, AGD84 to GDA94 (2) + GDA94 to GDA2020 (1), 1.01 m, Australia - AGD84
+
+    PROJ string:
+    +proj=pipeline +step +proj=axisswap +order=2,1 \
+                   +step +proj=unitconvert +xy_in=deg +xy_out=rad \
+                   +step +proj=push +v_3 \
+                   +step +proj=cart +ellps=aust_SA \
+                   +step +proj=helmert +x=-117.763 +y=-51.51 +z=139.061 \
+                                       +rx=-0.292 +ry=-0.443 +rz=-0.277 +s=-0.191 \
+                                       +convention=coordinate_frame \
+                   +step +proj=helmert +x=0.06155 +y=-0.01087 +z=-0.04019 \
+                                       +rx=-0.0394924 +ry=-0.0327221 +rz=-0.0328979 \
+                                       +s=-0.009994 +convention=coordinate_frame \
+                   +step +inv +proj=cart +ellps=GRS80 \
+                   +step +proj=pop +v_3 \
+                   +step +proj=unitconvert +xy_in=rad +xy_out=deg \
+                   +step +proj=axisswap +order=2,1
+
+    -------------------------------------
+    Operation n°3:
+
+    unknown id, AGD84 to GDA94 (5) + GDA94 to GDA2020 (2), 0.15 m, unknown domain of validity
+
+    PROJ string:
+    +proj=pipeline +step +proj=axisswap +order=2,1 \
+                   +step +proj=unitconvert +xy_in=deg +xy_out=rad \
+                   +step +proj=hgridshift +grids=National_84_02_07_01.gsb \
+                   +step +proj=hgridshift +grids=GDA94_GDA2020_conformal_and_distortion.gsb \
+                   +step +proj=unitconvert +xy_in=rad +xy_out=deg \
+                   +step +proj=axisswap +order=2,1
+
+    -------------------------------------
+    Operation n°4:
+
+    unknown id, AGD84 to GDA94 (5) + GDA94 to GDA2020 (3), 0.15 m, unknown domain of validity
+
+    PROJ string:
+    +proj=pipeline +step +proj=axisswap +order=2,1 \
+                   +step +proj=unitconvert +xy_in=deg +xy_out=rad \
+                   +step +proj=hgridshift +grids=National_84_02_07_01.gsb \
+                   +step +proj=hgridshift +grids=GDA94_GDA2020_conformal.gsb \
+                   +step +proj=unitconvert +xy_in=rad +xy_out=deg \
+                   +step +proj=axisswap +order=2,1
+
+One can see that the selected intermediate CRS that has been used is GDA94.
+This is a completely novel behaviour of PROJ 6 as opposed to the logic of PROJ.4 
+where datum transformations implied using EPSG:4326 / WGS 84 has the mandatory
+datum hub. PROJ 6 no longer hardcodes it as the mandatory datum hub, and relies
+on the database to find the appropriate hub(s).
+Actually, WGS 84 has been considered during the above lookup, because there are
+transformations between AGD84 and WGS 84 and WGS 84 and GDA2020. However those
+have been discarded in a step which we did not mention previously: just after
+the initial filtering of results and their sorting, there is a final filtering
+that is done. In the list of sorted results, given two operations A and B that
+have the same area of use, if B has an accuracy lower than A, and A does not use
+grids, or all the needed grids are available, then B is discarded.
+
+If one forces the datum hub to be considered to be EPSG:4326, ones gets:
+
+.. code-block:: shell
+
+    $ projinfo -s AGD84 -t GDA2020 --spatial-test intersects --pivot-crs EPSG:4326 -o PROJ
+
+    Candidate operations found: 2
+    -------------------------------------
+    Operation n°1:
+
+    unknown id, AGD84 to WGS 84 (7) + Inverse of GDA2020 to WGS 84 (2), 4 m, Australia - AGD84
+
+    PROJ string:
+    +proj=pipeline +step +proj=axisswap +order=2,1 \
+                   +step +proj=unitconvert +xy_in=deg +xy_out=rad \
+                   +step +proj=push +v_3 \
+                   +step +proj=cart +ellps=aust_SA \
+                   +step +proj=helmert +x=-117.763 +y=-51.51 +z=139.061 \
+                                       +rx=-0.292 +ry=-0.443 +rz=-0.277 \
+                                       +s=-0.191 +convention=coordinate_frame \
+                   +step +inv +proj=cart +ellps=GRS80 \
+                   +step +proj=pop +v_3 \
+                   +step +proj=unitconvert +xy_in=rad +xy_out=deg \
+                   +step +proj=axisswap +order=2,1
+
+    -------------------------------------
+    Operation n°2:
+
+    unknown id, AGD84 to WGS 84 (9) + Inverse of GDA2020 to WGS 84 (2), 4 m, Australia - AGD84
+
+    PROJ string:
+    +proj=pipeline +step +proj=axisswap +order=2,1 \
+                   +step +proj=unitconvert +xy_in=deg +xy_out=rad \
+                   +step +proj=hgridshift +grids=National_84_02_07_01.gsb \
+                   +step +proj=unitconvert +xy_in=rad +xy_out=deg \
+                   +step +proj=axisswap +order=2,1
+
+Those operations are less accurate, since WGS 84 is assumed to be equivalent to
+GDA2020 with an accuracy of 4 metre. This is an instance demonstrating that using
+WGS 84 as a hub systematically can be sub-optimal.
+
+There are still situations where the attempt to find a hub CRS does not work,
+because there is no such hub. This can occur for example when transforming from
+GDA94 to the latest realization at time of writing of WGS 84, WGS 84 (G1762).
+There are transformations between WGS 84 (G1762). Using the above described
+techniques, we would only find one non-ballpark operation taking the route:
+1. Conversion from GDA94 (geog2D) to GDA94 (geocentric): synthetized by PROJ
+2. Inverse of ITRF2008 to GDA94 (1): from EPSG
+3. Inverse of WGS 84 (G1762) to ITRF2008 (1): from EPSG
+4. Conversion from WGS 84 (G1762) (geocentric) to WGS 84 (G1762): synthetized by PROJ
+
+This is not bad, but the global validity area of use is "Australia - onshore and EEZ",
+whereas GDA94 has a larger area of use.
+There is another road that can be taken by going throug GDA2020 instead of ITRF2008.
+The GDA94 to GDA2020 transformations operate on the respective geographic CRS,
+whereas GDA2020 to WGS 84 (G1762) operate on the geocentric CRS. Consequently,
+GDA2020 cannot be identifier as a hub by a "simple" self-join SQL request on
+the coordinate operation table. This requires to do the join based on the datum
+referenced by the source and target CRS of each operation rather than the
+source and target CRS themselves. When there is a match, PROJ inserts the required
+conversions between geographic and geocentric CRS to have a consistent concatenated
+operation, like the following:
+1. GDA94 to GDA2020 (1): from EPSG
+2. Conversion from GDA2020 (geog2D) to GDA2020 (geocentric): synthetized by PROJ
+3. GDA2020 to WGS 84 (G1762) (1): frmo EPSG
+4. Conversion from WGS 84 (G1762) (geocentric) to WGS 84 (G1762) (geog2D): synthetized by PROJ
+
+Projected CRS to any target CRS
+---------------------------------------------------------------------------------
+
+This actually extends to any Derived CRS, whose Projected CRS is a well-known
+particular case. Such transformations are done in 2 steps:
+
+1. Use the inverse conversion of the derived CRS to its base CRS, typically an
+   inverse map projection.
+2. Find transformations from this base CRS to the target CRS. If the base CRS
+   is the target CRS, this step can be skipped.
+
+.. code-block:: shell
+
+    $ projinfo -s EPSG:32631 -t RGF93
+
+    Candidate operations found: 1
+    -------------------------------------
+    Operation n°1:
+
+    unknown id, Inverse of UTM zone 31N + Inverse of RGF93 to WGS 84 (1), 1 m, France
+
+    PROJ string:
+    +proj=pipeline +step +inv +proj=utm +zone=31 +ellps=WGS84 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1
+
+This is implemented by the ``createOperationsDerivedTo`` method
+
+For the symetric case, source CRS to a derived CRS, the above algorithm is applied
+by switching the source and target CRS, and then inverting the resulting operation(s).
+This is mostly a matter of avoiding to write very similar code twice. This logic
+is also applied to all below cases when considering the transformation between 2 different
+types of objects.
+
+.. _verttogeog:
+
+Vertical CRS to a Geographic CRS
+---------------------------------------------------------------------------------
+
+Such transformation is normally not meant as being used as standalone by PROJ
+users, but as an internal computation step of a Compound CRS to a target CRS.
+
+In cases where we are lucky, the PROJ database will have a transformation registered
+between those:
+
+.. code-block:: shell
+
+    $ projinfo -s "NAVD88 height" -t "NAD83(2011)" -o PROJ --spatial-test intersects
+    Candidate operations found: 11
+    -------------------------------------
+    Operation n°1:
+
+    INVERSE(DERIVED_FROM(EPSG)):9229, Inverse of NAD83(2011) to NAVD88 height (3), 0.015 m, USA - CONUS - onshore
+
+    PROJ string:
+    +proj=vgridshift +grids=g2018u0.gtx +multiplier=1
+
+
+But in cases where there is no match, the ``createOperationsVertToGeog`` method
+will be used to synthetize a ballpark vertical transformation, just taking care
+of unit changes, and axis reversal in case the vertical CRS was a depth rather than
+a height. Of course the results of such an operation are questionable, hence the
+ballpark qualifier and a unknown accuracy advertized for such an operation.
+
+Vertical CRS to a Vertical CRS
+---------------------------------------------------------------------------------
+
+Overall logic is similar to the above case. There might be direct operations in
+the PROJ database, involving grid transformations or simple offsets. The fallback
+case is to synthetize a ballpark transformation.
+
+This is implemented by the ``createOperationsVertToVert`` method
+
+.. code-block:: shell
+
+    $ projinfo -s "NGVD29 depth (ftUS)" -t "NAVD88 height" --spatial-test intersects -o PROJ
+
+    Candidate operations found: 3
+    -------------------------------------
+    Operation n°1:
+
+    unknown id, Inverse of NGVD29 height (ftUS) to NGVD29 depth (ftUS) + NGVD29 height (ftUS) to NGVD29 height (m) + NGVD29 height (m) to NAVD88 height (3), 0.02 m, USA - CONUS east of 89°W - onshore
+
+    PROJ string:
+    +proj=pipeline +step +proj=axisswap +order=1,2,-3 +step +proj=unitconvert +z_in=us-ft +z_out=m +step +proj=vgridshift +grids=vertcone.gtx +multiplier=0.001
+
+    -------------------------------------
+    Operation n°2:
+
+    unknown id, Inverse of NGVD29 height (ftUS) to NGVD29 depth (ftUS) + NGVD29 height (ftUS) to NGVD29 height (m) + NGVD29 height (m) to NAVD88 height (2), 0.02 m, USA - CONUS 89°W-107°W - onshore
+
+    PROJ string:
+    +proj=pipeline +step +proj=axisswap +order=1,2,-3 +step +proj=unitconvert +z_in=us-ft +z_out=m +step +proj=vgridshift +grids=vertconc.gtx +multiplier=0.001
+
+    -------------------------------------
+    Operation n°3:
+
+    unknown id, Inverse of NGVD29 height (ftUS) to NGVD29 depth (ftUS) + NGVD29 height (ftUS) to NGVD29 height (m) + NGVD29 height (m) to NAVD88 height (1), 0.02 m, USA - CONUS west of 107°W - onshore
+
+    PROJ string:
+    +proj=pipeline +step +proj=axisswap +order=1,2,-3 +step +proj=unitconvert +z_in=us-ft +z_out=m +step +proj=vgridshift +grids=vertconw.gtx +multiplier=0.001
+
+
+Compound CRS to a Geographic CRS
+---------------------------------------------------------------------------------
+
+A typical example of a Compound CRS is a CRS made of a geographic or projected CRS
+as the horizontal component, and a vertical CRS. E.g. "NAD83 + NAVD88 height"
+
+When the horizontal component of the compound source CRS is a projected CRS, we
+first look for the operation from this source CRS to another compound CRS made
+of the geographic CRS base of the projected CRS,
+like "NAD83 / California zone 1 (ftUS) + NAVD88 height" to "NAD83 + NAVD88 height",
+which ultimately goes to one of the above described case. Then we can consider
+the transformation from a compound CRS made of a geographic CRS to another geographic CRS.
+
+It first starts by the vertical transformations from the vertical CRS of the
+source compound CRS to the target geographic CRS, using the strategy detailed
+in :ref:`Vertical CRS to a Geographic CRS <verttogeog>`
+
+What we did not mention is that when there is not a transformation registered
+between the vertical CRS and the target geographic CRS, PROJ attempts to find
+transformations between that vertical CRS and any other geographic CRS. This is
+clearly an approximation.
+If the research of the vertical CRS to the target geographic CRS resulted in
+operations that use grids that are not available, as another approximation, we
+research operations from the vertical CRS to the source geographic CRS for the
+vertical component.
+
+Once we got those more or less accurate vertical transformations, we must consider
+the horizontal transformation(s). The algorithm iterates over all found vertical
+transformations and look for their target geographic CRS. This will be used as
+the interpolation CRS for horizontal transformations. PROJ will then look for
+available transformations from the source geographic CRS to the interpolation CRS
+and from the interpolation CRS to the target geographic CRS. There is then a
+3-level loop to create the final set of operations chaining together:
+
+- the horizontal transformation from the source geographic CRS to the interpolation CRS
+- the vertical transformation from the source vertical CRS to the interpolation CRS
+- the horizontal transformation from the interpolation CRS to the target geographic CRS.
+
+This is implemented by the ``createOperationsCompoundToGeog`` method
+
+Example:
+
+.. code-block:: shell
+
+    $ projinfo -s "NAD83(NSRS2007) + NAVD88 height" -t "WGS 84 (G1762)" --spatial-test intersects --summary
+
+    Candidate operations found: 21
+    unknown id, Inverse of NAD83(NSRS2007) to NAVD88 height (1) + NAD83(NSRS2007) to WGS 84 (1) + WGS 84 to WGS 84 (G1762), 3.05 m, USA - CONUS - onshore
+    unknown id, Inverse of NAD83(HARN) to NAD83(NSRS2007) (1) + Inverse of NAD83(HARN) to NAVD88 height (7) + NAD83(HARN) to WGS 84 (1) + WGS 84 to WGS 84 (G1762), 3.15 m, USA - CONUS south of 41°N, 95°W to 78°W - onshore
+    unknown id, Inverse of NAD83(HARN) to NAD83(NSRS2007) (1) + Inverse of NAD83(HARN) to NAVD88 height (7) + NAD83(HARN) to WGS 84 (3) + WGS 84 to WGS 84 (G1762), 3.15 m, USA - CONUS south of 41°N, 95°W to 78°W - onshore
+    unknown id, Inverse of NAD83(HARN) to NAD83(NSRS2007) (1) + Inverse of NAD83(HARN) to NAVD88 height (6) + NAD83(HARN) to WGS 84 (1) + WGS 84 to WGS 84 (G1762), 3.15 m, USA - CONUS south of 41°N, 112°W to 95°W - onshore
+    unknown id, Inverse of NAD83(HARN) to NAD83(NSRS2007) (1) + Inverse of NAD83(HARN) to NAVD88 height (6) + NAD83(HARN) to WGS 84 (3) + WGS 84 to WGS 84 (G1762), 3.15 m, USA - CONUS south of 41°N, 112°W to 95°W - onshore
+    unknown id, Inverse of NAD83(HARN) to NAD83(NSRS2007) (1) + Inverse of NAD83(HARN) to NAVD88 height (2) + NAD83(HARN) to WGS 84 (1) + WGS 84 to WGS 84 (G1762), 3.15 m, USA - CONUS north of 41°N, 112°W to 95°W
+    unknown id, Inverse of NAD83(HARN) to NAD83(NSRS2007) (1) + Inverse of NAD83(HARN) to NAVD88 height (2) + NAD83(HARN) to WGS 84 (3) + WGS 84 to WGS 84 (G1762), 3.15 m, USA - CONUS north of 41°N, 112°W to 95°W
+    unknown id, Inverse of NAD83(HARN) to NAD83(NSRS2007) (1) + Inverse of NAD83(HARN) to NAVD88 height (3) + NAD83(HARN) to WGS 84 (1) + WGS 84 to WGS 84 (G1762), 3.15 m, USA - CONUS north of 41°N, 95°W to 78°W
+    unknown id, Inverse of NAD83(HARN) to NAD83(NSRS2007) (1) + Inverse of NAD83(HARN) to NAVD88 height (3) + NAD83(HARN) to WGS 84 (3) + WGS 84 to WGS 84 (G1762), 3.15 m, USA - CONUS north of 41°N, 95°W to 78°W
+    unknown id, Inverse of NAD83(HARN) to NAD83(NSRS2007) (1) + Inverse of NAD83(HARN) to NAVD88 height (5) + NAD83(HARN) to WGS 84 (1) + WGS 84 to WGS 84 (G1762), 3.15 m, USA - CONUS south of 41°N, west of 112°W - onshore
+    unknown id, Inverse of NAD83(HARN) to NAD83(NSRS2007) (1) + Inverse of NAD83(HARN) to NAVD88 height (5) + NAD83(HARN) to WGS 84 (3) + WGS 84 to WGS 84 (G1762), 3.15 m, USA - CONUS south of 41°N, west of 112°W - onshore
+    unknown id, Inverse of NAD83(HARN) to NAD83(NSRS2007) (1) + Inverse of NAD83(HARN) to NAVD88 height (1) + NAD83(HARN) to WGS 84 (1) + WGS 84 to WGS 84 (G1762), 3.15 m, USA - CONUS north of 41°N, west of 112°W - onshore
+    unknown id, Inverse of NAD83(HARN) to NAD83(NSRS2007) (1) + Inverse of NAD83(HARN) to NAVD88 height (1) + NAD83(HARN) to WGS 84 (3) + WGS 84 to WGS 84 (G1762), 3.15 m, USA - CONUS north of 41°N, west of 112°W - onshore
+    unknown id, Inverse of NAD83(HARN) to NAD83(NSRS2007) (1) + Inverse of NAD83(HARN) to NAVD88 height (4) + NAD83(HARN) to WGS 84 (1) + WGS 84 to WGS 84 (G1762), 3.15 m, USA - CONUS north of 41°N, east of 78°W - onshore
+    unknown id, Inverse of NAD83(HARN) to NAD83(NSRS2007) (1) + Inverse of NAD83(HARN) to NAVD88 height (4) + NAD83(HARN) to WGS 84 (3) + WGS 84 to WGS 84 (G1762), 3.15 m, USA - CONUS north of 41°N, east of 78°W - onshore
+    unknown id, Inverse of NAD83(HARN) to NAD83(NSRS2007) (1) + Inverse of NAD83(HARN) to NAVD88 height (8) + NAD83(HARN) to WGS 84 (1) + WGS 84 to WGS 84 (G1762), 3.15 m, USA - CONUS south of 41°N, east of 78°W - onshore
+    unknown id, Inverse of NAD83(HARN) to NAD83(NSRS2007) (1) + Inverse of NAD83(HARN) to NAVD88 height (8) + NAD83(HARN) to WGS 84 (3) + WGS 84 to WGS 84 (G1762), 3.15 m, USA - CONUS south of 41°N, east of 78°W - onshore
+    unknown id, Ballpark geographic offset from NAD83(NSRS2007) to NAD83(FBN) + Inverse of NAD83(FBN) to NAVD88 height (1) + Ballpark geographic offset from NAD83(FBN) to WGS 84 (G1762), unknown accuracy, USA - CONUS - onshore, has ballpark transformation
+    unknown id, Ballpark geographic offset from NAD83(NSRS2007) to NAD83(2011) + Inverse of NAD83(2011) to NAVD88 height (3) + Ballpark geographic offset from NAD83(2011) to WGS 84 (G1762), unknown accuracy, USA - CONUS - onshore, has ballpark transformation
+    unknown id, Ballpark geographic offset from NAD83(NSRS2007) to NAD83(2011) + Inverse of NAD83(2011) to NAVD88 height (3) + Conversion from NAD83(2011) (geog2D) to NAD83(2011) (geocentric) + Inverse of ITRF2008 to NAD83(2011) (1) + Inverse of WGS 84 (G1762) to ITRF2008 (1) + Conversion from WGS 84 (G1762) (geocentric) to WGS 84 (G1762) (geog2D), unknown accuracy, USA - CONUS - onshore, has ballpark transformation
+    unknown id, NAD83(NSRS2007) to WGS 84 (1) + WGS 84 to WGS 84 (G1762) + Transformation from NAVD88 height to WGS 84 (G1762) (ballpark vertical transformation, without ellipsoid height to vertical height correction), unknown accuracy, USA - CONUS and Alaska; PRVI, has ballpark transformation
+
+
+CompoundCRS to CompoundCRS
+---------------------------------------------------------------------------------
+
+There is some similarity with the previous paragraph. We first research the
+vertical transformations between the vertical CRS. If such tranformation has
+a registered interpolation geographic CRS, then it is used. Otherwise we fallback
+to the geographic CRS of the source CRS.
+
+Finally, a 3-level loop to create the final set of operations chaining together:
+
+- the horizontal transformation from the source CRS to the interpolation CRS
+- the vertical transformation
+- the horizontal transformation from the interpolation CRS to the target CRS.
+
+This is implemented by the ``createOperationsCompoundToGeog`` method
+
+Example:
+
+.. code-block:: shell
+
+    $ projinfo -s "NAD27 + NGVD29 height (ftUS)" -t "NAD83 + NAVD88 height" --spatial-test intersects --summary
+
+    Candidate operations found: 20
+    unknown id, NGVD29 height (ftUS) to NAVD88 height (3) + NAD27 to NAD83 (1), 0.17 m, USA - CONUS east of 89°W - onshore
+    unknown id, NGVD29 height (ftUS) to NAVD88 height (2) + NAD27 to NAD83 (1), 0.17 m, USA - CONUS 89°W-107°W - onshore
+    unknown id, NGVD29 height (ftUS) to NAVD88 height (1) + NAD27 to NAD83 (1), 0.17 m, USA - CONUS west of 107°W - onshore
+    unknown id, NGVD29 height (ftUS) to NAVD88 height (3) + NAD27 to NAD83 (3), 1.02 m, unknown domain of validity
+    unknown id, NGVD29 height (ftUS) to NAVD88 height (2) + NAD27 to NAD83 (3), 1.02 m, unknown domain of validity
+    unknown id, NGVD29 height (ftUS) to NAVD88 height (1) + NAD27 to NAD83 (3), 1.02 m, unknown domain of validity
+    unknown id, NGVD29 height (ftUS) to NAVD88 height (3) + NAD27 to NAD83 (5), 1.02 m, unknown domain of validity, at least one grid missing
+    unknown id, NGVD29 height (ftUS) to NAVD88 height (3) + NAD27 to NAD83 (6), 1.52 m, unknown domain of validity, at least one grid missing
+    unknown id, NGVD29 height (ftUS) to NAVD88 height (2) + NAD27 to NAD83 (9), 1.52 m, unknown domain of validity, at least one grid missing
+    unknown id, NGVD29 height (ftUS) to NAVD88 height (1) + NAD27 to NAD83 (9), 1.52 m, unknown domain of validity, at least one grid missing
+    unknown id, NGVD29 height (ftUS) to NAVD88 height (3) + Ballpark geographic offset from NAD27 to NAD83, unknown accuracy, USA - CONUS east of 89°W - onshore, has ballpark transformation
+    unknown id, NGVD29 height (ftUS) to NAVD88 height (2) + Ballpark geographic offset from NAD27 to NAD83, unknown accuracy, USA - CONUS 89°W-107°W - onshore, has ballpark transformation
+    unknown id, NGVD29 height (ftUS) to NAVD88 height (1) + Ballpark geographic offset from NAD27 to NAD83, unknown accuracy, USA - CONUS west of 107°W - onshore, has ballpark transformation
+    unknown id, Transformation from NGVD29 height (ftUS) to NAVD88 height (ballpark vertical transformation) + NAD27 to NAD83 (1), unknown accuracy, USA - CONUS including EEZ, has ballpark transformation
+    unknown id, Transformation from NGVD29 height (ftUS) to NAVD88 height (ballpark vertical transformation) + NAD27 to NAD83 (3), unknown accuracy, Canada, has ballpark transformation
+    unknown id, Transformation from NGVD29 height (ftUS) to NAVD88 height (ballpark vertical transformation) + NAD27 to NAD83 (4), unknown accuracy, Canada - NAD27, has ballpark transformation
+    unknown id, Transformation from NGVD29 height (ftUS) to NAVD88 height (ballpark vertical transformation) + NAD27 to NAD83 (5), unknown accuracy, Canada - Quebec, has ballpark transformation, at least one grid missing
+    unknown id, Transformation from NGVD29 height (ftUS) to NAVD88 height (ballpark vertical transformation) + NAD27 to NAD83 (6), unknown accuracy, Canada - Quebec, has ballpark transformation, at least one grid missing
+    unknown id, Transformation from NGVD29 height (ftUS) to NAVD88 height (ballpark vertical transformation) + NAD27 to NAD83 (9), unknown accuracy, Canada - Saskatchewan, has ballpark transformation, at least one grid missing
+    unknown id, Transformation from NGVD29 height (ftUS) to NAVD88 height (ballpark vertical transformation) + Ballpark geographic offset from NAD27 to NAD83, unknown accuracy, World, has ballpark transformation
+
+
+When the source or target CRS is a BoundCRS
+---------------------------------------------------------------------------------
+
+The BoundCRS concept is an hybrid concept where a CRS is linked to a transformation
+from it to a hub CRS, typically WGS 84. This is a long-time practice in PROJ.4
+strings with the ``+towgs84``, ``+nadgrids`` and ``+geoidgrids`` keywords, or the
+``TOWGS84[]`` node of WKT 1. When encountering those attributes when parsing
+a CRS string, PROJ will create a BoundCRS object capturing this transformation.
+A BoundCRS object can also be provided with a WKT2 string, and in that case with
+a hub CRS being potentially different from WGS 84.
+
+Let's consider the case of a transformation between a BoundCRS
+("+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000
++ellps=airy +towgs84=446.448,-125.157,542.06,0.15,0.247,0.842,-20.489 +units=m"
+which used to be the PROJ.4 definition of "OSGB 1936 / British National Grid")
+and a target Geographic CRS, ETRS89.
+
+We apply the following steps:
+
+- transform from the base of the source CRS (that is the CRS wrapped by BoundCRS,
+  here a ProjectedCRS) to the geographic CRS of this base CRS
+- apply the transformation of the BoundCRS to go from the geographic CRS of this base CRS
+  to the hub CRS of the BoundCRS, in that instance WGS 84.
+- apply a transformation from the hub CRS to the target CRS.
+
+This is implemented by the ``createOperationsBoundToGeog`` method
+
+Example:
+
+.. code-block:: shell
+
+    $ projinfo -s "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.06,0.15,0.247,0.842,-20.489 +units=m +type=crs" -t ETRS89 -o PROJ
+
+    Candidate operations found: 1
+    -------------------------------------
+    Operation n°1:
+
+    unknown id, Inverse of unknown + Transformation from unknown to WGS84 + Inverse of ETRS89 to WGS 84 (1), unknown accuracy, Europe - ETRS89
+
+    PROJ string:
+    +proj=pipeline +step +inv +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +step +proj=push +v_3 +step +proj=cart +ellps=airy +step +proj=helmert +x=446.448 +y=-125.157 +z=542.06 +rx=0.15 +ry=0.247 +rz=0.842 +s=-20.489 +convention=position_vector +step +inv +proj=cart +ellps=GRS80 +step +proj=pop +v_3 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1
+
+There are other situations with BoundCRS, involving vertical transformations,
+or transforming to other objects than a geographic CRS, but the curious reader
+will have to inspect the code for the actual gory details.