[postgis-tickets] [PostGIS] #4817: ST_Transform fails when both +geoidgrids and +nadgrids are present

[PostGIS]

#4817: ST_Transform fails when both +geoidgrids and +nadgrids are present
----------------------+---------------------------
 Reporter:  aberenyi  |      Owner:  pramsey
     Type:  defect    |     Status:  new
 Priority:  medium    |  Milestone:  PostGIS 3.1.0
Component:  postgis   |    Version:  3.0.x
 Keywords:            |
----------------------+---------------------------
 The following works just fine using `cs2cs`.

 `$ echo 16.582941335 47.710445103 291.7692| cs2cs -f "%.3f"
 +init=epsg:4326 +to +proj=somerc +lat_0=47.14439372222222
 +lon_0=19.04857177777778 +k_0=0.99993 +x_0=650000 +y_0=200000 +ellps=GRS67
 +nadgrids=etrs2eov_notowgs.gsb +geoidgrids=geoid_eht2014.gtx +units=m
 +no_defs`

 However, the same transformation fails in PostGIS.

 `$ psql -d db -c "SELECT
 ST_AsEWKT(ST_Transform(ST_GeomFromText('POINT(16.582941335 47.710445103
 291.7692)', 4326), '+proj=somerc +lat_0=47.14439372222222
 +lon_0=19.04857177777778 +k_0=0.99993 +x_0=650000 +y_0=200000 +ellps=GRS67
 +nadgrids=etrs2eov_notowgs.gsb +geoidgrids=geoid_eht2014.gtx +units=m
 +no_defs'))"`

 `ERROR:  proj_crs_is_swapped: proj_crs_get_coordinate_system returned NULL
 CONTEXT:  SQL function "st_transform" statement 1`

 Interestingly, if I omit either `+geoidgrids` or `+nadgrids`
 `ST_Transform` doesn't throw an error, but the coordinates are off.

 `$ psql -d db -c "SELECT
 ST_AsEWKT(ST_Transform(ST_GeomFromText('POINT(16.582941335 47.710445103
 291.7692)', 4326), '+proj=somerc +lat_0=47.14439372222222
 +lon_0=19.04857177777778 +k_0=0.99993 +x_0=650000 +y_0=200000 +ellps=GRS67
 +geoidgrids=geoid_eht2014.gtx +units=m +no_defs'))"`

 `POINT(465007.88411299663 265848.07307208166 246.45166349798842)`
 (formatting omitted for brewity)

 `$ psql -d db -c "SELECT
 ST_AsEWKT(ST_Transform(ST_GeomFromText('POINT(16.582941335 47.710445103
 291.7692)', 4326), '+proj=somerc +lat_0=47.14439372222222
 +lon_0=19.04857177777778 +k_0=0.99993 +x_0=650000 +y_0=200000 +ellps=GRS67
 +nadgrids=etrs2eov_notowgs.gsb +units=m +no_defs'))"`

 `POINT(465092.51654231607 265877.1140839074 291.7692)`

 This suggests that I can work around the problem by using a 2-step
 approach, but I'd like to aviod that if possible.

 What am I doing wrong?

 Env:
 - I'm using the PostGIS docker conainter
 https://hub.docker.com/r/postgis/postgis

 `$  psql -d db -c "SELECT PostGIS_Full_Version();"`
 `POSTGIS="3.1.0dev 50b1e70" [EXTENSION] PGSQL="120" GEOS="3.9.0dev-
 CAPI-1.14.0" PROJ="7.2.0" LIBXML="2.9.4" LIBJSON="0.12.1"
 LIBPROTOBUF="1.3.1" WAGYU="0.5.0 (Internal)"`

 - The geiod model and the grid is available here:
 https://github.com/OSGeoLabBp/eov2etrs

-- 
Ticket URL: <https://trac.osgeo.org/postgis/ticket/4817>
PostGIS <http://trac.osgeo.org/postgis/>
The PostGIS Trac is used for bug, enhancement & task tracking, a user and developer wiki, and a view into the subversion code repository of PostGIS project.