[GIS] Geotools: unexpected difference between CRS.decode(“EPSG:4326”) and DefaultGeographicCRS.WGS84

coordinate systemgeotoolswgs84

I wanted to use geotools (15.3) to reproject data from EPSG:4326 to a local projection. To see if everything works as expected I printed the length of two test linestrings. The expected result is around 200 meters.

However I noticed that the result varies greatly depending on how I define the source CoordinateReferenceSystem. I thought these two CRS should be identical, however they are not:

  1. DefaultGeographicCRS.WGS84 -> works as expected
  2. CRS.decode("EPSG:4326") -> completely wrong results

For DefaultGeographicCRS.WGS84 the defined CRS and the (correct) lengths are:

GEOGCS["WGS84(DD)", 
    DATUM["WGS84", 
        SPHEROID["WGS84", 6378137.0, 298.257223563]], 
    PRIMEM["Greenwich", 0.0], 
    UNIT["degree", 0.017453292519943295], 
    AXIS["Geodetic longitude", EAST], 
    AXIS["Geodetic latitude", NORTH]]
northSouth: 199.70672773686056m
eastWest: 200.33865746495295m

For CRS.decode("EPSG:4326") the defined CRS and the (completely wrong) lengths are:

GEOGCS["WGS 84", 
    DATUM["World Geodetic System 1984", 
        SPHEROID["WGS 84", 6378137.0, 298.257223563, AUTHORITY["EPSG","7030"]], 
        AUTHORITY["EPSG","6326"]], 
    PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]], 
    UNIT["degree", 0.017453292519943295], 
    AXIS["Geodetic latitude", NORTH], 
    AXIS["Geodetic longitude", EAST], 
    AUTHORITY["EPSG","4326"]]
northSouth: 222.76855610197927m
eastWest: 346.2445236833851m

How is this possible or what am I doing wrong? The result is the same with geotools 14.0.

import org.geotools.factory.Hints;
import org.geotools.geometry.jts.JTS;
import org.geotools.geometry.jts.JTSFactoryFinder;
import org.geotools.geometry.jts.WKTReader2;
import org.geotools.referencing.CRS;
import org.geotools.referencing.crs.DefaultGeographicCRS;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
import org.opengis.referencing.operation.MathTransform;

import com.vividsolutions.jts.geom.Geometry;
import com.vividsolutions.jts.geom.GeometryFactory;

public class TransformFailMain {

    public static void main(String[] args) throws Exception {
        CoordinateReferenceSystem sourceCRS = CRS.decode("EPSG:4326");
        // CoordinateReferenceSystem sourceCRS = DefaultGeographicCRS.WGS84;
        CoordinateReferenceSystem targetCRS = CRS.decode("EPSG:31256");

        Hints hints = new Hints(Hints.CRS, sourceCRS);
        GeometryFactory geometryFactory = JTSFactoryFinder.getGeometryFactory(hints);
        WKTReader2 reader = new WKTReader2(geometryFactory);

        // two lineStrings in Vienna, both around 200m long
        Geometry linestringNorthSouth = reader.read("LINESTRING(16.368831 48.201796,16.368831 48.2)");
        Geometry linestringEastWest = reader.read("LINESTRING(16.368831 48.2,16.371526 48.2)");

        MathTransform transform = CRS.findMathTransform(sourceCRS, targetCRS);
        System.out.println(sourceCRS);
        System.out.println("northSouth: " + JTS.transform(linestringNorthSouth, transform).getLength() + "m");
        System.out.println("eastWest: " + JTS.transform(linestringEastWest, transform).getLength() + "m");
    }

}

Best Answer

I stumbled over a common pitfall.

  • for EPSG:4326 the coordinate order should be latitude, longitude
  • the javdoc for DefaultGeographicCRS.WGS84 state that the coordinate order should be longitude latitude

So you should choose wisely which CRS you use when e.g. parsing external WKT strings as I did.

More in-depth explanations of why this is the case are found in https://stackoverflow.com/a/13579921/1648538 and http://docs.geotools.org/latest/userguide/library/referencing/order.html. (Thank you for the hint, iant!)

Good to know: the order of lat/lon can be changed in geotools:

Hints hints = new Hints(Hints.FORCE_LONGITUDE_FIRST_AXIS_ORDER, Boolean.TRUE);
CRSAuthorityFactory factory = ReferencingFactoryFinder.getCRSAuthorityFactory("EPSG", hints);
CoordinateReferenceSystem crs = factory.createCoordinateReferenceSystem("EPSG:4326");
Related Question