[PROJ] Transforming from epsg:4326 to epsg:31468 produces a 2 meter offset

[Even Rouault]

projinfo is often the answer to such questions. Datum transformations 
*are* generally needed. What is generally no longer required in PROJ 6 
era is to use the explict +nadgrids / +towgs84 terms. PROJ will find 
which ones are needed, but the user may still have to arbitrate because 
there's not always a single possibility.

So the most precise transformation uses the BETA2007 grid, or you have 2 
potential Helmert-based transformations depending on the location

$ projinfo -s EPSG:4326 -t EPSG:31468 --spatial-test intersects -o PROJ
Candidate operations found: 3
-------------------------------------
Operation No. 1:

unknown id, Inverse of DHDN to WGS 84 (4) + 3-degree Gauss-Kruger zone 
4, 1 m, Germany - onshore - states of Baden-Wurtemberg, Bayern, Berlin, 
Brandenburg, Bremen, Hamburg, Hessen, Mecklenburg-Vorpommern, 
Niedersachsen, Nordrhein-Westfalen, Rheinland-Pfalz, Saarland, Sachsen, 
Sachsen-Anhalt, Schleswig-Holstein, Thuringen.

PROJ string:
+proj=pipeline
   +step +proj=axisswap +order=2,1
   +step +proj=unitconvert +xy_in=deg +xy_out=rad
   +step +inv +proj=hgridshift +grids=de_adv_BETA2007.tif
   +step +proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0 
+ellps=bessel
   +step +proj=axisswap +order=2,1

-------------------------------------
Operation No. 2:

unknown id, Inverse of DHDN to WGS 84 (2) + 3-degree Gauss-Kruger zone 
4, 3 m, Germany - states of former West Germany onshore - 
Baden-Wurtemberg, Bayern, Bremen, Hamburg, Hessen, Niedersachsen, 
Nordrhein-Westfalen, Rheinland-Pfalz, Saarland, Schleswig-Holstein.

PROJ string:
+proj=pipeline
   +step +proj=axisswap +order=2,1
   +step +proj=unitconvert +xy_in=deg +xy_out=rad
   +step +proj=push +v_3
   +step +proj=cart +ellps=WGS84
   +step +inv +proj=helmert +x=598.1 +y=73.7 +z=418.2 +rx=0.202 +ry=0.045
         +rz=-2.455 +s=6.7 +convention=position_vector
   +step +inv +proj=cart +ellps=bessel
   +step +proj=pop +v_3
   +step +proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0 
+ellps=bessel
   +step +proj=axisswap +order=2,1

-------------------------------------
Operation No. 3:

unknown id, Inverse of DHDN to WGS 84 (3) + 3-degree Gauss-Kruger zone 
4, 2 m, Germany - states of former East Germany - Berlin, Brandenburg; 
Mecklenburg-Vorpommern; Sachsen; Sachsen-Anhalt; Thuringen.

PROJ string:
+proj=pipeline
   +step +proj=axisswap +order=2,1
   +step +proj=unitconvert +xy_in=deg +xy_out=rad
   +step +proj=push +v_3
   +step +proj=cart +ellps=WGS84
   +step +inv +proj=helmert +x=612.4 +y=77 +z=440.2 +rx=-0.054 +ry=0.057
         +rz=-2.797 +s=2.55 +convention=position_vector
   +step +inv +proj=cart +ellps=bessel
   +step +proj=pop +v_3
   +step +proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0 
+ellps=bessel
   +step +proj=axisswap +order=2,1

Le 15/04/2021 à 14:51, Alan Snow a écrit :
> Original question: https://github.com/pyproj4/pyproj/discussions/814 
> <https://github.com/pyproj4/pyproj/discussions/814>
>
> What he is asking is if TOWGS84 parameters are needed for datum 
> transformations.
>
> Thanks 👍
>
> On Thu, Apr 15, 2021, 2:32 AM XCTrails <info at xctrails.org 
> <mailto:info at xctrails.org>> wrote:
>
>     Hi,
>
>     I have input polygon in epsg:31468 which I need to transform to
>     epsg:4326
>
>     My current code boils down to this:
>
>     project = pyproj.Transformer.from_crs("EPSG:31468","EPSG:4326",
>     always_xy=True)
>     poly = transform(project.transform, poly)
>
>     The problem however is that the resulting polygons have a 1-2
>     meter offset according to visual inspection in QGIS. On reddit,
>     someone suggested that I might also need a datum transformation,
>     on the pyproj4 issue tracker someone said, that datum
>     transformations are no longer needed in PROJ 6+ and that I should
>     ask for clarification here. So here I am ;-)
>
>     Can anyone shed any light on why I get these offsets and what I
>     can do to fix this?
>     There's no offset if I manually reproject this in QGIS and compare
>     it with the original.
>
>     Thanks,
>     Armin
>     _______________________________________________
>     PROJ mailing list
>     PROJ at lists.osgeo.org <mailto:PROJ at lists.osgeo.org>
>     https://lists.osgeo.org/mailman/listinfo/proj
>     <https://lists.osgeo.org/mailman/listinfo/proj>
>
>
> _______________________________________________
> PROJ mailing list
> PROJ at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/proj

-- 
http://www.spatialys.com
My software is free, but my time generally not.

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20210415/d6d1e1a4/attachment-0001.html>