Height not transformed from 3D promoted EPSG:4277 to EPSG:4979 if uk_os_OSTN15_NTv2_OSGBtoETRS.tif present

This is very similar to https://github.com/OSGeo/PROJ/issues/2225, but this time involving a horizontal grid and no compound coordinate reference systems.

Example of problem

Without the grid we get this:

echo 51.1788899 -1.8263548 147.885247288534 | cs2cs EPSG:4979 "$(projinfo EPSG:4277 --3d -o PROJJSON -q --single-line)" -f "%.12f"
51.178367122924 -1.824954404965 99.962219635025

and when uk_os_OSTN15_NTv2_OSGBtoETRS.tif can be found we get this:

echo 51.1788899 -1.8263548 147.885247288534 | cs2cs EPSG:4979 "$(projinfo EPSG:4277 --3d -o PROJJSON -q --single-line)" -f "%.12f"
51.178372136265 -1.824964570568 147.885247288534

Problem description

When the horizontal grid is used, the height is not transformed, but it should undergo the same Helmert transform that happens when no grid is available.

Expected Output

When the grid is found:

echo 51.1788899 -1.8263548 147.885247288534 | cs2cs EPSG:4979 "$(projinfo EPSG:4277 --3d -o PROJJSON -q --single-line)" -f "%.12f"
51.178372136265 -1.824964570568 99.962219635025

Environment Information

PROJ version 7.1.1
Ubuntu 18.04

Installation method

Compiled from source

@RogerLott May I ask you your opinion on the topic raised here ? That is a user wants to transform from OSGB 36 (EPSG:4277), promoted to a Geographic3D CRS, to WGS 84 Geographic3D CRS (EPSG:4979).

