[GIS] Unable to transform from EPSG:3857 to UTM in Geotools

coordinate systemgeotools

I'm using Geotools to perform some transformations between coordinate reference systems. I created a proof of concept taking data from Google Earth, which uses ESPG:4326, and converting it to UTM.

The first problem I encountered was identifying the UTM tile to work on. I think this gets solved using the following reference:
CRS.decode(String.format("AUTO2:42001,%s,%s", longitude, latitude ), true). Where longitud and latitude are the coordinates of the area I'm working on.

I was able to apply transformations between these two systems:

CoordinateReferenceSystem wgs84= CRS.decode("EPSG:4326", true);
CoordinateReferenceSystem utm=CRS.decode(String.format("AUTO2:42001,%s,%s", -3.691406,  40.403431 ), true);
MathTransform toMeters= CRS.findMathTransform(wgs84, utm);
MathTransform toDegrees= CRS.findMathTransform(utm, wgs84);

com.vividsolutions.jts.geom.GeometryFactory jtsGf= JTSFactoryFinder.getGeometryFactory();

Coordinate[] coordinates = new Coordinate[5];
coordinates[0]=new Coordinate(-3.679182821543855,40.43597099928637);
coordinates[1]=new Coordinate(-3.677717792288855,40.43591469226429);
coordinates[2]=new Coordinate(-3.677605067311292,40.43726977694883);
coordinates[3]=new Coordinate(-3.67914084885543,40.4373645084199);
coordinates[4]=new Coordinate(-3.679182821543855,40.43597099928637);

Geometry mRank = JTS.transform(jtsGf.createPolygon(jtsGf.createLinearRing(coordinates), null), toMeters);

When I try to run the same code using ESPG:3857 instead of ESPG:4362 I get the following error when trying to apply the transformation:

org.geotools.referencing.operation.projection.ProjectionException: The transform result may be 0 meters away from the expected position. Are you sure that the input coordinates are inside this map projection area of validity? The point is located 3°00.0'E away from the central meridian and 0°00.0'N away from the latitude of origin. The projection is "Transverse_Mercator".
at org.geotools.referencing.operation.projection.MapProjection.checkReciprocal(MapProjection.java:708)
at org.geotools.referencing.operation.projection.MapProjection.transform(MapProjection.java:903)
at org.geotools.referencing.operation.projection.MapProjection.transform(MapProjection.java:938)
at org.geotools.referencing.operation.transform.ConcatenatedTransformDirect.transform(ConcatenatedTransformDirect.java:81)
at org.geotools.geometry.jts.DefaultCoordinateSequenceTransformer.transform(DefaultCoordinateSequenceTransformer.java:123)
at org.geotools.geometry.jts.GeometryCoordinateSequenceTransformer.projectCoordinateSequence(GeometryCoordinateSequenceTransformer.java:295)
at org.geotools.geometry.jts.GeometryCoordinateSequenceTransformer.transformStraightLineString(GeometryCoordinateSequenceTransformer.java:237)
at org.geotools.geometry.jts.GeometryCoordinateSequenceTransformer.transformLineString(GeometryCoordinateSequenceTransformer.java:216)
at org.geotools.geometry.jts.GeometryCoordinateSequenceTransformer.transformPolygon(GeometryCoordinateSequenceTransformer.java:304)
at org.geotools.geometry.jts.GeometryCoordinateSequenceTransformer.transform(GeometryCoordinateSequenceTransformer.java:170)
at org.geotools.geometry.jts.JTS.transform(JTS.java:442)

The most weird part is ... transform result may be 0 meters away from the expected position..., which seems the perfect situation for a transformation.

Any clues on what is wrong?

EDIT: added failing code, stack trace, dependencies, UTM projection and transformation dumps

The failing code snippet

        CoordinateReferenceSystem wgs84 = CRS.decode("EPSG:3857 ", true);
    //systema utm para las coordenadas de madrid
    CoordinateReferenceSystem utm = CRS.decode(String.format("AUTO2:42001,%s,%s", -3.691406, 40.403431), true);
    MathTransform toMeters = CRS.findMathTransform(wgs84, utm);
    MathTransform toDegrees = CRS.findMathTransform(utm, wgs84);

    com.vividsolutions.jts.geom.GeometryFactory jtsGf = JTSFactoryFinder.getGeometryFactory();

    Coordinate[] coordinates = new Coordinate[5];
    coordinates[0] = new Coordinate(-3.679182821543855, 40.43597099928637);
    coordinates[1] = new Coordinate(-3.677717792288855, 40.43591469226429);
    coordinates[2] = new Coordinate(-3.677605067311292, 40.43726977694883);
    coordinates[3] = new Coordinate(-3.67914084885543, 40.4373645084199);
    coordinates[4] = new Coordinate(-3.679182821543855, 40.43597099928637);

    Geometry mRank = JTS.transform(jtsGf.createPolygon(jtsGf.createLinearRing(coordinates), null), toMeters);

