Partially-conformal variations of the Standard Molodensky datum transformation

Abstract:

Standard Molodensky is a recognised method of transforming coordinates between geodetic datums. Although less accurate than some other methods, it has the merit of being direct. That is to say it can be applied to geodetic coordinates, without involving Cartesian coordinates that give rise to difficulties in computing latitude. This paper considers the use of Standard Molodensky when at least one of the datums is 2D+1D in nature, meaning that that horizontal and vertical positions are obtained by different methods. This was generally the case before 3D positioning by satellites and is a widespread characteristic of local datums that are still used. The 2D+1D property weakens the argument for 3D conformality, and invites the possibility that different translation parameters might be used for horizontal and vertical shifts. The possibility of including a Z-rotation as a 7th parameter is also considered. Besides being ideal for those who favour the simplicity of Standard Molodensky, the variations introduced in this paper offer significant improvements in accuracy such as error reductions of 75%, 69% and 99% in the three selected case studies.

Keywords:
Standard Molodensky; datum transformations; ordinary least-squares; partial conformality; inverse transformations

1. Introduction

A geodetic datum is a reference frame to facilitate precise positioning. It must define (either directly or indirectly) a reference ellipsoid of a given size and shape, its centre, the direction of its polar axis, and the zero meridian.

Geodetic datums have developed over time for different countries and regions. When dealing with geodetic coordinates, there is a variety of reference ellipsoids, some chosen to fit the local geoid, others chosen arbitrarily. The advent of satellites has given rise to the concept of best-fitting geocentric reference ellipsoid. The version adopted by the International Union of Geodesy and Geophysics (IUGG) in 1979 was the Geodetic Reference System 1980 (GRS80) described in Moritz (2000).

GRS80 includes a reference ellipsoid for dynamic datums/reference frames based on the International Terrestrial Reference System. The corresponding reference frames are named after the year in which the frame started (eg ITRF92, ITRF2008, ITRF2014). They allow station positions to reflect the effect of plate tectonic motion, which is done by positional functions of time elapsed since the start time. Sometimes a dynamic datum gives rise to a static realisation at a particular epoch within a few years of the start time; an example is Geocentric Datum of Australia (GDA94) which is a static realisation of ITRF92 at epoch 1994.0.

The World Geodetic System 1984 (WGS84) has itself become a series of dynamic or time-varying datums, with each re-definition designed to be consistent to a few centimetres with the latest ITRF. Before the full implementation of NAVSTAR Global Positioning System (GPS) and other global navigation satellite systems (GNSS), 3D positioning of points on the earth’s surface was problematic. In particular, datums were ‘2D+1D’ in their use, in the sense that geodetic latitudes and longitudes were obtained in an entirely different manner to heights. Originally orthometric heights (obtained by levelling) were used directly; later, ellipsoidal heights were used, converted from orthometric heights by the application of a geoid model. For a time, this was even true in ITRFs or realisations of ITRS.

Transformation of coordinates between geodetic datums, generally referred to as datum transformations, are an ongoing tool for position determination. A growing number of countries have sophisticated models in accessible software that provide a high degree of accuracy, eg Brazil’s ProGriD (Santos et al. 2012). Such models, however, require a lot of observed data. Many countries (for example, Brunei, Ghana, Italy and Sudan) are not at that stage.

The range of possible models that are used as datum transformations is considerable. Varga, Grgić and Bašić (2017) describes and compares several direct approaches applied to Croatian data. Grgić, Varga and Bašić (2016) describes and compares 12 distortion-modelling methods applied to the same data, creating a grid that enables bilinear interpolation. A comprehensive comparison of transformation methods is beyond the scope of this paper. However, it is worth noting that multi-parameter models derived from small datasets could be introducing characteristics that reflect noise rather than physical reality. In some cases, notably when polynomials are involved, this can cause instability. Limited data and the drawbacks of over-elaborate models are factors which explain the continued use of simple transformations, in addition to ease of computation.

The notation used for ellipsoidal quantities is near-universal (a, e ² , f, b, v, ρ for semi-major axis, eccentricity squared, flattening, semi-minor axis, prime-vertical radius of curvature and meridional radius of curvature). The mathematical relationships can be found, for example, in Deakin (2004). Ellipsoidal quantities are quoted with subscripts s and t when it is necessary to distinguish between the ellipsoids of the source and target datums. In the context of source-to-target transformations, the prefix Δ is shorthand for ‘target minus source’; eg Δa = a _t - a _s . References to ‘linear units’ usually indicate metres and can be assumed to be such in this paper.

One of the simplest types of datum transformation is the one with just 3 shift parameters, namely ΔX, ΔY, ΔZ, applied to global Cartesian coordinates. Derived sets of 3 parameters for transforming local geodetic datums to the first realisation of WGS84 can be found in NIMA (2004). The differences between the respective ellipsoids in semi-major axes and flattening, ∆a and ∆f, are not normally counted as unknown parameters to be estimated from field data as they are derived from the definition of source and target reference ellipsoids. However, they need to be included in the model. Conceptually, this is a conformal transformation. To be precise, it is a translation since there is no change in scale or in the direction of the axes. The rigorous implementation of this transformation consists of 3 stages:

• conversion of geodetic coordinates to Cartesian coordinates in the source datum;

• application of the shifts ΔX, ΔY and ΔZ to the Cartesian coordinates;

• conversion of Cartesian coordinates to geodetic coordinates in the target datum.

The final stage is more difficult than the first, the problem being the computation of latitude. The traditional method involves iteration and is described, with the rest of the 3-stage algorithm, in Ruffhead (2016). Research around the world has produced a wide variety of non-iterative methods, many of which are listed in Featherstone and Claessens (2008).

The most direct methods for computing the 3-parameter datum transformation are those associated with Mikhail Sergeevich Molodensky (1909-1991). Changes in latitude, longitude and ellipsoidal height are computed as linear combinations of ΔX, ΔY, ΔZ, Δa and Δf. The version originating directly from Molodensky, Eremeev and Yurkina (1962) consists of the ‘Standard Molodensky’ formulae, given in subsection 2.1 below. This method is an approximate implementation of the 3-parameter conformal transformation. A computation process to make it an exact implementation was proposed in Ruffhead (2016).

