3

I am using GeoTools to transform geometries from Krovak projection (EPSG:5514) to WGS84. I am following GeoTools library documentation - http://docs.geotools.org/stable/userguide/library/referencing/crs.html but it fails on "No transform for classification "Krovak (North Orientated)"" error which I am unable to resolve.

String wkbRep = "010100000000000000f6a826c10000000090432dc1"; //point POINT (-742523 -958920) for testing in WKB
byte[] wkb = WKBReader.hexToBytes(wkbRep);
WKBReader wkbReader = new WKBReader();
Geometry sourceGeometry = wkbReader.read(wkb);
System.out.println(sourceGeometry.toString());

CoordinateReferenceSystem sourceCRS = CRS.decode("EPSG:5514");
CoordinateReferenceSystem targetCRS = CRS.decode("EPSG:4326");
MathTransform transform = CRS.findMathTransform(sourceCRS, targetCRS, true);
Geometry targetGeometry = JTS.transform(sourceGeometry, transform);
System.out.println(targetGeometry.toString());

This generates error:

Exception in thread "main" org.opengis.referencing.NoSuchIdentifierException: No transform for classification "Krovak (North Orientated)".
at org.geotools.referencing.operation.DefaultMathTransformFactory.lambda$getProvider$1(DefaultMathTransformFactory.java:267)
at java.util.Optional.orElseThrow(Optional.java:290)
at org.geotools.referencing.operation.DefaultMathTransformFactory.getProvider(DefaultMathTransformFactory.java:264)
at org.geotools.referencing.operation.DefaultMathTransformFactory.getDefaultParameters(DefaultMathTransformFactory.java:297)
at org.geotools.referencing.factory.epsg.DirectEpsgFactory.createCoordinateOperation(DirectEpsgFactory.java:2880)
at org.geotools.referencing.factory.BufferedAuthorityFactory.createCoordinateOperation(BufferedAuthorityFactory.java:968)
at org.geotools.referencing.factory.epsg.DirectEpsgFactory.createCoordinateReferenceSystem(DirectEpsgFactory.java:2149)
at org.geotools.referencing.factory.BufferedAuthorityFactory.createCoordinateReferenceSystem(BufferedAuthorityFactory.java:731)
at org.geotools.referencing.factory.AuthorityFactoryAdapter.createCoordinateReferenceSystem(AuthorityFactoryAdapter.java:779)
at org.geotools.referencing.factory.FallbackAuthorityFactory.createCoordinateReferenceSystem(FallbackAuthorityFactory.java:624)
at org.geotools.referencing.factory.AuthorityFactoryAdapter.createCoordinateReferenceSystem(AuthorityFactoryAdapter.java:779)
at org.geotools.referencing.factory.ThreadedAuthorityFactory.createCoordinateReferenceSystem(ThreadedAuthorityFactory.java:635)
at org.geotools.referencing.DefaultAuthorityFactory.createCoordinateReferenceSystem(DefaultAuthorityFactory.java:176)
at org.geotools.referencing.CRS.decode(CRS.java:533)
at org.geotools.referencing.CRS.decode(CRS.java:447)
at cz.tmobile.bigdata.internal.geotools.Transform.main(Transform.java:24)

I also tried to define the EPSG:5514 by WKT from EPSG.io:

"PROJCS[\"S-JTSK / Krovak East North\",GEOGCS[\"S-JTSK\",DATUM[\"System_Jednotne_Trigonometricke_Site_Katastralni\",SPHEROID[\"Bessel 1841\",6377397.155,299.1528128,AUTHORITY[\"EPSG\",\"7004\"]],TOWGS84[485.0,169.5,483.8,7.786,4.398,4.103,0],AUTHORITY[\"EPSG\",\"6156\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4156\"]],PROJECTION[\"Krovak\"],PARAMETER[\"latitude_of_center\",49.5],PARAMETER[\"longitude_of_center\",24.83333333333333],PARAMETER[\"azimuth\",30.28813972222222],PARAMETER[\"pseudo_standard_parallel_1\",78.5],PARAMETER[\"scale_factor\",0.9999],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"X\",EAST],AXIS[\"Y\",NORTH],AUTHORITY[\"EPSG\",\"5514\"]]"

And changing the code to:

