GDAL 3.0.1: coordinate transformations are inconsistent for different CRS types #1972

and-marsh · 2019-11-01T14:13:22Z

Expected behavior and actual behavior.

Hello,

I'm sorry to disturb you again, but there is something that I just can't understand with GDAL coordinate transformations.

I tried to force GDAL 3 to transform height for any system. It's important for us to make behaviour predictable in all cases and we hope you can make it a little more clear for us and other developers adopting GDAL 3.

All the tests have been performed using `GDAL 3.0.1` and `Proj4 6.2.0`:

I divided the tests into several cases:

1. 3D -> 3D

EPSG:4979 -> IGNF:RGNCGEO

Transform CRSs from registry:

./gdaltransform -s_srs EPSG:4979 -t_srs IGNF:RGNCGEO
160 -20
160 -20.0005291658479 -240.343272392638

Then get proj description of this CRS and try define target coordinate system with it:

./testepsg IGNF:
...
PROJ.4 rendering of [IGNF:RGNCGEO] = +proj=longlat +ellps=intl +no_defs

./gdaltransform -s_srs EPSG:4979 -t_srs '+proj=longlat +ellps=intl +no_defs' 
160 -20
160 -20 0

Height wasn't transformed, and second coordinate differs from example above. Let's add towgs84=0,0,0:

./gdaltransform -s_srs EPSG:4979 -t_srs '+proj=longlat +ellps=intl +no_defs +towgs84=0,0,0'
160 -20
160 -20.0005291658479 0

Now second coordinate was transformed correctly, but height still wasn't. Then add vunits=m to proj string:

./gdaltransform -s_srs EPSG:4979 -t_srs '+proj=longlat +ellps=intl +no_defs +towgs84=0,0,0 +vunits=m'
160 -20
160 -20.0005291658479 -240.343272392638

What are my conclusions based on this test?

when transformation is 3D -> 3D change in height above the datum is taken into account
proj4 string doesn't completely describe CRS unlike WKT. Transformations using CRS from registry and from proj4 string may be different, is it correct?
if +towgs84 params are not specified, it would be more correct to set it to default value +towgs84=0,0,0
to transform height using proj4 string it's necessary to specify +towgs84 and +vunits params GDAL 3.0.1: coordinate vertical transformation without geoid #1946 (comment)

Ok, next case.

EPSG:4979 -> EPSG:4985

The same flow as above:

./gdaltransform -s_srs EPSG:4979 -t_srs EPSG:4985
10 10
9.99984611111111 9.99995931437372 0

For some reason height wasn't transformed even if target CRS is defined from registry.

Now use proj string to force height transform:

./testepsg EPSG:4985
...
PROJ.4 rendering of [EPSG:4985] = +proj=longlat +ellps=WGS72 +towgs84=0,0,4.5,0,0,0.554,0.2263 +no_defs

./gdaltransform -s_srs EPSG:4979 -t_srs '+proj=longlat +ellps=WGS72 +towgs84=0,0,4.5,0,0,0.554,0.2263 +no_defs'
10 10
9.99984611111111 9.99995931437372 0

./gdaltransform -s_srs EPSG:4979 -t_srs '+proj=longlat +ellps=WGS72 +towgs84=0,0,4.5,0,0,0.554,0.2263 +no_defs +vunits=m'
10 10
9.99984611111111 9.99995931437372 -0.230804002843797

Here still 3D -> 3D transformation, but change in height above the datum is NOT taken into account, until it isn't forced by +vunits in proj4 string.

Is this behaviour expected?

2. 3D -> 2D

EPSG:4979 -> EPSG:4475

Here all works as expected:

./gdaltransform -s_srs EPSG:4979 -t_srs EPSG:4475
45 -13
44.9978863935678 -12.9974687407008 0

./testepsg EPSG:4475
...
PROJ.4 rendering of [EPSG:4475] = +proj=longlat +ellps=intl +towgs84=-381.788,-57.501,-256.673,0,0,0,0 +no_defs

./gdaltransform -s_srs EPSG:4979 -t_srs '+proj=longlat +ellps=intl +towgs84=-381.788,-57.501,-256.673,0,0,0,0 +no_defs'
45 -13
44.9978863935678 -12.9974687407008 0

