Detail Explanation of Coordinate Transformation Procedures Used in Geodesy and Geoinformatics on the Example of the Territory of the Republic of Croatia

: With strong technology improvement, especially GNSS, coordinate transformation has become an essential tool in everyday practice. For this reason, it is beneficial to have a good understanding of the mathematical process of coordinate transformation, as well as the method for calculating transformation parameters. The purpose of this article is to provide a detailed explanation of the coordinate transformation procedure and the calculation of transformation parameters for seven transformation models: 3-parameter, 5-parameter standard and abridged, 7-parameter, 8-parameter, 9-parameter, and 12-parameter. For each transformation model, a basic transformation equation, a procedure for calculating transformation parameters, inverse expression for reversible transformation, as well as characteristics of each model are provided. As an addition, for each model, a numerical example of calculating transformation parameters for the territory of the Republic of Croatia is provided. Additionally, for the example of the 7-parameter transformation, importance of reversible transformation, rotation convention, as well as the form of the rotation matrix used during transformation are demonstrated.


Introduction
Coordinate transformation represents a change of coordinates from one coordinate reference system to another coordinate reference system based on a different datum through a one-to-one relationship.There are two important terms in this definition-coordinate system and datum.The coordinate system is a set of mathematical rules for specifying how coordinates are to be assigned to points, whereas datum represents a parameter or set of parameters that define the position of the origin, the scale, and the orientation of a coordinate reference system.The two most commonly used coordinate systems in geodesy are the Cartesian coordinate system, which gives the position of points relative to n multiple axis, and the ellipsoidal coordinate system, in which position is specified by geodetic latitude, geodetic longitude, and ellipsoidal height.Coordinate reference systems that use at least two independent coordinate reference systems are named compound coordinate reference systems [1].
With respect to scale along the coordinate axis, coordinate transformations can be divided in two main groups: affine and similarity transformations.Similarity transformations, which preserve angles, are characterized by unified scale along all coordinate system's axis, whereas this is not the case for affine transformation parameters, which do not preserve angles.Coordinate transformation procedures are also distinguished by the number of transformation parameters.In geodesy and geoinformatics, often used coordinate transformation procedures, due to the number of parameters, are: 3-parameter [2,3], 5-parameter [2,4], 7-parameter [2,[5][6][7]], 8-parameter [8], 9-parameter [8], and 12-parameter [9].Transformation parameters for all mentioned transformation procedures can be obtained from identical points in two different coordinate reference frames.The increase in the number of transformation parameters allows more errors linked to the realization of the coordinate system to be eliminated, but with a consequence of very complicated mathematical expressions coordinate transformation procedures.
Many studies (e.g., [10][11][12][13][14]) considered different coordinate transformation procedures focusing on results but less on detail mathematical background, like derivation of numerical expressions.Therefore, this article focuses on a detail explanation of the aforementioned coordinate transformation procedures (mathematical derivations of transformation and inverse transformation equations along with numerical examples) for the example of the territory of the Republic of Croatia, utilizing coordinate transformation between the inherited coordinate reference system of the former Austrian-Hungarian Monarchy called Croatian State Coordinate System (cro.Hrvatski državni koordinatni sustav-HDKS) EPSG:3096 [15], with height component in Croatian Height Reference System 1875 (cro.Hrvatski visinski referentni sustav 1875-HVRS1875) EPSG:5195, and the official coordinate reference system-Croatian Terrestrial Reference System 1996 (cro.Hrvatski terestrički referentni sustav 1996-HTRS96) EPSG:4731 [16], with height component in Croatian Height Reference System 1971 (cro.Hrvatski visinski referentni sustav 1971-HVRS71) EPSG:5610, based on a set of 5194 identical points.Out of those identical 5194 points, 5000 were used for the calculation of the transformation parameters (Figure 1a), whereas the remaining 194 points were used for the external accuracy assessment.Points used for the external accuracy assessment were chosen in such a way that they are spread along the entire territory of the Republic of Croatia (Figure 1b).The height transformation between HVRS1785 and HVRS71 was excluded from this research in a way that all identical points in HDKS and HTRS96 have the same normal-orthometric height expressed in HVRS71.Transformation between HVRS1875 and HVRS71 is a part of the Croatian height transformation model [17], later included in the Croatian official HDKS/HTRS96 grid transformation model T7D [18].Therefore, all further mathematical expressions that utilize height component h denote normal-orthometric height H.

Inherited and Official Horizontal and Vertical Reference Systems of the Republic of Croatia
As a former member of Austrian-Hungarian Monarchy, the Republic of Croatia historically inherited the horizontal (2D) coordinate reference system (in 1901) called Croatian State Coordinate System (cro.Hrvatski državni koordinatni sustav-HDKS), Hermannskögel or HR1901.HDKS was realized by astro-geodetic measurements of the former Military Geographic Institute (MGI) from Vienna in a 1st-order triangulation network, adjusted in seven separate blocks (2nd, 5th, 6th, 22nd, 23rd, 24th, and 25th) [16].The origin (so-called fundamental point) of HDKS was situated at Hermannskögel hill, near Vienna with orientation to the nearby hill Hundsheimer Berg.HDKS utilizes a non geocentric geodetic datum that is defined by five parameters: two for surface definition that is a mathematical approximation of the earth shape-ellipsoid-and three positional parameters-coordinates of the fundamental point.For HDKS, the Bessel 1841 rotation ellipsoid was adopted (a = 6,377,397.15500m, 1/f = 299.15281886).At the fundamental point, surfaces of the selected ellipsoid and geoid coincide, which causes vertical deflection components as well as geoid undulation to be zero.Due to the coincidence of the ellipsoid and geoid at the fundamental point, vertical to geoid and normal to ellipsoid coincide as well, so ellipsoidal and orthometric heights are equal [19].
Following modern European trends, as well as the development of satellite positioning systems (primarily GPS and GLONASS), the poor documentation of inherited reference system inhomogeneity of HDKS (due to adjustment of separate block) was more than a good reason to redefine the national horizontal as well as height (and gravimetric) reference systems [21].
On 4 August 2004, the Government of the Republic of Croatia brought the Decree on establishing new official geodetic data and map projections of Republic of Croatia [16].According to the Decree, the Republic of Croatia adopted a new horizontal reference system called Croatian Terrestrial Reference System 1996 (cro.Hrvatski terestrički referentni sustav 1996-HTRS96), based on the European Terrestrial Reference System 1989 (ETRS89) with GRS80 level-ellipsoid (a = 6,378,137.00000m, 1/f = 298.25722210).The reference frame of HTRS96 is represented with the network of 78 1st-order trigonometric points (a part of the former Austrian-Hungarian 1st-order trigonometric network over the Croatian territory), with coordinates determined within combined solutions of GPS measuring campaigns: CROSLO'94 (ITRF96, e1994.4),CROREF'95 (ITRF96, e1995.7),and CROREF'96 (ITRF96, e1996.7) in mutual epochs of measurements 1995.55 with transformation to ETRS89 using EUREF Memo specifications [22].
Based on the same decree [16], Croatia adopted a new vertical reference system as well, called Croatian Height Reference System 1971 (cro.Hrvatski visinski referentni sustav 1971-HVRS71).HVRS71 was realized with 2nd High Accuracy Levelling (cro.II.nivelman visoke točnosti-IINVT), which was stretched over the territory of former Yugoslavia.Mean seal level was determined by measurements at five tide gauges on the eastern shore of the Adriatic Sea: Koper (Slovenia), Rovinj, Bakar, Split, and Dubrovnik for a period of 18.6 years.HVRS71 as well as HVRS1875 utilize normal-orthometric heights.Therefore, official CCRS of the Republic of Croatia includes horizontal coordinates in HTRS96 and a vertical component in HVRS71.

Review of Previous Studies
The first scientific studies of coordinate transformation procedures between different terrestrial reference frames appeared in the 1960s [2,5,6], so scientists like M. S. Molodensky (1909Molodensky ( -1991)), M. Bursa (1929Bursa ( -2019)), and H. Wolf (1910Wolf ( -1994) can be considered as pioneers of coordinate transformation.Development and usage of global navigation satellite systems (GNSS), primarly NAVSTAR GPS in the 1980s [23,24], for solving everyday geodetic tasks required linking national coordinate reference frames to global ones [3,7] and pointed out that, often, inherited (historical) coordinate reference systems are inappropriate for the application of new technologies.Therefore, finding the appropriate coordinate transformation procedure between those coordinate reference systems is essential [10][11][12][13]25].
For the territory of the Republic of Croatia, coordinate transformation procedure also became important with the usage of NAVSTAR GPS as a measuring technique in the late 1980s and early 1990s within first regional geodynamic (CRODYN) and EUREF GPS measuring campaigns [26][27][28].Later, two different transformation models were established to ensure the coordinate transformation procedure between HDKS and ETRS89/HTRS96, called DAT_ABMO and, currently, the official one, T7D, both of which are based on a 7-parameter transformation [14,18,21,[29][30][31][32][33][34].
The previous studies [15,18,19,27,29,35] have shown that transformation parameters between the inherited (HDKS) and new coordinate system (HTRS96) have been calculated multiple times since the establishment of the Republic of Croatia.The transformation parameters for the territory of the Republic of Croatia were first calculated in the year 2000 as part of the "Geodetic and Cartographic Datum of Republic of Croatia" project based on 120 identical points in HDKS and HTRS96 coordinate systems [19,29].In 2002, a new set of transformation parameters was calculated based on 241 identical points in two coordinate systems as part of the "Unique adjustment of 10 km GPS network of Republic of Croatia" project [35], resulting in improved 3D accuracy compared to the previous effort in 2000.The next set of transformation parameters was calculated in 2006 as part of the "Development of unique transformation model-HTRS96/HDKS" project, where a total of 1780 identical points in the inherited and new coordinate systems of the Republic of Croatia were used in the calculation, resulting in the transformation model T7D2006 [15].In order to improve the T7D2006 transformation model, a project, "New geoid model of Republic of Croatia and improvement of T7D transformation model", was initiated in 2009, which resulted in the new T7D2009 transformation model [18].Unique transformation parameters were calculated for this model based on 5200 points distributed across the entire territory of the Republic of Croatia.
The purpose of all the mentioned projects was to calculate transformation parameters for coordinate transformation in the territory of the Republic of Croatia, without the need for detailed explanations of the calculation procedure.For this reason, the aim of this article is to provide a detailed explanation of the procedure for calculating transformation parameters of all types of transformations, with transformation parameters given as solutions of numerical examples.

3-Parameter Transformation
This transformation (often called Molodensky 3-parameter transformation) belongs to the group of similarity transformations and is based only on translations along three coordinate axes (Figure 2).General expression for the 3-parameter transformation is [3] or, in matrix form, where X o is the position vector in the output coordinate system, t is the translation vector from the input to the output coordinate system, and X i is the position vector in the input coordinate system, all in meters.

Transformation Parameters
In the 3D case, three observation equations are needed for the transformation of a single point.Observation equations are linear as follows: where exponent i denotes ordinal number of the point (i = 1, 2, 3, . . ., n).
Design block matrix and observation vector, following Equation (3), are Least squares adjustment method using Equation (4) gives where P is a unit weight matrix.The solution of Equation ( 5) gives three translation parameters that correspond to the arithmetic mean of all coordinate differences of identical points in both systems: are coordinates in the output coordinate system, X i , Y i , Z i are coordinates in the input coordinate system, and k is the number of points.

Inverse Transformation
Using Equation (1), Cartesian rectangular coordinates in the output coordinate system can be calculated based on known Cartesian rectangular coordinates in the input coordinate system and transformation parameters.Reversion to the original coordinates in the input coordinate system is ensured by using the same transformation parameters obtained from Equation ( 6), so Equation (1) needs to be adapted by expressing X i on the left side of the equation:

Calculation of Transformation Parameters
Using Equations ( 4) and ( 5), transformation parameters for 3-parameter transformations are calculated, and the procedure of calculation is shown in Appendix A.1.Finally, transformation parameters for 3-parameter transformation from HDKS to HTRS96 calculated based on 5000 identical points are shown in Table 1.This transformation is also called the Molodensky transformation because Molodensky conducted differential equations that connect two ellipsoidal coordinate systems [2].Different approaches to those equations resulted in two types of this transformation [4]:

•
Standard 5-parameter Molodensky transformation, which represents the transformation between two ellipsoidal coordinate systems, including translations along three coordinate axes as well as changes of the semi-major axis and flattening between two ellipsoids (Figure 3) [4]; • Abridged 5-parameter Molodensky transformation, which represents simplified expressions of standard Molodensky trasnformation where ellipsoidal heights of points in the input coordinate system are unknown [4].Standard and abridged 5-parameter Molodensky transformations represent the approximation of a 3-parameter transformation, so they can be considered similarity transformations.An advantage of these transformations over 3-parameter transformations is the avoidance of conversion between ellipsoidal and Cartesian coordinates, and errors caused by approximation are negligible compared to the uncertainty of a 3-parameter transformation [36].

Transformation Parameters
Varga et al. [14] define the parameters of a 5-parameter transformation by taking three translation parameters from the 3-parameter transformation while obtaining the difference in ellipsoid parameters by subtracting the theoretical parameters given in definitions of the two ellipsoids.However, this method of determining transformation parameters assumes that the theoretical ellipsoid parameters and those realized in reality are completely identical, which is not the most accurate assumption, especially in the case of the Bessel 1841 ellipsoid, which is defined based on measurements in the local area.In this article, transformation parameters are calculated using the procedure proposed by Deakin [4], which is based on transformation equations.This approach eliminates the potential difference between theoretical and realized parameters because the parameters are computed directly form the mathematical expressions of transformation.
The first step in this procedure is to express the conversion from ellipsoidal to Cartesian 3D coordinates [37][38][39]: where N is the radius of the prime vertical plane of the ellipsoid, e 2 is the square of the eccentricity of the ellipsoid, and φ, λ, h are ellipsoidal coordinates.For N and e 2 , the following relations apply: where f is the flattening of the ellipsoid, which is function of semi-major (a) and semi-minor (b) axes of the ellipsoid.
Substituting Equations ( 9) and (10) into Equation ( 8) results in Cartesian X, Y, Z coordinates expressed as functions of a, f, φ, λ, and h: Using the total differential theorem [4,40], small changes of Cartesian X, Y, Z coordinates can be linked to small changes in a, f, φ, λ, and h, yielding mathematical expressions for the standard Molodensky transformation: where all δ represent differences of the values between two ellipsoids (e.g., δa = a o − a i , . . .).Equation ( 12) can be written in matrix form [4]: where matrices B and J are shown in Equations (A1) and (A2) in Appendix A.2. Changes in ellipsoidal coordinates can be extracted on the left side of Equation ( 13): The inverse of the matrix J cannot be determined in a simple way, so, it is necessary to use an adjoint matrix.The procedure of the calculation of the inverse matrix is described in the Appendix A.2.
After introducing Equations (A1) and (A4) to ( 14) and rearranging, changes in ellipsoidal coordinates are expressed in matrix form as: where M is the radius of curvature in the meridian: Matrices in the Equation ( 15) can be mutually multiplied, whereupon relations Finally, equations for the standard Molodensky 5-parameter transformation are [4] Abridged Molodensky transformation equations do not contain parameter h, so they can be obtained by setting the value h = 0 and simplifying coefficients next to δa and δ f in Equation (17).Using some approximations, final expressions for abridged Molodensky transformation are [4] δφ δh = cos φ cos λt x + cos φ sin λt y + sin φt z + (−1 + f sin 2 φ)δa + a i sin 2 φδ f Based on Equation (17), observation equations of the standard Molodensky transformation are formed as follows: According to Equation ( 19), design matrix and observation vector are According to the observation equations of the standard Molodensky transformation in Equation ( 19), based on Equation ( 18), the observation equations for the abridged Molodensky transformation follow: 21), design matrix and observation vector are In Equation ( 20), angular and linear values enter the adjustment, so it is necessary to homogenize the observations [41].Homogenization can be obtained by forming weight matrix P that defines relative relations between accuracy of geodetic latitude, longitude, and ellipsoidal height: where K is a favorably chosen coefficient, s φ is the position accuracy for geodetic latitude in angular units, s λ is the position accuracy for geodetic longitude in angular units, and s h is the accuracy of ellipsoidal height in linear units.
If the position accuracy is defined in linear units in the direction of geodetic latitude, longitude, and ellipsoidal height or is unknown and assumed to be the same amount in all directions, the accuracy assessment needs to be converted from linear to angular units using a so-called metric matrix [42].The metric matrix represents the relationship between differentials of curvilinear coordinates and their projections into a local geodetic coordinate system (ENU).Differentials can also be considered as free vectors [42], and by substituting differentials with vectors of standard deviations, the relationship between position accuracy of the point in linear and angular units is given by: where S N , S E , S U are position accuracies in the local geodetic coordinate system.Introducing Equation ( 24) in (23), gives the weight matrix: After definition of the design matrix, observation vector, and weight matrix, transformation parameters of the 5-parameter transformation can be calculated by the least squares adjustment method.

Inverse Transformation
Conducting the standard and abridged inverse transformation is not possible in a direct way because inverse transformation depends on the values φ i and λ i , which are unknown in the output coordinate system.In order to perform the inverse transformation, Ruffhead [43] gives the next solution:

•
In Equation ( 17) or (18), include negative values of the transformation parameters (t x , t y , t z , δa, δ f ) and coordinates of the point in the output coordinate system to obtain coordinates in the input coordinate system: • Obtained values of the transformed coordinates in the input coordinate system transform back to the output coordinate system: • Calculate the difference between transformed coordinates in the previous step and input ones and then subtract the difference from the coordinates obtained in the first step: Precision of this procedure is limited; to achieve a better precision, an iterative procedure is recommended.

Calculation of Transformation Parameters
Using Equations ( 20), (23), and (5), transformation parameters for 5-parameter transformations (standard and abridged) are calculated, and the calculation procedure is shown in Appendix A.2. Finally, transformation parameters for standard and abridged Molodensky 5-parameter transformations from HDKS to HTRS96 based on 5000 identical points are shown in Tables 2 and 3, respectively.From the calculated transformation parameters, it is evident that there are differences compared to the parameters provided by Varga et al. [14].Considering the magnitudes of these differences, it can be concluded that these discrepancies are not caused by the choice of identical points but have arisen due to different approaches in computing transformation parameters.However, it should be emphasized that the difference that has occurred in the ellipsoid parameters (da and df ) is not of significant magnitude because the translation parameters have also changed compared to those from the 3-parameter transformation.