If we look at transformations available in EPSG, we have available a grid-based transformation through the OSTN15_NTv2_OSGBtoETRS.gsb (EPSG:7710 "OSGB 1936 to WGS 84 (9)", but this is 2D only. And we have also a Helmert-based one (like EPSG:1314 "OSGB 1936 to WGS 84 (6)"), that could possibly be used to operate on 3D coordinates.

Would it make sense (as suggested by the reporter of this issue) for a software to "hybrid" them, that is use the grid-based transformation for transforming the longitude and latitude (since it is the most precise for horizontal coordinates), and the Helmert-one to transform the ellipsoidal height, and combine the results ? As OSGB 36 has no official corresponding 3D CRS, this may be a bit questionable. I was going to add: "For the sake of the argument, let's say there would be a corresponding official 3D CRS", but browsing a bit through EPSG transformations, there are not so many situations where we have a source and target CRS that admit both a 2D and 3D geographic CRS, and have grid-based transformation and Helmert transformation. The only instance I could find is for GDA94 to GDA2020, but GDA94_GDA2020_conformal.gsb is documented to give identical results (for the 2D part) to Helmert transformation EPSG:8048 "GDA94 to GDA2020 (1)". So in that case, the "hybridation" would be equivalent to doing the Helmert transformation in 3D.

@RogerLott May I ask you your opinion on the topic raised here ? That is a user wants to transform from OSGB 36 (EPSG:4277), promoted to a Geographic3D CRS, to WGS 84 Geographic3D CRS (EPSG:4979).

If we look at transformations available in EPSG, we have available a grid-based transformation through the OSTN15_NTv2_OSGBtoETRS.gsb (EPSG:7710 "OSGB 1936 to WGS 84 (9)", but this is 2D only. And we have also a Helmert-based one (like EPSG:1314 "OSGB 1936 to WGS 84 (6)"), that could possibly be used to operate on 3D coordinates.

Would it make sense (as suggested by the reporter of this issue) for a software to "hybrid" them, that is use the grid-based transformation for transforming the longitude and latitude (since it is the most precise for horizontal coordinates), and the Helmert-one to transform the ellipsoidal height, and combine the results ?

I'm no expert in geodesy and these hybrids have always looked debatable to me, but on the other hand, something similar is done for geoids where the interpolation CRS is not the same as the target CRS of the transformation. I would argue that this request is as correct or wrong as using the Helmert transform for the horizontal part and the Geoid for the vertical, which PROJ would attempt for example, when transforming between EPSG:4326 and EPSG:7405 having the OSGM15 grid available but not the horizontal one (suboptimal, but possible). Anyway, in this particular case, I would say that if the vertical component is an important part of the deliverables of a user, most probably they would be using ODN for the vertical CRS

The difference between the situation in this issue and the one in #2225 is that in the later, to do "WGS 84 + EGM 96 height" to "CH1903+ promoted as 3D", one can chain 2 operations: "WGS 84 + EGM 96 height to WGS 84 3D" and "WGS 84 3D to CH1903+". So in ISO 19111 concepts, which PROJ (mostly) follows, this is a concatenated operation. In the case of that ticket, this is not really chaining existing operations but creating a new hybrid one, that would be a PROJ specific Transformation. It would give the following pipeline that back ups the longitude latitude, does the Helmert-based operation (doing the geocentric conversions), restores the original longitude and latitude and finally does the grid-based operation

echo 51.1788899 -1.8263548 147.885247 | PROJ_NETWORK=ON src/cct -d 8 +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=push +v_1 +v_2 \
  +step +proj=cart +ellps=WGS84 \
  +step +inv +proj=helmert +x=446.448 +y=-125.157 +z=542.06 +rx=0.15 +ry=0.247 \
        +rz=0.842 +s=-20.489 +convention=position_vector \
  +step +inv +proj=cart +ellps=airy \
  +step +proj=pop +v_1 +v_2 \
  +step +inv +proj=hgridshift +grids=uk_os_OSTN15_NTv2_OSGBtoETRS.tif \
  +step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m \
  +step +proj=axisswap +order=2,1

outputs 51.17837214 -1.82496457 99.96221935

(to answer my own observation above in https://github.com/OSGeo/PROJ/issues/2384#issuecomment-712317221 about some geographic 2D CRS not having official 3D counterparts, the Swiss example is an interesting example. Despite CH1903+ being 2D only in EPSG, in Swisstopo doc https://www.swisstopo.admin.ch/content/swisstopo-internet/en/home/meta/search.download/swisstopo-internet/en/documents/geo-documents/refsys_e.pdf it is used also as a Geographic 3D CRS. For example in paragraph 7.1 page 19 where there are ellipsoidal heights in CH1903+)

I see the problem of this transformation not being expressible as a pipeline concatenating existing transformation from the EPSG catalog.

That said, I would like to understand if my assumption about EPSG:4979 to EPSG:7405 doing a similar thing is wrong or not. With projinfo -s EPSG:4979 -t EPSG:4277+5701 -o PROJ (which is a simpler version of going to EPSG:7405 for that matter) I can see that there's no transformation that would use both the geoid and the Helmert transform. However, if I do

echo 51.1788899 -1.8263548 147.885247288534 | cs2cs EPSG:4979 EPSG:4277+5701 -f "%.12f"
51.178367118505 -1.824954393177 100.368278693984

I get a horizontal transformation that I don't know from where it comes because adding PROJ_DEBUG=2 I get:

Inverse of ETRS89 to WGS 84 (1) + ETRS89 to ODN height (2) + ETRS89 to WGS 84 (1) + Inverse of OSGB 1936 to WGS 84 (6)

but I can't find that transformation in the output from projinfo -s EPSG:4979 -t EPSG:4277+5701 -o PROJ --spatial-test intersects

That said, I would like to understand if my assumption about EPSG:4979 to EPSG:7405 doing a similar thing is wrong or not.

That's a different topic that the case you submitted in this ticket. EPSG:4979 to EPSG:7405 is Geographic3D to CompoundCRS transformation. A different code path in PROJ than Geographic3D to Geographic3D. With current master (and 7.1 branch as well), I do get the same result as you for cs2cs. But PROJ_DEBUG=3 shows "Using coordinate operation Inverse of ETRS89 to WGS 84 (1) + ETRS89 to ODN height (2) + Inverse of OSGB 1936 to ETRS89 (2)" , and that's exactly what projinfo -s EPSG:4979 -t EPSG:4277+5701 -o PROJ -spatial-test intersects shows as first operation.

OSGB 1936 (EPSG::4277) is a horizontal 2D CRS. In principle, no height can be described as part of a geographic 2D CRS (I will return to this below). "3D promoted" makes no sense to me. Where did the input 'height' come from? Is it a measurement or is it a coordinate? If it is a measurement then it has no place in a coordinate transformation: we should be discussing the pipeline from one 2D CRS to some other 2D CRS (in this case OSGB 1936 to WGS 84 2D [EPSG::4326]), with the measurement being an attribute of the location regardless of the horizontal CRS in which the position is expressed. If it is a coordinate, which vertical CRS? The coordinate triple is then a compound CRS.

Geographic 3D <> Compound transformations are straightforward when the horizontal component of the compound CRS is the geographic 2D subset of the geographic 3D CRS. They are reversible. But they are 3D transformations and require a geoid model to process the change from ellipsoidal height in the geographic 3D CRS to a gravity-related height. If the horizontal component of the compound CRS and the geographic 3D CRS are associated with different geodetic datums, things become complicated at best and often geodetically, but not mathematically, impossible. To retain integrity of the height in the transformation pipeline a necessary intermediate CRS is the 3D geographic CRS which is a superset of the geographic 2D CRS forming the horizontal component of the compound CRS. In the case of an old 'classical' CRS such as OSGB 1936 this may not exist. Furthermore the necessary 3D transformations may or may not be available. When only a 2D transformation is available the 2D to 3D step in a geog2D (CRS A) <> geog3D (CRS A) <> geocentric (CRS A) <> geocentric (CRS B) <> geog3D (CRS B) <> geog2D (CRS B) pipeline is indeterninate.

So how does the use of a Helmert transformation in the gecentric domain fit into a pipeline geog2D (CRS A) <> geog3D (CRS A) <> geocentric (CRS A) <> geocentric (CRS B) <> geog3D (CRS B) <> geog2D (CRS B) ? From a mathematical perspective the transformation in the geocentric domain is very clearly 3D. From the geodetic perspective this may not be true. It will depend upon whether geoid height has been properly accounted for during the determination of the transformation. In many cases the supposedly 3D transformation has been derived such that can only be accurately used for transformation of horizontal coordinates. Most of the time this is poorly documented. But whenever the source or target CRS has only a 2D form, it should be assumed to be a transformation using a 3D method in the 2D domain.

I suggest that attempting to do a transformation from 2D EPSG::4277 to 3D EPSG::4979 is entirely inappropriate. The only bug I see is in permitting it in the pipeline.

[As an aside to this issue as it may not be relevant (because the 'height' is unclear), EPSG will publishing reversible geog3D to compound transformations in the near future. But only for cases where the geog3D and compound horizontal are associated with the same geodetic datum. After some months of discussion and deliberation we decided that a case including geodetic datum change had too many uncertainties and failure cases to be generally viable, although certain specific instances might be].

Thanks for your comments @RogerLott

"3D promoted" makes no sense to me.

This is the act of constructing a Geographic 3D CRS from a Geographic 2D CRS by adding an ellipsoidal height. projinfo EPSG:4277 --3d -o WKT2:2019 outputs

WKT2:2019 string:
GEOGCRS["OSGB 1936",
    DATUM["OSGB 1936",
        ELLIPSOID["Airy 1830",6377563.396,299.3249646,
            LENGTHUNIT["metre",1]],
        ID["EPSG",6277]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8901]],
    CS[ellipsoidal,3],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]],
        AXIS["ellipsoidal height (h)",up,
            ORDER[3],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
    REMARK["Promoted to 3D from EPSG:4277"]]

This may not be relevant/appropriate for OSGB 1936, but as the user explictly asked for that promotion to 3D, PROJ honors this request. There is also the interesting case I mentioned in a later comment about CH1903+. CH1903+ only exists as a Geographic 2D CRS in EPSG, but in Swisstopo documentation, they definitely mention test points in CH1903+ with an ellipsoidal height, so in the case, allowing CH1903+ to be promoted as a Geographic 3D CRS seems to be legitimate.

Thank you for the explanation of the Proj behaviour. When I said "3D promoted makes no sense to me", it was not that I did not understand that a height was being added to a 2D CRS, but that in the specific circumstances of this issue, involving OSGB 1936, there is no geodetic logic in this action. Nobody has ever determined an ellipsoidal height for OSGB 1936 - it just does not exist in geodesy. For some old datums the the geoid height is defined at the origin point of the datum, but because of the way this datum was defined (fixed to mid-19th century 2D coordinates of 12 points across the British Isles) this does not apply here. There is no right answer. And it is difficult to say what a best answer might be. Just because software can be made to do it does not mean it is appropriate to do so. I can only reiterate that in my view attempting to do a transformation from 2D EPSG::4277 to 3D EPSG::4979 is entirely inappropriate. Appending either a measurement or a gravity-related height to a geog2D CRS does make sense.

@rouault said:

There is also the interesting case I mentioned in a later comment about CH1903+. CH1903+ only exists as a Geographic 2D CRS in EPSG, but in Swisstopo documentation, they definitely mention test points in CH1903+ with an ellipsoidal height, so in the case, allowing CH1903+ to be promoted as a Geographic 3D CRS seems to be legitimate.

If I understand correctly, CH1903+ is equivalent to LV95, i.e. a modern, 3D, GNSS based system - at least to some extent traceable to ITRF. So the 2D version really is a "demoted" version of the 3D thing. So in that case it makes plenty of sense to quote an ellipsoidal height in connection with the latitude/longitude-pair.

But as @RogerLott says, for the classical 2D OSGB 1936 datum, ellipsoidal heights make no (geodetic-geophysical) sense.

Assuming, however, that your OSGB coordinate is situated directly on the ellipsoid (i.e. h=0), or in other words, that the original 20th century coordinate has been reduced to the ellipsoid, it may make some (geodetic-computational) sense to carry on in 3D: At least, the ETRF-height obtained will give you some indication of the error induced by reducing to the modern ellipsoid (by demotion), vs. the original reduction (by necessity) to the Airy ellipsoid with orientation defined by the OSGB36 triangulation.

In other words: That which is without meaning may make sense anyway. IIRC, Messrs. Poder and Engsager shared some thoughts on this in their Some Conformal Mappings... (1998).

Thanks for your insightful replies @RogerLott

Geographic 3D <> Compound transformations are straightforward when the horizontal component of the compound CRS is the geographic 2D subset of the geographic 3D CRS. They are reversible. But they are 3D transformations and require a geoid model to process the change from ellipsoidal height in the geographic 3D CRS to a gravity-related height. If the horizontal component of the compound CRS and the geographic 3D CRS are associated with different geodetic datums, things become complicated at best and often geodetically, but not mathematically, impossible. To retain integrity of the height in the transformation pipeline a necessary intermediate CRS is the 3D geographic CRS which is a superset of the geographic 2D CRS forming the horizontal component of the compound CRS.

This is good to know, because when I started learning about this I was scratching my head about why some software allows to combine coordinate systems like CH1903+ with geoid-powered vertical CRSs that use a different underlying geographic CRS like EGM96 and that didn't make any sense to me. This confirms my suspicion. I wonder if there's any good way to know which is the base 2D geographic CRS that a given vertical CRS can be compound with without having to implement complex heuristics that query the EPSG registry and the available grid files, because otherwise is very easy to mix and match CRS that don't make sense.

In fact, the transformation that I was requesting falls in this category of transformations that make little or no sense and users should avoid. I reserve a question regarding this to the end.

In the case of an old 'classical' CRS such as OSGB 1936 this may not exist.

Thanks for pointing this out, this is an important detail that I was naively neglecting. If there's a grid it's obviously because OSGB 1936 is pre-GNSS and the triangulated network contains tensions. That said, guessing that the 12 key benchmarks of the network were positioned using astronomical observations and are still based on the Airy 1830 ellipsoid, would it be incorrect to say that these points have an ellipsoidal height? Or is this CRS defined by a single control point and all the other benchmark coordinates were obtained by traditional triangulation on the projected CRS space? (which means that the grid was obtained by first unprojecting EPSG:27700 coordinates to EPSG:4277)

So how does the use of a Helmert transformation in the gecentric domain fit into a pipeline geog2D (CRS A) <> geog3D (CRS A) <> geocentric (CRS A) <> geocentric (CRS B) <> geog3D (CRS B) <> geog2D (CRS B) ? From a mathematical perspective the transformation in the geocentric domain is very clearly 3D. From the geodetic perspective this may not be true. It will depend upon whether geoid height has been properly accounted for during the determination of the transformation. In many cases the supposedly 3D transformation has been derived such that can only be accurately used for transformation of horizontal coordinates. Most of the time this is poorly documented. But whenever the source or target CRS has only a 2D form, it should be assumed to be a transformation using a 3D method in the 2D domain.

I see. This answers a question that I had written after reading the previous paragraph and then removed. Let me rephrase it in my own words to see if I understood it correctly. In this particular case, I get that OSGB 1936 to ETRS is the only transformation available in the EPSG catalog, and that one is purely 2D to 2D as it uses a grid, no question about extending it to 3D because it's ill-defined. Then, in those instances where a 2D-2D transformation is based on a Helmert transformation (and there's no 3D one), it implies that the 2D transformation is only adjusted for the "horizontal" surface defined by the region of interest or even by just a few control points and there's no guarantee that points above (or below) that surface will project on the correct 2D location when the Z coordinate is discarded. I still assume that between geographic 3D CRSs that are basically ellipsoidal (e.g. CHTRS95 and CH1903+), a Helmert transformation has good accuracy on altitude as well, is this correct?

