Incorrect unit conversion on affine transformation

[ktheijs]

What is the bug?

I get different results between the cs2cs app and the c# bindings when trying to transform a coordinate to a simple derived CRS. This only happens for coordinate systems that are not in foot.

I have prepared a basic example where I use a derived CRS to offset the X coordinate by 1000ftUS.

The result when using cs2cs is what I would expect. The c# version first seems to do an unnecessary unit conversion.

Input: 2300000.00 2000000.00 0.0
Output cs2cs: 2301000.00 2000000.00 0.00
Output c#: 702041.4020828041, 609601.2192024384, 0

Steps to reproduce the issue

Using C# bindings:

SpatialReference sourceCs = new SpatialReference("");
sourceCs.ImportFromEPSG(6453);
string targetCsWkt = @"DERIVEDPROJCRS[""Ground for NAD83(2011) / Idaho West (ftUS)"",BASEPROJCRS[""NAD83(2011) / Idaho West (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 Idaho West zone (US Survey feet)"",METHOD[""Transverse Mercator"",ID[""EPSG"",9807]],PARAMETER[""Latitude of natural origin"",41.6666666666667,ANGLEUNIT[""degree"",0.0174532925199433],ID[""EPSG"",8801]],PARAMETER[""Longitude of natural origin"",-115.75,ANGLEUNIT[""degree"",0.0174532925199433],ID[""EPSG"",8802]],PARAMETER[""Scale factor at natural origin"",0.999933333,SCALEUNIT[""unity"",1],ID[""EPSG"",8805]],PARAMETER[""False easting"",2624666.667,LENGTHUNIT[""US survey foot"",0.304800609601219],ID[""EPSG"",8806]],PARAMETER[""False northing"",0,LENGTHUNIT[""US survey foot"",0.304800609601219],ID[""EPSG"",8807]]],CS[Cartesian,2],AXIS[""easting (X)"",east,ORDER[1],LENGTHUNIT[""US survey foot"",0.304800609601219]],AXIS[""northing (Y)"",north,ORDER[2],LENGTHUNIT[""US survey foot"",0.304800609601219]],USAGE[SCOPE[""Engineering survey, topographic mapping.""],AREA[""United States (USA) - Idaho - counties of Ada; Adams; Benewah; Boise; Bonner; Boundary; Canyon; Clearwater; Elmore; Gem; Idaho; Kootenai; Latah; Lewis; Nez Perce; Owyhee; Payette; Shoshone; Valley; Washington.""],BBOX[41.99,-117.24,49.01,-114.32]],ID[""EPSG"",6453]],DERIVINGCONVERSION[""Grid to ground"",METHOD[""Similarity transformation"",ID[""EPSG"",9621]],PARAMETER[""Ordinate 1 of evaluation point in target CRS"",1000,LENGTHUNIT[""US survey foot"",0.304800609601219],ID[""EPSG"",8621]],PARAMETER[""Ordinate 2 of evaluation point in target CRS"",0,LENGTHUNIT[""US survey foot"",0.304800609601219],ID[""EPSG"",8622]],PARAMETER[""Scale factor for source CRS axes"",1,SCALEUNIT[""unity"",1],ID[""EPSG"",1061]],PARAMETER[""Rotation angle of source CRS axes"", 0,ANGLEUNIT[""degree"",0.0],ID[""EPSG"",8614]]],CS[Cartesian,2],AXIS[""easting (X)"",east,ORDER[1],LENGTHUNIT[""US survey foot"",0.304800609601219]],AXIS[""northing (Y)"",north,ORDER[2],LENGTHUNIT[""US survey foot"",0.304800609601219]]]";
SpatialReference targetCs= new SpatialReference(targetCsWkt);
CoordinateTransformation transform = new CoordinateTransformation(sourceCs, targetCs);
double[] argOut = new double[3];
t.TransformPoint(argOut, 2300000, 2000000, 0);
// Result = {702041.4020828041, 609601.2192024384, 0}

Using cs2cs in command line:

