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

Thu Apr 15 07:53:30 PDT 2021

I don’t understand why you are using EPSG:4326 (WGS84) when a 2 metre offset is not acceptable for you. EPSG:4326 is a general code for a collection (datum ensemble) of more precise defined datums (WGS84 realisations) with their own codes. The differences within the datums of the ensemble prevent a transformation to EPSG:4326 that is more precise than 2 metres. So if you what a higher accuracy, you will have to use a CRS based on another datum, e.g. EPSG:9000 (ITRF2014).

Regards, Jochem

From: PROJ <proj-bounces at lists.osgeo.org> On Behalf Of XCTrails
Sent: donderdag 15 april 2021 16:26
To: Even Rouault <even.rouault at spatialys.com>
Cc: PROJ <proj at lists.osgeo.org>
Subject: Re: [PROJ] Transforming from epsg:4326 to epsg:31468 produces a 2 meter offset

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<mailto: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<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

_______________________________________________

PROJ mailing list

PROJ at lists.osgeo.org<mailto: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.
_______________________________________________
PROJ mailing list
PROJ at lists.osgeo.org<mailto:PROJ at lists.osgeo.org>
https://lists.osgeo.org/mailman/listinfo/proj

Disclaimer:
De inhoud van dit bericht is uitsluitend bestemd voor geadresseerde.
Gebruik van de inhoud van dit bericht door anderen zonder toestemming van het Kadaster
is onrechtmatig. Mocht dit bericht ten onrechte bij u terecht komen, dan verzoeken wij u
dit direct te melden aan de verzender en het bericht te vernietigen.
Aan de inhoud van dit bericht kunnen geen rechten worden ontleend.

Disclaimer:
The content of this message is meant to be received by the addressee only.
Use of the content of this message by anyone other than the addressee without the consent
of the Kadaster is unlawful. If you have received this message, but are not the addressee,
please contact the sender immediately and destroy the message.
No rights can be derived from the content of this message.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20210415/0a9a3c3d/attachment-0001.html>