The question arises, however, whether 3D conformal transformations are appropriate for datums that are 2D+1D in character. This paper considers the scope for improving accuracy in the Standard Molodensky algorithm in cases where at least one of the datums is 2D+1D. The conformality requirement is applied to the horizontal components with one set of translation parameters, a different set being used to transform ellipsoidal heights.

A further aspect of investigation is the extent to which accuracy is further improved by a 7th parameter - rotation about the Z axis. This paper demonstrates that it is easily derived and easily applied.

2. Existing methods

It has been noted that there are a wide variety of transformation methods. The ones described here are Standard Molodensky formulae and other forms of Molodensky formulae. The standard version is the starting point for the 6-parameter variations proposed in this paper. The other forms are mentioned for completeness, and for the simple modification that allows a Z-rotation parameter.

2.1 Standard Molodensky formulae

Datum transformations based on three shift parameters ΔX, ΔY, ΔZ are often computed (approximately) by Standard Molodensky formulae, documented in Molodensky, Eremeev and Yurkina (1962). They are explicitly recommended by NIMA (2004) for most applications involving local geodetic datums and WGS84.

Although NIMA (2004) is a much-quoted source of specific Standard Molodensky transformations to WGS84, it should not be regarded as definitive. The shift parameters were computed from data that preceded GPS and other GNSS. Also, ‘WGS84’ in that context refers to the first realisation in 1987.

A detailed derivation of Standard Molodensky formulae can be found in Deakin (2004). Δϕ and Δλ are given below in radians, but can easily be converted to arc-seconds or degrees. The equation for Δh is given in linear units.

$\Delta \varphi =[-\Delta X\sin {\varphi }_s\cos {\lambda }_s-\Delta Y\sin {\varphi }_s\sin {\lambda }_s+\Delta Z\cos {\varphi }_s+\Delta a({\nu }_s{e}_s^2\sin {\varphi }_s\cos {\varphi }_s)/a_s+\Delta f({{\rho }_{s}a}_s/b_s+{\nu }_s{b}_s/a_s\left)\sin {\varphi }_s\cos {\varphi }_s\right]/({\rho }_s+h_s)$ (1)

$\Delta \lambda =(-\Delta X\sin {\lambda }_s+\Delta Y\cos {\lambda }_s)/\left[\left({\nu }_s+h_s\right)\cos {\varphi }_s\right]$ (2)

$\Delta h=\Delta X\cos {\varphi }_s\cos {\lambda }_s+\Delta Y\cos {\varphi }_s\sin {\lambda }_s+\Delta Z\sin {\varphi }_s-\Delta a(a_s/{\nu }_s)+\Delta f(b_s/a_s){\nu }_s{\sin }^2{\varphi }_s$ (3)

Equations (1) and (2) can be converted from radians into the same linear units used in (11):

$({\rho }_s+h_s)\Delta \varphi =-\Delta X\sin {\varphi }_s\cos {\lambda }_s-\Delta Y\sin {\varphi }_s\sin {\lambda }_s+\Delta Z\cos {\varphi }_s+\Delta a\left({\nu }_s{e}_s^2\sin {\varphi }_s\cos {\varphi }_s\right)/a_s+\Delta f({\rho }_s{a}_s/b_s+{\nu }_s{b}_s/a_s)\sin {\varphi }_s\cos {\varphi }_s$ (4)

$\left({\nu }_s+h_s\right)\Delta \lambda \cos {\varphi }_s=-\Delta X\sin {\lambda }_s+\Delta Y\cos {\lambda }_s$ (5)

Equations (4), (5) and (3) can be regarded as the equations of eneralizat geodetic shifts given that they are in terms of ϕ, λ and h, but with the same linear units.

2.2 Other forms of Molodensky formulae

It should be noted that Molodensky, Eremeev and Yurkina (1962) contains implementations of more general conformal transformations. The additional parameters are a change in scale ΔS and small rotations R _X , R _Y , R _Z in the global Cartesian axes. Details are given in Soler (1976). In practice these have rarely been used, probably because the resulting generalisations of equations (1)-(3) are far more complicated. It is worth noting, however, that the rotation term R _Z is easy to accommodate: it only needs to be added to the right-hand side of (2) (assuming it complies with the position-vector sign convention).

Molodensky, Eremeev and Yurkina (1962) also has a generalization which varies Δa and Δf. This is, of course, outside the scope of transformations between ellipsoids of known size and shape.

There are simplified versions of the standard form that were not proposed by Molodensky. The best known of these is the set of ‘Abridged Molodensky’ formulae, which do not use ellipsoidal heights. The formulae are given in (6)-(8), in radians, radians and linear units respectively.

$\Delta \varphi =[-\Delta X\sin {\varphi }_s\cos {\lambda }_s-\Delta Y\sin {\varphi }_s\sin {\lambda }_s+\Delta Z\cos {\varphi }_s+(a_s\Delta f+f_s\Delta a\left)\sin {2\varphi }_s\right]/{\rho }_s$ (6)

$\Delta \lambda =(-\Delta X\sin {\lambda }_s+\Delta Y\cos {\mathrm{\lambda }}_{\mathrm{s}})/\left({\nu }_s\cos {\varphi }_s\right)$ (7)

$\Delta h=\Delta X\cos {\varphi }_s\cos {\lambda }_s+\Delta Y\cos {\varphi }_s\sin {\lambda }_s+\Delta Z\sin {\varphi }_s+(a_s\Delta f+f_s\Delta a){\sin }^2{\varphi }_s-\Delta a$ (8)

