BelgianLambert_Ellipsoidal_3D -> BelgianLambert_OstendHeight transformation should probably be implemented by doing BelgianLambert_Ellipsoidal_3D -> WGS 84 3D and WGS84 3D -> BelgianLambert_OstendHeight

[ktheijs]

What is the bug?

I would expect to get the same output (within error margin) when I do both transformations:
EPSG:4326+Ellipsoidal -> EPSG:31370+5710
EPSG:4326+Ellipsoidal -> EPSG:31370+Ellipsoidal -> EPSG:31370+5710

However for different ways of formatting the WKT, I get different results.
When making a simple derived version of EPSG:31370 I also get yet different results.
I am not sure if there are multiple bugs causing this or if it all has the same underlying cause.

In the case of transforming between two compound systems with vertical CRS
VERTCRS["Ellipsoid (Metre)",VDATUM["Ellipsoid"],CS[vertical,1],AXIS["ellipsoidal height (h)",up,LENGTHUNIT["metre",1]], it seems that no transformation of the Z value is done. Does the code assume it is the same vertical CRS? Even though the WKT of the vertical CRS is the same it should be treated differently because EPSG:31370 uses a different ellipsoid from EPSG:4326.

Steps to reproduce the issue

In the following zip file you will find the test class in C# and the WKT definitions as separate files for ease of use outside of C#.

GdalSupportCase.zip

In the output below I use:

• WGS84 = EPSG:4326+Ellipsoidal
• LambertEllipsoidal = EPSG:31370+Ellipsoidal
• LambertOstend = EPSG:31370+5710

Using 3D Ellipsoidal--------------------------------------
WGS84 input point:
5.73340278,50.60861389,299.16

WGS84 -> LambertOstend:
246588.1241007616,145103.98252358567,254.9427094303551

WGS84 -> LambertEllipsoidal:
246588.20546968083,145103.98057698458,255.04905407130718

WGS84 -> LambertEllipsoidal -> LambertOstend:
246588.2054696809,145103.98057698365,210.83176220202247

Using Compound Ellipsoidal--------------------------------
WGS84 input point:
5.73340278,50.60861389,299.16

WGS84 -> LambertOstend:
246588.1241007616,145103.98252358567,254.9427094303551

WGS84 -> LambertEllipsoidal:
246588.20546968083,145103.98057698458,299.16

WGS84 -> LambertEllipsoidal -> LambertOstend:
246588.2054696809,145103.98057698365,254.94270813071532

Using Derived---------------------------------------------
WGS84 input point:
5.73340278,50.60861389,299.16

WGS84 -> LambertOstend:
246598.20546968083,145103.98057698458,210.8317622020225

WGS84 -> LambertEllipsoidal:
246598.20546968083,145103.98057698458,299.16

WGS84 -> LambertEllipsoidal -> LambertOstend:
246598.2054696809,145103.98057698365,254.94270813071532

Versions and provenance

GDAL 3.10.0, released 2024/11/01
PROJ 9.5.0

As released on https://www.gisinternals.com/query.html?content=filelist&file=release-1930-x64-gdal-3-10-0-mapserver-8-2-2.zip

Additional context

No response

[jratike80]

Have you checked and compared with projinfo what Proj plans to do with different inputs/outputs? For example:

projinfo -s EPSG:4979 -t @BelgianLambert_Ellipsoidal_3D.txt -o proj
Operation No. 1:

unknown id, Inverse of BD72 to WGS 84 (3) + Belgian Lambert 72, 1 m, Belgium - onshore.

PROJ string:
+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=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=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

projinfo -s EPSG:4979 -t @BelgianLambert_Ellipsoidal_Compound.txt -o proj

unknown id, Inverse of Transformation from Ellipsoid (Metre) to WGS 84 (ballpark vertical transformation, without ellipsoid height to vertical height correction) + Inverse of Ballpark geographic offset from BD72 to WGS 84 + Belgian Lambert 72, unknown accuracy, World, has ballpark transformation

PROJ string:
+proj=pipeline
  +step +proj=axisswap +order=2,1
  +step +proj=unitconvert +xy_in=deg +xy_out=rad
  +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

[ktheijs]

Thanks for the tip.

Problem 1: When transforming between ellipsoidal coordinate systems with a compound WKT, it does seem to ignore the Z value. (see the push and pop of the Z)

C:\OSGeo4W>projinfo -s @Wgs84_Ellipsoidal_Compound.txt -t @BelgianLambert_Ellipsoidal_Compound.txt -o proj
Candidate operations found: 3
-------------------------------------
Operation No. 1:

unknown id, Inverse of BD72 to WGS 84 (3) + Belgian Lambert 72, 1 m, Belgium - onshore.

PROJ string:
+proj=pipeline
  +step +proj=axisswap +order=2,1
  +step +proj=unitconvert +xy_in=deg +xy_out=rad
  +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

As opposed to doing the same thing but with a 3D WKT representation.

C:\OSGeo4W>projinfo -s @Wgs84_Ellipsoidal_3D.txt -t @BelgianLambert_Ellipsoidal_3D.txt -o proj
Candidate operations found: 3
-------------------------------------
Operation No. 1:

unknown id, Inverse of BD72 to WGS 84 (3) + Belgian Lambert 72, 1 m, Belgium - onshore.

PROJ string:
+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=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=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

Problem 2: When transforming Wgs84_Ellipsoidal_Compound -> BelgianDerived_OstendHeight a different pipeline is used then when transforming Wgs84_Ellipsoidal_Compound -> BelgianLambert_OstendHeight. I would expect that the pipeline is the same for both, only the derived one would have an additional affine step at the end.

C:\OSGeo4W>projinfo -s @Wgs84_Ellipsoidal_Compound.txt -t @BelgianDerived_OstendHeight.txt -o proj
Candidate operations found: 12
-------------------------------------
Operation No. 1:

unknown id, Conversion from WGS 84 to WGS 84 + Inverse of BD72 to WGS 84 (3) + ETRS89 to Ostend height (1) using BD72 to ETRS89 (3) + Belgian Lambert 72 + Deriving conversion, 1.04 m, Belgium - onshore.

PROJ string:
+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=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=push +v_1 +v_2
  +step +proj=hgridshift +grids=be_ign_bd72lb72_etrs89lb08.tif +omit_inv
  +step +inv +proj=vgridshift +grids=be_ign_hBG18.tif +multiplier=1
  +step +inv +proj=hgridshift +grids=be_ign_bd72lb72_etrs89lb08.tif +omit_fwd
  +step +proj=pop +v_1 +v_2
  +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
  +step +proj=affine +xoff=10 +s11=1 +s12=0 +yoff=0 +s21=-0 +s22=1

C:\OSGeo4W>projinfo -s @Wgs84_Ellipsoidal_Compound.txt -t @BelgianLambert_OstendHeight.txt -o proj
Candidate operations found: 3
-------------------------------------
Operation No. 1:

unknown id, Conversion from WGS 84 to WGS 84 + 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.

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

The derived flow actually seems to make more sense to me. So would 210.8317622020225 be the correct height? This is the same result as when I do Wgs84_Ellipsoidal_3D -> BelgianLambert_Ellipsoidal_3D -> BelgianLambert_OstendHeight. It seems like the direct transformation from WGS84 to Belgian Lambert with Ostend Height doesn't work correctly. Or is it the other way around?

[rouault]

Transferred as a PROJ ticket.
I don't think there's any PROJ real bug here. The issue stands more with the fact that your _Ellipsoidal_Compound.txt files are not really conformant with the ISO19111 / OGC Abstract Topic 2 model that PROJ uses. An ellipsoidal height is not a valid vertical CRS according to it. So Wgs84_Ellipsoidal_3D.txt should be used instead of Wgs84_Ellipsoidal_Compound.txt. And similarly BelgianLambert_Ellipsoidal_3D.txt should be used instead of BelgianLambert_Ellipsoidal_Compound.txt

[jratike80]

Do the problematic WKT strings perhaps come from some ESRI software, see Esri/projection-engine-db-doc#47?

[ktheijs]

@rouault , I understood that PROJ made an excemption for this kind of VERTCRS. How else can you create a derived VERTCRS based on ellipsoidal heights?

A PROJCRS with 3 axes is also not conformant with the standard but still supported by PROJ, correct?
If so, then we still have an issue. Because Wgs84_Ellipsoidal_3D -> BelgianLambert_OstendHeight should be the same as Wgs84_Ellipsoidal_3D -> BelgianLambert_Ellipsoidal_3D -> BelgianLambert_OstendHeight but it is not.

[jratike80]

This is also discussed in #3075.

[ktheijs]

Thanks for the reference @jratike80 . That confirms that this is a bug then.

I'm also still not sure which elevation is correct now for Ostend Height. 254.94 or 210.83.

Wgs84_Ellipsoidal_3D -> BelgianLambert_OstendHeight -> Z=254.94
Wgs84_Ellipsoidal_3D -> BelgianLambert_Ellipsoidal_3D -> BelgianLambert_OstendHeight -> Z=210.83
Wgs84_Ellipsoidal_3D -> BelgianDerived_OstendHeight -> Z=210.83

[jratike80]

That confirms that this is a bug then.

That is probably opinion based, but I do not really understand these things (I am a GDAL user, not a Proj developer). Have you been reading the standard https://docs.ogc.org/is/18-010r11/18-010r11.pdf and considered how you would define your coordinate systems as WKT2?
I do not understand what you mean by "A PROJCRS with 3 axes is also not conformant with the standard but still supported by PROJ, correct?" Can you give some example?

[ktheijs]

I wasn't aware that a PROJCRS with 3 axes is conformant. But even better that it is. This means it should be supported by PROJ too, right?

The WKT definitions as in the files Wgs84_Ellipsoidal_3D.txt, BelgianLambert_OstendHeight.txt, BelgianLambert_Ellipsoidal_3D.txt and BelgianDerived_OstendHeight.txt should be conformant then. You can find them in the zip file attached to the original post. I don't want to copy them here as they are pretty long.

Since these definitions should be conformant, I don't see how the different transformations can get such different Z values.

[rouault]

Wgs84_Ellipsoidal_3D -> BelgianLambert_OstendHeight -> Z=254.94

I would assume that one to be correct from the PROJ pipeline. returned by projinfo

Wgs84_Ellipsoidal_3D -> BelgianLambert_Ellipsoidal_3D -> BelgianLambert_OstendHeight -> Z=210.83

That one is probably dubious because I'm not sure ellipsoidal heights are really defined in BD72. PROJ uses a Helmert transformation to do WGS 84 --> BD72 and it uses it in 3D to change the ellipsoidal heights. But then BelgianLambert_Ellipsoidal_3D -> BelgianLambert_OstendHeight is likely even more dubious because there's no known transformation between BD72 ellipsoidal heights and Ostend Height, so PROJ (likely wrongly) uses BD72 ellipsoidal heights as if they were ETRS89 ellipsoidal heights to use the be_ign_hBG18.tif grid

[ktheijs]

PROJ seems to follow the defined transformations when possible, but there is also a bunch of custom logic to try and makes sense of a transformation without formal definition. What is PROJ's filosophy regarding transformations that are not formally defined? Is this specific scenario something that PROJ should be able to do or not according to it's filosophy? It is hard to create applications on top of PROJ if we cannot know what PROJ will transform correctly and what not.

You also did not mention the Wgs84_Ellipsoidal_3D -> BelgianDerived_OstendHeight case. When deriving coordinate systems they typically are custom coordinate systems, so they don't have any formal definition. So can I not expect that a transformation to a derived CRS (based on a well defined CRS) should work?

[rouault]

but there is also a bunch of custom logic to try and makes sense of a transformation without formal definition. What is PROJ's filosophy regarding transformations that are not formally defined?

Not sure what you mean by "formally defined". There isn't much philosophy, but a set of heuristics more or less good depending on situations, possibly with error implementations at times. A simplified version of what is implemented is in https://proj.org/en/stable/operations/operations_computation.html, but it doesn't reflect everything that is implemented. I'm now thinking that BelgianLambert_Ellipsoidal_3D -> BelgianLambert_OstendHeight transformation should probably be implemented by doing BelgianLambert_Ellipsoidal_3D -> WGS 84 3D and WGS84 3D -> BelgianLambert_OstendHeight rather than what it does currently.

You also did not mention the Wgs84_Ellipsoidal_3D -> BelgianDerived_OstendHeight case

The pipeline of that one looks fine

[ktheijs]

If Wgs84_Ellipsoidal_3D -> BelgianLambert_OstendHeight -> Z=254.94 is correct. How can Wgs84_Ellipsoidal_3D -> BelgianDerived_OstendHeight -> Z=210.83 also be correct?

The derived horizontal system only moves the X position 10m to the East. It should not impact the elevation.

[rouault]

It should not impact the elevation.

it does not for me:

$ echo 50.60861389 5.73340278 299.16 | PROJ_NETWORK=ON bin/cs2cs "$(cat /home/even/gdal/gdal/build_cmake/GdalSupportCase/Wgs84_Ellipsoidal_3D.txt)" "$(cat /home/even/gdal/gdal/build_cmake/GdalSupportCase/BelgianLambert_OstendHeight.txt)"
246588.12	145103.98 254.94

$ echo 50.60861389 5.73340278 299.16 | PROJ_NETWORK=ON bin/cs2cs "$(cat /home/even/gdal/gdal/build_cmake/GdalSupportCase/Wgs84_Ellipsoidal_3D.txt)" "$(cat /home/even/gdal/gdal/build_cmake/GdalSupportCase/BelgianDerived_OstendHeight.txt)"
246598.12	145103.98 254.94

[ktheijs]

You are right. I was doing Wgs84_Ellipsoidal_Compound -> BelgianDerived_OstendHeight instead of Wgs84_Ellipsoidal_3D -> BelgianDerived_OstendHeight.

As you've mentioned before, Wgs84_Ellipsoidal_Compound is not a valid CRS definition.