./gdaltransform -s_srs EPSG:4979 -t_srs '+proj=longlat +ellps=intl +towgs84=-381.788,-57.501,-256.673,0,0,0,0 +no_defs +vunits=m'
45 -13
44.9978863935678 -12.9974687407008 -1.45678852126002

Conclusions:

for 2D systems change in height above the datum is NOT taken into account, until it isn't forced by +vunits in proj4 string

3. 3D -> PROJCS

EPSG:4979 -> EPSG:2189

The most unclear case, first try to use CRS from registry:

./gdaltransform -s_srs EPSG:4979 -t_srs EPSG:2189
-28 39
413304.964602352 4317337.75826103 -58.6802441990003

Height is transformed despite the fact that target CRS is horizontal. Try to use proj string:

./testepsg EPSG:2189
...
PROJ.4 rendering of [EPSG:2189] = +proj=utm +zone=26 +ellps=intl +towgs84=-104,167,-38,0,0,0,0 +units=m +no_defs

./gdaltransform -s_srs EPSG:4979 -t_srs '+proj=utm +zone=26 +ellps=intl +towgs84=-104,167,-38,0,0,0,0 +units=m +no_defs'
-28 39
413304.964602352 4317337.75826103 0

As expected height isn't transformed because vunits isn't specified

./gdaltransform -s_srs EPSG:4979 -t_srs '+proj=utm +zone=26 +ellps=intl +towgs84=-104,167,-38,0,0,0,0 +units=m +no_defs +vunits=m'
-28 39
413304.964602352 4317337.75826103 0

But after adding vunits height still isn't transformed!

Conclusions:

For some reason 3D -> PROJCS transforms height. But why if PROJCS is horizontal CRS?
Adding +towgs84 and +vunits doesn't force height transformation for PROJCS. How can i force height transformation for projected CRS?

I didn't find any reasonable explanations for such differences in behaviour for different CRS types. Am I approaching this correctly or I should be looking for a different workflow in general?

My questions

Using a system from registry differs from a system defined with the same parameters manually(for example by extracting the proj4 string). Is this behaviour correct?
In example 1.2(EPSG:4979 -> EPSG:4985) the height doesn't get transformed even when converting coordinates between two 3D systems from the registry. Is there a reason for this?
3D -> PROJCS transformation transforms height despite PROJCS being a horizontal CRS.
Adding +towgs84 and +vunits to proj4 description of CRS doesn't force height transformation for PROJCS. How can I force height transformation for projected CRS?
Is it correct to add default parameters towgs84=0,0,0 if they aren't specified? Will it affect transformation?

I will be happy to create a guide for coordinate transformations with GDAL >= 3 to help other developers with the same sort of questions, but I don't understand this enough yet. Would you be interested in this?

Thank you!

Operating system

macOS Catalina 10.15

GDAL version and provenance

GDAL 3.0.1, PROJ 6.2.0

The text was updated successfully, but these errors were encountered:

rouault · 2019-11-01T14:44:50Z

Yes this is a complicated topic, furthermore complicated by the fact that PROJ 6.x/7dev is still in development/tuning in that area.

The fact that "gdaltransform -s_srs EPSG:4979 -t_srs EPSG:2189" causes a change in height is a PROJ bug. The datum shift applied is "Azores Central 1948 to WGS 84 (1)" EPSG:1886. The method for it is "Geocentric translations (geog2D domain)", so PROJ is wrongly extending it into the 3D domain, but it shouldn't alter the Z value since this tranformation isn't intended for that (at least, as much as I understand the EPSG guidance notes). Can you file a ticket in PROJ github tracker about that one ?

The fact that you specify 3D systems as input and output doesn't necessarily mean that PROJ can and will transform the Z value. For that it must find a transformation in its database. If it doesn't find one, then it will leave the Z untouched. It can also not be able to do horizontal datum transform if it has not the available data. If you want to have a deeper understanding of what is going on when investigating a transformation, I suggest you use the projinfo -s {srs_def} -t {srs_def} command that comes with PROJ 6.

Regarding EPSG:4979 to EPSG:4985 (WGS72 3D), the available transformation in EPSG is EPSG:1237 which is a "Position Vector transformation (geog2D domain)", so valid for 2D only.

It is expected that using PROJ strings and EPSG code doesn't always lead to the same results. PROJ strings are a lossy way of encoding a CRS as now stated in https://gdal.org/doxygen/classOGRSpatialReference.html#a271b3de4caf844135b0c61e634860f2b