The error stack trace:

org.geotools.referencing.operation.projection.ProjectionException: The transform result may be 0 meters away from the expected position. Are you sure that the input coordinates are inside this map projection area of validity? The point is located 3°00.0'E away from the central meridian and 0°00.0'N away from the latitude of origin. The projection is "Transverse_Mercator".

at org.geotools.referencing.operation.projection.MapProjection.checkReciprocal(MapProjection.java:708)
at org.geotools.referencing.operation.projection.MapProjection.transform(MapProjection.java:903)
at org.geotools.referencing.operation.projection.MapProjection.transform(MapProjection.java:938)
at org.geotools.referencing.operation.transform.ConcatenatedTransformDirect.transform(ConcatenatedTransformDirect.java:81)
at org.geotools.geometry.jts.DefaultCoordinateSequenceTransformer.transform(DefaultCoordinateSequenceTransformer.java:123)
at org.geotools.geometry.jts.GeometryCoordinateSequenceTransformer.projectCoordinateSequence(GeometryCoordinateSequenceTransformer.java:295)
at org.geotools.geometry.jts.GeometryCoordinateSequenceTransformer.transformStraightLineString(GeometryCoordinateSequenceTransformer.java:237)
at org.geotools.geometry.jts.GeometryCoordinateSequenceTransformer.transformLineString(GeometryCoordinateSequenceTransformer.java:216)
at org.geotools.geometry.jts.GeometryCoordinateSequenceTransformer.transformPolygon(GeometryCoordinateSequenceTransformer.java:304)
at org.geotools.geometry.jts.GeometryCoordinateSequenceTransformer.transform(GeometryCoordinateSequenceTransformer.java:170)
at org.geotools.geometry.jts.JTS.transform(JTS.java:442)
at com.byteflair.mercury.api.checkin.CheckinServiceTest.checkThatICanTransformAndBufferPolygon(CheckinServiceTest.java:236)
at sun.reflect.NativeMethodAccessorImpl.invoke0(Native Method)
at sun.reflect.NativeMethodAccessorImpl.invoke(NativeMethodAccessorImpl.java:62)
at sun.reflect.DelegatingMethodAccessorImpl.invoke(DelegatingMethodAccessorImpl.java:43)
at java.lang.reflect.Method.invoke(Method.java:483)
at org.junit.runners.model.FrameworkMethod$1.runReflectiveCall(FrameworkMethod.java:50)
at org.junit.internal.runners.model.ReflectiveCallable.run(ReflectiveCallable.java:12)
at org.junit.runners.model.FrameworkMethod.invokeExplosively(FrameworkMethod.java:47)
at org.junit.internal.runners.statements.InvokeMethod.evaluate(InvokeMethod.java:17)
at org.junit.runners.ParentRunner.runLeaf(ParentRunner.java:325)
at org.junit.runners.BlockJUnit4ClassRunner.runChild(BlockJUnit4ClassRunner.java:78)
at org.junit.runners.BlockJUnit4ClassRunner.runChild(BlockJUnit4ClassRunner.java:57)
at org.junit.runners.ParentRunner$3.run(ParentRunner.java:290)
at org.junit.runners.ParentRunner$1.schedule(ParentRunner.java:71)
at org.junit.runners.ParentRunner.runChildren(ParentRunner.java:288)
at org.junit.runners.ParentRunner.access$000(ParentRunner.java:58)
at org.junit.runners.ParentRunner$2.evaluate(ParentRunner.java:268)
at org.junit.internal.runners.statements.RunBefores.evaluate(RunBefores.java:26)
at org.junit.runners.ParentRunner.run(ParentRunner.java:363)
at org.junit.runner.JUnitCore.run(JUnitCore.java:137)
at com.intellij.junit4.JUnit4IdeaTestRunner.startRunnerWithArgs(JUnit4IdeaTestRunner.java:117)
at com.intellij.rt.execution.junit.JUnitStarter.prepareStreamsAndStart(JUnitStarter.java:234)
at com.intellij.rt.execution.junit.JUnitStarter.main(JUnitStarter.java:74)

