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

gmh04 commented 5 years ago

On epsg.io https://epsg.io/transform#s_srs=27700&t_srs=3857&x=327420.9886680&y=690284.5471100 the BNG easting/northing coodinates 327420.988668 690284.547110 reproject correctly to -352695.04 7578309.23. If I attempt to reproject the same coordinates using proj4j (0.1.0) with the following code:

Double easting = 327420.988668;
Double northing = 690284.54711;
CoordinateReferenceSystem BNG = crsFactory.createFromParameters("EPSG:27700", "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.06,0.15,0.247,0.842,-20.489 +units=m +no_defs");
CoordinateReferenceSystem MERCATOR = crsFactory.createFromName("EPSG:3857");
CoordinateTransform trans = ctFactory.createTransform(BNG, MERCATOR);
ProjCoordinate in = new ProjCoordinate(easting, northing);
ProjCoordinate out = new ProjCoordinate();
trans.transform(in, out);
log.info(String.format("%s : %s", out.x, out.y));

I get -352695.04030562507 : 7542788.088927755. The Y coord is around 30 miles south of the correct position. If I use version 1.0.0 of proj4j, I get -352535.43434034986 : 7578324.217405743, which is much closer but still around 200 meters out.

Incidently if I reproject first to EPSG:4326 then to EPSG:3857, I get the correct result.

gmh04 commented 5 years ago

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 commented 5 years ago

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 commented 4 years ago

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 commented 3 years ago

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 commented 3 years ago

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.

locationtech / proj4j