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> > 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 is unrelated to EPSG. The only official and authoritative source of EPSG definitions is http://epsg-registry.org/ (or alternatively the databases downloaded from http://www.epsg.org). 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 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> > 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
