proj.ProjTrans from WGS84 WKT to custom WKT produces incorrect results

Example of problem

I have a custom WKT string and I want to transfrom coordinates from WKS84 to that custom system. I am using C# wrapper with proj version 7.1.0 and the following code. The transformed coordinates are different from what NetTopologySuite and also our proprietary projection engine produce. My best guess is that proj does not perform datum transformation between the WGS84 ellipsoid and the sphere defined in the custom coordinate system. Is that so? Why? And what should I do to arrive at the right transformation result?

public void PointTransformationIsOk() 
        {
            var sourceWkt = "GEOGCS[\"WGS 84\", DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\", 6378137, 298.257223563, AUTHORITY[\"EPSG\", \"7030\"]], AUTHORITY[\"EPSG\", \"6326\"]], PRIMEM[\"Greenwich\", 0, AUTHORITY[\"EPSG\", \"8901\"]], UNIT[\"degree\", 0.0174532925199433, AUTHORITY[\"EPSG\", \"9122\"]], AUTHORITY[\"EPSG\", \"4326\"]]";
            var targetWkt = "PROJCS[\"Custom Projection\",GEOGCS[\"Unused\",DATUM[\"User defined\",SPHEROID[\"Sphere(Radius = 6370000)\",6370000,0]],PRIMEM[\"Greenwich\",0],UNIT[\"Degree\",0.0174532925199433]],PROJECTION[\"Lambert_Conformal_Conic_2SP\"],PARAMETER[\"False_Easting\",0],PARAMETER[\"False_Northing\",0],PARAMETER[\"Central_Meridian\",125],PARAMETER[\"Standard_Parallel_1\",30],PARAMETER[\"Standard_Parallel_2\",60],PARAMETER[\"Latitude_Of_Origin\",21.3000335693359],UNIT[\"Meter\",1]]"; ;

            var x = 120.80981;
            var y = 22.05211;
            var z = 0.0;

            var ctx = proj.ProjContextCreate();
            var crs2Crs = proj.ProjCreateCrsToCrs(ctx, sourceWkt, targetWkt, null);
            var coords = proj.ProjCoord(x, y, z);
            var result = proj.ProjTrans(crs2Crs, PJ_DIRECTION.PJ_FWD, coords);

            var xp = result.Enu.E;
            var yp = result.Enu.N;
            var zp = result.Enu.U;

            // Represented as WKT strings:
            // PROJ.............: POINT Z(-450415.3359880044  99264.43806309221 0)
            // NetTopologySuite.: POINT Z(-451226.5228159     83777.59511600249 0)
            // Proprietary......: POINT Z(-451226.52281590283 83777.59511600249 0)
        }

Problem description

We get almost the same results from NetTopologySuite and from our Proprietary projection engine, but PROJ produces different results. We would like to reconcile all three libraries.

Expected Output

POINT Z(-451226.5228159 83777.59511600249 0)

Environment Information

PROJ version 7.1.0
Operation System Information, Windows 7

Installation method

cmake and a c# wrapper

PROJ behaves as designed here. The issue is that you're transforming between the WGS 84 datum and a user defined datum, that has no expressed relationship with WGS 84. Hence PROJ will skip any datum transformation. To get the result of the other software, you'd need to add a TOWGS84[0,0,0] child node to the DATUM[] of your target WKT

$ echo 22.05211 120.80981 0 | cs2cs "GEOGCS[\"WGS 84\", DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\", 6378137, 298.257223563, AUTHORITY[\"EPSG\", \"7030\"]], AUTHORITY[\"EPSG\", \"6326\"]], PRIMEM[\"Greenwich\", 0, AUTHORITY[\"EPSG\", \"8901\"]], UNIT[\"degree\", 0.0174532925199433, AUTHORITY[\"EPSG\", \"9122\"]], AUTHORITY[\"EPSG\", \"4326\"]]"   "PROJCS[\"Custom Projection\",GEOGCS[\"Unused\",DATUM[\"User defined\",SPHEROID[\"Sphere(Radius = 6370000)\",6370000,0],TOWGS84[0,0,0]],PRIMEM[\"Greenwich\",0],UNIT[\"Degree\",0.0174532925199433]],PROJECTION[\"Lambert_Conformal_Conic_2SP\"],PARAMETER[\"False_Easting\",0],PARAMETER[\"False_Northing\",0],PARAMETER[\"Central_Meridian\",125],PARAMETER[\"Standard_Parallel_1\",30],PARAMETER[\"Standard_Parallel_2\",60],PARAMETER[\"Latitude_Of_Origin\",21.3000335693359],UNIT[\"Meter\",1]]"
-451226.52  83777.60 0.00

OSGeo / PROJ