Ellipsoidal terrain correction based on multi-cylindrical equal-area map projection of the reference ellipsoid

An operational algorithm for computation of terrain correction (or local gravity field modeling) based on application of closed-form solution of the Newton integral in terms of Cartesian coordinates in multi-cylindrical equal-area map projection of the reference ellipsoid is presented. Multi-cylindrical equal-area map projection of the reference ellipsoid has been derived and is described in detail for the first time. Ellipsoidal mass elements with various sizes on the surface of the reference ellipsoid are selected and the gravitational potential and vector of gravitational intensity (i.e. gravitational acceleration) of the mass elements are computed via numerical solution of the Newton integral in terms of geodetic coordinates {,,h}. Four base- edge points of the ellipsoidal mass elements are transformed into a multi-cylindrical equal-area map projection surface to build Cartesian mass elements by associating the height of the corresponding ellipsoidal mass elements to the transformed area elements. Using the closed-form solution of the Newton integral in terms of Cartesian coordinates, the gravitational potential and vector of gravitational intensity of the transformed Cartesian mass elements are computed and compared with those of the numerical solution of the Newton integral for the ellipsoidal mass elements in terms of geodetic coordinates. Numerical tests indicate that the difference between the two computations, i.e. numerical solution of the Newton integral for ellipsoidal mass elements in terms of geodetic coordinates and closed-form solution of the Newton integral in terms of Cartesian coordinates, in a multi-cylindrical equal-area map projection, is less than 1.6×10^–8 m²/s² for a mass element with a cross section area of 10×10 m and a height of 10,000 m. For a mass element with a cross section area of 1×1 km and a height of 10,000 m the difference is less than 1.5×10^–4m²/s². Since 1.5× 10^–4 m²/s² is equivalent to 1.5×10^–5m in the vertical direction, it can be concluded that a method for terrain correction (or local gravity field modeling) based on closed-form solution of the Newton integral in terms of Cartesian coordinates of a multi-cylindrical equal-area map projection of the reference ellipsoid has been developed which has the accuracy of terrain correction (or local gravity field modeling) based on the Newton integral in terms of ellipsoidal coordinates.Acknowledgments. This research has been financially supported by the University of Tehran based on grant number 621/4/859. This support is gratefully acknowledged. The authors are also grateful for the comments and corrections made to the initial version of the paper by Dr. S. Petrovic from GFZ Potsdam and the other two anonymous reviewers. Their comments helped to improve the structure of the paper significantly.     

Hyperboloidal coordinates: transformations and applications in special constructions

The vector-based algorithms for biaxial and triaxial ellipsoidal coordinates presented by Feltens (J Geod 82:493–504, 2008; 83:129–137, 2009) have been extended to hyperboloids of one sheet. For the backward transformation from Cartesian to hyperboloidal coordinates, of two iterative process candidates one was identified to be well suited. It turned out that a careful selection of the center of curvature is essential for the establishment of a stable and reliable iteration process. In addition, for zero hyperboloidal heights a closed solution is presented. The hyperboloid algorithms are again based on simple formulae and have been successfully tested for various theoretical hyperboloids. The paper concludes with a practical application example on a cooling tower construction.     

Vector methods to compute azimuth, elevation, ellipsoidal normal, and the Cartesian ( X , Y , Z ) to geodetic ( φ , λ , h ) transformation

Vector-based algorithms for the computation of azimuth, elevation and the ellipsoidal normal unit vector from 3D Cartesian coordinates are presented. As a by-product, the formulae for the ellipsoidal normal vector can also be used to iteratively transform rectangular Cartesian coordinates (X, Y, Z) into geodetic coordinates (φ, λ, h) for a height range from −5600 km to 10⁸ km. Comparisons with existing methods indicate that the new transformation can compete with them.     

Transforming ellipsoidal heights and geoid undulations between different geodetic reference frames

Transforming height information that refers to an ellipsoidal Earth reference model, such as the geometric heights determined from GPS measurements or the geoid undulations obtained by a gravimetric geoid solution, from one geodetic reference frame (GRF) to another is an important task whose proper implementation is crucial for many geodetic, surveying and mapping applications. This paper presents the required methodology to deal with the above problem when we are given the Helmert transformation parameters that link the underlying Cartesian coordinate systems to which an Earth reference ellipsoid is attached. The main emphasis is on the effect of GRF spatial scale differences in coordinate transformations involving reference ellipsoids, for the particular case of heights. Since every three-dimensional Cartesian coordinate system ‘gauges’ an attached ellipsoid according to its own accessible scale, there will exist a supplementary contribution from the scale variation between the involved GRFs on the relative size of their attached reference ellipsoids. Neglecting such a scale-induced indirect effect corrupts the values for the curvilinear geodetic coordinates obtained from a similarity transformation model, and meter-level apparent offsets can be introduced in the transformed heights. The paper explains the above issues in detail and presents the necessary mathematical framework for their treatment. An erratum to this article can be found at     

