Vertical datum transformation based on PROJ

[rhugonnet]

Hi all,
I've been trying to apply vertical transformation of data using pyproj. For example, converting heights from the WGS84 ellipsoid to the EGM96 geoid.
Using the EGM .gtx files of PROJ here, I get a good result using the old pyproj.Proj method with init and geoidgrids keyword argument (example inspired by https://gis.stackexchange.com/questions/340392/vertical-datum-transformation-using-pyproj):

import pyproj
lat = 43.70012234
lng = -79.41629234
z = 100
#WGS84 datum ellipsoid height
ellipsoid=pyproj.Proj(init="EPSG:4326")
#EGM96 geoid in Chile, we expect about 30 m difference
geoid=pyproj.Proj(init="EPSG:4326", geoidgrids='/home/atom/downloads/egm96_15.gtx')
print(pyproj.transform(ellipsoid, geoid, lng, lat, z))

for which I get:

:8: DeprecationWarning: This function is deprecated. See: https://pyproj4.github.io/pyproj/stable/gotchas.html#upgrading-to-pyproj-2-from-pyproj-1
(-79.41629234, 43.70012234, 137.28983978090645)

I looked for another way to do this with current Pyproj methods but couldn't find anything or any issue referenced on this.
I'm thus wondering:

1. Is there a non-deprecated method for performing the same vertical grid referencing + transformation?
2. Would it be a possibility to add an integrated download/checksum method in pyproj for automated vertical datum extraction, as exists in PROJ? Right now a lot of people do this manually in their own Python repos (e.g., https://github.com/adehecq/geoutils/tree/master/geoutils/EGM96) so there is quite a need in the community for this. I might not be aware of all existing solutions, but I think this belongs in Pyproj :)

Thanks a lot in advance!

[snowman2]

Is there a non-deprecated method for performing the same vertical grid referencing + transformation?

• https://gis.stackexchange.com/a/352300/144357
• https://pyproj4.github.io/pyproj/stable/gotchas.html#upgrading-to-pyproj-2-from-pyproj-1

Would it be a possibility to add an integrated download/checksum method in pyproj for automated vertical datum extraction, as exists in PROJ?

• For the download part: https://pyproj4.github.io/pyproj/stable/transformation_grids.html

automated vertical datum extraction, as exists in PROJ

If you are referring to the code with GDAL in it, then no. I am strongly opposed to adding GDAL as a dependency to pyproj due to the complexity it adds. However, I do like the concept of a centralizing the code somewhere for that. Maybe creating another package specifically for that purpose (pyproj-grids?). Potentially useful libraries to include in that package could include rioxarray and pyresample.

Another alternative would be to raise an issue upstream (https://github.com/osgeo/proj) to see if they would expose an API for reading grids through PROJ. If they do that, I would be more supportive of adding the logic directly to pyproj.

[snowman2]

Hopefully your question has been answered (I am assuming that is the case).

[rhugonnet]

Yes, thanks a lot! :)
And sorry for not getting back to you sooner, busy months.

[snowman907]

@rhugonnet how did you implement the vertical transformation?

Doing something like this obviously deosen't make sense:
transformer=Transformer.from_crs("epsg:4326","epsg:26934",geoidgrids=r'C:\Users\Eyal\Downloads\geoid09_ak.gtx') transformer.transform(lat_d,lon_d,abs_alt)

Did you need to write your own vertical transformation: https://gis.stackexchange.com/questions/352277/including-geoidgrids-when-initializing-projection-via-epsg/352300#352300??

[snowman2]

This should do it:

transformer=Transformer.from_crs(
    "epsg:4326",
    CRS(init="epsg:26934", geoidgrids=r'C:\Users\Eyal\Downloads\geoid09_ak.gtx')
)
 transformer.transform(lat_d,lon_d,abs_alt)

[snowman907]

thanks @snowman2 ,
so just to be clear I did:
from pyproj.crs import CRS transformer=Transformer.from_crs("epsg:4326",CRS(init="epsg:26934",geoidgrids=r'C:\Users\Eyal\Downloads\geoid09_ak.gtx')) transformer.transform(lat_d,lon_d,abs_alt)

this worked well. But I get the note that init is deprecated. So I tried without init
transformer=Transformer.from_crs("epsg:4326",CRS("epsg:26934",geoidgrids=r'C:\Users\Eyal\Downloads\geoid09_ak.gtx')) transformer.transform(lat_d,lon_d,abs_alt)
And got the following response:
`CRSError Traceback (most recent call last)
in
----> 1 transformer=Transformer.from_crs("epsg:4326",CRS("epsg:26934",geoidgrids=r'C:\Users\Eyal\Downloads\geoid09_ak.gtx'))
2 transformer.transform(lat_d,lon_d,abs_alt)

c:\users\eyal\appdata\local\programs\python\python37\lib\site-packages\pyproj\crs\crs.py in init(self, projparams, **kwargs)
294 projstring = _prepare_from_string(" ".join((projstring, projkwargs)))
295
--> 296 super().init(projstring)
297
298 @staticmethod

pyproj_crs.pyx in pyproj._crs._CRS.init()

CRSError: Invalid projection: epsg:26934 +geoidgrids=C:\Users\Eyal\Downloads\geoid09_ak.gtx +type=crs: (Internal Proj Error: proj_create: unrecognized format / unknown name)`
Is there yet a way to not use the deprecated way?
thanks again

[snowman2]

Is there yet a way to not use the deprecated way?

https://gis.stackexchange.com/a/352300/144357

[rhugonnet]

That's a lot of snowmans in the same conversation! ⛄

@snowman907 I kept the init= syntax for now in my code. If the method becomes deprecated in the future, it would certainly be good to have a one-liner to input a geoidgrid file implemented in pyproj.
It wasn't immediately clear to me (from pyproj documentation) that all available grids available there (https://github.com/OSGeo/PROJ-data) could be used directly if you have the package proj-data installed. For your grid:

crs_init = pyproj.CRS(4326)
# to remove warning
with warnings.catch_warnings():
    warnings.filterwarnings("ignore", module="pyproj")
        # your grid file is available here: https://github.com/OSGeo/PROJ-data/tree/master/us_noaa, filename works 
        # only if your pyproj.datadir is defined properly (good chance it's OK by default) and if you have proj-data
        crs_dest = pyproj.Proj(init="EPSG:26934", geoidgrids='us_noaa_geoid09_ak.tif').crs
transform = Transformer.from_crs(crs_init,crs_dest)

Along with a team of colleagues, we're in the process of writing a DEM package to put together (robustly and consistently, if possible) all the DEM tools we've been developing for the past ~10 years. It's early stages still, but here's the code for vertical grid shifting within a DEM class: https://github.com/GlacioHack/xdem/blob/main/xdem/dem.py#L116
@snowman2 We've also been writing so many things that now resemble rasterio (GDAL wrapper) and rioxarray (link GDAL -netCDF/xarray), it's time to stop reinventing the wheel and contribute to the ongoing efforts (even if it takes more time up front)!

[snowman907]

Thanks, @rhugonnet for your insight and @snowman2 for your hard work developing pyproj and responding patiently to a large number of issues with the same question.

Yay for snowmans lol