The fact of adding +vunits=m in gdaltransform -s_srs EPSG:4979 -t_srs '+proj=longlat +ellps=WGS72 +towgs84=0,0,4.5,0,0,0.554,0.2263 +no_defs +vunits=m' will make PROJ understand that the +towgs84 transformations is valid in the 3D domain (which is wrong given that the original EPSG transformation is valid only for horizontal)

…rt geog2D transformation to a geog3D Fixes for example EPSG:4979 to EPSG:2189, as raised in OSGeo/gdal#1972 (comment)

rouault · 2019-11-03T14:36:04Z

The fact that "gdaltransform -s_srs EPSG:4979 -t_srs EPSG:2189" causes a change in height is a PROJ bug. [...] Can you file a ticket in PROJ github tracker about that one ?

No longer needed. Fixed in PROJ per OSGeo/PROJ#1711

…rt geog2D transformation to a geog3D Fixes for example EPSG:4979 to EPSG:2189, as raised in OSGeo/gdal#1972 (comment)

and-marsh · 2019-11-05T12:30:24Z

Hello,

Thank you for your help!

Regarding to your answer, please, correct me if i'm wrong:

It is better to rely on GDAL in coordinate transformations since a forcing of transformation with proj4-string can provide wrong results
proj4-string representation is deprecated and it is better to use AUTORITY or WKT

But there are some features that i was able to use only with proj4-string:

Using vertical geoid shifts for datums for Europe and North America, during transformation from WGS84. GDAL 3.0.1: coordinate transformation with vertical shift #1836 (comment) .
Adding custom geoids to custom CRS or CRS from registry

Is there any workaround for this cases? Can i define geoid file for CRS using WKT or some functions provided by GDAL?

rouault · 2019-11-05T12:39:22Z

* It is better to rely on GDAL in coordinate transformations since a forcing of transformation with `proj4`-string can provide wrong results

* `proj4`-string representation is deprecated and it is better to use `AUTORITY` or `WKT`

Yes

Using vertical geoid shifts for datums for Europe and North America, during transformation from WGS84. #1836 (comment) .

If the vertical CRS is known and there is a transformation in the EPSG dataset, then it will be used

Adding custom geoids to custom CRS or CRS from registry

In PROJ master, the capability to add a GEOIDMODEL in WKT has been added per OSGeo/PROJ#1710. It was also possible to use the BOUNDCRS mechanism that gives a transformation between a VertCRS and a GeographicCRS using a grid.

e.g.

COMPOUNDCRS["unknown",
    GEOGCRS["unknown",
        DATUM["North American Datum 1983",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6269]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8901]],
        CS[ellipsoidal,2],
            AXIS["longitude",east,
                ORDER[1],
                ANGLEUNIT["degree",0.0174532925199433,
                    ID["EPSG",9122]]],
            AXIS["latitude",north,
                ORDER[2],
                ANGLEUNIT["degree",0.0174532925199433,
                    ID["EPSG",9122]]]],
    BOUNDCRS[
        SOURCECRS[
            VERTCRS["unknown",
                VDATUM["unknown"],
                CS[vertical,1],
                    AXIS["gravity-related height (H)",up,
                        LENGTHUNIT["metre",1,
                            ID["EPSG",9001]]]]],
        TARGETCRS[
            GEOGCRS["WGS 84",
                DATUM["World Geodetic System 1984",
                    ELLIPSOID["WGS 84",6378137,298.257223563,
                        LENGTHUNIT["metre",1]]],
                PRIMEM["Greenwich",0,
                    ANGLEUNIT["degree",0.0174532925199433]],
                CS[ellipsoidal,3],
                    AXIS["latitude",north,
                        ORDER[1],
                        ANGLEUNIT["degree",0.0174532925199433]],
                    AXIS["longitude",east,
                        ORDER[2],
                        ANGLEUNIT["degree",0.0174532925199433]],
                    AXIS["ellipsoidal height",up,
                        ORDER[3],
                        LENGTHUNIT["metre",1]],
                ID["EPSG",4979]]],
        ABRIDGEDTRANSFORMATION["unknown to WGS84 ellipsoidal height",
            METHOD["GravityRelatedHeight to Geographic3D"],
            PARAMETERFILE["Geoid (height correction) model file","@my.gtx",
                ID["EPSG",8666]]]]]