Vector method to compute the Cartesian (X, Y, Z) to geodetic ({\phi}, λ, h) transformation on a triaxial ellipsoid

The vector-based algorithm to transform Cartesian (X, Y, Z ) into geodetic coordinates (, λ, h) presented by Feltens (J Geod, 2007, doi:) has been extended for triaxial ellipsoids. The extended algorithm is again based on simple formulae and has successfully been tested for the Earth and other celestial bodies and for a wide range of positive and negative ellipsoidal heights.     

Helmert's projection of a ground point onto the rotational reference ellipsoid in topocentric cartesian coordinates

Summary In order to derive the ellipsoidal height of a point P_t, on the physical surface of the earth, and the direction of ellipsoidal normal through P_t, we presente here an iterative procedure rapidly convergent to compute, in a topocentric Cartesian system, the coordinates of Helmert's projection of the ground point P_t onto the reference ellipsoid of revolution .We derive as well the cofactor matrix of total vector of the topocentric coordinates of the above ground point and of its Helmert's projection so that to compute the variance of ellipsoidal height.     

Alignment of geodetic and satellite coordinate systems to the average terrestrial system

The parallelism of geodetic and satellite systems to the average terrestrial system is examined, under the assumption that a geodetic system is a fixed framework invariant with respect to geodetic network adjustment. In this case a geodetic system is rotated with respect to the average terrestrial system only about the ellipsoid normal of the initial point. The method is demonstrated using coordinates and covariance matrices for BC-4 and SECOR satellite tracking stations computed by Mueller and his co-workers. It is shown that the NAD geodetic system is scaled significantly larger than the satellite systems; the SECOR satellite systems have significant Z-rotations with respect to the average terrestrial system; and the ETH geodetic system may have a significant rotation with respect to the average terrestrial system.     

A symbolic analysis of Vermeille and Borkowski polynomials for transforming 3D Cartesian to geodetic coordinates

Closed form solutions for transforming 3D Cartesian to geodetic coordinates reduce the problem to finding the real solutions of the fourth degree latitude equation or variations of it. By using symbolic tools (Sturm–Habicht coefficients and subresultants mainly) we study the methods (and polynomials) proposed by Vermeille and Borkowski to solve this problem. For Vermeille approach, the region where it cannot be applied is completely characterized. For those points it is shown how to transform 3D Cartesian to geodetic coordinates and a new method for solving Vermeille equation for those cases not yet covered is introduced. Concerning Borkowski’s approach, the symbolic analysis produces a complete characterization of the singular cases (i.e. where multiple roots appear).     

Accurate algorithms to transform geocentric to geodetic coordinates

The problem of the transformation is reduced to solving of the equation $$2 sin (\psi - \Omega ) = c sin 2 \psi ,$$ where Ω = arctg[bz/(ar)], c = (a²?b²)/[(ar)²]^1/2 a andb are the semi-axes of the reference ellisoid, andz andr are the polar and equatorial, respectively, components of the position vector in the Cartesian system of coordinates. Then, the geodetic latitude is found as ?=arctg [(a/b tg ψ)], and the height above the ellipsoid as h = (r?a cos ψ)cos ψ + (z?b sin ψ)sin ψ. Two accurate closed solutions are proposed of which one is approximative in nature and the other is exact. They are shown to be superior to others, found in literature and in practice, in both or either accuracy and/or simplicity.     

GPS测量坐标转换实用性问题的分析

针对ＧＰＳ测量坐标转换方法中存在的问题,提出了强制符合平面四参数法和多项式拟合法,这两种方法能够有效的克服高程系统以及椭球参数不一致造成的误差,比较适合于工程自由坐标之间的转换;同时本文给出了基于“全球大地水准面的几何中心地球质心相重合”这一假设之上的莫洛金斯坐标转换法,该法不需要联测公共点即可将ＷＧＳ－８４坐标转换成本地局部坐标。上述几种方法减少了额外联测的工作量,提高了ＧＰＳ的使用效率。     

Transformation from Cartesian to Geodetic Coordinates Accelerated by Halley’s Method

