Converting GPS altitude from WGS84 to a specific EPSG #1073

shayannikoohemat · 2022-05-06T18:37:08Z

shayannikoohemat
May 6, 2022

I am trying to convert my GPS data from WGS84 EPSG:4326 to EPSG:6190 which is [Belge 1972 / Belgian Lambert 72 + Ostend height](https://epsg.io/6190).
I can convert the x and y successfully but no chance for z conversion. It returns the same GPS coordinates. One possibility is that the projection doesn't support height transformation but on ESPG website epsg:6190 is plus Ostend Height.

I tried it with ETRS89 instead of WGRS84 and also to convert Projections to .to_3d() but the same results.

How can I fix this? Is it the problem of PROJ or I am missing an argument in my transformation ?

Code Sample:

import pyproj

# longitude, latitude, altitude
points = [(3.747301, 50.415045, 88.763399),
          (3.747329, 50.415019, 88.866587),
          (3.747358, 50.414992, 88.960016)]

WGS84 = pyproj.CRS(4326).to_3d()
ETRS89 = pyproj.CRS(4937).to_3d()
ostend_height = pyproj.CRS(8370).to_3d()
target_epsg = pyproj.CRS(6190).to_3d()

transformer = pyproj.Transformer.from_crs(WGS84, target_epsg, always_xy=True)

for p in transformer.itransform(points):
    print('{:.3f} {:.3f} {:.3f}'.format(*p))

# result:
"""
105832.874 122868.785 88.763
105834.840 122865.876 88.867
105836.876 122862.856 88.960
"""

Expected result

Z should be in the range of 46.0 meters not 88.0.

Environment settings:

pyproj info:
    pyproj: 3.3.0
      PROJ: 8.2.0
  data dir: /home/usr/VENV/venv39/lib/python3.9/site-packages/pyproj/proj_dir/share/proj
user_data_dir: /home/shayan/.local/share/proj
PROJ DATA (recommended version): 1.8
PROJ Database: 1.2
EPSG Database: v10.038 [2021-10-21]
ESRI Database: ArcMap 12.8 [2021-05-06]
IGNF Database: 3.1.0 [2019-05-24]

System:
    python: 3.9.5 (default, Nov 23 2021, 15:27:38)  [GCC 9.3.0]
executable: /home/shayan/VENV/venv39/bin/python
   machine: Linux-5.4.0-109-generic-x86_64-with-glibc2.31

Answered by snowman2

May 11, 2022

Seems to be working:

>>> from pyproj import Transformer
>>> tr = Transformer.from_crs("EPSG:4979", "EPSG:6190")
>>> tr.transform(50.415045, 3.747301, 88.763399)
(105832.77837274286, 122868.77880230732, 46.2297680595955)
>>> tr = Transformer.from_crs("EPSG:4979", "EPSG:6190", always_xy=True)
>>> tr.transform(3.747301, 50.415045, 88.763399)
(105832.77837274286, 122868.77880230732, 46.2297680595955)
>>> import pyproj
>>> pyproj.show_versions()
pyproj info:
    pyproj: 3.3.2.dev0
      PROJ: 9.0.0
  data dir:~/miniconda/envs/pyproj/share/proj
user_data_dir: ~/.local/share/proj
PROJ DATA (recommended version): 1.9
PROJ Database: 1.2
EPSG Database: v10.054 [2022-02-13]
ESRI Database: ArcMap 12.9 [

View full answer

snowman2 · 2022-05-11T14:20:37Z

snowman2
May 11, 2022
Maintainer

Doesn't seem to be an issue with PROJ:

echo 50.415045 3.747301 88.763399 | cs2cs EPSG:4979 EPSG:6190
105832.78	122868.78 46.23

Candidate operations found: 3
-------------------------------------
Operation No. 1:

unknown id, Inverse of ETRS89 to WGS 84 (1) + ETRS89 to Ostend height (1) + Inverse of BD72 to ETRS89 (3) + Belgian Lambert 72, 1.03 m, Belgium - onshore., at least one grid missing

PROJ string:
+proj=pipeline
  +step +proj=axisswap +order=2,1
  +step +proj=unitconvert +xy_in=deg +xy_out=rad
  +step +inv +proj=vgridshift +grids=be_ign_hBG18.tif +multiplier=1
  +step +inv +proj=hgridshift +grids=be_ign_bd72lb72_etrs89lb08.tif
  +step +proj=lcc +lat_0=90 +lon_0=4.36748666666667 +lat_1=51.1666672333333
        +lat_2=49.8333339 +x_0=150000.013 +y_0=5400088.438 +ellps=intl

Grid be_ign_bd72lb72_etrs89lb08.tif needed but not found on the system. Can be obtained at https://cdn.proj.org/be_ign_bd72lb72_etrs89lb08.tif
Grid be_ign_hBG18.tif needed but not found on the system. Can be obtained at https://cdn.proj.org/be_ign_hBG18.tif

-------------------------------------
Operation No. 2:

unknown id, Inverse of ETRS89 to WGS 84 (1) + ETRS89 to Ostend height (1) + Inverse of BD72 to ETRS89 (2) + Belgian Lambert 72, 1.22 m, Belgium - onshore., at least one grid missing

PROJ string:
+proj=pipeline
  +step +proj=axisswap +order=2,1
  +step +proj=unitconvert +xy_in=deg +xy_out=rad
  +step +inv +proj=vgridshift +grids=be_ign_hBG18.tif +multiplier=1
  +step +proj=push +v_3
  +step +proj=cart +ellps=GRS80
  +step +inv +proj=helmert +x=-106.8686 +y=52.2978 +z=-103.7239 +rx=-0.3366
        +ry=0.457 +rz=-1.8422 +s=-1.2747 +convention=coordinate_frame
  +step +inv +proj=cart +ellps=intl
  +step +proj=pop +v_3
  +step +proj=lcc +lat_0=90 +lon_0=4.36748666666667 +lat_1=51.1666672333333
        +lat_2=49.8333339 +x_0=150000.013 +y_0=5400088.438 +ellps=intl

Grid be_ign_hBG18.tif needed but not found on the system. Can be obtained at https://cdn.proj.org/be_ign_hBG18.tif

-------------------------------------
Operation No. 3:

unknown id, Inverse of ETRS89 to WGS 84 (1) + ETRS89 to Ostend height (1) + Inverse of BD72 to ETRS89 (1) + Belgian Lambert 72, 2.02 m, Belgium - onshore., at least one grid missing

PROJ string:
+proj=pipeline
  +step +proj=axisswap +order=2,1
  +step +proj=unitconvert +xy_in=deg +xy_out=rad
  +step +inv +proj=vgridshift +grids=be_ign_hBG18.tif +multiplier=1
  +step +proj=push +v_3
  +step +proj=cart +ellps=GRS80
  +step +inv +proj=helmert +x=-99.1 +y=53.3 +z=-112.5 +rx=0.419 +ry=-0.83
        +rz=1.885 +s=-1 +convention=position_vector
  +step +inv +proj=cart +ellps=intl
  +step +proj=pop +v_3
  +step +proj=lcc +lat_0=90 +lon_0=4.36748666666667 +lat_1=51.1666672333333
        +lat_2=49.8333339 +x_0=150000.013 +y_0=5400088.438 +ellps=intl

Grid be_ign_hBG18.tif needed but not found on the system. Can be obtained at https://cdn.proj.org/be_ign_hBG18.tif

5 replies

shayannikoohemat May 11, 2022
Author

Thanks for the answer. I'm running this command echo 50.415045 3.747301 88.763399 | cs2cs EPSG:4979 EPSG:6190 in my terminal with export PROJ_NETWORK=ON get wrong results for Z: 105832.89 122868.78 88.76

shayannikoohemat May 11, 2022
Author

Is this the correct way to enable the PROJ_NETWORK ?!

snowman2 May 11, 2022
Maintainer

I am using PROJ 9, you may consider upgrading your version of PROJ.

shayannikoohemat May 11, 2022
Author

Yes, I tested your solutions in python as well, I am not getting the same results, so the issue should be from the PROJ version. Thanks.

I also noticed you are using epsg:4979 instead of '4326' for WGS84, is this because one is 3D and the other 2D?
Is it the same if I use ETRS89:4937? How should I make sure a Projection has 3D support?

snowman2 May 11, 2022
Maintainer

Just skipping a step:

>>> import pyproj
>>> pyproj.CRS(4326).to_3d().to_epsg()
4979

snowman2 · 2022-05-11T14:24:24Z

snowman2
May 11, 2022
Maintainer

It seems like the operation with missing grids are missing even though PROJ NETWORK is enabled:

>>> from pyproj.transformer import TransformerGroup
>>> tg = TransformerGroup("EPSG:4979", "EPSG:6190")
<TransformerGroup: best_available=True>
- transformers: 3
- unavailable_operations: 0
>>> print(tg.transformers[0].to_proj4(pretty=True))
+proj=pipeline
  +step +proj=axisswap +order=2,1
  +step +proj=unitconvert +xy_in=deg +z_in=m +xy_out=rad +z_out=m
  +step +proj=push +v_3
  +step +proj=cart +ellps=WGS84
  +step +inv +proj=helmert +x=-106.8686 +y=52.2978 +z=-103.7239 +rx=-0.3366
        +ry=0.457 +rz=-1.8422 +s=-1.2747 +convention=coordinate_frame
  +step +inv +proj=cart +ellps=intl
  +step +proj=pop +v_3
  +step +proj=lcc +lat_0=90 +lon_0=4.36748666666667 +lat_1=51.1666672333333
        +lat_2=49.8333339 +x_0=150000.013 +y_0=5400088.438 +ellps=intl

1 reply

snowman2 May 11, 2022
Maintainer

Related to #1072

snowman2 · 2022-05-11T16:07:24Z

snowman2
May 11, 2022
Maintainer

Seems to be working:

>>> from pyproj import Transformer
>>> tr = Transformer.from_crs("EPSG:4979", "EPSG:6190")
>>> tr.transform(50.415045, 3.747301, 88.763399)
(105832.77837274286, 122868.77880230732, 46.2297680595955)
>>> tr = Transformer.from_crs("EPSG:4979", "EPSG:6190", always_xy=True)
>>> tr.transform(3.747301, 50.415045, 88.763399)
(105832.77837274286, 122868.77880230732, 46.2297680595955)
>>> import pyproj
>>> pyproj.show_versions()
pyproj info:
    pyproj: 3.3.2.dev0
      PROJ: 9.0.0
  data dir:~/miniconda/envs/pyproj/share/proj
user_data_dir: ~/.local/share/proj
PROJ DATA (recommended version): 1.9
PROJ Database: 1.2
EPSG Database: v10.054 [2022-02-13]
ESRI Database: ArcMap 12.9 [2022-02-18]
IGNF Database: 3.1.0 [2019-05-24]

System:
    python: 3.10.4 | packaged by conda-forge | (main, Mar 24 2022, 17:39:04) [GCC 10.3.0]
executable: ~/miniconda/envs/pyproj/bin/python
   machine: Linux-5.13.0-40-generic-x86_64-with-glibc2.31

Python deps:
~/miniconda/envs/pyproj/lib/python3.10/site-packages/_distutils_hack/__init__.py:30: UserWarning: Setuptools is replacing distutils.
  warnings.warn("Setuptools is replacing distutils.")
   certifi: 2021.10.08
    Cython: 0.29.28
setuptools: 60.10.0
       pip: 22.0.4

I recommend ensuring PROJ_NETWORK=ON.

1 reply

shayannikoohemat May 11, 2022
Author

I also noticed it works with PROJ 8 if I use this transformation:

>>> tr = pyproj.Transformer.from_crs(pyproj.CRS(4979), pyproj.CRS(31370).to_3d(), always_xy=True)
>>> tr.transform(3.747301, 50.415045, 88.763399)
(105832.87382028469, 122868.78452432249, 45.737076096236706)

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Converting GPS altitude from WGS84 to a specific EPSG #1073

{{title}}

Replies: 3 comments 7 replies

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{editor}}'s edit

{{editor}}'s edit

Select a reply

Converting GPS altitude from WGS84 to a specific EPSG #1073

shayannikoohemat May 6, 2022

Code Sample:

Expected result

Environment settings:

Replies: 3 comments · 7 replies

snowman2 May 11, 2022 Maintainer

shayannikoohemat May 11, 2022 Author

shayannikoohemat May 11, 2022 Author

snowman2 May 11, 2022 Maintainer

shayannikoohemat May 11, 2022 Author

snowman2 May 11, 2022 Maintainer

snowman2 May 11, 2022 Maintainer

snowman2 May 11, 2022 Maintainer

snowman2 May 11, 2022 Maintainer

shayannikoohemat May 11, 2022 Author

shayannikoohemat
May 6, 2022

Replies: 3 comments 7 replies

snowman2
May 11, 2022
Maintainer

shayannikoohemat May 11, 2022
Author

shayannikoohemat May 11, 2022
Author

snowman2 May 11, 2022
Maintainer

shayannikoohemat May 11, 2022
Author

snowman2 May 11, 2022
Maintainer

snowman2
May 11, 2022
Maintainer

snowman2 May 11, 2022
Maintainer

snowman2
May 11, 2022
Maintainer

shayannikoohemat May 11, 2022
Author