String krovakWKT = "PROJCS[\"S-JTSK / Krovak East North\",GEOGCS[\"S-JTSK\",DATUM[\"System_Jednotne_Trigonometricke_Site_Katastralni\",SPHEROID[\"Bessel 1841\",6377397.155,299.1528128,AUTHORITY[\"EPSG\",\"7004\"]],TOWGS84[485.0,169.5,483.8,7.786,4.398,4.103,0],AUTHORITY[\"EPSG\",\"6156\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4156\"]],PROJECTION[\"Krovak\"],PARAMETER[\"latitude_of_center\",49.5],PARAMETER[\"longitude_of_center\",24.83333333333333],PARAMETER[\"azimuth\",30.28813972222222],PARAMETER[\"pseudo_standard_parallel_1\",78.5],PARAMETER[\"scale_factor\",0.9999],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"X\",EAST],AXIS[\"Y\",NORTH],AUTHORITY[\"EPSG\",\"5514\"]]";
CoordinateReferenceSystem sourceCRS = CRS.parseWKT(krovakWKT);
//rest is the same

This sadly produces point which is 215 m offset and with switched coordinates.

Edit 1: I uploaded printscreen from QGIS to show the offset -

QGIS printscreen

The expected WGS84 coordinates are 50.836672, 14.262571; given by the online transformation tool.

Edit 2: I found the way how to deal with the reverse order of returned coordinates here. I just added following line to the main method:

    System.setProperty("org.geotools.referencing.forceXY", "true");
  • Did you try any other Krovak projections available at epsg.io site? – TomazicM May 3 at 14:00
  • Thank you for commenting! Yes I tried several but without any luck. Correct me if I am wrong but there is no alternative to the EPSG5514, right? I am trying to transform data from geoportal.cuzk.cz/(S(odzarxqqjegl10qv2or3mh1r))/… – Petr Švec May 3 at 14:28
  • The 215 m offset sounds like a missing geographic/datum transformation. See EPSG:1622. – mkennedy May 3 at 16:20
  • @mkennedy, thank you for the suggestion. I tried the EPSG:1622 WKT but the output offset is much higher. I believe it is caused by different covered areas - only Czech Republic in case of EPSG:1622 as opposed to Czech Republic and Slovakia in case of EPSG:5514. – Petr Švec May 6 at 9:10
1

It appears that the version of the EPSG database that GeoTools is using doesn't recognise the EPSG:5514 code you are using. Currently, GeoTools is using v8.6 and there is an outstanding ticket to update it but no one has sponsored the work and no volunteer has found time to do it yet.

You don't specify what your expected answer is so I can't check but by adding

