[PROJ] Best practices for converting 2D CRS to 3D CRS

[Sean Lilley]

We occasionally receive datasets that have height values in feet and only
specify a 2D CRS. It seems that recent versions of proj/gdal (6.3.0/3.0.3)
promote to 3D by assuming the z-axis is meters rather than inferring the
height units from other components. Earlier versions would infer the units
as feet which was what we expected for those datasets.

I want to confirm if this change is intentional. I'm guessing that it is
and that promoting 2D to 3D is not a simple problem. Nevertheless does
anyone have a workaround that infers the expected vertical units from the
2D CRS most of the time?

I'm using the gdal C++ api and inferring the height units from the 2D CRS's
linear units. I haven't tried building this exact snippet so there may be
compiler errors, but it gives the general idea.

            int inputEpsg = 2263;
            int outputEpsg = 4979;
            double x = 1002759.79006824;
            double y = 220148.669988647;
            double z = 1797.1066;
            double expectedHeight = 547.759;

            OGRSpatialReference inputSpatialReference;
            inputSpatialReference.importFromEPSG(inputEpsg);

inputSpatialReference.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);

            // Guess the units of the vertical axis for 2D CRS
            double verticalToMeters = 1.0;
            int axesCount = inputSpatialReference.GetAxesCount();
            if (axesCount == 2)
            {
                verticalToMeters = inputSpatialReference.GetLinearUnits();
            }

            OGRSpatialReference outputSpatialReference;
            outputSpatialReference.importFromEPSG(outputEpsg);

outputSpatialReference.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);

            OGRCoordinateTransformation* transformation
= OGRCreateCoordinateTransformation(&inputSpatialReference,
&outputSpatialReference);
            transformation->Transform(1, &x, &y, &z);

            z *= verticalToMeters;
            assert(z == expectedHeight);

Is there a cleaner way for me to achieve what I want to do in the gdal C++
API? Is inferring the height units from the 2D CRS linear units generally
an ok approach?

Thanks,
Sean
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20200228/afd5fb1c/attachment.html>