The Abridged formulae are shown alongside those of Standard Molodensky in Seppelin (1974), which says the Molodensky formulae had been ‘previously selected for use with WGS66’. DMA (1987, section 7.2.4.3.2) noted that ‘the Abridged Molodensky Datum Transformation Formulas have been used more extensively that the Standard Formulas’ but on grounds of accuracy firmly recommends the latter in ‘any new software involving the Molodensky Datum Transformation Formulas’. Despite this advice, Abridged Molodensky is still widely used (see, for example, IOGP 2020 and Yun, Lee and Song 2020). Varga, Grgić and Bašić (2017) found that both Abridged and Standard Molodensky fitted Croatian horizontal data as precisely as the 3-stage 3-parameter conformal transformation (which they call “Molod-3p”).

3. New methods

This section proposes new variations on the Standard Molodensky formulae. It also describes how to compute the inverse formulae. The particular least-squares method chosen for this study to optimise the parameters is also described.

3.1 New variations on Standard Molodensky formulae

Partially-conformal variation (PCV) of Standard Molodensky is a new method with up to 7 unknown parameters, not counting Δa or Δf. It is applied to geodetic coordinates. A convenient abbreviation for the method is SM-PCV. The full set of parameters consists of:

• translation parameters ΔX _hor , ΔY _hor , ΔZ _hor for use in the horizontal-shift equations (1) and (2);

• rotation parameter R _Z (as per the position vector convention) for application to the longitude after the Standard Molodensky longitude formula;

• translation parameters ΔX _ver , ΔY _ver , ΔZ _ver for use in the vertical-shift equation (3).

These are not to be confused with the 7 parameters of the 3D conformal transformation (which include 3 rotation parameters and a scale change). If all 7 SM-PCV parameters are used, the transformation formulae are

$\Delta \varphi =[-\Delta X_{hor}\sin {\varphi }_s\cos {\lambda }_s-\Delta Y_{hor}\sin {\varphi }_s\sin {\lambda }_s+\Delta Z_{hor}\cos {\varphi }_s+\Delta a({\nu }_s{e}_s^2\sin {\varphi }_s\cos {\varphi }_s)/a_s+\Delta f({\rho }_s{a}_s/b_s+{\nu }_s{b}_s/a_s\left)\sin {\varphi }_s\cos {\varphi }_s\right]/({\rho }_s+h_s)$ (9)

$\Delta \lambda =R_Z+(-\Delta X_{hor}\sin {\lambda }_s+\Delta Y_{hor}\cos {\lambda }_s)/\left[\right({\nu }_s+h_s\left)\cos {\varphi }_s\right]$ (10)

$\Delta h=\Delta X_{ver}\cos {\varphi }_s\cos {\lambda }_s+\Delta Y_{ver}\cos {\varphi }_s\sin {\lambda }_s+\Delta Z_{ver}\sin {\varphi }_s-\Delta a(a_s/{\nu }_s)+\Delta f(b_s/a_s){\nu }_s{\sin }^2{\varphi }_s$ (11)

Equations (9) and (10) are in radians but can be converted into the same linear units used in (11). The resulting equations (12) and (13), when coupled with equation (11), describe homogenised geodetic shifts.

$({\rho }_s+h_s)\Delta \varphi =-\Delta X_{hor}\sin {\varphi }_s\cos {\lambda }_s-\Delta Y_{hor}\sin {\varphi }_s\sin {\lambda }_s+\Delta Z_{hor}\cos {\varphi }_s+\Delta a\left({\nu }_s{e}_s^2\sin {\varphi }_s\cos {\varphi }_s\right)/a_s+\Delta f({\rho }_s{a}_s/b_s+{\nu }_s{b}_s/a_s)\sin {\varphi }_s\cos {\varphi }_s$ (12)

$({\nu }_s+h_s)\Delta \lambda \cos {\varphi }_s={({\nu }_s+h_s)R}_Z\cos {\varphi }_s-\Delta X_{hor}\sin {\lambda }_s+\Delta Y_{hor}\cos {\lambda }_s$ (13)

In equations (10) and (13), R _Z is in radians. If R _Z is disregarded, meaning no allowance for Z-rotation, these equations become a 6-parameter variation.

Equations (12) and (13) have some similarity to the approach used by Molnár and Timár (2005). They scaled the longitude difference Δλ by ϕ _s , but the relative weight of their expressions for Δϕ and Δλcosϕ _s made no allowance for differences in radii of curvature. Moreover, their equations were for Abridged Molodensky rather than Standard Molodensky, and were entirely for the purpose of ptimizing the horizontal transformation without regard for ellipsoidal heights.

Equations (9) to (11) are insensitive to errors in the ellipsoidal heights. This can easily be shown by test computations: a change of 10 metres in h _s will have a negligible effect on Δϕ, Δh and Δh, although h _t will experience an equivalent change. This characteristic of Standard Molodensky is shared by the 3-stage 3-parameter method. For both methods, a rough approximation of h _s should not affect the horizontal accuracy of the transformation.

3.2 Computing the inverse transformation

Given the parameters for the forward transformation (ϕ _s , λ _s , h _s ) → (ϕ _t , λ _t , h _t ), there remains the issue of how to compute the reverse (or inverse) transformation (ϕ _t , λ _t , h _t ) → (ϕ _s , λ _s , h _s ).

The formulae (9)-(11) can be used in reverse, but a modification is needed to cover the 7-parameter variation when R _Z is non-zero. In the forward transformation, R _Z is applied to the longitude after the translations; therefore, in the reverse transformation, -R _Z must be applied to the longitude before application of the reverse translations. Intermediate longitudes λ _t ^U (where U stands for ‘unrotated’) are computed using

${\lambda }_t^{\mathrm{U}}={\lambda }_t-R_Z$ (14)

The ‘simple inverse’ (ϕ _s ^SI , λ _s ^SI , h _s ^SI ) is obtained by computing the following shifts and adding them to ϕ _t , λ _t , h _t . (${\varphi }_S^{SI}={\varphi }_t+{\Delta \varphi }^{SI}$, etc.) If h _t is not known, a dummy value should be assigned so that the resulting h _s ^SI represents the same point.

