<div dir="ltr"><div>Back to report some kind of success.</div><div><br></div><div>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.</div><div><br></div><div>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:</div><div><br></div><div>+proj=pipeline<br>+step +inv +proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0 +ellps=bessel<br>+step +proj=hgridshift +grids=BETA2007.gsb<br>+step +proj=unitconvert +xy_in=rad +xy_out=deg</div><div><br></div><div>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!<br></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">Am Do., 15. Apr. 2021 um 16:26 Uhr schrieb XCTrails <<a href="mailto:info@xctrails.org">info@xctrails.org</a>>:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div>Ok, I clearly need to research this more judging from your answers so far.</div><div><br></div><div>Maybe as background from where the problem originates<span></span>:</div><div><br></div><div>I'm getting shapes in
<span>EPSG:31468</span> 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).</div><div><br></div><div>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.<br></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">Am Do., 15. Apr. 2021 um 15:19 Uhr schrieb Even Rouault <<a href="mailto:even.rouault@spatialys.com" target="_blank">even.rouault@spatialys.com</a>>:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
<div>
<p>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.<br>
</p>
<p>So the most precise transformation uses the BETA2007 grid, or you
have 2 potential Helmert-based transformations depending on the
location<br>
</p>
<p>$ projinfo -s EPSG:4326 -t EPSG:31468 --spatial-test intersects
-o PROJ<br>
Candidate operations found: 3<br>
-------------------------------------<br>
Operation No. 1:<br>
<br>
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.<br>
<br>
PROJ string:<br>
+proj=pipeline<br>
+step +proj=axisswap +order=2,1<br>
+step +proj=unitconvert +xy_in=deg +xy_out=rad<br>
+step +inv +proj=hgridshift +grids=de_adv_BETA2007.tif<br>
+step +proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0
+ellps=bessel<br>
+step +proj=axisswap +order=2,1<br>
<br>
-------------------------------------<br>
Operation No. 2:<br>
<br>
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.<br>
<br>
PROJ string:<br>
+proj=pipeline<br>
+step +proj=axisswap +order=2,1<br>
+step +proj=unitconvert +xy_in=deg +xy_out=rad<br>
+step +proj=push +v_3<br>
+step +proj=cart +ellps=WGS84<br>
+step +inv +proj=helmert +x=598.1 +y=73.7 +z=418.2 +rx=0.202
+ry=0.045<br>
+rz=-2.455 +s=6.7 +convention=position_vector<br>
+step +inv +proj=cart +ellps=bessel<br>
+step +proj=pop +v_3<br>
+step +proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0
+ellps=bessel<br>
+step +proj=axisswap +order=2,1<br>
<br>
-------------------------------------<br>
Operation No. 3:<br>
<br>
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.<br>
<br>
PROJ string:<br>
+proj=pipeline<br>
+step +proj=axisswap +order=2,1<br>
+step +proj=unitconvert +xy_in=deg +xy_out=rad<br>
+step +proj=push +v_3<br>
+step +proj=cart +ellps=WGS84<br>
+step +inv +proj=helmert +x=612.4 +y=77 +z=440.2 +rx=-0.054
+ry=0.057<br>
+rz=-2.797 +s=2.55 +convention=position_vector<br>
+step +inv +proj=cart +ellps=bessel<br>
+step +proj=pop +v_3<br>
+step +proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0
+ellps=bessel<br>
+step +proj=axisswap +order=2,1<br>
<br>
</p>
<div>Le 15/04/2021 à 14:51, Alan Snow a
écrit :<br>
</div>
<blockquote type="cite">
<div dir="auto">Original question: <a href="https://github.com/pyproj4/pyproj/discussions/814" target="_blank">https://github.com/pyproj4/pyproj/discussions/814</a>
<div dir="auto"><br>
</div>
<div dir="auto">What he is asking is if TOWGS84 parameters are
needed for datum transformations.</div>
<div dir="auto"><br>
</div>
<div dir="auto">Thanks 👍</div>
</div>
<br>
<div class="gmail_quote">
<div dir="ltr" class="gmail_attr">On Thu, Apr 15, 2021, 2:32 AM
XCTrails <<a href="mailto:info@xctrails.org" target="_blank">info@xctrails.org</a>> wrote:<br>
</div>
<blockquote class="gmail_quote">
<div dir="ltr">
<div>Hi,</div>
<div><br>
</div>
<div>
I have input polygon in epsg:31468 which I need to
transform to epsg:4326 <br>
</div>
<div><br>
</div>
<div>My current code boils down to this:</div>
<div><br>
</div>
<div>project =
pyproj.Transformer.from_crs("EPSG:31468","EPSG:4326",
always_xy=True)</div>
<div>poly = transform(project.transform, poly)</div>
<div><br>
</div>
<div>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 ;-)</div>
<div><br>
</div>
<div>Can anyone shed any light on why I get these offsets
and what I can do to fix this?<br>
</div>
<div>There's no offset if I manually reproject this in QGIS
and compare it with the original.<br>
</div>
<div><br>
</div>
<div>Thanks,</div>
<div>Armin<br>
</div>
</div>
_______________________________________________<br>
PROJ mailing list<br>
<a href="mailto:PROJ@lists.osgeo.org" rel="noreferrer" target="_blank">PROJ@lists.osgeo.org</a><br>
<a href="https://lists.osgeo.org/mailman/listinfo/proj" rel="noreferrer noreferrer" target="_blank">https://lists.osgeo.org/mailman/listinfo/proj</a><br>
</blockquote>
</div>
<br>
<fieldset></fieldset>
<pre>_______________________________________________
PROJ mailing list
<a href="mailto:PROJ@lists.osgeo.org" target="_blank">PROJ@lists.osgeo.org</a>
<a href="https://lists.osgeo.org/mailman/listinfo/proj" target="_blank">https://lists.osgeo.org/mailman/listinfo/proj</a>
</pre>
</blockquote>
<pre cols="72">--
<a href="http://www.spatialys.com" target="_blank">http://www.spatialys.com</a>
My software is free, but my time generally not.</pre>
</div>
_______________________________________________<br>
PROJ mailing list<br>
<a href="mailto:PROJ@lists.osgeo.org" target="_blank">PROJ@lists.osgeo.org</a><br>
<a href="https://lists.osgeo.org/mailman/listinfo/proj" rel="noreferrer" target="_blank">https://lists.osgeo.org/mailman/listinfo/proj</a><br>
</blockquote></div>
</blockquote></div>