Developing low-cost automated tool for integrating maps with GNSS satellite positioning data

1 Introduction

The main issue in geodesy is (Hirvonen 1960): “Find the space coordinates of any point P at the physical surface S of the earth when a sufficient number of geodetic operations have been carried out along S.” Hence, the selection of a suitable frame to define the point P on spatial coordinates is critical. Nowadays, the world is witnessing a new era in data collection methods and visualizing the Earth’s physical features on a map due to the tremendous development in geospatial technologies. In practice, the position of each point on a map is determined by a particular reference system, which must be identified accurately to attain a high positional quality of spatial data (Snyder 1997, Kneissl et al. 2011). In contrast, most existing maps were created by classical surveys adopted on various local or regional datums that are non-geocentric and chosen to provide the best fit for the Earth’s shape in a specific geographic area. The advent of satellite positioning techniques such as the global navigation satellite system (GNSS) revolutionized position fixing methods in geodetic sciences. The GNSS observations are based on the world geodetic system 1984 (WGS84), and with the growing exchange of geographic data, both locally and globally, position information has to be given in both datums (Yang 1999, Kumar 1988). As a result, there is a necessity to establish a relationship between two reference frames to ensure the consistency of the coordinates (Vaníček and Steeves 1996).

Three-dimensional coordinate transformation is broadly applied for data fusion of various coordinate systems to build an interpolation surface with minimal distortions and estimate the coordinates of non-common locations (Chen and Hill 2005, Even-Tzur 2018). In general, several datum shift methods have been investigated in the literature to change between WGS84 and local datums that require a case-by-case evaluation due to deformation in terrestrial networks (Harvey 1986, Abd-Elmotaal 1994, Deakin and Leahy 1994, Newsome and Harvey 2003, Chuan and Yi 2014, Marx 2017, Yang et al. 2020, Yang et al. 2022). According to their simplicity and accuracy, 3D conformal models such as Helmert, Molodensky–Badekas, Bursa–Wolf, and polynomials (Bursa 1962, Molodensky et al. 1962, Wolf 1963, Krakiwsky and Thomson 1974, Kutoglu et al. 2002, Soycan 2005) are widely utilized in surveying, photogrammetry, geographic information systems, geodesy, LIDAR, and terrestrial laser scanning (Kashani 2006, Paffenholz and Bae 2012, Zeng et al. 2018).

The transformation of parameters from GNSS-derived coordinates to national terrestrial datum or vice versa include axes rotations, origin transitions, and scaling factors (Deakin 1998). They are defined using reference points in two sets of coordinates such that their spatial distribution and quantity impact the accuracy of results (Závoti and Kalmár 2016, Ioannidou and Pantazis 2020). The redundant stations, which are more than the minimum necessary number of common ones, are required to provide the best solution because the least-squares method can be employed to calculate the polynomial coefficients (Johnson and Faunt 1992, Challis 1995). Besides, the performance of the transformation algorithm is primarily evaluated by computing its parameters accuracy and the root mean square error of control points. This contribution aims to develop a low-cost spatial tool in the Microsoft Visual Studio environment to convert positional information to an alternative datum using 7-parameters similarity models and second-order conformal polynomial fit.

2 Materials and methods

Sustainable development goals (SDGs) as a roadmap to construct a better society for humankind must balance socio-economic and environmental components (GA 2015). Human capital formation, government policies, and technological progress are essential for carrying out these aims (Benhabib and Spiegel 2005, Banerjee and Roy 2014). A Geospatial database is a vital pillar in digital government transformation to underpin SDG attainment. In addition, implementing geodetic infrastructure based on an accurate and stable coordinate reference frame is significant and will provide a tie between scientists, policymakers, and geospatial society (Barbero et al. 2019). Previously, the geographic coordinates, including latitude φ, longitude λ, and orthometric or geoidal height H, have been established worldwide by classical surveying techniques, but at present, GNSS technology is applied to measure latitude and longitude and ellipsoidal height h on WGS84, as indicated in Figure 1. Hence, the most central issue in building geodatabases is to unify multi-sources spatial referenced data.

Figure 1

Geodetic and Cartesian coordinate system.

2.1 Case study description

The test site is situated in the Golden Triangle of Ghana, between latitudes 4°30′ N and 7°30′ N and longitude 3° W and 1° E (Figure 2), with a total area of 238,540 km² and a maximum elevation of 880 m (Higazi 2005). The national ellipsoid applied in Ghana for all geospatial purposes is the Accra datum based on the War Office 1926 spheroid, with European Petroleum Survey Group (EPSG) code of 2136, which has the following linear parameters: semi-major axis a = 6378300.00 m, semi-minor axis b = 6356751.689 m, and inverse flattening 1/f = 296.0 (Ayer and Fosu 2008). The conformal projection adopted in Ghana for mapping the geographic position of features is a transverse Mercator projection that covers an area by a single grid zone extending longitudinally 6°, from north to south, with a central meridian of λ _o = 1° W, the origin of latitude φ _o = is 4° 40′ N with scale factor K = 0.99975. False easting E _o = 274320.00 m and false northing N _o = 0.00 km are introduced in the mapping system to keep the positive value of coordinates (Thomas et al. 2000). The performance of the developed tool is verified using a dataset acquired from a previous study to build a numerical model across several control points (Ziggah et al. 2019) to change coordinates from the geocentric datum into the Ghana local reference frame.

Figure 2

Spatial distribution of GNSS control points in the study area.

2.2 Research methodology

The use of information processing technology to conduct critical functions in the mapping industry has been growing in recent years. The increasing diversity of computerized geographical data sources creates new challenges from the technical point of view that must be overcome to enable efficient integration. Indeed, if building a complete software package is required, available free spatial tools cannot be directly integrated since they are either built on old transformation techniques or lack the integrity of commercial software with a suitable user interface. Accordingly, this investigation endeavors to build a low-cost tool for harmonizing heterogeneous geospatial data sources within a coordinate reference system by providing the knowledge and framework required to develop a tool that can run on its own for coordinate transformation, while preserving that adaptability to be implemented in a complete software package. First, various 3D transformation models will be reviewed to highlight their characteristics and mathematical equations to be implemented in the developed tool. Also, the performance of the software tool will be evaluated to verify that it meets all specified requirements and positional quality. It aims to prove that all performance-degrading faults have been eradicated before the system is used. Finally, conclusions can be drawn based on the study results. Accordingly, Figure 3 presents a methodology flowchart used in this research.

Figure 3

Graphical depiction of the research methodology.

2.3 Similarity transformation model

A similarity transformation is a coordinate change that contains just a scale, rotation, and translation. The Bursa–Wolf model is a simplified form of 3D Helmert’s approach utilized to determine the datum transformation parameters between any two 3D reference systems. It defines a geometrical relationship between the source (X _S , Y _S , Z _S ) and target (X _T , Y _T , Z _T ) datums using common pairs of points whose coordinates are delivered in both the systems. A seven-parameter function includes three transition components (δ _X , δ _Y , δ _Z ) along the X, Y, and Z-axes, respectively, three axes rotations (R _X , R _Y , R _Z ), and a scale factor (K), as shown in Figure 4. The rotation angles are positive counterclockwise when considered from the axis’s positive end towards the center of the coordinate. The conversion expression of Bursa–Wolf is given in matrix form by equation (1) (Blewitt et al. 1992, Güneş and Demir 2021, Kalu et al. 2022).

(1) X T Y T Z T = X S Y S Z S + ( 1 + K ) 1 R Z − R Y − R Z 1 R X R Y − R X 1 X S Y S Z S + δ X δ Y δ Z .

Figure 4

Bursa–Wolf algorithm of the datum transformation.

The Bursa–Wolf transformation parameters are identified from a redundant set of n points by applying parametric least-squares adjustment to fit data to a mathematical surface. This approach is formulated as a Gauss–Markov model to produce a minimum difference between the observed and computed coordinates of control points (Chang 2015). Equation (2) shows the matrix form of parametric least-squares adjustment of the Bursa–Wolf method with its solution given in equation (3) (Mikhail and Gracie (1981), Ogundare (2018)).

(2) V 3 n × 1 + B 3 n × m ⋅ Δ m × 1 = f 3 n × 1 ,

(3) N m × m = B t m × 3 n . B 3 n × m , t m × 1 = B t m × 3 n . f 3 n × 1 , Δ = N − 1 t , V 3 n × 1 = f 3 n × 1 − B 3 n × m ⋅ Δ m × 1 , Q ΔΔ = N − 1 ,

(4) v position = v X S 2 + v Y S 2 + v Z S 2 . σ o = V t . V 3 n − m ,

where B is the numerical coefficient matrix of parameters, V is the vector of observational residuals, Δ is the vector of unknown parameters m, N is the square symmetric matrix, f is the vector of numerical constants, Q _ΔΔ is the cofactor matrix which describes the accuracy of the estimated parameters, v is the residual vector of each common point position, and σ _o is the posterior reference standard deviation.

(5) B 3 n × m = 0 − Z S 1 Y S 1 1 0 0 X S 1 Z S 1 0 − X S 1 0 1 0 Y S 1 − Y S 1 X S 1 0 0 0 1 Z S 1 … … … … … … … … … … … … … … … … … … 0 − Z S n Y S n 1 0 0 X S n Z S n 0 − X S n 0 1 0 Y S n − Y S n X S n 0 0 0 1 Z S n ,

(6) V t 3 n × 1 = v X S 1 v Y S 1 v Z S 1 ⋯ ⋯ ⋯ v X S n v Y S n v Z S n , f t 3 n × 1 = X T 1 − X S 1 Y T 1 − Y S 1 Z T 1 − Z S 1 ⋯ ⋯ ⋯ X T n − X S n Y T n − Y S n Z T n − Z S n , Δ t m × 1 = R X R Y R Z δ X δ Y δ Z K .

The Molodensky–Badekas model, as a similarity transformation, includes a centroid to release the strong correlation between the adjusted parameters when utilized across a limited Earth surface area (El-Mowafy et al. 2009). The average of coordinates of control points is employed to obtain the centroid coordinates (X _c , Y _c , Z _c ) for the Molodensky–Badekas algorithm, as shown in equation (7) (Kutoglu et al. 2002). Hence, the source system coordinates are shifted to the centroid according to equation (8). The general mathematical relationship of the model is described in equation (9). Figure 5 demonstrates the centroidal system whose origin is at a centroid, and its axes are parallel to the source one. The matrixial form of equation (9) can be written as mentioned previously in Bursa–Wolf for applying the least-square solution with a difference in the B matrix of numerical coefficients indicated in equation (10).

(7) X c Y c Z c t = ∑ i = 1 n X S i n ∑ i = 1 n Y S i n ∑ i = 1 n Z S i n t ,

(8) X ¯ Y ¯ Z ¯ t = X S − X c Y S − Y c Z S − Z c t ,

(9) X T Y T Z T = X S Y S Z S + ( 1 + K ) 1 R Z − R Y − R Z 1 R X R Y − R X 1 X ¯ Y ¯ Z ¯ + δ X δ Y δ Z ,

(10) B 3 n × m = 0 − Z ¯ 1 Y ¯ 1 1 0 0 X ¯ 1 Z ¯ 1 0 − X ¯ 1 0 1 0 Y ¯ 1 − Y ¯ 1 X ¯ 1 0 0 0 1 Z ¯ 1 … … … … … … … … … … … … … … … … … … 0 − Z ¯ n Y ¯ n 1 0 0 X ¯ n Z ¯ n 0 − X ¯ n 0 1 0 Y ¯ n − Y ¯ n X ¯ n 0 0 0 1 Z ¯ n .

Figure 5

Graphical representation of Molodensky–Badekas method.

3 Proposed software

Geospatial data represent real-world features in a suitable coordinate system. Hence, the coordinates must be converted to harmonize heterogeneous georeferenced data from different sources. This study built a low-cost spatial tool called “TopoGC” in the Microsoft Visual Studio environment to compare source and destination geodetic coordinates based on the process indicated in Figure 6. The three-parameter, five-parameter, Bursa–Wolf, Molodensky–Badekas, and second-degree polynomial transformation algorithms were addressed to be automated. By defining projection parameters, the converted coordinates can be mapped to a flat surface using an appropriate grid reference system. Figure 7 depicts the life cycle steps of the implemented software. The development process has six phases which start with requirement identification according to software system application and finish with evaluation and maintenance against the design principles. The logic plan’s detail is expressed in the design stage by conducting the flowchart and pseudocode algorithm and then coded using Visual C++ programming language. However, the TopoGC tool enables users to modify point coordinates inserted manually or from MS-Excel. The outputs can be represented, printed, and exported to MS-Excel. Indeed, the ActiveX control is constructed to share information between this application and MS-Excel using Visual Basic language.

Figure 6

Data transformation and mapping process in the developed tool.

Figure 7

Schematic diagram of developed tool life cycle.

3.1 Estimating transformation parameters

The TopoGC tool is a Windows-based application with a user-friendly interface for datum transformation and mapping operations. TopoGC toolset automatically changes three-dimensional geodetic coordinates from one geographic datum to another using the transformation definition dialog box, as shown in Figure 8. The mathematical model parameters are assigned after inserting the coordinates of control points in both reference frameworks. The spatial data coordinates can also be imported from an MS-Excel file reordered within the list view or saved in a text file. The report of calculation results appears in a new window, which includes the following elements:

• Type of transformation algorithm utilized in computation.

• Reference points coordinates.

• Transformation parameters with their accuracy.

• Residual of observations.

• Standard error of unit weight.

Figure 8

Dialog for calculating transformation parameters.

3.1.1 Conversion into a new reference system

The coordinates transformation to another surface is performed from the Final Coordinate box according to the calculated parameters. The coordinates should be input, and the conversion into a new system will be achieved, and the output data are displayed in the list view that may be printed, exported to MS-Excel, or projected using an appropriate mapping system. On the other hand, the moving from

3.1.2 Geodetic vs geocentric coordinates

As is well-known, the coordinates in a geodetic system are characterized as either geodetic or geocentric. The conversion between these two systems is one of the fundamental tasks in computational geodesy. Equation (18) directly converts the spheroidal coordinates of point B in Figure 1 into Cartesian coordinates, while the inverse solution is provided in equation (19) (Borkowski 1989, Ligas 2013). The developed tool has automated the relationship between Cartesian and geodetic coordinates on a biaxial ellipsoid.

(18) X = ( ν + h ) cos φ cos λ , Y = ( ν + h ) cos φ sin λ , Z = ( ν ( 1 − e 2 ) + h ) sin φ ,

(19) λ = cos − 1 X P f , φ = tan − 1 Z + e / 2 ⋅ b ⋅ sin 3 ω P f − e 2 ⋅ a ⋅ cos 3 ω , h = P f cos φ − ν ,

where e ² is the first eccentricity squared and ν is the radius of curvature in the prime vertical.

ω = tan − 1 Z ⋅ a P f ⋅ b , P f = ( X 2 + Y 2 ) , e 2 = 1 − b a 2 .

4 Software verification

Changing the coordinates between different reference systems is vital for fusing various spatial data sources. In this regard, the performance of the TopoGC application will be verified using data illustrated in Figure 2. However, the results of mathematical equations given in the Sections 2 and 3 were generated utilizing a developed tool and compared with others from a commercial software package (Geographic Calculator) and manually in an MS-Excel spreadsheet. Generally, the comparison results showed identical values to the manual approach and minor variations (in millimeters) with the commercial one. In practice, the research underlines the importance of combining the quantitative and qualitative approaches in implementing an accurate analysis of datum transformation algorithm choice. Figure 9 shows the flowchart adopted for 3D transformation coordinates. The parameters of the Bursa–Wolf and Molodensky–Badekas models with their accuracies were estimated and recorded in Table 1. The residuals of control points were also calculated, which are identical for both the systems. The obtained standard deviation values of unit weight show the accuracy of the conversion harmonization process for every dataset. Another validation of interpolated surface can be achieved by determining the coordinate differences at checkpoints between the observed and computed data from the derived parameters. Indeed, the Molodensky–Badekas transformation requires more information (the centroid coordinates) than the Bursa–Wolf method, which is why the Bursa–Wolf is more common. Table 2 indicates the obtained polynomial algorithm parameters and standard deviations, representing the model’s accuracy in fitting the control points. Generally, the finding’s reliability relies on the accuracy of the observations in both the systems and the proper distribution and coverage of data in the investigated area. Table 2 indicates the obtained polynomial algorithm parameters and standard deviations, representing the model’s accuracy in fitting the control points. Generally, the finding’s reliability relies on the accuracy of the observations in both the systems and the proper distribution and coverage of data in the investigated area.

Figure 9

Depiction for developing datum transformation process.

The findings of the above tables can be summarized in the following expression to facilitate a comparative analysis of the investigated models. These values reflect statistically significant results on a well-fit surface for the training dataset in the study area. Determining a suitable transformation method relates to the distortion at data points and spatial pattern distribution. However, comparing descriptive statistics between tested approaches indicates that the seven- and five-parameter models can predict coordinates with an accuracy better than others, while the standard errors in the polynomial model showed the highest trend. Additionally, another investigation was applied based on the visual inspection of residual values in Cartesian coordinates between known and extracted data from estimated surfaces.

Figure 10 confirms that the polynomial model resulted in the highest ΔX error among various addressed transformation approaches, whereas the three-parameter method had a comparable value of 67% of the points, while the five- and seven-parameter methods were 33%. The ΔY values in Figure 11 show that the five- and seven-parameter methods are superior to the others. At the same time, they had more significant distortions than the X-direction. On the other hand, the worst performance occurred in the polynomial. Similarly, the Z-coordinate differences in Figure 12 clarify that the five and seven parameter techniques also generated the best interpolator than the others. Besides, the capability of investigated methods was evaluated by computing the vector values of each testing dataset and representing them, as shown in Figure 13. The assessment results quality give insight into the ability of the five- and seven-parameter models compared to the considered ones, while polynomial ranks the last.

Figure 10

X-coordinate differences between actual and computed values at checkpoints.

Figure 11

Differences between observed and calculated dataset Y-coordinate values.

Figure 12

Z-coordinate differences between measured and reproduced values at checkpoints.

Figure 13

Interface to transform between geodetic and Cartesian coordinates.

Besides the numerical comparison of the three models provided in this study, other literature sources have investigated the suitability of these models, and their conclusions match those of this study (Featherstone 1997, Ziggah and Youjian 2013, Ziggah et al. 2019).

5 Conclusion

In the past, the national geodetic control networks have been constructed by accurate astronomical and geodetic observations on the non-geocentric datum. In recent years, the dynamic development of modern satellite positioning devices has brought new techniques to capture data based on a geocentric reference system. However, integrating multi-source geospatial data increases operational efficiency and improves the mapping industry. The diversity of reference systems was addressed by transforming positional information into an alternative datum using mathematical models in which the parameters set are defined from control points. This study considered this issue by discussing the theoretical aspect of the geodetic harmonization process of spatial datasets and developing a tool in the Microsoft Visual Studio environment using the least-squares approach to estimate the correct values of transformation parameters. The results indicated that the developed automated tool provides high precision in computing transformation models. Nevertheless, the concept of an ideal transformation algorithm faces a real challenge of being influenced by the distortion of terrestrial network positions with those acquired from GNSS receivers. The quantitative analysis of transformation algorithms showed that five- and seven-parameter models are the best fit surface for the common points, followed by the three-parameter and polynomial algorithms.