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

Thu Apr 15 07:26:23 PDT 2021

Ok, I clearly need to research this more judging from your answers so far.

Maybe as background from where the problem originates:

I'm getting shapes in EPSG:31468 which will subsequently be imported/mapped
to OpenstreetMap. Eventually, there will be updates to the EPSG:31468
shapes and I'm in the process of writing a script which should detect shape
changes/deletions/additions. In order to be able to geometrically compare
the 2 datasets, my current workaround is to perform the transformation from
EPSG:31468 to EPSG:4326 manually in QGIS, the diff then runs on 2 4326
geojson files and produces reasonable polygon matches with an intersection
over union approach. Ideally, I'd like to get rid of the manual
transformation step in QGIS (along the way understanding the underlying
issue).

So, if I understand this correctly: creating a pyproj transformation using
one of the mentioned PROJ strings/pipelines should yield better results?
Within the limits of WGS84 - which might by in the same range as my
initially observed offset.

Am Do., 15. Apr. 2021 um 15:19 Uhr schrieb Even Rouault <
even.rouault at spatialys.com>:

> 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
>
> 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> 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
>> https://lists.osgeo.org/mailman/listinfo/proj
>>
>
> _______________________________________________
> PROJ mailing listPROJ at lists.osgeo.orghttps://lists.osgeo.org/mailman/listinfo/proj
>
> -- http://www.spatialys.com
> My software is free, but my time generally not.
>
> _______________________________________________
> PROJ mailing list
> PROJ at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/proj
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20210415/7d08548a/attachment.html>