55514=PROJCS["S-JTSK / Krovak East North",GEOGCS["S-JTSK",DATUM["System_Jednotne_Trigonometricke_Site_Katastralni",SPHEROID["Bessel 1841",6377397.155,299.1528128,AUTHORITY["EPSG","7004"]],TOWGS84[572.213,85.334,461.94,-4.9732,-1.529,-5.2484,3.5378],AUTHORITY["EPSG","6156"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4156"]],PROJECTION["Krovak"],PARAMETER["latitude_of_center",49.5],PARAMETER["longitude_of_center",24.83333333333333],PARAMETER["azimuth",30.28813972222222],PARAMETER["pseudo_standard_parallel_1",78.5],PARAMETER["scale_factor",0.9999],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],AUTHORITY["EPSG","5514"]]

to my geotools.props file (in my resources folder). Note there are 3 variants of the projection given on epsg.io, with accuracies between 1 and 6m given.

I get the following answer to this code:

PROJCS["S-JTSK / Krovak East North", 
  GEOGCS["S-JTSK", 
    DATUM["System_Jednotne_Trigonometricke_Site_Katastralni", 
      SPHEROID["Bessel 1841", 6377397.155, 299.1528128, AUTHORITY["EPSG","7004"]], 
      TOWGS84[572.213, 85.334, 461.94, -4.9732, -1.529, -5.2484, 3.5378], 
      AUTHORITY["EPSG","6156"]], 
    PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]], 
    UNIT["degree", 0.017453292519943295], 
    AXIS["Longitude", EAST], 
    AXIS["Latitude", NORTH], 
    AUTHORITY["EPSG","4156"]], 
  PROJECTION["Krovak"], 
  PARAMETER["latitude_of_center", 49.50000000000001], 
  PARAMETER["longitude_of_center", 24.833333333333332], 
  PARAMETER["azimuth", 30.288139722222223], 
  PARAMETER["pseudo_standard_parallel_1", 78.5], 
  PARAMETER["scale_factor", 0.9999], 
  PARAMETER["false_easting", 0.0], 
  PARAMETER["false_northing", 0.0], 
  UNIT["m", 1.0], 
  AXIS["x", EAST], 
  AXIS["y", NORTH], 
  AUTHORITY["EPSG","5514"]]
POINT (-742523 -958920)
POINT (50.83599142455958 14.26583534258449)

code:

package com.ianturton.cookbook.projections;

import org.geotools.geometry.jts.JTS;
import org.geotools.referencing.CRS;
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.GeometryFactory;
import org.locationtech.jts.geom.Point;
import org.locationtech.jts.io.WKTReader;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
import org.opengis.referencing.operation.MathTransform;

public class ReprojectTest {

  public static void main(String[] args) throws Exception {
    GeometryFactory geometryFactory = new GeometryFactory();
    CoordinateReferenceSystem sourceCRS = CRS.decode("geotools:55514");
    System.out.println(sourceCRS);
    CoordinateReferenceSystem targetCRS = CRS.decode("EPSG:4326");
    WKTReader reader = new WKTReader(geometryFactory );
    Point point = (Point) reader.read("POINT (-742523 -958920)");

    System.out.println(point.toString());

    MathTransform mTrans = CRS.findMathTransform(sourceCRS, targetCRS);
    Geometry transformed = JTS.transform(point, mTrans);
    System.out.println(transformed.toText());
  }

}

Update

On further investigation I'm more confused than before:

a map of somewhere in the Czech Republic

The brown dot is your expected value, the purple one next to it is from gdaltransform (good) the pink one to the South East is GeoTools with the same parameters as the gdal one (3 param TOWGS84 589,76,480,0,0,0,0), this matches exactly the default 5514 answer, while the orange one to the north is the ESRI proj (with no bursa wolf params). Using the more accurate 7 param transforms with GeoTools leads to a small cluster around the SE dot.

So it is possible that the GeoTools Krovak projection code has a bug or possibly doesn't know how to handle the axis order flip used in EPSG:5514, but I'm a bit out of my depth now. I would raise an issue with some well known test points.

2nd Update

Working on a hunch that it might be a Bursa Wolf issue (most of my projection issues are) I tried removing them from the proj string (and setting lenient to true in transform) gives me an answer that matches yours (almost) of POINT (50.83667109256844 14.262570143324908).

5514-GDALNBW= PROJCS["S-JTSK / Krovak East North",GEOGCS["S-JTSK",DATUM["System_Jednotne_Trigonometricke_Site_Katastralni",SPHEROID["Bessel 1841",6377397.155,299.1528128,AUTHORITY["EPSG","7004"]]], PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4156"]],PROJECTION["Krovak"],PARAMETER["latitude_of_center",49.5],PARAMETER["longitude_of_center",24.83333333333333],PARAMETER["azimuth",30.28813972222222],PARAMETER["pseudo_standard_parallel_1",78.5],PARAMETER["scale_factor",0.9999],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],AUTHORITY["EPSG","5514"]] 

enter image description here

Green point is the new GeoTools one, purple is gdaltransform and your expected is mostly covered by the green one.

With any luck this will help whoever tries to fix it what is going on. Do you have any building outlines or something we can try reprojecting to test to see who is right?

  • Hello Ian, thank you for the answer, I tried your code. I was unable to find the geotools.props file so I used parseWKT method and got the same output point as you did. Sadly, the issue remains, the point is off by cca 250 m. – Petr Švec May 6 at 8:43
  • I uploaded printscreen from QGIS - imgur.com/xcMroa4 . The expected WGS84 coordinates are 50.836672, 14.262571, given by the transformation tool from martin.hinner.info/geo . – Petr Švec May 6 at 8:53
  • The second Update made it! I tested the WKT without the said transformation on couple of polygons and it behaves as expected. Offset is no more than a meter. Just one more question, how do I reverse the coordinates given by the geotools? Instead of WKT point "POINT (50.83667109256844 14.262570143324908)" how do I get "POINT (14.26257014332490850.83667109256844)"? Is it something in CRS WKT definition? – Petr Švec May 6 at 12:24
  • you can add System.setProperty("org.geotools.referencing.forceXY", "true"); see docs.geotools.org/latest/userguide/library/referencing/… – Ian Turton May 6 at 15:41

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.