${\Delta \varphi }^{\mathrm{S}\mathrm{I}}=[\Delta X_{hor}\sin {\varphi }_t\cos {\lambda }_t^{\mathrm{U}}+\Delta Y_{hor}\sin {\varphi }_t\sin {\lambda }_t^{\mathrm{U}}-\Delta Z_{hor}\cos {\varphi }_t-\Delta a({\nu }_t{e}_t^2\sin {\varphi }_t\cos {\varphi }_t)/a_t-\Delta f({\rho }_t{a}_t/b_t+{\nu }_t{b}_t/a_t\left)\sin {\varphi }_t\cos {\varphi }_t\right]/({\rho }_t+h_t)$ (15)

${\Delta \lambda }^{\mathrm{S}\mathrm{I}}=(\Delta X_{hor}\sin {\lambda }_t^{\mathrm{U}}-\Delta Y_{hor}\cos {\lambda }_t^{\mathrm{U}})/\left[\right({\nu }_t+h_t\left)\cos {\varphi }_t\right]$ (16)

${\Delta h}^{\mathrm{S}\mathrm{I}}=-\Delta X_{ver}\cos {\varphi }_t\cos {\lambda }_t^{\mathrm{U}}-\Delta Y_{ver}\cos {\varphi }_t\sin {\lambda }_t^{\mathrm{U}}-\Delta Z_{ver}\sin {\varphi }_t+\Delta a(a_t/{\nu }_t)-\Delta f(b_t/a_t){\nu }_t{\sin }^2{\varphi }_t$ (17)

If formulae (9)-(11) are applied to (ϕ _s ^SI , λ _s ^SI , h _s ^SI ), with those coordinates used to derive ρ _s and v _s , there will be a misclosure that can be as much as 0.050 metres in 3D distance. The misclosure (ϕ _t ^SMV - ϕt, λ _t ^SMV - λ _t , h _t ^SMV - h _t ) is a good indication of the error in the simple inverse. (SMV denotes the application of Standard Molodensky Variation to the simple inverse.)

If the user wishes a corrected inverse (CI) which is more exact, the misclosure can be applied in much the same way as in Sections 2.4, 2.5 and 2.13 of Ruffhead and Whiting (2020).

${\varphi }_s^{\mathrm{C}\mathrm{I}}={\varphi }_s^{\mathrm{S}\mathrm{I}}-\left({\varphi }_t^{\mathrm{S}\mathrm{M}\mathrm{V}}-{\varphi }_t\right)$ (18)

${\lambda }_s^{\mathrm{C}\mathrm{I}}={\lambda }_s^{\mathrm{S}\mathrm{I}}-\left({\lambda }_t^{\mathrm{S}\mathrm{M}\mathrm{V}}-{\lambda }_t\right)$ (19)

$h_s^{\mathrm{C}\mathrm{I}}=h_s^{\mathrm{S}\mathrm{I}}-\left(h_t^{\mathrm{S}\mathrm{M}\mathrm{V}}-h_t\right)$ (20)

In the 6-parameter variation (where R _Z is zero), both the forward and reverse transformations involve the same amount of computation as the 3-parameter Standard Molodensky. This is because ΔX is simply replaced by ΔX _hor or ΔX _ver , with similar replacements occurring for ΔX and ΔY. In the reverse transformation, λ _t ^U will be λ _t since reverse application of R _Z is not needed.

For the 7-parameter variation, the only extra computation in the forward transformation is the addition of R _Z in equation (10). The only extra computation in the reverse transformation is the subtraction of R _Zin equation (14).

3.3 Computation of parameters

For the purposes of this paper, least-squares optimisation means ‘ordinary’ least-squares (possibly weighted). It is worth noting that the alternative of ‘total’ least-squares has become increasingly popular for deriving datum transformations (See, for example, Fang 2014 and Tao, Hua and Feng 2020.) Papers that compare the two methods include Pan et al. (2015) and Laari, Ziggah and Annan (2016).

The traditional method of optimising parameters ΔX, ΔY and ΔZ for Standard Molodensky has been mean Cartesian shifts. See, for example, NIMA (2004, p.7-4). Mean Cartesian shifts are calculated from the Cartesian coordinates of the control points in the source and target datums. In terms of least-squares, these are the optimum parameters for the 3-parameter conformal transformation applied to (X _S , Y _S , Z _S ) which is not quite the same as the two Molodensky transformations.

The rigorous method for optimising ΔX, ΔY and ΔZ for Standard Molodensky is to use equations (4), (5) and (3) for each control point, with residual terms (which like the equations themselves will be in linear units). Using x to denote the vector containing the parameters (in this case ΔX, ΔY and ΔZ), the observation equations for the n control points will take the form

$\mathbf{A}\mathbf{x}+\mathbf{v}=\mathbf{b}$ (21)

The 3n ´ 3 design matrix A contains the coefficients of ΔX, ΔY and ΔZ for each equation and each point. Vector b contains the 3n ‘observed’ quantities arising from equations (4), (5) and (3). These are

$\left({\rho }_s+h_s\right)\Delta \varphi -\Delta a\left({\nu }_s{e}_s^2\sin {\varphi }_s\cos {\varphi }_s\right)/a_s-\Delta f({\rho }_s{a}_s/b_s+{\nu }_s{b}_s/a_s)\sin {\varphi }_s\cos {\varphi }_s$

$\left({\nu }_s+h_s\right)\Delta \lambda \cos {\varphi }_s$

$\Delta h+\Delta a(a_s/{\nu }_s)-\Delta f(b_s/a_s){\nu }_s{\sin }^2{\varphi }_s$

for each point. Vector v consists of the 3n residuals.

The quantity to be minimised can be expressed as v ^T Wv, where W is a square symmetric matrix of weights. Ideally, W should consist of variances and covariances, to reflect variations in data quality and correlation. In practice, W is frequently taken to be the identity matrix.