and-marsh · 2019-11-05T15:28:28Z

If the vertical CRS is known and there is a transformation in the EPSG dataset, then it will be used

Yes, you are right, it works out of box with PROJ 6.2.1 and GDAL 3.0.2(not sure that GDAL update is necessary for that), amazing! Thank you!

In PROJ master, the capability to add a GEOIDMODEL in WKT has been added per OSGeo/PROJ#1710.

Will it supported via GDAl-API or it is only PROJ business?

It was also possible to use the BOUNDCRS mechanism that gives a transformation between a VertCRS and a GeographicCRS using a grid

...
ABRIDGEDTRANSFORMATION["unknown to WGS84 ellipsoidal height",
            METHOD["GravityRelatedHeight to Geographic3D"],
            PARAMETERFILE["Geoid (height correction) model file","@my.gtx",
...

Should it be added manually or there is some functionality for this purpose? I didn't found any API to create or modify such BOUNDCRS in GDAL docs.

rouault · 2019-11-05T17:13:09Z

Will it supported via GDAl-API or it is only PROJ business?

Good question. Honestly this topic of custom grids is a non-trivial one when using a pure WKT based approach like PROJ does now, so we're a bit on the experiment edge for now I think.
You've several options: importFromProj4() with your +geoidgrids, but of course you loose the automated horizontal adjustments if you have a non-WGS84 datum, so you may also have to use your own +towgs84/+nadgrids. Or edit manually a WKT (or PROJJSON) to include the BoundCRS node. The downside of the BoundCRS approach that a user discovered is if you're VerticalCRS uses US Feet unit for example, then PROJ will think that the grid expects US Feet as input... . With the GEOIDGRID["PROJ my.gtx"], then it assumes that my.gtx expects metre and will do the appropriate vertical unit conversion

and-marsh · 2019-11-07T07:44:00Z

Ok, thank you, i am looking forward to new updates!

rouault · 2020-04-05T16:41:41Z

Closing as I think the discussion solved the issue

rouault added a commit to rouault/PROJ that referenced this issue Nov 3, 2019

createOperations(): in some circumstances we wrongly promoted a Helme…

5e5997a

…rt geog2D transformation to a geog3D Fixes for example EPSG:4979 to EPSG:2189, as raised in OSGeo/gdal#1972 (comment)

rouault mentioned this issue Nov 3, 2019

createOperations(): in some circumstances we wrongly promoted a Helmert geog2D transformation to a geog3D OSGeo/PROJ#1711

Merged

backporting bot pushed a commit to OSGeo/PROJ that referenced this issue Nov 3, 2019

createOperations(): in some circumstances we wrongly promoted a Helme…

97a468b

…rt geog2D transformation to a geog3D Fixes for example EPSG:4979 to EPSG:2189, as raised in OSGeo/gdal#1972 (comment)

rouault closed this as completed Apr 5, 2020

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

GDAL 3.0.1: coordinate transformations are inconsistent for different CRS types #1972

GDAL 3.0.1: coordinate transformations are inconsistent for different CRS types #1972

and-marsh commented Nov 1, 2019

rouault commented Nov 1, 2019

rouault commented Nov 3, 2019

and-marsh commented Nov 5, 2019

rouault commented Nov 5, 2019

and-marsh commented Nov 5, 2019

rouault commented Nov 5, 2019

and-marsh commented Nov 7, 2019

rouault commented Apr 5, 2020

GDAL 3.0.1: coordinate transformations are inconsistent for different CRS types #1972

GDAL 3.0.1: coordinate transformations are inconsistent for different CRS types #1972

Comments

and-marsh commented Nov 1, 2019

Expected behavior and actual behavior.

All the tests have been performed using GDAL 3.0.1 and Proj4 6.2.0:

1. 3D -> 3D

2. 3D -> 2D

3. 3D -> PROJCS

My questions

Operating system

GDAL version and provenance

rouault commented Nov 1, 2019

rouault commented Nov 3, 2019

and-marsh commented Nov 5, 2019

rouault commented Nov 5, 2019

and-marsh commented Nov 5, 2019

rouault commented Nov 5, 2019

and-marsh commented Nov 7, 2019

rouault commented Apr 5, 2020

All the tests have been performed using `GDAL 3.0.1` and `Proj4 6.2.0`: