Alternative transformation from Cartesian to geodetic coordinates by least squares for GPS georeferencing applications T. Soler a,n , J.Y. Han b , N.D. Weston a a Spatial Reference System Division, National Geodetic Survey, NOAA, Silver Spring, MD 20910, USA b Department of Civil Engineering, National Taiwan University, No. 1, Sec. 4, Roosevelt Road, Taipei 10617, Taiwan article info Article history: Received 25 August 2011 Received in revised form 27 October 2011 Accepted 28 October 2011 Available online 6 November 2011 Keywords: Cartesian coordinates Geodetic coordinates Inverse transformation Least squares abstract The inverse transformation of coordinates, from Cartesian to curvilinear geodetic, or symbolically (x,y,z)-(l,j,h) has been extensively researched in the geodetic literature. However, published formula- tions require that the application must be deterministically implemented point-by-point individually. Recently, and thanks to GPS technology, scientists have made available thousands of determinations of the coordinates (x,y,z) at a single point perhaps characterized by different observational circumstances such as date, length of occupation time, distance and geometric distribution of reference stations, etc. In this paper a least squares (LS) solution is introduced to determine a unique set of geodetic coordinates, with accompanying accuracy predictions all based on the given sets of individual (x,y,z) GPS-obtained values and their variance–covariance matrices. The (x,y,z) coordinates are used as pseudo-observations with their attached stochastic information in the LS process to simultaneously compute a unique set of (l,j,h) curvilinear geodetic coordinates from different observing scenarios. Published by Elsevier Ltd. 1. Introduction Many scientists have investigated the so-called, non-trivial, inverse transformation of coordinates from Cartesian (x,y,z) to curvilinear (orthogonal) geodetic coordinates (l,j,h). Both sets of coordinates are defined with respect to any arbitrary geodetic Cartesian reference frame and, in the case of geodetic coordinates, a complementary rotational ellipsoid with center at the origin of the Cartesian frame, its semi-minor axis coincident with the z-axis, and its semi-major axis on the equatorial plane defined by the x and y axes. The selected rotational ellipsoid is typically the GRS80 ellipsoid as adopted by the International Association of Geodesy (Moritz, 1992). For an exhaustive study of the many inverse trans- formations available to the user, the readers may consult Featherstone and Claessens (2008) and Awange et al. (2010, p. 157) where they will find a partial list of approaches to solve this specific transformation problem. A recent article by Shu and Li (2010) cites newly developed algorithms to compute geodetic coordinates not mentioned in any of the above mentioned references. For complete- ness, it should also be mentioned that the International Earth Rotation and Reference System Service (IERS) recommends the use of Fukushima’s (1999) iterative method. However, all these transfor- mation equations and algorithms were analytically developed to compute coordinates in a one-by-one point basis, that is, given the Cartesian coordinates of a point determine the equivalent curvilinear geodetic coordinates at the same point; therefore they are determi- nistic methods but not stochastic methods. The alternative method presented herein takes advantage of the large number of (x,y,z) determinations that we sometimes have on hand these days at a single particular point when processing GPS data. The main idea that we are proposing is to obtain the most accurate curvilinear geodetic coordinates obtain- able from the complete set of available GPS ‘‘detrended’’ time series (x,y,z) coordinates (e.g. Teza et al., 2010). Not only are the results going to be statistically meaningful, in a least squares (LS) sense, because as a byproduct the full variance–covariance (v–c) matrix for this uniquely derived triplet of coordinates can be determined. This statistical element is missing from the standard transformation formulas mentioned above. Other options, as for example, taking the weighted mean of all individually computed (x,y,z) values and finally transform them to (l,j,h) using any of the current methods will not be as complete, statistically speak- ing, as the procedure that will be outlined in this paper. If nothing else, because the v–c matrix of each individually processed GPS point (x,y,z) is known and taken into consideration in the LS solution. This information is not properly exploited when taking any other type of statistical sample averages. As an immediate practical application of this procedure, one may think of the calculation of a single unique set of geodetic coordinates for a point, referred to a predefined datum ellipsoid, determined from a set of original GPS-processed (x,y,z) solutions. The intention here is to get the ‘‘best’’ (l,j,h) coordinates for each Contents lists available at SciVerse ScienceDirect journal homepage: www.elsevier.com/locate/cageo Computers & Geosciences 0098-3004/$ - see front matter Published by Elsevier Ltd. doi:10.1016/j.cageo.2011.10.026 n Corresponding author. Tel.: þ1 301 7133205x157; fax: þ1 301 7134324. E-mail address: tom.soler@noaa.gov (T. Soler). Computers & Geosciences 42 (2012) 100–109