4

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");
    }

}
  • Check you switched the coordinates when you change from one to the other – Ian Turton Nov 29 '16 at 18:16
  • I checked it many times before posting here because I could not believe it. The posted output is from the Java program above (see the two definitions of CoordinateReferenceSystem sourceCRS - I commented out one for each run - the rest of the program stays the same) – evod Nov 29 '16 at 20:18
  • 4
    But one of the CRS is lat, lon and the other is lon, lat. – Ian Turton Nov 29 '16 at 20:22
5

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");
1

Your issue is due to the handling of axis order in the two coordinate systems, while for one of them the lines are in Austria, for the other they are in the Middle East and so get distorted when reprojected to EPSG:31256.

I fixed your example in this manner:

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

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

    Geometry linestringEastWest;
    Geometry linestringNorthSouth;
    for (CoordinateReferenceSystem s : sourceCRS) {
      if (CRS.getAxisOrder(s).equals(AxisOrder.EAST_NORTH)) {
        // two lineStrings in Vienna, both around 200m long
        linestringNorthSouth = reader.read("LINESTRING(16.368831 48.201796,16.368831 48.2)");
        linestringEastWest = reader.read("LINESTRING(16.368831 48.2,16.371526 48.2)");
      } else {
        linestringNorthSouth = reader.read("LINESTRING( 48.201796 16.368831, 48.2 16.368831)");
        linestringEastWest = reader.read("LINESTRING( 48.2 16.368831, 48.2 16.371526)");
      }
      MathTransform transform = CRS.findMathTransform(s, targetCRS);
      System.out.println(s);
      System.out.println("northSouth: " + JTS.transform(linestringNorthSouth, transform).getLength() + "m");
      System.out.println("eastWest: " + JTS.transform(linestringEastWest, transform).getLength() + "m");
    }
  }

Which gives the correct result of:

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: 199.70671652817677m
eastWest: 200.33864623436713m
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
  • This will work, however I think the solution I got from the geotools website (using Hints.FORCE_LONGITUDE_FIRST_AXIS_ORDER) is more feasible, since the WKT may come from an external source. Thanks for pointing out the source of the problem! – evod Dec 1 '16 at 14:47
  • if your WKT comes from an external source nothing will save you as you have no idea of their interpretation of the axis order for geographic points. – Ian Turton Dec 1 '16 at 14:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.