Geog3D to Geog2D+GravityRelatedHeight (PL txt)

Formula

Note: These formulas have been transcribed from EPSG Guidance Note #7-2. Users are encouraged to use that document rather than the text which follows as reference because limitations in the transcription will be avoided.

This method transforms coordinates between a geographic 3D CRS and a compound CRS consisting of separate geographic 2D and vertical CRSs. The transformation is three-dimensional and reversible.

A transformation using this operation method consists of two parts:
a) The transformation of coordinates from the geographic 3D to or from the geographic 2D CRS;
b) The transformation of ellipsoidal height to or from gravity-related height.

However, complexities arise if the source 3D and horizontal component of the target compound CRS are based on different geodetic datums. To account for any offset between and non-parallelism of the ellipsoid surfaces, the method is entirely dependent upon a suitable three-dimensional transformation existing between the geographic 3D and geographic 2D CRSs, and this may not always be the case. For this reason the general case is considered to be out of scope here and this method is restricted to circumstances where the source geographic 3D CRS and the target geographic 2D CRS are both based on the same geodetic datum, i.e. the horizontal component of the compound CRS is the two-dimensional subset of the geographic 3D CRS. This geographic 2D CRS is also the CRS used for interpolation of the geoid model file grid.

Forward transformation from geographic 3D to compound
Because of the restriction above, latitude and longitude of the point in the source CRS do not change as a result of the transformation; therefore φt = φs and λt = λs, where the subscripts 's' and ‘t’ indicates the coordinates in the source and target CRSs respectively.

The vertical component of the forward transformation is as for the Geographic3D to GravityRelatedHeight (gtx) method, code 9665. The method requires that the ‘Interpolation CRS’ - the horizontal CRS to which the grid containing the geoid heights (ζ) are referenced - to be specified. That will be the geographic 2D CRS which is both the horizontal subset of the geographic 3D CRS and the horizontal component of the compound CRS. The unit of measure in which the geoid heights are given is part of the file format definition.

The forward case is straightforward:
φt = φs
λt = λs
Then using (φ, λ) for interpolation for the geoid height at the point,
H = h – ζ

Reverse transformation from compound to geographic 3D
For the reverse case the height offset ζ is first interpolated using the geographic 2D coordinates (φ, λ).
Then:
φt = φs
λt = λs
h = H + ζ

Example

For coordinate transformation: ETRS89 to ETRS89 + NAP height.
The Interpolation CRS is ETRS89 geographic 2D, CRS code 4258.

A point in the geographic 3D coordinate reference system ETRS89, code 4937, with:
	latitude 		φ = 51°59'10.80033"N	=  +51.986333425°
	longitude 		λ =   4°37'48.72315"E 	=    +4.630200875°
	ellipsoidal height 	h = 36.7595m 

is converted to the geographic 2D coordinate reference system ETRS89, code 4258, as: 
	point 			φ  = +51.986333426°	λ = +4.630200875°

To find its NAP height, first obtain the geoid height (ζ) at each of the surrounding grid nodes:
	NW node 		φ  = +51.9875°		λ = +4.6200°		ζ = +43.5376m
	NE node 		φ  = +51.9875°		λ = +4.6400°		ζ = +43.5398m
	SE node 		φ  = +51.9750°		λ = +4.6400°		ζ = +43.5479m
	SW node 		φ  = +51.9750°		λ = +4.6200°		ζ = +43.5455m

Use bi-linear interpolation in ETRS89 to obtain the geoid height (ζ):
	point			φ  = +51.986333426°	λ = +4.630200875°	ζ = +43.5395m

Then the NAP height  	H = h  –  ζ = 36.7595 − 43.5395 = −6.7800m

The coordinates in compound CRS "ETRS89 + NAP height", code 9283, are: 
	latitude 		φ  = 51°59'10.80033"N
	longitude 		λ  =   4°37'48.72315"E
	NAP height 		H = −6.7800m 

For the reverse transformation to find the ellipsoidal ETRS89 height, first obtain the geoid height (ζ) at each of the surrounding grid nodes and use bi-linear interpolation in ETRS89 to obtain the geoid height (ζ) at the point, as for the forward computation above:
	point			φ  = +51.986333426°	λ = +4.630200875°	ζ = +43.5395m

Then the ellipsoidal ETRS89 height h = H  +  ζ = −6.7800 + 43.5395 = 36.7595m

The coordinates in geographic 3D coordinate reference system ETRS89, code 4937, are: 
	latitude 		φ = 51°59'10.80033"N	= +51.986333425°
	longitude 		λ  =   4°37'48.72315"E 	=   +4.630200875°
	ellipsoidal height 	h = 36.7595m

EPSG:1100