[PROJ] Possible error with +datum=potsdam and PROJ >= 5.0.0

classic Classic list List threaded Threaded
4 messages Options
Reply | Threaded
Open this post in threaded view
|

[PROJ] Possible error with +datum=potsdam and PROJ >= 5.0.0

Alexander Njemz
Hi,

we noticed differences in transforming geographic data from EPSG:31468
to WGS 84 using ST_Transform on different PostGIS installations and it
seems to come from different PROJ versions.

For example, we have a point from CAD drawings of the Nuremberg subway
network with the coordinates 4435877.4477166822 5475816.0547993957, so
we know it must be on the train tracks.

Converting this point using PROJ 4.9.3 gives us this point [1] with the
coordinates 11.1148165905971 49.4159011654427 and using PROJ 5.0.1
gives us this point [2] with the coordinates 11.1161576118633
49.4169712194789. However, updating the proj4text in the
spatial_ref_sys table for the srid 31468 from

+proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0 +datum=potsdam
+units=m +no_defs

to

+proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0 +ellps=bessel
on the PostGIS installation with version 5.0.1 of PROJ again gives
correct results.

PROJ version >= 5.0.0 seems to incorrectly handle conversions using
+datum=potsdam. From [3] one can download test data to verify correct
conversion from DHDN / GK to ETRS89 / UTM. Testing conversion of the
value in line 7 of the downloadable CSV file [4] the results for PROJ
version 4.9.3 and versions >= 5.0.0 differ. I have executed the
following command using PROJ versions 4.9.3, 5.0.0, 5.1.0, and 6.0.0

echo 4395886.918912 5819485.694352 | cs2cs +proj=tmerc +lat_0=0
+lon_0=12 +k=1 +x_0=4500000 +y_0=0 +datum=potsdam +units=m +no_defs +to
+init=epsg:25832

The output for 4.9.3 is:

599474.76 5817502.61 42.57

The output for the other versions is:

599559.03 5817663.04 0.00

Is this a bug in PROJ? And if so, what would be the recommended
workaround?

Kind Regards,


Alexander Njemz



[1] <a href="https://www.google.de/maps/place/49%C2%B024'57.2%22N+11%C2%B">https://www.google.de/maps/place/49%C2%B024'57.2%22N+11%C2%B
006'53.3%22E/@49.4158889,11.1126169,751m/data=!3m2!1e3!4b1!4m5!3m4!1s0x
0:0x0!8m2!3d49.4159012!4d11.1148166

[2] <a href="https://www.google.de/maps/place/49%C2%B025'01.1%22N+11%C2%B">https://www.google.de/maps/place/49%C2%B025'01.1%22N+11%C2%B
006'58.2%22E/@49.4169722,11.113978,750m/data=!3m2!1e3!4b1!4m5!3m4!1s0x0
:0x0!8m2!3d49.4169712!4d11.1161576

[3] http://crs.bkg.bund.de/crseu/crs/descrtrans/BeTA/de_dhdn2etr
s_beta.php

[4] http://crs.bkg.bund.de/crseu/crs/descrtrans/BeTA/BETA2007tes
tdaten.csv


--
------------------------------------
Alexander Njemz
Mondula GmbH

Geschäftsführer Michael Looft, Dr. Jan F. Ortmann
Amtsgericht Hamburg HRB 68092
Gerichtsstand Hamburg

Behringstrasse 28 a, H1, OG2,
22765 Hamburg
fon: ++49 (0)40 - 35 70 33 - 00
fon: ++49 (0)40 - 35 70 33 - 12

mailto:[hidden email]
mailto:[hidden email]
http://www.mondula.com
------------------------------------
Diese Email enthält vertrauliche und/oder rechtlich geschützte
Informationen. Wenn Sie nicht der richtige Adressat sind oder diese
Email irrtümlich erhalten haben, informieren Sie bitte sofort den
Absender und vernichten Sie diese Mail. Das unerlaubte Kopieren sowie
die unbefugte Weitergabe der Mail ist  nicht gestattet.
------------------------------------
This e-mail may contain confidential and/or privileged information. If
you are not the intended recipient (or have received this e-mail in
error) please notify the sender immediatley and destroy this e-mail.
Any unauthorised copying, disclosure or distribution of the material in
this e-mail is strictly forbidden.
------------------------------------
_______________________________________________
PROJ mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Possible error with +datum=potsdam and PROJ >= 5.0.0

Even Rouault-2
Alexander,

In PROJ 4.9.3 and earlier, +datum=potsdam expanded to
+ellps=bessel +towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.7

In PROJ 5.0.0 and later, +datum=postdam was changed to expand to
+ellps=bessel +nadgrids=@BETA2007.gsb

That is it will use the more precise BETA2007.gsb grid rather than the towgs84
clause, but this require the grid to be installed alongside other PROJ grids.
If not present, then no datum shift is done, hence the result you see.

So with PROJ 5.0.0 and later, you could get the same result as PROJ 4.9.3 if
you don't have the grid installed, by replacing "+datum=potsdam" with
"+ellps=bessel +towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.7"
(I don't understand why you'd get a correct result with just +ellps=bessel
without a +towgs84 or +nadgrids. That doesn't make sense to me)

With proj 5.2.0

$ echo "4395886.918912 5819485.694352" | cs2cs -f "%.4f" +proj=tmerc +lat_0=0
+lon_0=12 +k=1 +x_0=4500000 +y_0=0 +ellps=bessel
+towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.7 +to +init=epsg:25832
599474.7609 5817502.6096 42.5689

With proj 6, things will work transparently:

- with BETA2007.gdb installed, you get exactly the result from the .csv file:

$ echo "4395886.918912 5819485.694352" | cs2cs -f "%.4f" EPSG:31464 EPSG:25832
599474.9342 5817502.6270 0.0000

- without it installed, it fallbacks to the Helmert transformation

$ echo "4395886.918912 5819485.694352" | cs2cs -f "%.4f" EPSG:31464 EPSG:25832
599474.7609 5817502.6096 0.0000

Even


> Hi,
>
> we noticed differences in transforming geographic data from EPSG:31468
> to WGS 84 using ST_Transform on different PostGIS installations and it
> seems to come from different PROJ versions.
>
> For example, we have a point from CAD drawings of the Nuremberg subway
> network with the coordinates 4435877.4477166822 5475816.0547993957, so
> we know it must be on the train tracks.
>
> Converting this point using PROJ 4.9.3 gives us this point [1] with the
> coordinates 11.1148165905971 49.4159011654427 and using PROJ 5.0.1
> gives us this point [2] with the coordinates 11.1161576118633
> 49.4169712194789. However, updating the proj4text in the
> spatial_ref_sys table for the srid 31468 from
>
> +proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0 +datum=potsdam
> +units=m +no_defs
>
> to
>
> +proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0 +ellps=bessel
> on the PostGIS installation with version 5.0.1 of PROJ again gives
> correct results.
>
> PROJ version >= 5.0.0 seems to incorrectly handle conversions using
> +datum=potsdam. From [3] one can download test data to verify correct
> conversion from DHDN / GK to ETRS89 / UTM. Testing conversion of the
> value in line 7 of the downloadable CSV file [4] the results for PROJ
> version 4.9.3 and versions >= 5.0.0 differ. I have executed the
> following command using PROJ versions 4.9.3, 5.0.0, 5.1.0, and 6.0.0
>
> echo 4395886.918912 5819485.694352 | cs2cs +proj=tmerc +lat_0=0
> +lon_0=12 +k=1 +x_0=4500000 +y_0=0 +datum=potsdam +units=m +no_defs +to
> +init=epsg:25832
>
> The output for 4.9.3 is:
>
> 599474.76 5817502.61 42.57
>
> The output for the other versions is:
>
> 599559.03 5817663.04 0.00
>
> Is this a bug in PROJ? And if so, what would be the recommended
> workaround?
>
> Kind Regards,
>
>
> Alexander Njemz
>
>
>
> [1] <a href="https://www.google.de/maps/place/49%C2%B024'57.2%22N+11%C2%B">https://www.google.de/maps/place/49%C2%B024'57.2%22N+11%C2%B
> 006'53.3%22E/@49.4158889,11.1126169,751m/data=!3m2!1e3!4b1!4m5!3m4!1s0x
> 0:0x0!8m2!3d49.4159012!4d11.1148166
>
> [2] <a href="https://www.google.de/maps/place/49%C2%B025'01.1%22N+11%C2%B">https://www.google.de/maps/place/49%C2%B025'01.1%22N+11%C2%B
> 006'58.2%22E/@49.4169722,11.113978,750m/data=!3m2!1e3!4b1!4m5!3m4!1s0x0
>
> :0x0!8m2!3d49.4169712!4d11.1161576
>
> [3] http://crs.bkg.bund.de/crseu/crs/descrtrans/BeTA/de_dhdn2etr
> s_beta.php
>
> [4] http://crs.bkg.bund.de/crseu/crs/descrtrans/BeTA/BETA2007tes
> tdaten.csv


--
Spatialys - Geospatial professional services
http://www.spatialys.com
_______________________________________________
PROJ mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Possible error with +datum=potsdam and PROJ >= 5.0.0

Even Rouault-2
Actually, I realize there must be something fishy with your PROJ installation.
The BETA2007.gsb grid is part of the
https://download.osgeo.org/proj/proj-datumgrid-1.8.zip package which is
normally installed together with PROJ with usual binary distributions, so the
changed to expanding to +nadgrids=@BETA2007.gsb was reasonable.
If you compile PROJ from source, you must unzip its content in
${proj_install_prefix}/share/proj

> Alexander,
>
> In PROJ 4.9.3 and earlier, +datum=potsdam expanded to
> +ellps=bessel +towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.7
>
> In PROJ 5.0.0 and later, +datum=postdam was changed to expand to
> +ellps=bessel +nadgrids=@BETA2007.gsb
>
> That is it will use the more precise BETA2007.gsb grid rather than the
> towgs84 clause, but this require the grid to be installed alongside other
> PROJ grids. If not present, then no datum shift is done, hence the result
> you see.
>
> So with PROJ 5.0.0 and later, you could get the same result as PROJ 4.9.3 if
> you don't have the grid installed, by replacing "+datum=potsdam" with
> "+ellps=bessel +towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.7"
> (I don't understand why you'd get a correct result with just +ellps=bessel
> without a +towgs84 or +nadgrids. That doesn't make sense to me)
>
> With proj 5.2.0
>
> $ echo "4395886.918912 5819485.694352" | cs2cs -f "%.4f" +proj=tmerc
> +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0 +ellps=bessel
> +towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.7 +to +init=epsg:25832
> 599474.7609 5817502.6096 42.5689
>
> With proj 6, things will work transparently:
>
> - with BETA2007.gdb installed, you get exactly the result from the .csv
> file:
>
> $ echo "4395886.918912 5819485.694352" | cs2cs -f "%.4f" EPSG:31464
> EPSG:25832 599474.9342 5817502.6270 0.0000
>
> - without it installed, it fallbacks to the Helmert transformation
>
> $ echo "4395886.918912 5819485.694352" | cs2cs -f "%.4f" EPSG:31464
> EPSG:25832 599474.7609 5817502.6096 0.0000
>
> Even
>
> > Hi,
> >
> > we noticed differences in transforming geographic data from EPSG:31468
> > to WGS 84 using ST_Transform on different PostGIS installations and it
> > seems to come from different PROJ versions.
> >
> > For example, we have a point from CAD drawings of the Nuremberg subway
> > network with the coordinates 4435877.4477166822 5475816.0547993957, so
> > we know it must be on the train tracks.
> >
> > Converting this point using PROJ 4.9.3 gives us this point [1] with the
> > coordinates 11.1148165905971 49.4159011654427 and using PROJ 5.0.1
> > gives us this point [2] with the coordinates 11.1161576118633
> > 49.4169712194789. However, updating the proj4text in the
> > spatial_ref_sys table for the srid 31468 from
> >
> > +proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0 +datum=potsdam
> > +units=m +no_defs
> >
> > to
> >
> > +proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0 +ellps=bessel
> > on the PostGIS installation with version 5.0.1 of PROJ again gives
> > correct results.
> >
> > PROJ version >= 5.0.0 seems to incorrectly handle conversions using
> > +datum=potsdam. From [3] one can download test data to verify correct
> > conversion from DHDN / GK to ETRS89 / UTM. Testing conversion of the
> > value in line 7 of the downloadable CSV file [4] the results for PROJ
> > version 4.9.3 and versions >= 5.0.0 differ. I have executed the
> > following command using PROJ versions 4.9.3, 5.0.0, 5.1.0, and 6.0.0
> >
> > echo 4395886.918912 5819485.694352 | cs2cs +proj=tmerc +lat_0=0
> > +lon_0=12 +k=1 +x_0=4500000 +y_0=0 +datum=potsdam +units=m +no_defs +to
> > +init=epsg:25832
> >
> > The output for 4.9.3 is:
> >
> > 599474.76 5817502.61 42.57
> >
> > The output for the other versions is:
> >
> > 599559.03 5817663.04 0.00
> >
> > Is this a bug in PROJ? And if so, what would be the recommended
> > workaround?
> >
> > Kind Regards,
> >
> >
> > Alexander Njemz
> >
> >
> >
> > [1] <a href="https://www.google.de/maps/place/49%C2%B024'57.2%22N+11%C2%B">https://www.google.de/maps/place/49%C2%B024'57.2%22N+11%C2%B
> > 006'53.3%22E/@49.4158889,11.1126169,751m/data=!3m2!1e3!4b1!4m5!3m4!1s0x
> > 0:0x0!8m2!3d49.4159012!4d11.1148166
> >
> > [2] <a href="https://www.google.de/maps/place/49%C2%B025'01.1%22N+11%C2%B">https://www.google.de/maps/place/49%C2%B025'01.1%22N+11%C2%B
> > 006'58.2%22E/@49.4169722,11.113978,750m/data=!3m2!1e3!4b1!4m5!3m4!1s0x0
> >
> > :0x0!8m2!3d49.4169712!4d11.1161576
> >
> > [3] http://crs.bkg.bund.de/crseu/crs/descrtrans/BeTA/de_dhdn2etr
> > s_beta.php
> >
> > [4] http://crs.bkg.bund.de/crseu/crs/descrtrans/BeTA/BETA2007tes
> > tdaten.csv


--
Spatialys - Geospatial professional services
http://www.spatialys.com
_______________________________________________
PROJ mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Possible error with +datum=potsdam and PROJ >= 5.0.0

Alexander Njemz
On Tue, 2019-05-07 at 21:52 +0200, Even Rouault wrote:

> Actually, I realize there must be something fishy with your PROJ
> installation.
> The BETA2007.gsb grid is part of the
> https://download.osgeo.org/proj/proj-datumgrid-1.8.zip package which
> is
> normally installed together with PROJ with usual binary
> distributions, so the
> changed to expanding to +nadgrids=@BETA2007.gsb was reasonable.
> If you compile PROJ from source, you must unzip its content in
> ${proj_install_prefix}/share/proj

Ok, thank you. That makes a lot of sense. It is a little unfortunate
that there is no warning or anything when BETA2007.gsb is not
available.

I was using a PostGIS Docker image based on Alpine Linux and indeed,
the grid file is not there. I will try to contact the maintainer of the
PROJ package for Alpine Linux, and ask him if he could update the
package.

Thanks again.

Best,

Alexander

>
> > Alexander,
> >
> > In PROJ 4.9.3 and earlier, +datum=potsdam expanded to
> > +ellps=bessel +towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.7
> >
> > In PROJ 5.0.0 and later, +datum=postdam was changed to expand to
> > +ellps=bessel +nadgrids=@BETA2007.gsb
> >
> > That is it will use the more precise BETA2007.gsb grid rather than
> > the
> > towgs84 clause, but this require the grid to be installed alongside
> > other
> > PROJ grids. If not present, then no datum shift is done, hence the
> > result
> > you see.
> >
> > So with PROJ 5.0.0 and later, you could get the same result as PROJ
> > 4.9.3 if
> > you don't have the grid installed, by replacing "+datum=potsdam"
> > with
> > "+ellps=bessel +towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.7"
> > (I don't understand why you'd get a correct result with just
> > +ellps=bessel
> > without a +towgs84 or +nadgrids. That doesn't make sense to me)
> >
> > With proj 5.2.0
> >
> > $ echo "4395886.918912 5819485.694352" | cs2cs -f "%.4f"
> > +proj=tmerc
> > +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0 +ellps=bessel
> > +towgs84=598.1,73.7,418.2,0.202,0.045,-2.455,6.7 +to
> > +init=epsg:25832
> > 599474.7609 5817502.6096 42.5689
> >
> > With proj 6, things will work transparently:
> >
> > - with BETA2007.gdb installed, you get exactly the result from the
> > .csv
> > file:
> >
> > $ echo "4395886.918912 5819485.694352" | cs2cs -f "%.4f" EPSG:31464
> > EPSG:25832 599474.9342 5817502.6270 0.0000
> >
> > - without it installed, it fallbacks to the Helmert transformation
> >
> > $ echo "4395886.918912 5819485.694352" | cs2cs -f "%.4f" EPSG:31464
> > EPSG:25832 599474.7609 5817502.6096 0.0000
> >
> > Even
> >
> > > Hi,
> > >
> > > we noticed differences in transforming geographic data from
> > > EPSG:31468
> > > to WGS 84 using ST_Transform on different PostGIS installations
> > > and it
> > > seems to come from different PROJ versions.
> > >
> > > For example, we have a point from CAD drawings of the Nuremberg
> > > subway
> > > network with the coordinates 4435877.4477166822
> > > 5475816.0547993957, so
> > > we know it must be on the train tracks.
> > >
> > > Converting this point using PROJ 4.9.3 gives us this point [1]
> > > with the
> > > coordinates 11.1148165905971 49.4159011654427 and using PROJ
> > > 5.0.1
> > > gives us this point [2] with the coordinates 11.1161576118633
> > > 49.4169712194789. However, updating the proj4text in the
> > > spatial_ref_sys table for the srid 31468 from
> > >
> > > +proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0
> > > +datum=potsdam
> > > +units=m +no_defs
> > >
> > > to
> > >
> > > +proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0
> > > +ellps=bessel
> > > on the PostGIS installation with version 5.0.1 of PROJ again
> > > gives
> > > correct results.
> > >
> > > PROJ version >= 5.0.0 seems to incorrectly handle conversions
> > > using
> > > +datum=potsdam. From [3] one can download test data to verify
> > > correct
> > > conversion from DHDN / GK to ETRS89 / UTM. Testing conversion of
> > > the
> > > value in line 7 of the downloadable CSV file [4] the results for
> > > PROJ
> > > version 4.9.3 and versions >= 5.0.0 differ. I have executed the
> > > following command using PROJ versions 4.9.3, 5.0.0, 5.1.0, and
> > > 6.0.0
> > >
> > > echo 4395886.918912 5819485.694352 | cs2cs +proj=tmerc +lat_0=0
> > > +lon_0=12 +k=1 +x_0=4500000 +y_0=0 +datum=potsdam +units=m
> > > +no_defs +to
> > > +init=epsg:25832
> > >
> > > The output for 4.9.3 is:
> > >
> > > 599474.76 5817502.61 42.57
> > >
> > > The output for the other versions is:
> > >
> > > 599559.03 5817663.04 0.00
> > >
> > > Is this a bug in PROJ? And if so, what would be the recommended
> > > workaround?
> > >
> > > Kind Regards,
> > >
> > >
> > > Alexander Njemz
> > >
> > >
> > >
> > > [1] https://www.google.de/maps/place/49%C2%B024'57.2%22N+1
> > > 1%C2%B
> > > 006'53.3%22E/@49.4158889,11.1126169,751m/data=!3m2!1e3!4b1!4m5!3m
> > > 4!1s0x
> > > 0:0x0!8m2!3d49.4159012!4d11.1148166
> > >
> > > [2] https://www.google.de/maps/place/49%C2%B025'01.1%22N+1
> > > 1%C2%B
> > > 006'58.2%22E/@49.4169722,11.113978,750m/data=!3m2!1e3!4b1!4m5!3m4
> > > !1s0x0
> > >
> > > :0x0!8m2!3d49.4169712!4d11.1161576
> > >
> > > [3] http://crs.bkg.bund.de/crseu/crs/descrtrans/BeTA/de_dh
> > > dn2etr
> > > s_beta.php
> > >
> > > [4] http://crs.bkg.bund.de/crseu/crs/descrtrans/BeTA/BETA2
> > > 007tes
> > > tdaten.csv
>
>
--
------------------------------------
Alexander Njemz
Mondula GmbH

Geschäftsführer Michael Looft, Dr. Jan F. Ortmann
Amtsgericht Hamburg HRB 68092
Gerichtsstand Hamburg

Behringstrasse 28 a, H1, OG2,
22765 Hamburg
fon: ++49 (0)40 - 35 70 33 - 00
fon: ++49 (0)40 - 35 70 33 - 12

mailto:[hidden email]
mailto:[hidden email]
http://www.mondula.com
------------------------------------
Diese Email enthält vertrauliche und/oder rechtlich geschützte
Informationen. Wenn Sie nicht der richtige Adressat sind oder diese
Email irrtümlich erhalten haben, informieren Sie bitte sofort den
Absender und vernichten Sie diese Mail. Das unerlaubte Kopieren sowie
die unbefugte Weitergabe der Mail ist  nicht gestattet.
------------------------------------
This e-mail may contain confidential and/or privileged information. If
you are not the intended recipient (or have received this e-mail in
error) please notify the sender immediatley and destroy this e-mail.
Any unauthorised copying, disclosure or distribution of the material in
this e-mail is strictly forbidden.
------------------------------------
_______________________________________________
PROJ mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/proj