PROJ projection does not match with epsg.io/transform

[filipkral]

Describe the bug I'm getting unexpected results when projecting from British National Grid to Lambert Equal Area (on ETRS1989 datum).

To Reproduce Explore the values asserted in this code

using MaxRev.Gdal.Core;
using OSGeo.OSR;
using System;
using System.Threading.Tasks;
using Xunit;

namespace DHI.WaterData.Projections.Provider.Test.WKTS
{
    public class MatchWKTStringsTest
    {
        [Fact]
        public async Task TransformPointIsOk()
        {
            var x = 826158.063;
            var y = 2405844.125;

            var sourceWkt =
                "PROJCS[\"OSGB 1936 / British National Grid\",GEOGCS[\"OSGB 1936\",DATUM[\"OSGB_1936\",SPHEROID[\"Airy 1830\",6377563.396,299.3249646,AUTHORITY[\"EPSG\",\"7001\"]],TOWGS84[446.448,-125.157,542.06,0.15,0.247,0.842,-20.489],AUTHORITY[\"EPSG\",\"6277\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4277\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",49],PARAMETER[\"central_meridian\",-2],PARAMETER[\"scale_factor\",0.9996012717],PARAMETER[\"false_easting\",400000],PARAMETER[\"false_northing\",-100000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH],AUTHORITY[\"EPSG\",\"27700\"]]";
            var targetWkt =
                "PROJCS[\"ETRS89 / LAEA Europe\",GEOGCS[\"ETRS89\",DATUM[\"European_Terrestrial_Reference_System_1989\",SPHEROID[\"GRS 1980\",6378137,298.257222101,AUTHORITY[\"EPSG\",\"7019\"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY[\"EPSG\",\"6258\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4258\"]],PROJECTION[\"Lambert_Azimuthal_Equal_Area\"],PARAMETER[\"latitude_of_center\",52],PARAMETER[\"longitude_of_center\",10],PARAMETER[\"false_easting\",4321000],PARAMETER[\"false_northing\",3210000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AUTHORITY[\"EPSG\",\"3035\"]]";

            Proj6.Configure();

            var sourceSr = new SpatialReference(string.Empty);
            sourceSr.ImportFromWkt(ref sourceWkt);
            sourceSr.SetAxisMappingStrategy(AxisMappingStrategy.OAMS_TRADITIONAL_GIS_ORDER);

            var targetSr = new SpatialReference(string.Empty);
            targetSr.ImportFromWkt(ref targetWkt);
            targetSr.SetAxisMappingStrategy(AxisMappingStrategy.OAMS_TRADITIONAL_GIS_ORDER);

            var transformation = new CoordinateTransformation(sourceSr, targetSr);
            if (transformation == null)
                throw new ArgumentException("Invalid projection WKT");

            var projected = new double[3];
            transformation.TransformPoint(projected, x, y, 0.0);
            var px = projected[0];
            var py = projected[1];

            var ex = 4316131D;
            var ey = 5330926D;
            Assert.Equal(ex, px, 0);
            Assert.Equal(ey, py, 0);
        }
    }
}

Expected behavior coord	British n. grid	expected projected	actual projected
x	826158.063	4316131.0	4316331.549
y	2405844.125	5330926.0	5331101.982

The expected projected coordinates come from: https://epsg.io/transform#s_srs=27700&t_srs=3035&x=826158.0630000&y=2405844.1250000 The expected projected coordinates were also retuned from another, proprietary. transformation library I used.

Environment information:

• OS (version): Windows 10 and Windows 11.
• Package version (core): 3.2.0.320
• Package version (runtime): 3.2.0.250

Additional context This may be ultimately a PROJ issue, but I'm using it through MaxRev.gdal so this seems to be the right place to start. I've been researching about PROJ (and the changes in version 4, 6, and onwards), but nothing really explained this discrepancy. I suspect the datum transformation is not applied because when I changed the TOWGS parameters of the British National Grid WKT string to all zeros, I got exactly the same results. Also the SetAxisMappingStrategy calls have no impact in this case. It may also be that we are missing the transformation grids needed for this transformation. But how do we install grid shift files with MaxRev.gdal?

[MaxRev-Dev]

Not a bug. I'm building these packages from original sources using the release branch. Just checked via docker run and the result from test matches with output. So maybe, there is some error in the web implementation. Anyway, if you do not mind, I will add your test to this repo for coverage.

docker run --rm osgeo/proj /bin/bash -c "echo 826158.063 2405844.125 | cs2cs +init=epsg:27700 +to +init=epsg:3035"

Output (docker): 4316331.55 5331101.98 0.00

[filipkral]

OK, I took it with https://github.com/OSGeo/PROJ/issues/2955