My dependencies based on mvn dependency:list | grep geotools:

[INFO]    org.geotools:gt-referencing:jar:15-SNAPSHOT:compile
[INFO]    org.geotools:gt-main:jar:15-SNAPSHOT:compile
[INFO]    org.geotools:gt-api:jar:15-SNAPSHOT:compile
[INFO]    org.geotools:gt-metadata:jar:15-SNAPSHOT:compile
[INFO]    org.geotools:gt-epsg-hsql:jar:15-SNAPSHOT:compile
[INFO]    org.geotools:gt-opengis:jar:15-SNAPSHOT:compile

Also tried with version 14.2 with same results

The wgs84 projection:

PROJCS["WGS 84 / Pseudo-Mercator", 
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 longitude", EAST], 
AXIS["Geodetic latitude", NORTH], 
AUTHORITY["EPSG","4326"]], 
PROJECTION["Popular Visualisation Pseudo Mercator",AUTHORITY["EPSG","1024"]], 
PARAMETER["semi_minor", 6378137.0], 
PARAMETER["latitude_of_origin", 0.0], 
PARAMETER["central_meridian", 0.0], 
PARAMETER["scale_factor", 1.0], 
PARAMETER["false_easting", 0.0], 
PARAMETER["false_northing", 0.0], 
UNIT["m", 1.0], 
AXIS["Easting", EAST], 
AXIS["Northing", NORTH], 
AUTHORITY["EPSG","3857"]]

The utm projection:

PROJCS["WGS 84 / Auto UTM", 
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]], 
PROJECTION["Transverse_Mercator"], 
PARAMETER["central_meridian", -3.0], 
PARAMETER["latitude_of_origin", 0.0], 
PARAMETER["scale_factor", 0.9996], 
PARAMETER["false_easting", 500000.0], 
PARAMETER["false_northing", 0.0], 
UNIT["m", 1.0], 
AXIS["Easting", EAST], 
AXIS["Northing", NORTH]]

The toMeters transformation:

CONCAT_MT[INVERSE_MT[PARAM_MT["Transverse_Mercator", 
PARAMETER["semi_major", 6378137.0], 
PARAMETER["semi_minor", 6356752.314245179], 
PARAMETER["central_meridian", -3.0], 
PARAMETER["latitude_of_origin", 0.0], 
PARAMETER["scale_factor", 0.9996], 
PARAMETER["false_easting", 500000.0], 
PARAMETER["false_northing", 0.0]]], 
PARAM_MT["Popular Visualisation Pseudo Mercator", 
PARAMETER["semi_major", 6378137.0], 
PARAMETER["semi_minor", 6378137.0], 
PARAMETER["latitude_of_origin", 0.0], 
PARAMETER["central_meridian", 0.0], 
PARAMETER["scale_factor", 1.0], 
PARAMETER["false_easting", 0.0], 
PARAMETER["false_northing", 0.0]]]

Best Answer

I am unable to recreate this error, using the following code, which is basically coppied and pasted from the question:

CoordinateReferenceSystem wgs84= CRS.decode("EPSG:4326", true);
CoordinateReferenceSystem google = CRS.decode("EPSG:3857", true);
CoordinateReferenceSystem utm=CRS.decode(String.format("AUTO2:42001,%s,%s", -3.691406,  40.403431 ), true);
MathTransform toMeters= CRS.findMathTransform(wgs84, utm);
MathTransform toDegrees= CRS.findMathTransform(utm, wgs84);
MathTransform utmToGoogle = CRS.findMathTransform(utm, google);

com.vividsolutions.jts.geom.GeometryFactory jtsGf= JTSFactoryFinder.getGeometryFactory();

Coordinate[] coordinates = new Coordinate[5];
coordinates[0]=new Coordinate(-3.679182821543855,40.43597099928637);
coordinates[1]=new Coordinate(-3.677717792288855,40.43591469226429);
coordinates[2]=new Coordinate(-3.677605067311292,40.43726977694883);
coordinates[3]=new Coordinate(-3.67914084885543,40.4373645084199);
coordinates[4]=new Coordinate(-3.679182821543855,40.43597099928637);