The 7-Parameter Transformation
The 7-parameter transformation, or Helmert's transformation, is the most widespread and most used coordinate transformation in geodesy and geoinformatics and belongs to the group of similar transformations.It is based on three translations along and three rotations around three (X, Y, Z) axes as well as unique scale for all axes.In practice, two approximations of Helmert's transformation are used: Bursa-Wolf (Figure 4a) and Molodensky-Badekas (Figure 4b).Bursa-Wolf transformation rotation is performed around the origin of the coordinate system, and the general expression for transformation is [44] or, in matrix form, where X o is the position vector in the output coordinate system, t is the translation vector from the input to the output coordinate system, µ is scale, R is the rotation matrix from the input to the output coordinate system, and X i is the position vector in the input coordinate system.The translation vector as well as the position vectors are expressed in meter units, rotation matrix elements in radians, and scale can be expressed as where dµ is a scale change in ppm (parts per million, 10 −6 ) or ppb (parts per billion, 10 −9 ).Molodensky-Badekas transformation represents a generalization of the Bursa-Wolf transformation by defining the point around which the rotation is performed, called the centroid of the identical points and is given by expression where the X 0 , Y 0 , and Z 0 coordinates are defined by the expressions where X i , Y i , Z i are coordinates in the input coordinate system, and k is the number of points.
Translation and scale change in the 7-parameter transformation are unambiguously defined, whereas that is not the case with rotations.Knowing the rotation angles is not enough to carry out the rotation of the coordinate system.It is necessary to know the convention that defines the rotation matrix (direction and type of rotations) and the order of matrix multiplication.Rotation of the axes can be done according to coordinate frame rotation (CFR) convention, which implies a fixed position vector while the coordinate system rotates (Figure 5a) and according to position vector transformation (PVT) convention, which implies a fixed coordinate system while the position vector rotates (Figure 5b) [43].In Figure 5, angle α z represents the rotation angle around Z-axis, and angle θ represents the angle between the positive direction of the X-axis and the position vector in the input coordinate system.Based on these conventions, rotations can be observed in two ways: rotation of the output coordinate system around the input coordinate system, which corresponds to CFR convention, and rotation of the output coordinate system around the input coordinate system, which corresponds to PVT convention.In the exchange of transformation parameters, mandatory metadata is the rotation convention so the appropriate rotation matrix can be used in the process of coordinate transformation.In Equation (30), an abridged version of the rotation matrix is given.The complete rotation matrix can be obtained by multiplications of simple rotation matrices around three coordinate axes [8].
By rotating the output coordinate system around the input one (where the Z-axis is fixed), from Figure 5a, the following can be expressed: Terms for the output coordinate system can be written out by addition formulas for sine and cosine of angle differences: By including the first two terms in Equation ( 34) in (35), the relationship between two coordinate systems can be obtained (where Z i = Z o because the Z-axis is fixed): or, in matrix form, where R z represents the rotation matrix around the Z-axis.Other two rotation matrices (around the X and Y axes) can be obtained in the same way: The complete rotation matrix is a result of a multiplication of three rotation matrices from Equations ( 37) and (38).In the case when rotation angles are larger, the order of multiplication is very important, whereas for small angles, the order importance is not substantial.In practice, the system is assumed to rotate around the X, Y, and Z axes, respectively.In that case, the complete rotation matrix for the CFR convention is [8] R ZYX CFR =   cos α y cos α z sin α x sin α y cos α z + cos α x sin α z − cos α x sin α y cos α z + sinα x sin α z − cos α y sin α z − sin α x sin α y sin α z + cos α x cos α z cos α x sin α y sin α z + sin α x cos α z sin α y − sin α x cos α y cos α x cos α y   By applying the same procedure, from Figure 5b, the complete rotation matrix for the PVT convention can also be obtained by multiplying three rotation matrices: When transforming coordinates between two reference coordinate systems, where both have a globally defined geodetic datum, small rotation angles are expected.Therefore, the overall rotation matrix can be simplified by approximating its elements (sin α i = α i , cos α i = 1, sin α i sin α j = 0), so Equations ( 39) and ( 41) have a simpler form: To illustrate the importance of the rotation matrix form, a numerical transformation of a point from HDKS to HTRS96 based on a set of identical transformation parameters in three different cases for CFR convention-using a simple rotation matrix (Equation ( 42)), using the rotation matrix R XYZ , and using the rotation matrix R ZYX (Equation ( 39))-will be presented below.
Based on transformation parameters from HDKS to HTRS96: From the results above, it is evident that the choice of the rotation matrix during the transformation between global and local coordinate systems has an impact on the transformed coordinate, with a value of around 0.5-1 cm.Therefore, when acquiring transformation parameters, it is essential to clearly define which rotation matrix was used in the computation process so that the same matrix can be used in the transformation procedure.
When it comes to transformation between two global coordinate systems, the choice of the rotation matrix is not critical because the rotation angles are so small that they will not have a significant impact on the transformed coordinate.

Transformation Parameters
From the expressions for rotation matrix, it is clear that Equation ( 29) is not a linear function, and it is necessary to linearize it to form design matrix A and observation vector l [8].Due to complexity of linearization, it will be shown only for the case where the output coordinate system rotates around the input one (CFR convention); for the second case, similar procedure applies.
To linearize Equation ( 29), approximate values of rotation angles (α 0 x , α 0 y , α 0 z ) and corresponding corrections (δα x , δα y , δα z ) must be introduced.Definite values of rotation angles can be obtained by [8]: α z = α 0 z + δα z Based on approximate values of rotation angles, an approximate rotation matrix can be calculated.In abridged form, it can be expressed as [ Elements of the rotation matrix R are functions of rotation angles that can be linearized by approximated values of rotation angles [8]: where all derivatives from Equation ( 46) must be evaluated for approximate rotation angles.In matrix form, Equation ( 46) is as follows [8]: where matrices E, F, and G are shown in Appendix A.3.Scale factor can also be decomposed into two parts-approximate value of scale factor (µ 0 ) and corresponding correction (δµ)-so a definite value of scale factor can be calculated as [8] µ = µ 0 + δµ (48) Including terms from Equations ( 47) and (48) into Equation ( 29) yields [8] Because δα i and δµ have very small values, their mutual multiplication is negligible, so the observation equation in matrix form can be obtained [8]: Finally, based on the observation equation, design matrix A and observation vector l can be formed: It is important to know the order of the unknown parameters when forming design matrix A. In this case, the order of unknown parameters is t X , t Y , t Z , α X , α Y , α Z , and µ, so coefficients a ij in design matrix A have the following values: where index i represents the input coordinate frame, and exponent i represents the ordinal number of the point.

Inverse Transformation
To perform the inverse transformation for the 7-parameter transformation, Equation (29) needs to be modified by expressing the element X i on the left side of the equation: Scale factor µ is a scalar value, so µ −1 is very easy to calculate, but calculation of the inverse of the rotation matrix is more diffcult.To calculate the inverse of the rotation matrices in Equations ( 39) and ( 41), the rule of Cramer needs to be applied: where the determinant of the matrix can be calculated by the rule of Sarrus.
After calculation, the inverse rotation matrix for the CFR convention is cos α y cos α z − cos α y sin α z sin α y sin α x sin α y cos α z + cos α x sin α z cos α x cos α z − sin α x sin α y sin α z − sin α x cos α y sinα x sin α z − cos α x sin α y cos α z cos α x sin α y sin α z + sin α x cos α z cos α x cos α y   (54) and for the PVT convention is: cos α y cos α z cos α y sin α z − sin α y sin α x sin α y cos α z − cos α x sin α z sin α x sin α y sin α z + cos α x cos α z sin α x cos α y cos α x sin α y cos α z + sinα x sin α z cos α x sin α y sin α z − sin α x cos α z cos α x cos α y   (55) Based on the calculated inverses of the rotation matrices, it is evident that the rotation matrix is orthonormal, which means that R −1 = R T .
To illustrate the purpose of a reversible transformation, an example of transforming a point from HDKS to HTRS96 and back to HDKS based on two sets of parameters and the reversible transformation is shown.Based on 5000 identical points in the HDKS and HTRS96, it is possible to calculate two sets of parameters:

•
From HDKS to HTRS96: The procedure of reversible transformation (from HDKS to HTRS96 and back to HDKS) using two sets of parameters involves transformation from HDKS to HTRS96 based on the first set of parameters, and then transformation of the resulting point from HTRS96 back to HDKS based on the second set of parameters.For point P with the coordinates in HDKS:  By comparing these two points, it is evident that their coordinates are not entirely equal.In the second case, if only one set of transformation parameters is used (from HDKS to HTRS96) and the return to the HDKS coordinate system is performed by using the inverse expression (Equation ( 52)), the output point in HDKS has the following coordinates: General expression for the 8-parameter transformation is similar to the 7-parameter transformation, except the scale factor is not the same for all axes, so it can be written as [8] where S represents matrix of the scale change, and other elements are same as those in Equation ( 29).In Equation ( 57), the abridged form of the rotation matrix is given, whereas, for the calculation of transformation parameters, rotation matrix from the CFR convention (Equation ( 39)) are used.

Transformation Parameters
Content of the translation vector and rotation matrix were defined in the previous section.The difference compared to the 7-parameter transformation is the scale factor that, in this case, cannot be written as a scalar, but as a matrix because scale factor is not unique for all axes.By introducing approximate values of the scale factor and corresponding corrections, the scale matrix S can be written as [8] where M and P are Rotation matrix R was defined in Equation (47).Introducing Equations ( 47) and ( 58) into (56) yields The product of elements δα i and δµ i has negligible value, so Equation (59) can be simplified.From the simplified expression, observation equations can be formed and written [8]: Based on observation equations, design matrix A and observation vector l can be formed: It is also important to know what order of unknown parameters in the design matrix [8].In this case, order of the unknown parameters is t x , t y , t z , α x , α y , α z , µ xy , and µ z , so coefficients a ij in the design matrix have the following values [8]:

Inverse Transformation
To perform the inverse transformation, Equation (56) needs to be modified by expressing the element X i on the left side of the equation: where matrix R −1 is defined by Equation ( 54) or (55), and the inverse of the scale matrix can be easy calculated:

Calculation of Transformation Parameters
Using Equations ( 5) and (61), transformation parameters for the 8-parameter transformation are calculated, and the procedure is shown in Appendix A.4.Finally, parameters of the 8-parameter transformation from HDKS to HTRS96 based on 5000 identical points are shown in Table 5.

The 9-Parameter Transformation
Affine 9-parameter transformation consists of three translations and three rotations around three coordinate axes and three scale factors, which are different for all coordinate axes (Figure 7).Because of the non-unique scale factor, this transformation does not belong to the group of conformal (similarity) transformations.With the change of angles, there is also a change of position, direction, size, and shape of the transformed object.General expression for the 9-parameter transformation is [8] In Equation (65), the abridged form of the rotation matrix is written, and the complete form of the rotation matrix is given in Equation (39).

Transformation Parameters
For the linearization of the general expression, identical relations mentioned in 7parameter and 8-parameter transformation (Equations ( 44)-(47)) can be applied.The difference between 7-parameter and 8-parameter transformation is an unequal scale factor for all axes.By introducing approximate values of the scale factor (µ 0 x , µ 0 y , µ 0 z ) and corresponding corrections (δµ x , δµ y , δµ z ), the scale matrix S for the 9-parameter transformation can be written as [8] where Introducing Equations ( 47) and (66) in the general expression for the 9-parameter transformation (Equation (64)) yields As δα i and δµ i have very small values, their product will be negligible, so the observation equations can be written as [8] Based on the observation equations, design matrix A and observation vector l can be formed: It is also important to know what the order of unknown parameters in design matrix A is.In this case, the order of the unknown parameters is t x , t y , t z , α x , α y , α z , µ x , µ y , and µ z , so coefficients a ij in design matrix A have the following values:

Inverse Transformation
To perform the inverse transformation for the 9-parameter transformation, Equation (64) needs to be modified by expressing the element X i on the left side of the equation: where matrix R −1 is defined by Equation ( 54) or (55), and the inverse of the scale matrix can be easily calculated:

.3. Calculation of Transformation Parameters
Using Equations ( 5) and (69), transformation parameters for 9-parameter transformations are calculated, and the procedure is shown in Appendix A.5.Finally, parameters of the 9-parameter transformation from HDKS to HTRS96 based on 5000 identical points are shown in Table 6.With the change of the origin position, orientation of axes, and scale of coordinate system that exist in the 9-parameter transformation, the 12-parameter transformation also implies skew of coordinate axes [9].This implies that orthogonal axes of the input coordinate system are not transformed into orthogonal axes in the output coordinate system.The 12-parameter transformation is not often used in common geodetic use and GIS software due to its complexity and the insufficient accuracy improvement it offers compared to other types of transformations.However, the study of Wu et al. [46] shows that using 12-parameter transformation for the coordinate transformation of airborne LiDAR data and terrestrial LiDAR data have 2-3 times greater accuracy than 7-parameter transformation.Also, the transformation is part of different software programs in medical image computing for mapping points in one 3D image volume into their geometrically corresponding points in another, related 3D image volume [8].
General expression for 12-parameter transformation is [9] where X o represents the position vector in the output coordinate system, t is the translation vector, H is the skew matrix, S is the scale matrix, R is the rotation matrix, and X i is the position vector in the input coordinate system.Equation (72), written in matrix form, is

Transformation Parameters
To define a 12-parameter transformation, rotation matrix, scale matrix, and skew matrix are often combined in one matrix by mutual multiplication.Result of the multiplication is the matrix U, which can be written as [9] where parameters u ij along with the parameters t x , t y and t z represent 12 transformation parameters for the 12-parameter transformation.
Based on Equations ( 72) and (74), observation equations can be written [9]: Following Equation (75), design matrix A and observation vector l can be obtained as

Inverse Transformation
To perform the inverse transformation for the 12-parameter transformation, Equation (72) needs to be modified by expressing the element X i on the left side of the equation: where matrix R −1 is defined by Equations ( 54) or (55), inverse of the scale matrix is defined by Equation ( 71), and inverse of the skew matrix can be easily calculated: In a shorter form, Equation (77) can be written as

Calculation of Transformation Parameters
Using Equations ( 76) and ( 5), transformation parameters for 12-parameter transformations are calculated, and the procedure is shown in Appendix A.6.Finally, parameters of 12-parameter transformation from HDKS to HTRS96 based on 5000 identical points are shown in Table 7.

External Accuracy Assessment
After calculating the transformation parameters for all transformation models, an external accuracy assessment for each model based on 194 points in the territory of Croatia (Figure 1b) was performed.For each of the 194 points, the horizontal difference between measured coordinates in HTRS96 and corresponding coordinates in HTRS96 obtained by transformation from HDKS using the transformation parameters was calculated as follows: where dE is difference in the easting coordinate, and dN is the difference in the northing coordinate.
For the obtained horizontal differences, statistical indicators (minimum, maximum, average, and standard deviation) for all transformation models were calculated (Table 8).8 provide an average positional accuracy rating for the entire transformation area represented by a set of 194 control points.Based on these statistical indicators, it can be concluded that an increase in the number of parameters leads to a better transformation accuracy.Therefore, the 3-parameter and 5-parameter transformations yield the worst results, whereas the 9-parameter and 12-parameter transformations yield the best results.However, it should be noted that the 7-parameter and 8-parameter transformations do not yield significantly worse results than the 9-parameter transformation.Based on these statistical indicators, the accuracy rating for a specific transformation area of the Croatian territory cannot be observed.Therefore, a visual representation of the external accuracy assessment was created (Figure 8), allowing examination of areas of the Croatian territory with the best and worst results.

Discussion
Coordinate transformation is a very important and frequently used mathematical procedure in geodesy and geoinformatics when dealing with different coordinate systems.The aim of this article was to unify all coordinate transformations, explain the computation procedure of the transformation parameters, define the main transformation metadata, and explain the reversible transformation on the example of the Croatian territory.In Varga et al. [14], a similar issue was processed, but in this article, the computation procedure of transformation parameters was shown in a more detailed way.
In a coordinate transformation procedure, a set of transformation parameters is not sufficient to carry out the transformation unambiguously.Additional metadata about the rotation matrix, which is used to calculate transformation parameters, are also required.The first metadata is the convention that defines the rotation matrix.Although it is usual that CFR convention is used for transformation between global and local, and the PVT convention between two global coordinate systems, it is good to know which convention determines the parameters in case the user is unsure about the systems in question.The second important metadata is rotation matrix form, based on which the transformation parameters were calculated, so that the same rotation matrix can be used in transformation procedure.Using different forms of rotation matrices (simple rotation matrix, rotation matrix R XYZ , rotation matrix R ZYX ) can result in different output coordinates with errors at the centimeter level, as demonstrated in the numerical example of the 7-parameter transformation.This highlights the importance of this metadata.
Reversible transformation, in cases where it is necessary to switch from one system to another and then return to the original one, should be performed with a single set of parameters, wherein the general expression for the transformation must be inverted.In that way, upon returning to the original system, identical coordinates to the ones we started with will be obtained, whereas this will not be the case if two sets of parameters and the same expression for the transformation is used, as demonstrated by the numerical example in the case of a 7-parameter transformation.
Selecting a specific transformation model in practice largely depends on the expected results: more precisely, the level of accuracy that is aimed to achieve with the transformation.Models that use a smaller number of parameters typically yield results with lower accuracy, but their application is more straightforward.On the other hand, models using a larger number of parameters deliver higher accuracy results but are mathematically more complex.The application of 3-parameter and 5-parameter models is suitable when lower accuracy is acceptable and there is a need to transform a large amount of data.In practice, 7-parameter models are most commonly used because they effectively align reference frames at a global level and belong to conformal models while maintaining satisfactory accuracy.Theoretical and realized ellipsoid parameters differ, resulting in errors in the conversion of ellipsoidal coordinates to Cartesian, upon which transformation parameters are calculated.Due to the inability to define whether the minor and major semi-axes are proportionally different, the 8-parameter model is suitable because it assigns one scale factor to the minor semi-axis and another scale factor to the major semi-axis.
By calculating the transformation parameters of 3-parameter, 7-parameter, 8-parameter, 9-parameter, and 12-parameter transformations, similar results were obtained as in Varga et al. [14].The results are not exactly the same because of the different common points used for the computation of the transformation parameters.The main difference compared to Varga et al. [14] can be noticed in the 5-parameter Molodensky transformation, more precisely in the approach to calculate transformation parameters.In Varga et al. [14], parameters of the 5-parameter transformation are calculated in a way that values of δa i δ f are obtained from the definition of the Bessel1841 and GRS80 ellipsoids, whereas the translation parameters fully correspond with those from 3-parameter Molodensky transformations.In this article, transformation parameters are determined by the procedure from Deakin [4] based on standard and abridged Molodensky transformation formulas.In this procedure, all five transformation parameters are obtained by means of least squares estimation.Because of the amount of identical points and therefore the size of design matrix A and observation vector l, only the first block matrix of the calculation will be shown.
Taking into account that the first point is Matrix J can be written in more general form: J =   j 11 j 12 j 13 j 21 j 22 j 23 j 31 j 32 j 33   and, then, the inverse matrix is Determinant of matrix J can easily be obtained using the rule of Sarrus, and the adjoint matrix can be obtained in several steps.Each element j ij of matrix J has a minor m ij and a cofactor c ij .The minor, of all elements, is the determinant of elements of the matrix remaining after row i and column j are deleted.Cofactors of the matrix J (c ij = (−1) i+j m ij ) form cofactor matrix C, whose transpose is the adjoint matrix.Finally, the inverse matrix can be calculated from Equation (A3):   Taking into account that the first point is Because of the amount of points and therefore the size of design matrix A and observation vector l, only the first block matrix of the calculation will be shown.
Taking into account that the first point is Because of the amount of points and therefore the size of design matrix A and observation vector l, only the first block matrix of the calculation will be shown.
Taking into account that the first point is Because of the amount of points and therefore the size of design matrix A and observation vector l, only the first block matrix of the calculation will be shown.
Taking into account that the first point is

Figure 1 .
Figure 1.Common points within the territory of the Republic of Croatia: (a) points for calculation of transformation parameters; (b) points for external accuracy evaluation.


By including corresponding terms and solving partial derivatives, final expressions for matrices E, F, and G are Because of the amount of identical points and therefore the size of design matrix A and observation vector l, only the first block matrix of the calculation will be shown.