C:\Temp>cs2cs +init=epsg:6453 +to "DERIVEDPROJCRS[\"Ground for NAD83(2011) / Idaho West (ftUS)\",BASEPROJCRS[\"NAD83(2011) / Idaho West (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 Idaho West zone (US Survey feet)\",METHOD[\"Transverse Mercator\",ID[\"EPSG\",9807]],PARAMETER[\"Latitude of natural origin\",41.6666666666667,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8801]],PARAMETER[\"Longitude of natural origin\",-115.75,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8802]],PARAMETER[\"Scale factor at natural origin\",0.999933333,SCALEUNIT[\"unity\",1],ID[\"EPSG\",8805]],PARAMETER[\"False easting\",2624666.667,LENGTHUNIT[\"US survey foot\",0.304800609601219],ID[\"EPSG\",8806]],PARAMETER[\"False northing\",0,LENGTHUNIT[\"US survey foot\",0.304800609601219],ID[\"EPSG\",8807]]],CS[Cartesian,2],AXIS[\"easting (X)\",east,ORDER[1],LENGTHUNIT[\"US survey foot\",0.304800609601219]],AXIS[\"northing (Y)\",north,ORDER[2],LENGTHUNIT[\"US survey foot\",0.304800609601219]],USAGE[SCOPE[\"Engineering survey, topographic mapping.\"],AREA[\"United States (USA) - Idaho - counties of Ada; Adams; Benewah; Boise; Bonner; Boundary; Canyon; Clearwater; Elmore; Gem; Idaho; Kootenai; Latah; Lewis; Nez Perce; Owyhee; Payette; Shoshone; Valley; Washington.\"],BBOX[41.99,-117.24,49.01,-114.32]],ID[\"EPSG\",6453]],DERIVINGCONVERSION[\"Grid to ground\",METHOD[\"Similarity transformation\",ID[\"EPSG\",9621]],PARAMETER[\"Ordinate 1 of evaluation point in target CRS\",1000,LENGTHUNIT[\"US survey foot\",0.304800609601219],ID[\"EPSG\",8621]],PARAMETER[\"Ordinate 2 of evaluation point in target CRS\",0,LENGTHUNIT[\"US survey foot\",0.304800609601219],ID[\"EPSG\",8622]],PARAMETER[\"Scale factor for source CRS axes\",1,SCALEUNIT[\"unity\",1],ID[\"EPSG\",1061]],PARAMETER[\"Rotation angle of source CRS axes\", 0,ANGLEUNIT[\"degree\",0.0],ID[\"EPSG\",8614]]],CS[Cartesian,2],AXIS[\"easting (X)\",east,ORDER[1],LENGTHUNIT[\"US survey foot\",0.304800609601219]],AXIS[\"northing (Y)\",north,ORDER[2],LENGTHUNIT[\"US survey foot\",0.304800609601219]]]" --3d - E
2300000.00 2000000.00 0.0
2301000.00      2000000.00 0.00

For readability I repeat the WKT of the derived CRS:

DERIVEDPROJCRS["Ground for NAD83(2011) / Idaho West (ftUS)",
            BASEPROJCRS["NAD83(2011) / Idaho West (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 Idaho West zone (US Survey feet)",
					METHOD["Transverse Mercator",ID["EPSG",9807]],
					PARAMETER["Latitude of natural origin",41.6666666666667,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],
					PARAMETER["Longitude of natural origin",-115.75,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],
					PARAMETER["Scale factor at natural origin",0.999933333,SCALEUNIT["unity",1],ID["EPSG",8805]],
					PARAMETER["False easting",2624666.667,LENGTHUNIT["US survey foot",0.304800609601219],ID["EPSG",8806]],
					PARAMETER["False northing",0,LENGTHUNIT["US survey foot",0.304800609601219],ID["EPSG",8807]]],
				CS[Cartesian,2],
				AXIS["easting (X)",east,ORDER[1],LENGTHUNIT["US survey foot",0.304800609601219]],
				AXIS["northing (Y)",north,ORDER[2],LENGTHUNIT["US survey foot",0.304800609601219]],
				USAGE[SCOPE["Engineering survey, topographic mapping."],AREA["United States (USA) - Idaho - counties of Ada; Adams; Benewah; Boise; Bonner; Boundary; Canyon; Clearwater; Elmore; Gem; Idaho; Kootenai; Latah; Lewis; Nez Perce; Owyhee; Payette; Shoshone; Valley; Washington."],BBOX[41.99,-117.24,49.01,-114.32]],
				ID["EPSG",6453]],
            DERIVINGCONVERSION["Grid to ground",
            	METHOD["Similarity transformation",ID["EPSG",9621]],
            	PARAMETER["Ordinate 1 of evaluation point in target CRS",1000,LENGTHUNIT["US survey foot",0.304800609601219],ID["EPSG",8621]],
            	PARAMETER["Ordinate 2 of evaluation point in target CRS",0,LENGTHUNIT["US survey foot",0.304800609601219],ID["EPSG",8622]],
            	PARAMETER["Scale factor for source CRS axes",1,SCALEUNIT["unity",1],ID["EPSG",1061]],
            	PARAMETER["Rotation angle of source CRS axes", 0,ANGLEUNIT["degree",0.0],ID["EPSG",8614]]],
            CS[Cartesian,2],
            AXIS["easting (X)",east,
            	ORDER[1],
            	LENGTHUNIT["US survey foot",0.304800609601219]],
            AXIS["northing (Y)",north,
            	ORDER[2],
            	LENGTHUNIT["US survey foot",0.304800609601219]]]