As shown, for example, in Bomford (1980) and Cross (1983), the parameter vector x which minimises v ^T Wv is the solution of the normal equations A ^T WAx = A ^T Wb.

The author has derived Standard Molodensky parameters for several datasets (Great Britain, Western Australia, Sweden, Ghana) and found them to be within 0.025 metres of the mean datum shifts computed from Cartesian shifts. NIMA’s technically-incorrect approach is therefore excusable.

The main advantage of homogenised geodetic shifts is the separation of horizontal shifts from shifts in height. This allows one set of parameters ΔX _hor , ΔY _hor , ΔZ _hor to be used for the horizontal shift and another set ΔX _ver , ΔY _ver , ΔZ _ver to be used for the vertical shift. The former can be derived from the equations for Δϕ and Δλ, and the latter from the equations for Δh. This is the key to deriving the parameters in the partially-conformal variations of Standard Molodensky.

For the Standard Molodensky variations, there are 6 parameters to be optimised with the possibility of a systematic shift (R _Z ) in longitude as a 7th. The parameter vector x to be used in equation (21) will consist of ΔX _hor , ΔY _hor , ΔZ _hor , R _Z , ΔX _ver , ΔY _ver and ΔZ _ver . The elements of the 3n ´ 7 design matrix A will consist of the coefficients of those terms arising from the homogenised equations (12), (13) and (11). Vector b will contain the same 3n quantities as before.

4. Results

Each of the transformations considered involved at least one local geodetic datum. The datasets are relatively old (circa 1990s) and vary in size from 20 to 82 common points. They were selected on the basis of what was available to the author. Each set of coordinates is 2D+1D in nature, with orthometric heights converted to ellipsoidal heights via a geoid model. The latter cannot be identified in all cases.

The Swedish and Australian datasets have been used by Andrei (2006) and Awange, Bae and Claessens (2008) to derive datum transformations, although the method considered in their studies was affine transformations.

For this study, Standard Molodensky and the new variations were each derived by ordinary least squares optimisation using unit weights. No weighting information was available, and indeed both Andrei (2006) and Awange, Bae and Claessens (2008) assumed unit weights in their studies. As in those investigations, no points were set aside as testing points and root-mean-square (RMS) residuals were derived from all the data points.

The quality-of-fit of the SM-PCVs is compared primarily with Standard Molodensky to demonstrate a large improvement without increased computation. The parameters and RMSs are given to greater precision than would normally be expected given the uncertainty of the calculation (1 decimal place of a metre, 3 decimal places of an arc-second). This is purely to assist checks by software developers.

Comparisons are made with Bursa-Wolf. This is a simplification of the 7-parameter conformal transformation applied to Cartesian coordinates (see, for example, Ruffhead and Whiting 2020). It is indirect in the sense it requires conversions of source geodetic coordinates to Cartesians and Cartesians to target geodetic coordinates. There is no attempt to model two different types of conformality. The accuracy advantage of Bursa-Wolf over direct methods is greatly reduced when compared to SM-PCVs.

The author also programmed partially-conformal variations on Abridged Molodensky. Given the arguments in 1.3 against the continued use of Abridged Molodensky, the results are given in lesser detail, but sufficient to show that the PCVs give near-identical improvements.

4.1 Case study 1 (Western Australia)

This study involves datum transformations from Australian Geodetic Datum 1984 (AGD84) to Geocentric Datum of Australia 1994 (GDA94). The area of application is Western Australia. There are 82 data points known in both datums and their distribution is illustrated in Figure 1.

The points are the 82 stations of the STATEFIX GPS network, described in Agustan and Featherstone (2004) and more fully in Stewart et al. (1997). The coordinates were computed in 1996 by the Geodesy Group at Curtin University of Technology from observations provided by the Western Australian Department of Land Administration (later superseded by Landgate).

From the age of the data, it is clear that neither set of ellipsoidal heights could have come directly from GNSS surveys. Enquiries made to Landgate failed to reveal the precise methods by which orthometric heights were converted to AGD84 ellipsoidal heights, although it seems likely that the conversion used the 1971 astrogeodetic datum produced by the Division of National Mapping in Canberra. For the GDA94 heights, the geoid model used was AUSGeoid98.

The optimised parameters and the statistical properties of the fit are given in Table 1. In this instance, the extra translation parameters have a much bigger impact than the rotation parameter.

Figure 1:
Data points of case study 1 (Western Australia)

Table 1:
Statistics of optimised models for case study 1.

Parameters and properties	Std Mol (3P)	6P SM-PCV	7P SM-PCV
ΔX hor	-140.488 m	-129.902 m	-133.451 m
ΔY hor	-33.932 m	-52.064 m	-54.210 m
ΔZ hor	142.251 m	153.218 m	153.254 m
R Z			-0.151112²
ΔX ver	-140.488 m	-127.212 m	-127.212 m
ΔY ver	-33.932 m	-30.226 m	-30.226 m
ΔZ ver	142.251 m	134.653 m	134.653 m
Δϕ RMS	2.2917 m	0.5388 m	0.5052 m
Δλ RMS	1.7678 m	0.5926 m	0.5721 m
Δh RMS	1.5223 m	0.2760 m	0.2760 m
2D RMS	2.8943 m	0.8009 m (72.3%)	0.7632 m (73.6%)
3D RMS	3.2703 m	0.8472 m (74.1%)	0.8116 m (75.2%)

• Percentage values indicate improvement (reduction) with respect to Std Mol values.

The variations on Standard Molodensky have a similar accuracy advantage over the 3-stage 3-parameter conformal transformation, since the optimised form of the latter has a horizontal RMS of 2.8946 m and a 3D RMS of 3.2705 m. The variations also have the simplicity advantage because the geodetic coordinates do not need to be converted to and from Cartesian coordinates.

The variations on Standard Molodensky have a simplicity advantage over Bursa-Wolf and Molodensky-Badekas which are also 3-stage methods requiring conversion to and from Cartesian coordinates. They are, however, slightly more accurate than 7P SM-PCV since they have a horizontal RMS of 0.7142 m and a 3D RMS of 0.7583 m.

