Vertical part from Compound ignored transforming to DerivedProjected

[jjimenezshaw]

Example of problem

When I try to transform from a Compound CRS like EPSG:6318+6360 (vertical in US foot) to a DerivedProjected 3D CRS, it seems to be ignoring the vertical component of the compound CRS. However it is not the case when source is a geographic 3D (with vertical in US foot), or when destination is a projected CRS in 3D.

Compound to DerivedProjected 3D

$ projinfo -o proj -s EPSG:6318+6360 -t 'DERIVEDPROJCRS["Custom Site Calibrated 3D CRS",
    BASEPROJCRS["Transverse Mercator centered in area of interest",
        BASEGEOGCRS["NAD83(2011)",
            DATUM["NAD83 (National Spatial Reference System 2011)",
                ELLIPSOID["GRS 1980",6378137,298.257222101,
                    LENGTHUNIT["metre",1]]],
            PRIMEM["Greenwich",0,
                ANGLEUNIT["degree",0.0174532925199433]]],
        CONVERSION["Transverse Mercator",
            METHOD["Transverse Mercator",
                ID["EPSG",9807]],
            PARAMETER["Latitude of natural origin",41,
                ANGLEUNIT["degree",0.0174532925199433],
                ID["EPSG",8801]],
            PARAMETER["Longitude of natural origin",-73,
                ANGLEUNIT["degree",0.0174532925199433],
                ID["EPSG",8802]],
            PARAMETER["Scale factor at natural origin",1,
                SCALEUNIT["unity",1],
                ID["EPSG",8805]],
            PARAMETER["False easting",0,
                LENGTHUNIT["US survey foot",0.304800609601219],
                ID["EPSG",8806]],
            PARAMETER["False northing",0,
                LENGTHUNIT["US survey foot",0.304800609601219],
                ID["EPSG",8807]]]],
    DERIVINGCONVERSION["Affine transformation as PROJ-based",
        METHOD["PROJ-based operation method: +proj=pipeline  +step +proj=unitconvert +xy_in=m +xy_out=us-ft +z_in=m +z_out=us-ft +step +proj=affine +xoff=5 +yoff=52 +zoff=10 +step +proj=unitconvert +xy_in=us-ft +xy_out=m +z_in=us-ft +z_out=m"]],
    CS[Cartesian,3],
        AXIS["(Y)",north,
            ORDER[1],
            LENGTHUNIT["US survey foot",0.304800609601219,
                ID["EPSG",9003]]],
        AXIS["(X)",east,
            ORDER[2],
            LENGTHUNIT["US survey foot",0.304800609601219,
                ID["EPSG",9003]]],
        AXIS["(Z)",up,
            ORDER[3],
            LENGTHUNIT["US survey foot",0.304800609601219,
                ID["EPSG",9003]]]]'
Candidate operations found: 1
-------------------------------------
Operation No. 1:

unknown id, Null geographic offset from NAD83(2011) to NAD83(2011) + Transverse Mercator + Affine transformation as PROJ-based, 0 m, World

PROJ string:
+proj=pipeline
  +step +proj=axisswap +order=2,1
  +step +proj=unitconvert +xy_in=deg +xy_out=rad
  +step +proj=tmerc +lat_0=41 +lon_0=-73 +k=1 +x_0=0 +y_0=0 +ellps=GRS80
  +step +proj=unitconvert +xy_in=m +xy_out=us-ft +z_in=m +z_out=us-ft
  +step +proj=affine +xoff=5 +yoff=52 +zoff=10
  +step +proj=unitconvert +xy_in=us-ft +xy_out=m +z_in=us-ft +z_out=m
  +step +proj=unitconvert +xy_in=m +z_in=m +xy_out=us-ft +z_out=us-ft
  +step +proj=axisswap +order=2,1

Geographic 3D (vert in US foot) to DerivedProjected

(see that source vertical units are considered)