By using Halley’s third-order formula to find the root of a non-linear equation, we develop a new iterative procedure to solve an irrational form of the “latitude equation”, the equation to determine the geodetic latitude for given Cartesian coordinates. With a limit to one iteration, starting from zero height, and minimizing the number of divisions by means of the rational form representation of Halley’s formula, we obtain a new non-iterative method to transform Cartesian coordinates to geodetic ones. The new method is sufficiently precise in the sense that the maximum error of the latitude and the relative height is less than 6 micro-arcseconds for the range of height, −10 km ≤ h ≤ 30,000 km. The new method is around 50% faster than our previous method, roughly twice as fast as the well-known Bowring’s method, and much faster than the recently developed methods of Borkowski, Laskowski, Lin and Wang, Jones, Pollard, and Vermeille.     

The Bruns formula in three dimensions

The Bruns formula is generalized to three dimensions with the derivation of equations expressing the height anomaly vector or the geoid undulation vector as a function of the disturbing gravity potential and its spatial derivatives. It is shown that the usual scalar Bruns formula provides not the separation along the normal to the reference ellipsoid but the component of the relevant spatial separation along the local direction of normal gravity. The above results which hold for any type of normal potential are specialized for the usual Somigliana-Pizzetti normal field so that the components of the geoid undulation vector are expressed as functions of the parameters of the reference ellipsoid, the disturbing potential and its spatial derivatives with respect to three types of curvilinear coordinates, ellipsoidal, geodetic and spherical. Finally the components of the geoid undulation vector are related to the deflections of the vertical in a spherical approximation.     

An alternative algebraic algorithm to transform Cartesian to geodetic coordinates

The algorithm to transform from 3D Cartesian to geodetic coordinates is obtained by solving the equation of the Lagrange parameter. Numerical experiments show that geodetic height can be recovered to 0.5 mm precision over the range from −6×10⁶ to 10¹⁰ m. Electronic Supplementary Material: Supplementary material is available in the online version of this article at     

The Meissl scheme for the geodetic ellipsoid

We present a variant of the Meissl scheme to relate surface spherical harmonic coefficients of the disturbing potential of the Earth’s gravity field on the surface of the geodetic ellipsoid to surface spherical harmonic coefficients of its first- and second-order normal derivatives on the same or any other ellipsoid. It extends the original (spherical) Meissl scheme, which only holds for harmonic coefficients computed from geodetic data on a sphere. In our scheme, a vector of solid spherical harmonic coefficients of one quantity is transformed into spherical harmonic coefficients of another quantity by pre-multiplication with a transformation matrix. This matrix is diagonal for transformations between spheres, but block-diagonal for transformations involving the ellipsoid. The computation of the transformation matrix involves an inversion if the original coefficients are defined on the ellipsoid. This inversion can be performed accurately and efficiently (i.e., without regularisation) for transformation among different gravity field quantities on the same ellipsoid, due to diagonal dominance of the matrices. However, transformations from the ellipsoid to another surface can only be performed accurately and efficiently for coefficients up to degree and order 520 due to numerical instabilities in the inversion.     

空间直角坐标计算大地坐标的抛物线逼近法

采用抛物线逼近法求解大地纬度和大地高,先计算空间点在椭球面上的子午面坐标,然后求解点的大地纬度和大地高。     

Harmonic analysis of the Earth's gravitational field by means of semi-continuous ephemerides of a low Earth orbiting GPS-tracked satellite. Case study: CHAMP

An algorithm for the determination of the spherical harmonic coefficients of the terrestrial gravitational field representation from the analysis of a kinematic orbit solution of a low earth orbiting GPS-tracked satellite is presented and examined. A gain in accuracy is expected since the kinematic orbit of a LEO satellite can nowadays be determined with very high precision, in the range of a few centimeters. In particular, advantage is taken of Newton's Law of Motion, which balances the acceleration vector with respect to an inertial frame of reference (IRF) and the gradient of the gravitational potential. By means of triple differences, and in particular higher-order differences (seven-point scheme, nine-point scheme), based upon Newton's interpolation formula, the local acceleration vector is estimated from relative GPS position time series. The gradient of the gravitational potential is conventionally given in a body-fixed frame of reference (BRF) where it is nearly time independent or stationary. Accordingly, the gradient of the gravitational potential has to be transformed from spherical BRF to Cartesian IRF. Such a transformation is possible by differentiating the gravitational potential, given as a spherical harmonics series expansion, with respect to Cartesian coordinates by means of the chain rule, and expressing zero- and first-order Ferrer's associated Legendre functions in terms of Cartesian coordinates. Subsequently, the BRF Cartesian coordinates are transformed into IRF Cartesian coordinates by means of the polar motion matrix, the precession–nutation matrices and the Greenwich sidereal time angle (GAST). In such a way a spherical harmonic representation of the terrestrial gravitational field intensity with respect to an IRF is achieved. Numerical tests of a resulting Gauss–Markov model document not only the quality and the high resolution of such a space gravity spectroscopy, but also the problems resulting from noise amplification in the acceleration determination process.     

