How to use GeoTools MathTransform?

Question

I need to transform points from EPSG:4312 to WGS84 using GeoTools Java library. But I'm not sure if I use it correctly, also I'm not sure it provides correct results.

See the following code sample:

    @Test
    public void testSingle4312ToWGS84() throws FactoryException, TransformException {
        double latitude = 29.0;
        double longitude = -99.0;

        CoordinateReferenceSystem sourceCrs = CRS.decode("EPSG:4312");
        CoordinateReferenceSystem targetCrs = CRS.decode("EPSG:4326");
        MathTransform engine = CRS.findMathTransform(sourceCrs, targetCrs, false);

        DirectPosition source = new DirectPosition2D(sourceCrs, longitude, latitude);
        DirectPosition target = new DirectPosition2D(targetCrs);
        engine.transform(source, target);

        System.out.println("x,y=" + target.getCoordinate()[0] + "," + target.getCoordinate()[1]);
    }

When running above code I get the following results:

x,y=-81.00483829765083,-150.99434404015307

Which is wrong, because -150 is not a valid latitude.

When I created the source object, I initialized X with longitude and Y with latitude, following the rationale that actually X is longitude and Y is latitude (according to the geographic community). However, if I reverse them:

DirectPosition source = new DirectPosition2D(sourceCrs, latitude, longitude);

I now get a result that looks more real:

x,y=29.003866857343915,-98.99222530180596

Of course I have to interpret X as latitude and Y as longitude now.

So, question #1, should we go with the inverse rational treating X as latitude and Y as longitude?

Next, I wanted to validate the result. First, trying to convert the same point using an Oracle SQL SDO function:

select sdo_cs.transform (
MDSYS.SDO_GEOMETRY(2001, 4312, MDSYS.SDO_POINT_TYPE(-99, 29, NULL), NULL, NULL),
8307
) from dual;

provides the following result:

MDSYS.SDO_GEOMETRY(2001, 8307, MDSYS.SDO_POINT_TYPE(-98.9924748682774, 29.0036029261273, NULL), NULL, NULL)

The converted point is similar with the one produced by GeoTools, but there is a difference of ~40 meters.

Next, I tried the transformation offered by this website, which provided the same results as Oracle. I also used Luciad AIXM Viewer, a graphical representation tool for GML, and this also displays my point at the same position computed by Oracle.

So, I have 3 independent tools providing the same result point which is 40 meters away from the one computed by GeoTools. Question #2, why is Geotools failing? Is Geotools reliable for CRS transformations?

Answer 1

For your first issue with axis order you have fallen for a common beginner's problem of assuming that you know what the axis order of a projection is, and in the case of ESPG:4326 that it is fixed.

For example, the following code:

CoordinateReferenceSystem targetCrs = CRS.decode("EPSG:4326");

System.out.println("EPSG:4326 - " + CRS.getAxisOrder(targetCrs));
System.out.println("WGS84 - " + CRS.getAxisOrder(DefaultGeographicCRS.WGS84));

gives this output:

EPSG:4326 - NORTH_EAST
WGS84 - EAST_NORTH

But if you add the line Hints.putSystemDefault(Hints.FORCE_LONGITUDE_FIRST_AXIS_ORDER, Boolean.TRUE); the output becomes:

EPSG:4326 - EAST_NORTH
WGS84 - EAST_NORTH

So, in general you should never rely on a "known" axis order in your code, instead using something like:

double latitude = 29.0;
double longitude = -99.0;

CoordinateReferenceSystem sourceCrs = CRS.decode("EPSG:4312");
CoordinateReferenceSystem targetCrs = CRS.decode("EPSG:4326");

MathTransform engine = CRS.findMathTransform(sourceCrs, targetCrs, false);

DirectPosition2D source;
if (CRS.getAxisOrder(sourceCrs).equals(org.geotools.referencing.CRS.AxisOrder.EAST_NORTH)) {
  source = new DirectPosition2D(sourceCrs, longitude, latitude);
} else {
  source = new DirectPosition2D(sourceCrs, latitude, longitude);
}
DirectPosition target = new DirectPosition2D(targetCrs);
engine.transform(source, target);

if (CRS.getAxisOrder(targetCrs).equals(org.geotools.referencing.CRS.AxisOrder.EAST_NORTH)) {
  System.out.println("lon,lat=" + target.getCoordinate()[0] + "," + target.getCoordinate()[1]);
} else {
  System.out.println("lon,lat=" + target.getCoordinate()[1] + "," + target.getCoordinate()[0]);
}

As for the accuracy of the transform that would depend on which referencing module you imported and which CRS definition Oracle are using. With the gt-epsg-hsql module loaded I see the same ToWGS84[601.705, 84.263, 485.227, 4.7354, -1.3145, -5.393, -2.3887] matrix as the most accurate transform provided in the EPSG registry. If there is an NTv2 transformation available then adding that to your project should improve precision.