Geoid loaded but not recognized in custom EPSG crs

[WaltDixie]

The vertical datum "Ostend Height" is predefined in the following CRS predefined in Metashape:

Code: [Select]

crs=Metashape.CoordinateSystem('EPSG::8370')
crs_equivalent=Metashape.CoordinateSystem('EPSG::3812+5710')
The vertical datum EPSG code is 5110, the associated horizontal & vertical CS EPSG code is 3812/5710. So the above two crs are equivalent.

I downloaded the geoid from https://github.com/OSGeo/PROJ-data/tree/master/be_ign
I adapted the TIFFTAG_IMAGEDESCRIPTION so that MEtashape recognizes the geoid as a vertical datum with epsg code 5110.
I'm able to do a coordinate transformation from lonlat towards this system.

Code: [Select]

crsref=Metashape.CoordinateSystem('EPSG::4326')
coord = Metashape.Vector([2,50,10]) #lonlatz
crsref.transform(coord,crsref,crs) #works
crsref.transform(coord,crsref,crs_equivalent) #works, same output

However, if I define a custom crs with the same vertical crs EPSG code (5710), running a similar transformation returns the error "Vertical datum missing".

Code: [Select]

crscustom=Metashape.CoordinateSystem('EPSG::32631+5710')
crsref.transform(coord,crsref,crscustom) #fails, reports 'Vertical datum missing'

My expectation is that it would not matter which horizontal CRS is used for the vertical datum to be recognized.
I can't detect any difference in the wkt string.
So, can you explain what is missing here?

[WaltDixie]

Anybody that can help with this issue?

[Paulo]

Hello WaltDixie,

I guess since the Ostend height Vertical datum is defined referred to ETRS89 geographic system, you must transform to a projected CS with ETRS89 geographic CS and thus instead of 'EPSG::32631+5710' or 'WGS84/ UTM zone 31N + Ostend height', you should use 'EPSG::25831+5710'  or  'ETRS89 / UTM zone 31N + Ostend height' and it works as in transforming point 5, 51, 70 (lon,lat,h) from WGS84 to ETRS89/ UTM 31 + Ostend height:

Code: [Select]

wgs84 = ps.CoordinateSystem("EPSG::4326")
crs1 = ps.CoordinateSystem("EPSG::32631+5710")
crs2 = ps.CoordinateSystem("EPSG::25831+5710")

wgs84,crs1,crs2
Out[10]: 2023-03-22 12:33:48
2023-03-22 12:33:48 (<CoordinateSystem 'WGS 84 (EPSG::4326)'>,
2023-03-22 12:33:48  <CoordinateSystem 'WGS 84 / UTM zone 31N + Ostend height'>,
2023-03-22 12:33:48  <CoordinateSystem 'ETRS89 / UTM zone 31N + Ostend height'>)

crs1.transform(ps.Vector((5,51,70)),wgs84,crs1)
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Input In [11], in <cell line: 1>()
----> 1 crs1.transform(ps.Vector((5,51,70)),wgs84,crs1)
RuntimeError: Vertical datum missing

crs1.transform(ps.Vector((5,51,70)),wgs84,crs2)
Out[12]: 2023-03-22 12:38:29 Vector([640333.2963369302, 5651728.682650693, 26.972396499421023])
Of course practically ETRS89 and WGS84 are the same but Ostend is defined with ETRS89 so all transformations to Ostend ht must use a ETRS89 based geo system....