$ projinfo -o proj -s 'GEOGCRS["NAD83(2011)",
     DATUM["NAD83 (National Spatial Reference System 2011)",
         ELLIPSOID["GRS 1980",6378137,298.257222101,
             LENGTHUNIT["metre",1]]],
     PRIMEM["Greenwich",0,
         ANGLEUNIT["degree",0.0174532925199433]],
     CS[ellipsoidal,3],
         AXIS["geodetic latitude (Lat)",north,
             ORDER[1],
             ANGLEUNIT["degree",0.0174532925199433]],
         AXIS["geodetic longitude (Lon)",east,
             ORDER[2],
             ANGLEUNIT["degree",0.0174532925199433]],
         AXIS["ellipsoidal height (h)",up,
             ORDER[3],
             LENGTHUNIT["US survey foot",0.304800609601219,
                 ID["EPSG",9003]]]]' -t 'DERIVEDPROJCRS["Custom Site Calibrated 3D CRS",
     BASEPROJCRS["Transverse Mercator centered in area of interest",
         BASEGEOGCRS["NAD83(2011)",
             DATUM["NAD83 (National Spatial Reference System 2011)",
                 ELLIPSOID["GRS 1980",6378137,298.257222101,
                     LENGTHUNIT["metre",1]]],
             PRIMEM["Greenwich",0,
                 ANGLEUNIT["degree",0.0174532925199433]]],
         CONVERSION["Transverse Mercator",
             METHOD["Transverse Mercator",
                 ID["EPSG",9807]],
             PARAMETER["Latitude of natural origin",41,
                 ANGLEUNIT["degree",0.0174532925199433],
                 ID["EPSG",8801]],
             PARAMETER["Longitude of natural origin",-73,
                 ANGLEUNIT["degree",0.0174532925199433],
                 ID["EPSG",8802]],
             PARAMETER["Scale factor at natural origin",1,
                 SCALEUNIT["unity",1],
                 ID["EPSG",8805]],
             PARAMETER["False easting",0,
                 LENGTHUNIT["US survey foot",0.304800609601219],
                 ID["EPSG",8806]],
             PARAMETER["False northing",0,
                 LENGTHUNIT["US survey foot",0.304800609601219],
                 ID["EPSG",8807]]]],
     DERIVINGCONVERSION["Affine transformation as PROJ-based",
         METHOD["PROJ-based operation method: +proj=pipeline  +step +proj=unitconvert +xy_in=m +xy_out=us-ft +z_in=m +z_out=us-ft +step +proj=affine +xoff=5 +yoff=52 +zoff=10 +step +proj=unitconvert +xy_in=us-ft +xy_out=m +z_in=us-ft +z_out=m"]],
     CS[Cartesian,3],
         AXIS["(Y)",north,
             ORDER[1],
             LENGTHUNIT["US survey foot",0.304800609601219,
                 ID["EPSG",9003]]],
         AXIS["(X)",east,
             ORDER[2],
             LENGTHUNIT["US survey foot",0.304800609601219,
                 ID["EPSG",9003]]],
         AXIS["(Z)",up,
             ORDER[3],
             LENGTHUNIT["US survey foot",0.304800609601219,
                 ID["EPSG",9003]]]]'
Candidate operations found: 1
-------------------------------------
Operation No. 1:

unknown id, Inverse of Transformation from NAD83(2011) to NAD83(2011) + Transverse Mercator + Affine transformation as PROJ-based, 0 m, World

PROJ string:
+proj=pipeline
  +step +proj=axisswap +order=2,1
  +step +proj=unitconvert +xy_in=deg +z_in=us-ft +xy_out=rad +z_out=m
  +step +proj=tmerc +lat_0=41 +lon_0=-73 +k=1 +x_0=0 +y_0=0 +ellps=GRS80
  +step +proj=unitconvert +xy_in=m +xy_out=us-ft +z_in=m +z_out=us-ft
  +step +proj=affine +xoff=5 +yoff=52 +zoff=10
  +step +proj=unitconvert +xy_in=us-ft +xy_out=m +z_in=us-ft +z_out=m
  +step +proj=unitconvert +xy_in=m +z_in=m +xy_out=us-ft +z_out=us-ft
  +step +proj=axisswap +order=2,1

Compound to Projected 3D

$ projinfo -o proj -s EPSG:6318+6360 -t 'PROJCRS["NAD83(2011) / Wyoming East (ftUS)",
     BASEGEOGCRS["NAD83(2011)",
         DATUM["NAD83 (National Spatial Reference System 2011)",
             ELLIPSOID["GRS 1980",6378137,298.257222101,
                 LENGTHUNIT["metre",1]]],
         PRIMEM["Greenwich",0,
             ANGLEUNIT["degree",0.0174532925199433]],
         ID["EPSG",6318]],
     CONVERSION["SPCS83 Wyoming East zone (US Survey feet)",
         METHOD["Transverse Mercator",
             ID["EPSG",9807]],
         PARAMETER["Latitude of natural origin",40.5,
             ANGLEUNIT["degree",0.0174532925199433],
             ID["EPSG",8801]],
         PARAMETER["Longitude of natural origin",-105.166666666667,
             ANGLEUNIT["degree",0.0174532925199433],
             ID["EPSG",8802]],
         PARAMETER["Scale factor at natural origin",0.9999375,
             SCALEUNIT["unity",1],
             ID["EPSG",8805]],
         PARAMETER["False easting",656166.6667,
             LENGTHUNIT["US survey foot",0.304800609601219],
             ID["EPSG",8806]],
         PARAMETER["False northing",0,
             LENGTHUNIT["US survey foot",0.304800609601219],
             ID["EPSG",8807]]],
     CS[Cartesian,3],
         AXIS["easting (X)",east,
             ORDER[1],
             LENGTHUNIT["US survey foot",0.304800609601219]],
         AXIS["northing (Y)",north,
             ORDER[2],
             LENGTHUNIT["US survey foot",0.304800609601219]],
         AXIS["up (Z)",up,
             ORDER[3],
             LENGTHUNIT["US survey foot",0.304800609601219]]]'

Candidate operations found: 2
-------------------------------------
Operation No. 1:

unknown id, Conversion from NAVD88 height (ftUS) to NAVD88 height + Inverse of NAD83(2011) to NAVD88 height (3) + SPCS83 Wyoming East zone (US Survey feet), 0.015 m, United States (USA) - CONUS onshore - Alabama; Arizona; Arkansas; California; Colorado; Connecticut; Delaware; Florida; Georgia; Idaho; Illinois; Indiana; Iowa; Kansas; Kentucky; Louisiana; Maine; Maryland; Massachusetts; Michigan; Minnesota; Mississippi; Missouri; Montana; Nebraska; Nevada; New Hampshire; New Jersey; New Mexico; New York; North Carolina; North Dakota; Ohio; Oklahoma; Oregon; Pennsylvania; Rhode Island; South Carolina; South Dakota; Tennessee; Texas; Utah; Vermont; Virginia; Washington; West Virginia; Wisconsin; Wyoming., at least one grid missing

PROJ string:
+proj=pipeline
  +step +proj=axisswap +order=2,1
  +step +proj=unitconvert +xy_in=deg +z_in=us-ft +xy_out=rad +z_out=m
  +step +proj=vgridshift +grids=us_noaa_g2018u0.tif +multiplier=1
  +step +proj=tmerc +lat_0=40.5 +lon_0=-105.166666666667 +k=0.9999375
        +x_0=200000.00001016 +y_0=0 +ellps=GRS80
  +step +proj=unitconvert +xy_in=m +z_in=m +xy_out=us-ft +z_out=us-ft

Grid us_noaa_g2018u0.tif needed but not found on the system. Can be obtained at https://cdn.proj.org/us_noaa_g2018u0.tif

-------------------------------------
Operation No. 2:

unknown id, Transformation from NAVD88 height (ftUS) to NAD83(2011) (ballpark vertical transformation, without ellipsoid height to vertical height correction) + SPCS83 Wyoming East zone (US Survey feet), unknown accuracy, World, has ballpark transformation

PROJ string:
+proj=pipeline
  +step +proj=axisswap +order=2,1
  +step +proj=unitconvert +z_in=us-ft +z_out=m
  +step +proj=unitconvert +xy_in=deg +z_in=m +xy_out=rad +z_out=m
  +step +proj=tmerc +lat_0=40.5 +lon_0=-105.166666666667 +k=0.9999375
        +x_0=200000.00001016 +y_0=0 +ellps=GRS80
  +step +proj=unitconvert +xy_in=m +z_in=m +xy_out=us-ft +z_out=us-ft

Problem description

The transformation should consider the vertical component of the compound, including the units management.

Expected Output

If I guess it correctly, something like this (convert z from us-ft to m, and apply vgridshift [in operation 1])