Geometry mRank = JTS.transform(jtsGf.createPolygon(jtsGf.createLinearRing(coordinates), null), utmToGoogle);
System.out.println("google"+mRank);
mRank = JTS.transform(jtsGf.createPolygon(jtsGf.createLinearRing(coordinates), null), toDegrees);
System.out.println("wgs84:"+mRank);

I get the following answer and no errors (and they both line up if I paste them into QGIS):

google POLYGON ((-833646.8252261352 40.59902002387106, -833646.8237650455 40.59896348970921, -833646.8236526242 40.600324038505704, -833646.825184276 40.60041915196819...
wgs84:POLYGON ((-7.488776846583527 0.0003647072020764, -7.488776833458335 0.0003647066942288, -7.488776832448438 0.0003647189162468, -7.488776846207499 0.0003647197706586,...

Update 1 mvn dependency:list | grep geotools gives:

[INFO]    org.geotools:gt-api:jar:15-SNAPSHOT:compile
[INFO]    org.geotools:gt-coverage:jar:15-SNAPSHOT:compile
[INFO]    org.geotools:gt-cql:jar:15-SNAPSHOT:compile
[INFO]    org.geotools:gt-data:jar:15-SNAPSHOT:compile
[INFO]    org.geotools:gt-epsg-hsql:jar:15-SNAPSHOT:compile
[INFO]    org.geotools:gt-epsg-wkt:jar:15-SNAPSHOT:compile
[INFO]    org.geotools:gt-grid:jar:15-SNAPSHOT:compile
[INFO]    org.geotools:gt-main:jar:15-SNAPSHOT:compile
[INFO]    org.geotools:gt-metadata:jar:15-SNAPSHOT:compile
[INFO]    org.geotools:gt-opengis:jar:15-SNAPSHOT:compile
[INFO]    org.geotools:gt-referencing:jar:15-SNAPSHOT:compile
[INFO]    org.geotools:gt-render:jar:15-SNAPSHOT:compile
[INFO]    org.geotools:gt-shapefile:jar:15-SNAPSHOT:compile
[INFO]    org.geotools:gt-swing:jar:15-SNAPSHOT:compile

Update2

I modified the code to try more forward and backward conversions:

MathTransform toMeters= CRS.findMathTransform(wgs84, utm);
MathTransform toDegrees= CRS.findMathTransform(utm, wgs84);
MathTransform utmToGoogle = CRS.findMathTransform(utm, google);
MathTransform wgsToGoogle = CRS.findMathTransform(wgs84, google);
MathTransform googleToUTM = CRS.findMathTransform(google,utm);
 //[.....create poly as before ...]
Polygon wgsPoly = jtsGf.createPolygon(jtsGf.createLinearRing(coordinates), null);
System.out.println("wgs:"+wgsPoly);

Geometry utmPoly = JTS.transform(wgsPoly, toMeters);
System.out.println("utm:"+utmPoly);
Geometry googlePoly1 = JTS.transform(utmPoly, utmToGoogle);
System.out.println("google1:"+googlePoly1);
Geometry googlePoly2 = JTS.transform(wgsPoly, wgsToGoogle);
System.out.println("google2:"+googlePoly2);
Geometry utmPoly2 = JTS.transform(googlePoly1, googleToUTM);
System.out.println("utm2"+utmPoly2);

this gives the following:

wgs:POLYGON ((-3.679182821543855 40.43597099928637, -3.677717792288855 40.43591469226429, -3.677605067311292 40.43726977694883...
utm:POLYGON ((442395.32424547715 4476369.013946666, 442519.53321862826 4476361.809516019, 442530.24786807405 4476512.15038715...
google1:POLYGON ((-409564.75822962367 4929500.130541016, -409401.67191895976 4929491.895343309, -409389.1234318578 4929690.085509476,...
google2:POLYGON ((-409564.75822962145 4929500.130541016, -409401.67191895755 4929491.895343313, -409389.1234318556 4929690.085509479...
utm2POLYGON ((442395.32424547547 4476369.013946667, 442519.5332186266 4476361.809516017, 442530.24786807236 4476512.150387148,

which seem to agree to within expected tolerances.

Best Answer

Related Solutions

[GIS] QGIS 1.8.0 – Define shapefile as UTM 15 NAD83 Feet

[GIS] ArcGIS Desktop, shapefile and EPSG code

Related Question