Dear Martin, Thank you for your answer!
According to http://epsg.io/1676, epsg.io uses the same parameters for the transformation. As you suggested, I tested the two steps separately and the 4cm difference occurs in the second step: 1. EPSG:4326 to EPSG:4150 From (46.218, 9.584) -> Apache SIS (46.21921455327953, 9.585262954192986) -> epsg.io<http://epsg.io/transform#s_srs=4326&t_srs=4150&x=9.5840000&y=46.2180000> (46.2192145, 9.585263) 1. EPSG:4150 to EPSG:2056 From (46.21921455327953, 9.585262954192986) -> Apache SIS (2765526.380321678, 1120768.5597209763) -> epsg.io<http://epsg.io/transform#s_srs=4150&t_srs=2056&x=9.5852630&y=46.2192146> (2765526.34, 1120768.56) I also performed the same transformation in DotSpatial with the following parameters: Authority null string AuthorityCode 0 int AuxiliarySphereType NotSpecified DotSpatial.Projections.AuxiliarySphereType CentralMeridian 7.4395833333333332 double? FalseEasting 2600000 double? FalseNorthing 1200000 double? Geoc false bool GeographicInfo {DotSpatial.Projections.GeographicInfo} DotSpatial.Projections.GeographicInfo Datum {DotSpatial.Projections.Datum} DotSpatial.Projections.Datum DatumType Param3 DotSpatial.Projections.DatumType Description null string NadGrids null string[] Name null string Proj4DatumName null string Spheroid {DotSpatial.Projections.Spheroid} DotSpatial.Projections.Spheroid Code "BR" string EquatorialRadius 6377397.155 double InverseFlattening 299.15281280000335 double KnownEllipsoid Bessel_1841 DotSpatial.Projections.Proj4Ellipsoid Name "Bessel_1841" string PolarRadius 6356078.9628181886 double Static members Non-Public members ToWGS84 {double[3]} double[] [0] 674.374 double [1] 15.056 double [2] 405.346 double Static members Non-Public members Meridian {DotSpatial.Projections.Meridian} DotSpatial.Projections.Meridian Code 8901 int Longitude 0 double Name "Greenwich" string pm null string Non-Public members Name null string Unit {DotSpatial.Projections.AngularUnit} DotSpatial.Projections.AngularUnit Name "Degree" string Radians 0.017453292519943295 double Non-Public members Non-Public members IsGeocentric false bool IsLatLon false bool IsSouth false bool IsValid true bool Lam0 0.12984522414316146 double Lam1 0.819474068676122 double Lam2 0 double LatitudeOfOrigin 46.952405555555572 double? LongitudeOfCenter null double? M null double? Name null string NoDefs true bool Over false bool Phi0 0.819474068676122 double Phi1 0.819474068676122 double Phi2 0 double ScaleFactor 1 double StandardParallel1 46.952405555555572 double? StandardParallel2 null double? Transform {DotSpatial.Projections.Transforms.SwissObliqueMercator} DotSpatial.Projections.Transforms.ITransform {DotSpatial.Projections.Transforms.SwissObliqueMercator} A 6377397.155 double B 6356078.9628181886 double E 0.08169683122252705 double Es 0.00667437223180207 double FromMeter 1 double IsElliptical true bool K0 1 double Lam0 0.12984522414316146 double Name "Swiss_Oblique_Mercator" string OneEs 0.99332562776819788 double Phi0 0.819474068676122 double Proj4Aliases null string[] Proj4Name "somerc" string ROneEs 1.0067192187991747 double Ra 1.5680378306312334E-07 double ToMeter 1 double X0 2600000 double Y0 1200000 double The results obtained in DotSpatial are similar to the ones from epsg.io, containing the 4cm difference to the Apache SIS results. Could you please give me an advice on how to reduce/eliminate this difference? Thank you! Best regards, Andreea From: Martin Desruisseaux <[email protected]> Sent: Wednesday, 17 October 2018 15:31 To: PLOCON Andreea <[email protected]> Cc: [email protected] Subject: Re: Hotine Oblique Mercator (variant B) Hello Andreea, and welcome! Le 17/10/2018 à 14:08, PLOCON Andreea a écrit : I am currently using the Apache SIS implementation of Hotine Oblique Mercator (variant B) together with the EPSG dataset 9.5.4<epsg-registry.org> to transform coordinates from WGS84 to EPSG:2056. For example, (46.218, 9.584) -> (2765526.380321678, 1120768.5597209763). I compared the transformed coordinates with the results obtained from epsg.io<http://epsg.io/transform#s_srs=4326&t_srs=2056&x=9.5840000&y=46.2180000><https://emea01.safelinks.protection.outlook.com/?url=http%3A%2F%2Fepsg.io%2Ftransform%23s_srs%3D4326%26t_srs%3D2056%26x%3D9.5840000%26y%3D46.2180000&data=02%7C01%7Candreea.plocon%40hexagon.com%7C206ff618d1834d86bf5f08d63434bca9%7C1b16ab3eb8f64fe39f3e2db7fe549f6a%7C0%7C0%7C636753798449057084&sdata=oszXtTj0tAY%2BajadRUvQx7lmCKiXrjB8LogTrJ1EZhw%3D&reserved=0> and there is a difference of about 4 cm in the x coordinate. Before to discuss this difference, a would like to write a note of caution. Maybe you are already aware of this issue, but I would like to mention it for the benefit of the list. Despite its name, http://epsg.io<https://emea01.safelinks.protection.outlook.com/?url=http%3A%2F%2Fepsg.io&data=02%7C01%7Candreea.plocon%40hexagon.com%7C206ff618d1834d86bf5f08d63434bca9%7C1b16ab3eb8f64fe39f3e2db7fe549f6a%7C0%7C0%7C636753798449067089&sdata=OL4FU1nCXokfQW3PASPnShVMC%2FonmSRwQ09cUwPYr3Q%3D&reserved=0> is unrelated to EPSG. The only official and authoritative source of EPSG definitions is http://epsg-registry.org/<https://emea01.safelinks.protection.outlook.com/?url=http%3A%2F%2Fepsg-registry.org%2F&data=02%7C01%7Candreea.plocon%40hexagon.com%7C206ff618d1834d86bf5f08d63434bca9%7C1b16ab3eb8f64fe39f3e2db7fe549f6a%7C0%7C0%7C636753798449067089&sdata=bz3xOtaP3slsZeXRVSqbYtzwwihQtXRgN%2FXJfw3kZG8%3D&reserved=0> (or alternatively the databases downloaded from http://www.epsg.org<https://emea01.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.epsg.org&data=02%7C01%7Candreea.plocon%40hexagon.com%7C206ff618d1834d86bf5f08d63434bca9%7C1b16ab3eb8f64fe39f3e2db7fe549f6a%7C0%7C0%7C636753798449077098&sdata=0ifrqwlra6qNDxXJ3f0U8YG7waHpXVlda6ai2%2BUss5c%3D&reserved=0>). The definitions found on epsg.io differ from authoritative definitions for most geographic CRS (different axis order) and some other cases like Google projection EPSG:3857 (epsg.io declares a wrong projection method, which result in an error of about 30 km at U.K. latitude). I have been told that if EPSG was not a non-profit committee, that may have taken legal action asking epsg.io to change their name (disclaimer: I'm not talking in the name of anyone on EPSG side). The epsg.io site is nevertheless useful for checking which definitions are effectively used by various popular software. The transformation from EPSG:4326 to EPSG:2056 involves a datum shift. There is various methods for applying datum shift, depending on which information is available, with different precision. In order to explain the 4 cm difference observed between Apache SIS and epsg.io, we need to know which transformation methods were used by each library. In Apache SIS case, this information can be inspected by System.out.println(operation) where 'operation' is the result of the call to CRS.findOperation(…). The output is as below (omitting some lines for brevity): ConcatenatedOperation["Geographic2D[“WGS 84”] ⟶ ProjectedCRS[“CH1903+ / LV95”]", SourceCRS[GeodeticCRS["WGS 84", ...], TargetCRS[ProjectedCRS["CH1903+ / LV95", ...], CoordinateOperationStep["Inverse operation", Method["Geocentric translations (geog2D domain)"], Parameter["X-axis translation", -674.374, Unit["metre", 1]], Parameter["Y-axis translation", -15.056, Unit["metre", 1]], Parameter["Z-axis translation", -405.346, Unit["metre", 1]], OperationAccuracy[1.0]], CoordinateOperationStep["Swiss Oblique Mercator 1995", Method["Hotine Oblique Mercator (variant B)"], Parameter["Latitude of projection centre", 46.952405555555565, Unit["degree", 0.017453292519943295]], Parameter["Longitude of projection centre", 7.439583333333332, Unit["degree", 0.017453292519943295]], ...etc...], Area["Liechtenstein; Switzerland."], BBox[45.82, 5.96, 47.81, 10.49]] The key information is that this transformation is executed as a concatenation of two steps: 1. A geocentric translation of (-674.374, -15.056, -405.346) metres. 2. The Oblique Mercator projection. The parameters for step 1 were fetched from the EPSG database, with opposite sign because we apply the transformation in opposite direction (from WGS84 to CH1903+): http://epsg-registry.org/?display=entity&urn=urn:ogc:def:coordinateOperation:EPSG::1676<https://emea01.safelinks.protection.outlook.com/?url=http%3A%2F%2Fepsg-registry.org%2F%3Fdisplay%3Dentity%26urn%3Durn%3Aogc%3Adef%3AcoordinateOperation%3AEPSG%3A%3A1676&data=02%7C01%7Candreea.plocon%40hexagon.com%7C206ff618d1834d86bf5f08d63434bca9%7C1b16ab3eb8f64fe39f3e2db7fe549f6a%7C0%7C0%7C636753798449077098&sdata=sB3MGHzLfKh3nhPbmC71fgvrPJKirUP6Gc%2BD3vXOYVQ%3D&reserved=0> For comparing with epsg.io, we would need to know which parameters they used. I'm not aware of an easy way to do that with Proj.4. Note however that this situation may improve with the ongoing "gdalbarn" effort. In the main time, one possible approach may be to test those two steps separately, by testing the transformation from EPSG:4326 to EPSG:4150 before to test the conversion from EPSG:4150 to EPSG:2056. I looked at the implementation of ObliqueMercator<https://github.com/apache/sis/blob/geoapi-4.0/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ObliqueMercator.java><https://emea01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fapache%2Fsis%2Fblob%2Fgeoapi-4.0%2Fcore%2Fsis-referencing%2Fsrc%2Fmain%2Fjava%2Forg%2Fapache%2Fsis%2Freferencing%2Foperation%2Fprojection%2FObliqueMercator.java&data=02%7C01%7Candreea.plocon%40hexagon.com%7C206ff618d1834d86bf5f08d63434bca9%7C1b16ab3eb8f64fe39f3e2db7fe549f6a%7C0%7C0%7C636753798449087103&sdata=ksuw7i30fytbTcxGRcMehBHTI9ty%2BVCrVhqp98Na50Q%3D&reserved=0> and compared the formulas to the ones mentioned in IOGP Publication 373-7-2 – Geomatics Guidance Note number 7, part 2 – August 2018 (Coordinate Conversions & Transformations including Formulas). Could you please tell me why the methods ObliqueMercator.transform and ObliqueMercator.inverseTransform do not consider the differences between variant A and variant B described in the mentioned document (see Page 60-61 Forward case – computation of u and Reverse case – computation of v’ and u’). The equations given by EPSG guidance notes involve two kinds of variables: * Variables that are fully determined by the projection parameters and do not depend on the latitude and longitude of the points to project. * Variables that depend, directly or indirectly, on the latitude and longitude of the point to project. The first category of variables can be computed once for ever in the map projection constructor. For performance reason we try to keep in the transform and inverseTransform methods only the variables that depend, directly or indirectly, on the values of the coordinates to project. When doing this separation, we see that all the logic that depends on whether we are using variant A or variant B appears in the constructor. The same phenomenon is observed with most EPSG projections we have implemented so far (Mercator, Lambert, etc.): those variants are just different ways to specify the parameters for initializing the projection. After that point, the "kernel" is the same for all variants. Could this cause the differences between the coordinates transformed with Apache SIS and the ones obtained from espg.io which I mentioned before? I think that difference is datum shift methods is more likely. This can be verified by testing the two steps separately, as described above. Regards, Martin