I suggest that attempting to do a transformation from 2D EPSG::4277 to 3D EPSG::4979 is entirely inappropriate. The only bug I see is in permitting it in the pipeline.

My original concern about this particular case was not whether it made sense (which now it's clear that it doesn't) to use OSGB 1936 in 3D, but whether in an application that deals with point clouds, which are inherently 3D, anybody would try to use this CRS. I need to transform some data from WGS84 3D to the point cloud CRS, and for the UK, my opinion is that using OSGB 1936 / British National Grid + ODN height should always be the preference, for which transformations to/from WGS84 are well defined, but I don't know if this the well established practice or not.

This "bug" report was triggered when we realized that there was something incorrect when trying to convert 3D WGS84 coordinates to this extension of OSGB 1936 to 3D. What you are stating is that essentially, this shouldn't be done. However, in a practical and general sense I find it quite limiting as there are many CRS for which there's no 3D equivalent (Even gave the example of CH1903+). But at the same time, it also means that transformations to/from 3D WGS84 involving a CRS registered as 2D in EPSG don't necessarily make sense (but again, they may make sense in same cases as Even pointed out).

By the way, what's the situation with NAD83 (EPSG:4269), which unlike more modern realizations doesn't have a 3D equivalent?

Coming back to OSGB 1936, is it reasonable to assume that no professional surveyor will ever try to use OSGB 1936 projected or geographic (EPSG:27700 or EPSG:4277) with input data that is essentially 3D (even if the deliverables needed are simply 2D) and they will stick to the compound with ODN height (EPSG:7405)?

@rouault after all this discussion my impression is that this ticket can be closed, right?

OSGeo / PROJ