Reprojecting British National Grid (27700) coords to Mercator (3857) using Proj4j about proj4j

[gmh04]

Note I also raised this on stackexchange: https://gis.stackexchange.com/questions/326430/reprojecting-british-national-grid-27700-coords-to-mercator-3857-using-proj4

[kervel]

same issue when converting from belgian lambert72 (31370) to 3857. here the error is 50 meters.
i worked around the issue by doing a double conversion: when doing 31370-->3857 it is inaccurate but when i do 31370-->wgs84-->3857 its much better.

also see this: https://gis.stackexchange.com/questions/221552/proj4j-not-precise-for-epsg3857-transformations if i use geotrellis proj4j with the coordinates mentioned there i get the same results as the guy asking the question.

[onlymate]

Hello everyone,

I investigated the same problem.

Description

I've tried to transform points from EPSG:2056 to EPSG:3857.
The point in EPSG:3857 had an offset to the original position.
A workaround for me was to transform EPSG:2056 to EPSG:4326 and then to EPSG:3857.

Results

This is the result of transforming EPSG:2056 -> EPSG:3857
http://epsg.io/map#srs=3857&x=829151.230272&y=5933820.784235&z=15&layer=streets

This is the result of transforming EPSG:4326 -> EPSG:3857
http://epsg.io/map#srs=3857&x=829045.228513&y=5933605.151438&z=15&layer=streets

This is the result of transforming EPSG:2056 -> EPSG:4326 -> EPSG:3857
http://epsg.io/map#srs=3857&x=829045.229627&y=5933605.149807&z=15&layer=streets

Code Example

import org.locationtech.proj4j.BasicCoordinateTransform;
import org.locationtech.proj4j.CRSFactory;
import org.locationtech.proj4j.CoordinateReferenceSystem;
import org.locationtech.proj4j.CoordinateTransform;
import org.locationtech.proj4j.ProjCoordinate;

public class App 
{
	/**
	 * This example shows the investigated bug
	 * 
	 * If I transform a point from EPSG:2056 -> EPSG:3857 there is a little offset in the result
	 * If I transform a point from EPSG:4326 -> EPSG:3857 it's fine
	 * If I transform a point from EPSG:2056 -> EPSG:4326 -> EPSG:3857 the offset is minimal
	 */
    public static void main( String[] args )
    {
    	CRSFactory crsFactory = new CRSFactory();
        CoordinateReferenceSystem epsg2056 = crsFactory.createFromName("EPSG:2056");
        CoordinateReferenceSystem epsg4326 = crsFactory.createFromName("EPSG:4326");
        CoordinateReferenceSystem epsg3857 = crsFactory.createFromName("EPSG:3857");
        
        CoordinateTransform transform2056To3857 = new BasicCoordinateTransform(epsg2056, epsg3857);
        CoordinateTransform transform2056To4326 = new BasicCoordinateTransform(epsg2056, epsg4326);
        CoordinateTransform transform4326To3857 = new BasicCoordinateTransform(epsg4326, epsg3857);
        
        // pSrcBern2056 is the coordinate where I investigated the bug (Bern, Switzerland)
        // http://epsg.io/map#srs=2056&x=2600670.52&y=1199667.32&z=15&layer=streets
        ProjCoordinate pSrcBern2056 = new ProjCoordinate(2600670.52, 1199667.32);
        
        // pSrcBern4326 is a simple reference coordinate (same position in Bern, Switzerland)
        // https://epsg.io/map#srs=4326&x=7.447440&y=46.948090&z=15&layer=streets
        ProjCoordinate pSrcBern4326 = new ProjCoordinate(7.44744,46.94809);
        
        ProjCoordinate pDstBern = transform2056To3857.transform(pSrcBern2056, new ProjCoordinate());
        ProjCoordinate pDstBern2 = transform4326To3857.transform(pSrcBern4326, new ProjCoordinate());
        
        ProjCoordinate p2056To4326 = transform2056To4326.transform(pSrcBern2056, new ProjCoordinate());
        ProjCoordinate pDstBern3 = transform4326To3857.transform(p2056To4326, new ProjCoordinate());
        
        // This is the result EPSG:2056 -> EPSG:3857 (BUG!)
        // http://epsg.io/map#srs=3857&x=829151.230272&y=5933820.784235&z=15&layer=streets
        System.out.println(pDstBern.toString()); // ProjCoordinate[829151.2300676685 5933820.784982888 NaN]

        // This is the result EPSG:4326 -> EPSG:3857
        // http://epsg.io/map#srs=3857&x=829045.228513&y=5933605.151438&z=15&layer=streets
        System.out.println(pDstBern2.toString()); // ProjCoordinate[829045.2285134573 5933605.151437666 NaN]
        
        // This is the result EPSG:2056 -> EPSG:4326 -> EPSG:3857
        // http://epsg.io/map#srs=3857&x=829045.229627&y=5933605.149807&z=15&layer=streets
        System.out.println(pDstBern3.toString()); // ProjCoordinate[829045.2301110257 5933605.150096844 NaN]
    }
}

[bosborn]

I looked into this a little bit. At this place, EPSG:3857 has an unknown datum causing the transformation to short cut. I tried some changes that allowed the method to continue if either datum was TYPE_3PARAM or TYPE_7PARAM. This resulted in correct longitude values but still incorrect latitude values. Tried comparing to proj4 and proj4js but couldn't figure out how to fix in the algorithm itself. proj4js is getting around this issue by transforming to WGS84 first.

[bosborn]

Fixed by pull request #75 #76 #77

For reference and documenting attempting to solve this through fixing the algorithm....

When changing the logic here to continue when the source or target datum transform type is TYPE_3PARAM or TYPE_7PARAM, the longitude calculated is correct here. However the latitude calculated here is incorrect.

Example

Source: EPSG:27700
Source towgs84: 446.448,-125.157,542.06,0.15,0.247,0.842,-20.489
Target: EPSG:3857
Coordinate: 327420.988668, 690284.547110

Transformation Step	Actual Value	Expected Value (when differs)
EPSG:27700	327420.988668, 690284.547110
Geodetic	-0.055272477580890766, 0.9791275748904039
Geocentric	3559822.9034467353, -196960.84735299222
WGS84	3560203.529194202, -197071.2700198092
Geodetic	-0.05529750149700846, 0.9760128672308176	-0.05529750149700847, 0.9791262638543461
Projected	-0.05529750149700846, 1.182600513116566	-0.05529750149700846, 1.1881697155477464
EPSG:3857	-352695.04030562507, 7542788.088927755	-352695.04030562507, 7578309.225014557

Reverse latitude transformations also fail in online coordinate converters from EPSG:3857 WGS84 values
https://epsg.io/transform#s_srs=4978&t_srs=3857&x=3551955.425120555&y=-196614.70501398295
Also tried here and here.