Least-squares models were also derived for Abridged Molodensky and the corresponding 6- and 7-parameter variations. The 6-parameter variation reduces the 3D RMS distance residual by 74.1% and the horizontal RMS distance residual by 72.3%. For the 7-parameter variation, the corresponding reductions are 75.2% and 73.6%.

4.2 Case study 2 (Great Britain)

This study involves datum transformations from Ordnance Survey Great Britain 1936 (OSGB 36) to WGS84. The area of application is Great Britain and there are 44 data points known in both datums. Their distribution is illustrated in Figure 2. The combined dataset was received by a colleague of the author circa 2001, so it is not recent. The ‘WGS84’ designation is clearly an early realisation of the datum.

Advice from the UK Defence Geographic Centre was that the points are Doppler stations and that the WGS84 coordinates came from Doppler precise positioning, not GPS. Geoid models must have been used to convert orthometric heights to ellipsoidal heights. (It was not possible to establish which of the many possibilities those models were.) This suggests that the ellipsoidal heights were much less accurate than the latitudes and longitudes.

Figure 2:
Data points of case study 2 (Great Britain)

The optimised parameters and the statistical properties of the fit are given in Table 2. In this instance, the extra translation parameters have a bigger impact than the rotation parameter, but the latter makes a significant contribution.

Table 2:
Statistics of optimised models for case study 2.

Parameters and properties	Std Mol (3P)	6P SM-PCV	7P SM-PCV
ΔX hor	376.414 m	453.370 m	452.520 m
ΔY hor	-111.291 m	-114.524 m	-134.223 m
ΔZ hor	431.660 m	538.810 m	538.793 m
R Z			1.091748²
ΔX ver	376.414 m	369.571 m	369.571 m
ΔY ver	-111.291 m	-156.683 m	-156.683 m
ΔZ ver	431.660 m	434.664 m	434.664 m
Δϕ RMS	7.5257 m	1.9818 m	1.6032 m
Δλ RMS	2.7466 m	1.9220 m	1.6039 m
Δh RMS	1.5964 m	1.0872 m	1.0872 m
2D RMS	8.0112 m	2.7608 m (65.5%)	2.2678 m (71.7%)
3D RMS	8.1687m	2.9671 m (63.7%)	2.5149 m (69.2%)

• Percentage values indicate improvement (reduction) with respect to Std Mol values.

The variations on Standard Molodensky have a similar accuracy advantage over the 3-stage 3-parameter conformal transformation, since the optimised form of the latter has a horizontal RMS of 8.0146 m and a 3D RMS of 8.1720 m. The variations also have the simplicity advantage because the geodetic coordinates do not need to be converted to and from Cartesian coordinates.

The variations on Standard Molodensky have a simplicity advantage over Bursa-Wolf and Molodensky-Badekas which are also 3-stage methods requiring conversion to and from Cartesian coordinates. They do, however, broadly match 7P SM-PCV for accuracy since they have a horizontal RMS of 2.2522 m and a 3D RMS of 2.5196 m.

Least-squares models were also derived for Abridged Molodensky and the corresponding 6- and 7-parameter variations. The 6-parameter variation reduces the 3D RMS distance residual by 63.6% and the horizontal RMS distance residual by 65.5%. For the 7-parameter variation, the corresponding reductions are 69.2% and 71.7%.

4.3 Case study 3 (Sweden)

This study involves datum transformations from SWEREF93 (The Swedish realisation of ETRS89) to a local reference coordinate system designated RT90/RH70 by Andrei (2006). The latter is ‘a mixture of the Swedish triangulation network RT90 and the 2nd Swedish precise levelling network RH70’. The geoid model associated with RT90 is RN02 (Engberg and Lilje, 2001). There are 20 data points known in both datums. Their distribution is illustrated in Figure 3.

Figure 3:
Data points of case study 3 (Sweden)

The geocentric Cartesian coordinates are given in Table 4.1 of Andrei (2006), but with columns in the wrong order. For a correct listing, see Table 5 of Amiri-Simkooei (2018).

The optimised parameters and the statistical properties of the fit are given in Table 3. In this instance, the rotation parameter had a much greater overall impact than the 3 extra translation parameters, although the latter caused a dramatic reduction in the ∆h RMS residual.

The variations on Standard Molodensky have a similar accuracy advantage over the 3-stage 3-parameter conformal transformation, since the optimised form of the latter has a horizontal RMS of 12.6148 m and a 3D RMS of 13.9134 m. The variations also have the simplicity advantage because the geodetic coordinates do not need to be converted to and from Cartesian coordinates.

The variations on Standard Molodensky have a simplicity advantage over Bursa-Wolf and Molodensky-Badekas which are also 3-stage methods requiring conversion to and from Cartesian coordinates. They are slightly more accurate than 7P SM-PCV, since they have a horizontal RMS of 0.1296 m and a 3D RMS of 0.1796 m, but the massive advantage they have over Standard Molodensky is largely eliminated.

Table 3:
Statistics of optimised models for case study 3.

Parameters and properties	Std Mol (3P)	6P SM-PCV	7P SM-PCV
ΔX hor	-498.396 m	-502.211 m	-471.993 m
ΔY hor	36.640 m	35.655 m	-66.133 m
ΔZ hor	-563.431 m	-569.957 m	-569.643 m
R Z			7.134725²
ΔX ver	-498.396 m	-416.328 m	-416.328 m
ΔY ver	36.640 m	-99.283 m	-99.283 m
ΔZ ver	-563.431 m	-585.556 m	-585.556 m
Δϕ RMS	5.1664 m	4.8139 m	0.1115 m
Δλ RMS	11.5042 m	11.6407 m	0.1057 m
Δh RMS	5.8693 m	0.1293 m	0.1293 m
2D RMS	12.6111 m	12.5968 m (0.1%)	0.1536 m (98.8%)
3D RMS	13.9100 m	12.5975 m (9.4%)	0.2008 m (98.6%)

