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

Sat Apr 17 10:07:40 PDT 2021

Back to report some kind of success.

First of all, none of the pipeline strings produced by projinfo here work
in my case. I either get a Pipeline: Bad step definition: inv (failed to
load datum shift file) error (although I have a copy of de_adv_BETA2007.tif
in my pyproj.datadir.get_data_dir()), or I get a null geometry when I use
BETA2007.gsb instead.

Before diving into yet another rathole and with now knowing that I can
initialize a pyproj.Transformer from a pipeline string, I looked again at
QGIS - which actually shows me a pipeline string when loading the
EPSG:31468 into a EPSG:4326 project - claiming that this transformation has
a 1m accuracy:

+proj=pipeline
+step +inv +proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0
+ellps=bessel
+step +proj=hgridshift +grids=BETA2007.gsb
+step +proj=unitconvert +xy_in=rad +xy_out=deg

With that at last, I get identical results between transforming the shapes
in QGIS or in Python. I know I'm still only scratching the surface of all
the details involved here, but I think as a CS and not a GIS guy, I can
live with that for now ;-) So thanks for pushing me into the right
direction!

Am Do., 15. Apr. 2021 um 16:26 Uhr schrieb XCTrails <info at xctrails.org>:

> 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/20210417/3a1b043b/attachment.html>