Incorrect coordinate conversion results from ProjNET when converting across horizontal datums

Question

Been fighting this for almost a week now, so I'm hoping that someone here has some insight or can tell me the obvious thing I'm missing. I am trying to add programmatic coordinate conversion on the end of one of the products I've developed in C#. My boss recommended using ProjNET, so I got it from NuGet and implemented the code without much trouble.

At the testing stage, however, things broke. Through trial and error, I found that my results were perfect if the client stayed within the same datum, but any cross-datum conversions were incorrect. Our use case is primarily cross-datum conversions.

I did a bunch of searching and found that ProjNET requires both the source coordinate system and target coordinate system to include the ToWGS84 parameter: https://github.com/NetTopologySuite/ProjNet4GeoAPI/wiki/Projecting-points-from-one-coordinate-system-to-another

Ironically (and tragically), the ProjNET.SRID library on NuGet, which translates SRIDs into WKTs, does not provide the ToWGS84 parameter in their WKTs. My boss and I found another source to translate SRIDs that includes the ToWGS84 parameter, but the results are still incorrect, possibly because the source is old and the WKTs might be outdated.

Here's the core of the code that's giving me trouble, using Nevada West NAD27 and Nevada West NAD83, both in US Survey feet:

ProjNet.CoordinateSystems.CoordinateSystemFactory csFactory = new CoordinateSystemFactory();
ProjNet.CoordinateSystems.CoordinateSystem sourceCoordSys = csFactory.CreateFromWkt("PROJCS[\"NAD27 / Nevada West\",  GEOGCS[\"NAD27\",    DATUM[\"North American Datum 1927\",      SPHEROID[\"Clarke 1866\", 6378206.4, 294.9786982138982, AUTHORITY[\"EPSG\",\"7008\"]],      TOWGS84[-4.2, 135.4, 181.9, 0.0, 0.0, 0.0, 0.0],      AUTHORITY[\"EPSG\",\"6267\"]],    PRIMEM[\"Greenwich\", 0.0, AUTHORITY[\"EPSG\",\"8901\"]],    UNIT[\"degree\", 0.017453292519943295],    AXIS[\"Geodetic longitude\", EAST],    AXIS[\"Geodetic latitude\", NORTH],    AUTHORITY[\"EPSG\",\"4267\"]],  PROJECTION[\"Transverse Mercator\"],  PARAMETER[\"central_meridian\", -118.58333333333333],  PARAMETER[\"latitude_of_origin\", 34.75],  PARAMETER[\"scale_factor\", 0.9999],  PARAMETER[\"false_easting\", 500000.0],  PARAMETER[\"false_northing\", 0.0],  UNIT[\"foot_survey_us\", 0.30480060960121924],  AXIS[\"Easting\", EAST],  AXIS[\"Northing\", NORTH],  AUTHORITY[\"EPSG\",\"32009\"]]");
ProjNet.CoordinateSystems.CoordinateSystem targetCoordSys = csFactory.CreateFromWkt("PROJCS[\"NAD83 / Nevada West (ft US)\",  GEOGCS[\"NAD83\",    DATUM[\"North American Datum 1983\",      SPHEROID[\"GRS 1980\", 6378137.0, 298.257222101, AUTHORITY[\"EPSG\",\"7019\"]],      TOWGS84[1.0, 1.0, -1.0, 0.0, 0.0, 0.0, 0.0],      AUTHORITY[\"EPSG\",\"6269\"]],    PRIMEM[\"Greenwich\", 0.0, AUTHORITY[\"EPSG\",\"8901\"]],    UNIT[\"degree\", 0.017453292519943295],    AXIS[\"Geodetic longitude\", EAST],    AXIS[\"Geodetic latitude\", NORTH],    AUTHORITY[\"EPSG\",\"4269\"]],  PROJECTION[\"Transverse Mercator\"],  PARAMETER[\"central_meridian\", -118.58333333333333],  PARAMETER[\"latitude_of_origin\", 34.75],  PARAMETER[\"scale_factor\", 0.9999],  PARAMETER[\"false_easting\", 2624666.6667],  PARAMETER[\"false_northing\", 13123333.3333],  UNIT[\"foot_survey_us\", 0.30480060960121924],  AXIS[\"Easting\", EAST],  AXIS[\"Northing\", NORTH],  AUTHORITY[\"EPSG\",\"3423\"]]");

CoordinateTransformationFactory transformFactory = new CoordinateTransformationFactory();
ICoordinateTransformation coordTransform = transformFactory.CreateFromCoordinateSystems(sourceCoordSys, targetCoordSys);

double XTransfer, YTransfer;
<For Loop that iterates through resultTable>
{
    (XTransfer, YTransfer) = coordTransform.MathTransform.Transform(
        (double)resultTable.Rows[i][eastingColName],
        (double)resultTable.Rows[i][northingColName]);

    resultTable.Rows[i][eastingColName] = XTransfer;
    resultTable.Rows[i][northingColName] = YTransfer;
}

And here's one of the test points I've been trying to correct:

Source SRID = 32009
Easting = 214491.354
Northing = 679690.231

And the results for target SRID = 3423:

From https://mygeodata.cloud/cs2cs/

Easting = 2338881.88508
Northing = 13803025.9945

From my code:

Easting = 2338932.15719927, difference of 50.272 US Survey feet
Northing = 13803002.2532437, difference of -23.741 US Survey feet

As an aside, we are not devoted to ProjNET, so I tried implementing the above in DotSpatial. I had the exact same experience of breaking down in cross-datum transformations. They have the solution to use grid shift files, but those didn't appear to have any effect once loaded. If anyone has insight into those, I'll gladly take it.

This is a similar question to here, but I cannot comment there due to the newness of my account, I obviously do not have an answer, and even the official strings from EPSG (recommended in a comment) do not contain the ToWGS84 parameter.

Answer 1

Unfortunately I found an answer that is unsatisfying, technically. I was never able to get the cross-datum transforms to be consistently accurate to less than 5 US feet using ProjNET, DotSpatial, or GDAL. GDAL got the closest, with a maximum of ~9 US feet deviation (that I observed), but never less than 1 US foot deviation.

Ultimately I ended up using a paid conversion library: https://effectiveobjects.com/prolatnet5

that is accurate to the thousandths of a US foot according to umpire I was using. It uses an odd classification system, but it automatically handles the decisions of when to use WGS84 vs Grid Shift, plus the implementation of both.

If you can't afford ProLat and are okay with minor relative inaccuracy, GDAL ended up being relatively easy to use once I found the right way to set it up: https://stackoverflow.com/a/37577459/1655707

It might take some tinkering even after installing from NuGET.

The technically satisfying answer would have been to figure out how to get the free libraries to correctly and automatically use their conversions that they appear to support, but ProLat works and I've wasted too much company time on this already.