• Percentage values indicate improvement (reduction) with respect to Std Mol values.

It must be admitted that the 7-parameter variation is an improvement on Standard Molodensky because a transformation limited to translation parameters is highly unsuitable for the Swedish dataset. The rotation R _Z between zero meridians is a very significant part of the transformation from SWEREF93 to RT90/RH70.

Least-squares models were also derived for Abridged Molodensky and the corresponding 6- and 7-parameter variations. The 6-parameter variation reduces the 3D RMS distance residual by 9.4% and the horizontal RMS distance residual by 0.1%. For the 7-parameter variation, the corresponding reductions are 98.6% and 98.8%.

5. Discussion

While the place of Standard Molodensky in the field of datum transformations is declining, its appeal as a simple direct method has led to its continued use, particular in countries where high-accuracy surface-fitting models are not yet available. It is therefore worth examining variations that preserve its simplicity whilst improving its accuracy.

In all three case studies, the 7-parameter variations improve the accuracy of Standard Molodensky by more than two-thirds, in terms of the reduction in residuals. The reasons for this vary from the alternative translation terms in Cases 1 and 2 to the rotation term in Case 3. The improvement in the 6-parameter variation over Standard Molodensky is substantial in Cases 1 and 2, but not in Case 3.

While these Standard Molodensky variations are not necessarily more accurate than the 3-stage methods which use Cartesian coordinates (eg Bursa-Wolf), they are simpler to apply and are an important addition to the direct methods. The application involves no more computation than the basic Standard Molodensky, apart from the addition of a constant term in the 7-parameter case.

The horizontal-accuracy comparisons have particular importance in cases where ellipsoidal heights are ignored (either because of unavailability or suspected inaccuracy). The 6-parameter SM-PCV then becomes 2D Standard Molodensky using horizontally-derived ΔX, ΔY, ΔZ as opposed to mean Cartesian shifts; the reductions in RMS horizontal distance residual were 72.3% and 65.5% in the first two case studies, although negligible in the third where the main source of error has not been addressed. The 7-parameter SM-PCV becomes rotationally-modified 2D Standard Molodensky using horizontally-derived ΔX, ΔY, ΔZ; the reduction in RMS horizontal distance residual is a massive 98.8% in the third case study.

Direct SM-PCV methods are worth investigating for transformations between local datums and modern realisations of global datums, where the case for full 3D conformality does not apply. The possibility of a change in Z-rotation is worth pursuing not only for its possible significance but also because it is easy to include in the observation equations and even easier to apply.

The primary conclusion can be stated as follows. On the evidence of the case studies in this paper, those who advocate the Standard Molodensky transformation for its simplicity should consider the partially-conformal variations as an equally-simple but more accurate alternative. The key limitation is that SM-PCV methods are only appropriate when at least one of the datums is of type 2D+1D.

As with any transformation method derived from a single application to the entire set of data points, the RMS of the residuals from SM-PCV methods will not give a true indication of their accuracy as a predictor. If such information is required within the region covered by the data points, there is always the option of cross-validation, a process described succinctly in Berrar (2018).

ACKNOWLEDGMENTS

The author thanks Professor Joseph Awange of Curtin University for providing control data for Western Australia, and Ordnance Survey for providing control data for Great Britain. The author also thanks the reviewers for their constructive comments.

REFERENCES

• Agustan and Featherstone, W. E. 2004. Reprocessing the Western Australian Statefix Geodetic GPS Network Using Commercial Software Paper presented at the 3rd FIG Regional Conference, Jakarta, Indonesia, October 3-7, 2004. 14 pp.
• Amiri-Simkooei, A.R. 2018. Parameter estimation in 3D affine and similarity transformation: implementation of variance component estimation. Journal of Geodesy, 92 (11), pp.1285-1297. https://doi.org/10.1007/s00190-018-1119-1
  » https://doi.org/10.1007/s00190-018-1119-1
• Andrei, C. O. 2006. 3D affine coordinate transformations MSc thesis. School of Architecture and the Built Environment. Royal Institute of Technology (KTH), Stockholm, Sweden.
• Awange, J.L., Bae, K.-H., and Claessens, S.J. 2008. Procrustean solution of the 9-parameter transformation problem. Earth, Planets and Space, 60 (6), pp.529-537. https://doi.org/10.1186/BF03353115
  » https://doi.org/10.1186/BF03353115
• Berrar, D. 2019. Cross-Validation. In: Encyclopedia of Bioinformatics and Computational Biology Elsevier, pp.542-545. https://doi.org/10.1016/B978-0-12-809633-8.20349-X
  » https://doi.org/10.1016/B978-0-12-809633-8.20349-X
• Bomford, G. 1980. Least Squares: Parametric method. In: Geodesy (fourth edition). Oxford, Clarendon Press. pp.707-711.
• Cross, P. A. 1983. The combined least squares process. In: Advanced Least Squares Applied to Position-fixing Working Paper No. 6, Department of Land Surveying, North East London Polytechnic. Ch. 3.
• Deakin, R. E. 2004. The Standard and Abridged Moldensky [sic] Coordinate Transformation Formulae Department of Mathematical and Geospatial Sciences, RMIT University, Melbourne, Australia, 22 pp.
• DMA, 1987. Supplement to Department of Defense World Geodetic System 1984 Technical Report: Part I - Methods, Techniques, and Data Used in WGS84 Development, DMA TR 8350.2-A, National Imagery and Mapping Agency, Washington, WA, USA.
• Engberg, L. E. and Lilje, M. 2001. A New Co-ordinate System for Sweden TS5.5 Reference Frame in Practice, FIG XXII International Congress Washington, D.C. USA, April 19-26 2002.
• Fang, X. 2014. A total least squares solution for geodetic datum transformations. Acta Geodaetica et Geophysica, 49 (2), pp.189-207. https://doi.org/10.1007/s40328-014-0046-8
  » https://doi.org/10.1007/s40328-014-0046-8