Operation No. 1:
PROJ string:
+proj=pipeline
  +step +proj=axisswap +order=2,1
  +step +proj=unitconvert +xy_in=deg +z_in=us-ft +xy_out=rad +z_out=m
  +step +proj=vgridshift +grids=us_noaa_g2018u0.tif +multiplier=1
  +step +proj=tmerc +lat_0=41 +lon_0=-73 +k=1 +x_0=0 +y_0=0 +ellps=GRS80
  +step +proj=unitconvert +xy_in=m +xy_out=us-ft +z_in=m +z_out=us-ft
  +step +proj=affine +xoff=5 +yoff=52 +zoff=10
  +step +proj=unitconvert +xy_in=us-ft +xy_out=m +z_in=us-ft +z_out=m
  +step +proj=unitconvert +xy_in=m +z_in=m +xy_out=us-ft +z_out=us-ft
  +step +proj=axisswap +order=2,1

Operation No. 2:
PROJ string:
+proj=pipeline
  +step +proj=axisswap +order=2,1
  +step +proj=unitconvert +z_in=us-ft +z_out=m
  +step +proj=unitconvert +xy_in=deg +z_in=m +xy_out=rad +z_out=m
  +step +proj=tmerc +lat_0=41 +lon_0=-73 +k=1 +x_0=0 +y_0=0 +ellps=GRS80
  +step +proj=unitconvert +xy_in=m +xy_out=us-ft +z_in=m +z_out=us-ft
  +step +proj=affine +xoff=5 +yoff=52 +zoff=10
  +step +proj=unitconvert +xy_in=us-ft +xy_out=m +z_in=us-ft +z_out=m
  +step +proj=unitconvert +xy_in=m +z_in=m +xy_out=us-ft +z_out=us-ft
  +step +proj=axisswap +order=2,1

Environment Information

• PROJ version (proj 9.1.0)
• Operation System Information: Ubuntu 20.04

Installation method

• Compiled from source

[jjimenezshaw]

After debugging into the code, I got some clues, and maybe a solution.

In the method CoordinateOperationFactory::Private::createOperationsCompoundToGeog gets to the code

        // Only do a vertical transformation if the target CRS is 3D.
        const auto dstSingle = dynamic_cast<crs::SingleCRS *>(targetCRS.get());
        if (dstSingle &&
            dstSingle->coordinateSystem()->axisList().size() == 2) {

where tartegCRS is the base of the derived crs, in the example PROJCRS["Transverse Mercator centered in area of interest" ... that is 2D, in spite of its "parent" crs is 3D: DERIVEDPROJCRS["Custom Site Calibrated 3D CRS"

Because it is 2D it is not doing the vertical component operation.

Trying to fix it, I modified the method CoordinateOperationFactory::Private::createOperationsDerivedTo doing:

    crs::CRSNNPtr baseCRS = derivedSrc->baseCRS();
    if (derivedSrc->coordinateSystem()->axisList().size() == 3 && 
        derivedSrc->baseCRS()->coordinateSystem()->axisList().size() == 2)
    {
        baseCRS = derivedSrc->baseCRS()->promoteTo3D(std::string(), nullptr);
    }
    auto opsSecond =
        createOperations(baseCRS, targetCRS, context);

However, a few lines later, when it calls createComputeMetadata, says "Inconsistent chaining of CRS in operations". True, the target of step 0 is 2D, while the source of step 1 is 3D.

So I think the solution is to create the baseCRS directly as 3D.
One option is in the method WKTParser::Private::buildDerivedProjectedCRS, there I add at the end this:

DerivedProjectedCRSNNPtr
WKTParser::Private::buildDerivedProjectedCRS(const WKTNodeNNPtr &node) {
    ...
    if (baseProjCRS->coordinateSystem()->axisList().size() == 2 &&
        cs->axisList().size() == 3) {
        ProjectedCRSPtr aux = std::dynamic_pointer_cast<ProjectedCRS>(
            baseProjCRS->promoteTo3D(std::string(), nullptr).as_nullable());
        if (aux)
            baseProjCRS = NN_NO_CHECK(aux);
    }
    return DerivedProjectedCRS::create(buildProperties(node), baseProjCRS,
                                       conversion, cs);
}

A second option is doing the promoteTo3D directly inside DerivedProjectedCRS::create, promoting to 3D in the same conditions (if the cs is 3D).

The first option would affect only WKT inputs (are there other parsers I should check?), while the second affects every case, not giving a choice to a potential C++ user, but also preventing the problem in the transformation.

Which solution (if any) makes more sense?

[jjimenezshaw]

Ups, I didn't see that @rouault had already implemented a solution (I do not refresh GitHub page enough, and I got no emails this time). Good to see that we reached to the same conclusion :)

[jjimenezshaw]

Thanks @rouault !