orbisgis / cts

Projection library in Java
GNU Lesser General Public License v3.0
49 stars 15 forks source link

Why result of transformation is not stable ? #158

Open mcherb opened 1 month ago

mcherb commented 1 month ago

Maybe I use this lib badly but why the results of the same conversion is not the same between different runs.

@Test
    void test_crs() throws IllegalCoordinateException, CoordinateOperationException, CRSException {
        for (int j = 0; j < 10; j++) {

            CRSFactory cRSFactory = new CRSFactory();

            RegistryManager registryManager = cRSFactory.getRegistryManager();
            registryManager.addRegistry(new EPSGRegistry());

            CoordinateReferenceSystem sourceSystem = cRSFactory.getCRS("EPSG:2154");
            CoordinateReferenceSystem destinationSystem = cRSFactory.getCRS("EPSG:4326");

            CoordinateOperation coordinateOperation = CoordinateOperationFactory.createCoordinateOperations((GeodeticCRS) sourceSystem, (GeodeticCRS) destinationSystem).iterator().next();

            double[] result = coordinateOperation.transform(new double[]{615599.0, 6840122.0});
            System.out.println("X == " + result[0] + "    Y == " + result[1]);
        }
    }

When running the code above it shows me this logs :

X == 1.8540873946173166 Y == 48.65575744955268 X == 1.8540873946173166 Y == 48.65575744955268 X == 1.8540873946173166 Y == 48.65575744955268 X == 1.8540873946173166 Y == 48.65575744955268 X == 1.8540873946173166 Y == 48.65575744955268 X == 1.8540873946173166 Y == 48.65575744955268 X == 1.854087394617317 Y == 48.65575745048813 X == 1.854087394617317 Y == 48.65575745048813 X == 1.8540873946173166 Y == 48.65575744955268 X == 1.854087394617317 Y == 48.65575745048813

Is there a way to fixe this behavior ?

ebocher commented 3 weeks ago

This is due to a rounding precision. You have a 10-6 precision from meters to lat/lon that seems to be enough

mcherb commented 3 weeks ago

The bug seems to be produced by the initialization phase.

When I change the code like below I get a stable result

    @Test
    void test_crs() throws IllegalCoordinateException, CoordinateOperationException, CRSException {

        CRSFactory cRSFactory = new CRSFactory();

        RegistryManager registryManager = cRSFactory.getRegistryManager();
        registryManager.addRegistry(new EPSGRegistry());
        CoordinateReferenceSystem sourceSystem = cRSFactory.getCRS("EPSG:2154");
        CoordinateReferenceSystem destinationSystem = cRSFactory.getCRS("EPSG:4326");

        for (int j = 0; j < 10; j++) {

            CoordinateOperation coordinateOperation = CoordinateOperationFactory.createCoordinateOperations((GeodeticCRS) sourceSystem, (GeodeticCRS) destinationSystem).iterator().next();

            double[] result = coordinateOperation.transform(new double[]{615599.0, 6840122.0});
            System.out.println("X == " + result[0] + "    Y == " + result[1]);
        }
    }

X == 1.8540873946173166 Y == 48.65575744955268 X == 1.8540873946173166 Y == 48.65575744955268 X == 1.8540873946173166 Y == 48.65575744955268 X == 1.8540873946173166 Y == 48.65575744955268 X == 1.8540873946173166 Y == 48.65575744955268 X == 1.8540873946173166 Y == 48.65575744955268 X == 1.8540873946173166 Y == 48.65575744955268 X == 1.8540873946173166 Y == 48.65575744955268 X == 1.8540873946173166 Y == 48.65575744955268 X == 1.8540873946173166 Y == 48.65575744955268

mcherb commented 3 weeks ago

@ebocher

This is due to a rounding precision. You have a 10-6 precision from meters to lat/lon that seems to be enough

But what explains the difference ?

mukoki commented 20 hours ago

Hi, Weird behaviour ! I can give some explanations. Will be more difficult to fix it properly. The normal way to go from 2154 to WGS84 is : Lambert 93 -> Geographic (RGF93) (= Geographic WGS 84) Here, we have recognized that the target datum (WGS884) is compatible with the source datum (RGF93), and that there is no transformation to go from the former to the second. If we don't recognize those datum are "compatible", the full transformation becomes : Lambert 93 -> Geographic (RGF93) -> Cartesian RGF93 -> Cartesian WGS84 -> Geographic (WGS84) It is almost the same as there is no transformation from RGF93 to WGS84, but as they have slightly different ellipsoids, at the end, coordinates may be a bit different. The first time, the first coordinate operation is used and it is correct. At the same time, CTS builds a new transformation from RGF93 to WGS84 using the full path and records it as a possible path from Lambert93 to WGS84. The method you used to choose an appropriate coordinate operation returns an unordered set of possible paths from source crs to target crs. The first time, there is only one possible path (i.e. coordinate operation), but after that, the operation following the full path through the different datum has been added as a new possible path. Why the first element of the returned set is not always the same ? I can't explain. But there is no guarantee about element order in a set. To make the choice more deterministic, there is two helper methods in CoordinateOperationFactory : getMostPrecise and getMostPrecise3DTransformation.

And of course, computing the CoordinateOperation out of the loop is a good thing to do !

Hope it helps