Versions and provenance

I'm using GDAL 3.8.4 (with proj 9.3.0) downloaded from https://www.gisinternals.com/release.php.

Additional context

No response

[jratike80]

Projinfo prints a warning, I do not know what it means. At least changing CS into ID is not a right fix.

Warning: Parsing error : syntax error, unexpected CS, expecting ID. Error occurred around: t",0.304800609601219],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["easting (X)",east,

[ktheijs]

@jratike80 The warning is caused by the CS definition in the BASEPROJCRS. You could remove that part but it results in the same issue.

[rouault]

email sent to crs.swg@lists.ogc.org:

[CRS.SWG] WKT: DERIVEDPROJCRS should allow CS for BASEPROJCRS
Hi,

While analyzing /~https://github.com/OSGeo/gdal/issues/9732, I'm realizing that there's a flow in the grammar of DERIVEDPROJCRS, at least when the BASEPROJCRS isn't an identified object. Currently it doesn't allow the BASEPROJCRS to have a CS node. But this leads to the following issue.

Let's consider the below WKT definition that complies with the current grammar,

DERIVEDPROJCRS["Ground for NAD83(2011) / Idaho West (ftUS)",
            BASEPROJCRS["unnamed BASEPROJCRS",
                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 Idaho West zone (US Survey feet)",
                    METHOD["Transverse Mercator",ID["EPSG",9807]],
                    PARAMETER["Latitude of natural origin",41.6666666666667,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],
                    PARAMETER["Longitude of natural origin",-115.75,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],
                    PARAMETER["Scale factor at natural origin",0.999933333,SCALEUNIT["unity",1],ID["EPSG",8805]],
                    PARAMETER["False easting",2624666.667,LENGTHUNIT["US survey foot",0.304800609601219],ID["EPSG",8806]],
                    PARAMETER["False northing",0,LENGTHUNIT["US survey foot",0.304800609601219],ID["EPSG",8807]]],
                ],
            DERIVINGCONVERSION["Grid to ground",
                METHOD["Similarity transformation",ID["EPSG",9621]],
                PARAMETER["Ordinate 1 of evaluation point in target CRS",1000,LENGTHUNIT["US survey foot",0.304800609601219],ID["EPSG",8621]],
                PARAMETER["Ordinate 2 of evaluation point in target CRS",0,LENGTHUNIT["US survey foot",0.304800609601219],ID["EPSG",8622]],
                PARAMETER["Scale factor for source CRS axes",1,SCALEUNIT["unity",1],ID["EPSG",1061]],
                PARAMETER["Rotation angle of source CRS axes", 0,ANGLEUNIT["degree",0.0],ID["EPSG",8614]]],
            CS[Cartesian,2],
            AXIS["easting (X)",east,
                ORDER[1],
                LENGTHUNIT["US survey foot",0.304800609601219]],
            AXIS["northing (Y)",north,
                ORDER[2],
                LENGTHUNIT["US survey foot",0.304800609601219]]]

Here, I've on purpose renamed the BASEPROJCRS to "unnamed BASEPROJCRS" to demonstrate the case where we can't do a database lookup to retrieve the CS of the BASEPROJCRS (in the original report, the BASEPROJCRS is EPSG:6453 "NAD83(2011) / Idaho West (ftUS)" that has Us survey foot axis.)

Now let's say a user asks to transform from EPSG:6453 to the above DERIVEDPROJCRS.  The logic an engine like PROJ applies is to chain a transformation from EPSG:6453 to the baseprojCRS of the DERIVEDPROJCRS, and then apply the DerivingConversion of the DERIVEDPROJCRS. All good, but when the CS[] of BASEPROJCRS isn't defined, then this BASEPROJCRS isn't a fully well defined Projected CRS. In particular due to the absence of the CS[], you may wrongly infer a metre based CS, and then this will cause the transformation from EPSG:6453 to BASEPROJCRS to have a unwanted wrong us-ft to metre conversion.

So we should probably allow CS[] as a valid subnode of BASEPROJCRS[] to avoid any ambiguity  (by the way, while exploring the issue, I see that I don't have this issue with PROJJSON as the derived_projected_crs grammar just links to the standard projected_crs definition for the base_crs, avoiding any "optimization" like done in WKT)

Even 

Working on workaround on PROJ (and potentially GDAL) side