• Featherstone, W.E. and Claessens, S.J. 2008. Closed-form transformation between geodetic and ellipsoidal coordinates. Studia Geophysica et Geodaetica, 52 (1), pp.1-18. https://doi.org/10.1007/s11200-008-0002-6
  » https://doi.org/10.1007/s11200-008-0002-6
• Grgić, M., Varga, M., and Bašić, T. 2016. Empirical research of interpolation methods in distortion modeling for the coordinate transformation between local and global geodetic datums. Journal of Surveying Engineering, 142 (2), 05015004. https://doi.org/10.1061/(ASCE)SU.1943-5428.0000154
  » https://doi.org/10.1061/(ASCE)SU.1943-5428.0000154
• IOGP, 2020, Geomatics Guidance Note 7, part 2: Coordinate Conversions & Transformations including Formulas. Report 373-7-2. International Association of Oil & Gas Producers (IOGP) Report 373-7-2. Available at: <Available at: https://www.iogp.org/wp-content/uploads/2019/09/373-07-02.pdf > [Accessed 21 May 2022].
  » https://www.iogp.org/wp-content/uploads/2019/09/373-07-02.pdf
• Laari, P. B. Ziggah, Y. Y. and Annan, R. F. 2016. Determination of 3D transformation parameters for the Ghana geodetic reference network using ordinary least squares and total least squares techniques. International Journal of Geomatics and Geosciences, 7 (3), pp.245-261.
• Molnár, G. and Timár, G. 2005. Determination of the parameters of the abridging Molodensky formulae providing the best horizontal fit. Geophysical Research Abstracts, 7 (01018).
• Molodensky, M. S. Eremeev, V. F. and Yurkina, M. I. 1962. Methods for Study of the External Gravitational Field and Figure of the Earth Translated from Russian. Israel Programme for Scientific Translations.
• Moritz, H. 2000. Geodetic Reference System 1980. Journal of Geodesy , 74 (1), pp.128-133. https://doi.org/10.1007/s001900050278
  » https://doi.org/10.1007/s001900050278
• NIMA, 2004. Department of Defense World Geodetic System 1984: Its definition and relationships with local geodetic systems 3rd ed. Technical Report no. 8350.2, National Imagery and Mapping Agency (NIMA), Washington, WA, USA. (Date on cover page is 2000, but 2004 is date of enclosed amendment.) Available at: <Available at: https://gis-lab.info/docs/nima-tr8350.2-wgs84fin.pdf > [Accessed 15 April 2022].
  » https://gis-lab.info/docs/nima-tr8350.2-wgs84fin.pdf
• Pan, G., Zhou, Y., Sun, H., and Guo, W. 2015. Linear observation based total least squares. Survey Review, 47 (340), pp.18-27. https://doi.org/10.1179/1752270614Y.0000000090
  » https://doi.org/10.1179/1752270614Y.0000000090
• Ruffhead, A.C. 2016. The SMITSWAM method of datum transformations consisting of Standard Molodensky in two stages with applied misclosures. Survey Review , 48 (350), pp.376-384. https://doi.org/10.1080/00396265.2016.1191748
  » https://doi.org/10.1080/00396265.2016.1191748
• Ruffhead, A. C. and Whiting, B. M. 2020. Introduction to geodetic datum transformations and their reversibility Surveying Working Paper No 01/2020, School of Architecture, Computing and Engineering, University of East London. 30 pp.
• Santos, M.F., Santos, M.C., Oliveira, L.C., Costa, S.A., Azevedo, J.B., and Galo, M. 2012. ProGriD: The transformation package for the adoption of SIRGAS2000 in Brazil. In: S. Kenyon, M.C. Pacino, and U. Marti, eds. Geodesy for Planet Earth International Association of Geodesy Symposia, 136, Springer, Berlin, Heidelberg, pp.869-875. https://doi.org/10.1007/978-3-642-20338-1_109
  » https://doi.org/10.1007/978-3-642-20338-1_109
• Seppelin, T. O. 1974. The Department of Defense World Geodetic System 1972. Presented at the International Symposium on Problems Relating to the Redefinition of North American Geodetic Networks, Fredericton, New Brunswick, Canada, 20-25 May 1974. [Prepared by the World Geodetic System Committee, DMA, Washington, DC.]
• Soler, T. 1976. On differential transformations between cartesian and curvilinear (geodetic) coordinates. Department of Geodetic Science Report No. 236, Ohio State University. 89 pp.
• Stewart, M. P. Ding, X. L. Tsakiri, M. and Featherstone, W. E. 1997. 1996 STATEFIX Project, Final Contract Report to the Department of Land Administration, School of Spatial Sciences, Curtin University of Technology. 100 pp.
• Tao, W., Hua, X., and Feng, S. 2020. Performing 3D similarity transformation with large rotation angles using constrained multivariate total least squares. Boletim de Ciencias Geodesicas, 26 (4), e2020021. https://doi.org/10.1590/s1982-21702020000400021
  » https://doi.org/10.1590/s1982-21702020000400021
• Varga, M., Grgić, M., and Bašić, T. 2017. Empirical comparison of the Geodetic Coordinate Transformation Models: a case study of Croatia. Survey Review , 49 (352), pp.15-27. https://doi.org/10.1080/00396265.2015.1104092
  » https://doi.org/10.1080/00396265.2015.1104092
• Yun, S., Lee, H., and Song, J. 2020. Studies on derivation of appropriate geodetic system transformation schemes for spatial data. Journal of the Korean Society of Surveying, Geodesy, Photogrammetry and Cartography, 38 (6), pp.561-571 [in Korean with English abstract]. https://doi.org/10.7848/KSGPC.2020.38.6.561
  » https://doi.org/10.7848/KSGPC.2020.38.6.561

Publication Dates

• Publication in this collection
  12 Dec 2022
• Date of issue
  2022

History

• Received
  03 Jan 2022
• Accepted
  11 Sept 2022