Accurate area determination in the cadaster: case study of Slovenia

This paper discusses methodological problems of accurate area determination in the cadaster. The paper contrasts the ambiguous legal definition of the parcel boundary and parcel area in relation to the theoretically well-defined geodetic parcel boundary and the geodetic parcel area on the reference ellipsoid. To align with the real world, parcel area must account for terrain elevation. Various approximate methods for area determination which can be used in the cadaster are tested. A highly accurate method for parcel area computation is proposed, based on an equal-area projection. Considering the geodetic parcel area as a reference, the achievable accuracy of different methods is evaluated. For this analysis, the coordinates of the parcel boundary points are treated as error-free. Finally, the relevance of various systematic errors is discussed in relation to the statistical uncertainty of the parcel area, which could be gained by an a real-time kinematic GNSS survey. A case study is presented for the territory of Slovenia, its georeferencing rules, land demarcation pattern, and characteristics of its topography. Based on the results of this study, some general recommendations for the parcel area determination are given.     

斜轴变形椭球高斯投影方法

针对东西跨度较大的线路,借助最小二乘法建立斜轴参考椭球,以减小高斯投影横坐标;通过坐标系转换理论,推导出测区在各坐标系下的空间直角坐标,进而确定测区相对于斜轴参考椭球上的大地坐标;利用椭球变换法建立斜轴变形椭球以减小因高程引起的投影变形。以某铁路线为例,可知"斜轴变形椭球高斯投影方法"可大大减小投影后横轴方向分量,避免高斯投影分带现象,同时有效减小高程及其引起的投影变形。该方法数学模型严谨、运算过程清晰,便于编制相关软件,可投入工程使用。     

W0: an estimate in the Finnish Height Datum N60, epoch 1993.4, from twenty-five GPS points of the Baltic Sea Level Project

Based upon a data set of 25 points of the Baltic Sea Level Project, second campaign 1993.4, which are close to mareographic stations, described by (1) GPS derived Cartesian coordinates in the World Geodetic Reference System 1984 and (2) orthometric heights in the Finnish Height Datum N60, epoch 1993.4, we have computed the primary geodetic parameter W ₀(1993.4) for the epoch 1993.4 according to the following model. The Cartesian coordinates of the GPS stations have been converted into spheroidal coordinates. The gravity potential as the additive decomposition of the gravitational potential and the centrifugal potential has been computed for any GPS station in spheroidal coordinates, namely for a global spheroidal model of the gravitational potential field. For a global set of spheroidal harmonic coefficients a transformation of spherical harmonic coefficients into spheroidal harmonic coefficients has been implemented and applied to the global spherical model OSU 91A up to degree/order 360/360. The gravity potential with respect to a global spheroidal model of degree/order 360/360 has been finally transformed by means of the orthometric heights of the GPS stations with respect to the Finnish Height Datum N60, epoch 1993.4, in terms of the spheroidal “free-air” potential reduction in order to produce the spheroidal W ₀(1993.4) value. As a mean of those 25 W ₀(1993.4) data as well as a root mean square error estimation we computed W ₀(1993.4)=(6 263 685.58 ± 0.36) kgal × m. Finally a comparison of different W ₀ data with respect to a spherical harmonic global model and spheroidal harmonic global model of Somigliana-Pizetti type (level ellipsoid as a reference, degree/order 2/0) according to The Geodesist's Handbook 1992 has been made. Received: 7 November 1996 / Accepted: 27 March 1997     

Interpolation of Point Values From Isoline Maps

A Bidirectional Hermitian Spline (BHS) method for the estimation of point values from isoline maps is presented and compared to three other methods. Hermitian splines are used and first derivatives are estimated by either Akima's method or by a clamped cubic spline, if Akima's method returns a zero first derivative. Every desired point value is interpolated twice, once by each of two orthogonally-directed splines. The two spline estimates are then averaged using the error formula for Hermitian splines. In addition, a periodic Hermitian spline is constructed around the study-area perimeter (representing a cross-sectional profile of the edge) to damp undesirable edge effects. Point values can be estimated from small-scale isoline maps drawn in spherical coordinates or from large-scale isoline maps drawn in Cartesian coordinates.