Transforming Coordinates from One Conformal Grid to Another : Lauf’s Method
Background
There is nothing new in the mathematics at the core of this article. It appeared, however, that the explanation and means to perform the mathematics to achieve such a coordinate transformation was now unavailable for historic data. This article seeks to rectify both of these issues. Furthermore, as transforming historic coordinate positions such as those on the Transverse Mercator (TM), yard grid of the Australian R502 map series to equivalent GPS metric coordinates was complex, an acceptable transformation solution, as advanced by Lauf, is evaluated.
In 1961, a paper titled Conformal Transformations from One Map Projection to Another, using Divided Difference Interpolation by Gordon Bertram Lauf (1914-1984) was published in the Bulletin Geodesique (International Union of Geodesy and Geophysics – IUGG). The paper was based on some six years of work which was initially published as Lauf’s Doctor of Science thesis in 1951. Subsequently in 1983 this methodology was published in Geodesy and Map Projections, commissioned by the then Royal Melbourne Institute of Technology, and in the National Mapping Council’s 1986, Special Publication 10, The Australian Geodetic Datum – Technical Manual. Lauf’s Method, as it became known, permitted the conformal transformation of coordinates from any conformal or orthomorphic map projection to any other conformal map projection.
Lauf’s Method
The usefulness of the Lauf Method lay in its ability to undertake the coordinate transformation without any knowledge of the basis of the grids being used. For instance, as stated by Lauf, the method does not depend upon a knowledge of the latitudes of the standard parallels, the longitude of the central meridian, the scale factor nor upon the orientation of the two systems. It does not even require a knowledge of the unit of measurement used nor even of the type of projection employed. Lauf warned, however, it must be emphasised that in the strictest sense, if no information of any kind about the two projections is available, no unique solution of the problem is possible. Any arbitrarily chosen set of coordinates for the result (within limits) could be shown to be consistent with some projection, even if it were purely fanciful.
Lauf (1983) mentioned the 1940s work of Vladimir Kirilov Hristov (also written as Wl. K. Hristow, 1902-1979), founder and first Director of the Central Laboratory of Geodesy at the Bulgarian Academy of Sciences (BAS). Hristov developed a method for transforming Gauss Conformal Coordinates from one grid zone to the adjacent zone, providing that both projections were orthomorphic or conformal transformations from the same spheroid and that the rectangular coordinates of several points common to both systems were also known. Wilhelm Jordan and Otto Eggert published in their 1941 Handbuch der Vermessungskunde (Surveyors Handbook), a method of transformation based largely on Hristov's method. The Jordan-Eggert method was then adopted for use on the Australian Map Grid (AMG66) and appropriate tables computed. Formulae, tables and proforma calculation sheets were given in the 1972 Australian Map Grid - Technical Manual, Special Publication No.7, section 5.12. Jordan’s Handbook of Geodesy by the US Army Map Service provided an English translation of the Jordan-Eggert work and a useful insight into the theory behind Lauf’s work.
Lauf described his method as a bivariate Interpolation into a table of complex numbers. Simplistically, with several common control points having known grid coordinates in each projection these coordinate pairs were converted to complex numbers and their divided differences taken. This process set up the transformation parameters which could then be used to convert the grid coordinates of any number of points between the projections.
The example can be thought of as having an old map grid and a new map grid that when constructed used different conformal projections to represent that portion of the Earth’s surface. Each projection, however, retained the angles between directions but distances could be distorted; Mercator and Lambert were two such projections. A common point in both map grids had the coordinate pair y and x in the old grid and Y and X in the new grid (the y and Y axis was north south and x and X axis was east west, by definition). The grid values of y,x and Y,X in each projection came from a function of ψ and λ, where ψ was the isometric latitude and λ the longitude of a point (the isometric latitude was a nonlinear function of the geodetic latitude proportional to the meridian distance and the perpendicular distance from the axis of revolution to the point concerned; it was a dimensionless quantity). Thus, the coordinate pairs in each grid were a function of the same point’s latitude and longitude on the Earth’s surface, represented in complex number form as :
y + ix = f (ψ + iλ); and
Y + iX = F (ψ + iλ).
By moving from the real plane to the complex (or Argand) plane, any coordinate pair say a,b in the real plane may be thought of as being a + ib in the complex plane, any polynomial associated with this value could be solved even if the square root of a negative number occurred. Complex numbers also had an advantage over ordinary vectors in that they obeyed the usual commutative, associative and distributive laws of arithmetic. In representing a coordinate pair as a complex number the i had no value and the term +i could be considered as a separator between the two values of the complex pair. Also i squared equalled -1 when complex numbers were algebraically combined and for brevity, z = (a + ib). (The development of complex numbers is credited to 16th century Italian Gerolamo Cardano (1501-1576), a most influential mathematician of the Renaissance, who proved that having a negative term inside a square root could lead to the solution of an equation. Up until then, it was thought to be impossible to find the square root of a negative number. Later, in the 18th century, mathematician Carl Friedrich Gauss (1777-1856) consolidated Cardano’s work.)
Thus from the theory of functions of a complex variable, as the coordinate pair Y and X in the new grid and the coordinate pair y and x in the old grid came from the same geographical point having latitude ψ and longitude λ, then the coordinate pair Y and X in the new grid were a function of the coordinate pair y and x in the old grid or :
Y + iX = f(y+ ix) alternatively Z = f(z).
With the coordinates of four control points common to both projected grids converted to complex pairs, a table could then be constructed as shown below.
|
Common point |
Real coordinate pair |
Complex coordinate pair |
Divided Differences of the Complex Numbers |
||||
|
Old grid |
New grid |
Old grid |
New grid |
First |
Second |
Third |
|
|
P1 |
y1,x1 |
Y1,X1 |
z1 = y1+ix1 |
Z1 = Y1+iX1 |
|
|
|
|
|
|
|
|
|
ΛZ11 = (Z2-Z1)/(z2-z1) |
|
|
|
P2 |
y2,x2 |
Y2,X2 |
z2 = y2+ix2 |
Z2 = Y2+iX2 |
|
ΛZ21 = (ΛZ12- ΛZ11)/(z3-z1) |
|
|
|
|
|
|
|
ΛZ12 = (Z3-Z2)/(z3-z2) |
|
ΛZ31= (ΛZ22- ΛZ21)/(z4-z1) |
|
P3 |
y3,x3 |
Y3,X3 |
z3 = y3+ix3 |
Z3 = Y3+iX3 |
|
ΛZ22 = (ΛZ13- ΛZ12)/(z4-z2) |
|
|
|
|
|
|
|
ΛZ13 = (Z4-Z3)/(z4-z3) |
|
|
|
P4 |
y4,x4 |
Y4,X4 |
z4 = y4+ix4 |
Z4 = Y4+iX4 |
|
|
|
This table now contained the transformation parameters from which any coordinates in the old grid (zo) could be converted to coordinates in the new grid (Zn) by :
Zn = Z1 + (zo-z1)ΛZ11 + (zo-z1) (zo-z2)ΛZ21 + (zo-z1) (zo-z2) (zo-z3)ΛZ31
alternatively, and as means of a check for hand calculation :
Zn = Z4 + (zo-z4)ΛZ13 + (zo-z4) (zo-z3)ΛZ22 + (zo-z4) (zo-z3) (zo-z2)ΛZ31.
Completing the table and deriving a final result required the subtraction, division and multiplication of complex numbers. These processes could be accomplished using the basic operations as described in Annex A. When calculators and computers became available these machines often had subroutines to handle complex number operations which sped up the data processing. Additionally, in the earliest days of these machines the use of complex numbers had the benefit in requiring reduced storage and programming. A coordinate pair could be stored as a complex number in one location instead of two locations as for real numbers and a complex number could be manipulated in one step instead of operating individually on each real number of the coordinate pair. Microsoft EXCEL today handles the task with ease. The basic algebraic operations may be programmed into EXCEL but EXCEL also has inbuilt complex number functions. Lauf’s Method using EXCEL’s complex number functions was utilised in the analysis that follows.
A Practical Use of Lauf’s Method
Australia’s first national map coverage was the mainly uncontoured R502 series. This map series comprised 540 printed map sheets at a scale of 1: 250 000 essentially controlled horizontally by astrofixes and with spot heights determined by barometric heighting techniques. Detail was plotted from 1: 50 000 scale aerial photography using a Transverse Mercator (TM) projection and a corresponding yard grid (Clarke 1858 spheroid) which covered the whole of Australia in 8 zones. Each zone was 5 degrees of longitude wide with a half degree of common overlap. This was a simple projection with no provision for a scale factor and each zone’s true origin at 34 degrees south latitude. False origins were 400 000 yards west and 800 000 yards south of the true origin. This projection was used for the R502 map series until replaced in 1966 when the Australian Map Grid (AMG) was adopted. After adoption of the AMG, almost half the R502 maps were overprinted with the 10 000 metre AMG grid in cyan. A diagram of the 8 zone, Transverse Mercator projection parameters may be viewed via this link.
Following the completion and adjustment of the Australian Geodetic Survey, Australia adopted the Australian National Spheroid (ANS) and as its survey and mapping coordinate system the Australian Geodetic Datum (AGD66). The accompanying National Topographic Map Series (NTMS) of 1: 250 000 scale maps utilised the metric Universal Transverse Mercator (UTM) projection as the basis for the Australian Map Grid (AMG). The AMG was in metres with a central scale factor of 0.9996 superimposed on all projected distances. The zones were 6 degrees of longitude wide with 0.5 degrees of overlap, numbered from 47 (central meridian at 99°E) to 58 (central meridian at 165°E), the true origin of each zone being the intersection of its central meridian longitude with the equator. False origins were 500 000 metres west and 10 000 000 metres south of the true origin.
To later conform to global positioning technology the Geocentric Datum of Australia or GDA94 was established along with the Map Grid of Australia (MGA94). MGA retained the same basic parameters as the AMG. A diagram of the 8 zone, Universal Transverse Mercator (UTM) projection over mainland Australia may be viewed via this link.
Below is a summary of Australian datums as published in the Intergovernmental Committee on Surveying and Mapping, Geocentric Datum of Australia 2020 Technical Manual.
Summary of the parameters of Australian datums.
|
Datum |
Geographic Coordinate Set |
Grid Coordinates |
Ellipsoid /Spheroid |
Semi-major axis (m) |
Inverse Flattening |
|
Clarke 1858 |
Sydney Observatory |
TM |
Clarke 1858 |
6378293.645 |
294.26 |
|
AGD66 |
AGD66 |
AMG66 |
ANS |
6378160.000 |
298.25 |
|
GDA94 |
GDA94 |
MGA94 |
GRS80 |
6378137.000 |
298.257222101 |
The grid coordinate shift between Clarke 1858, in yards, and either of the other two metric systems was such that the only way to convert from Clarke 1858 grid coordinates was to first convert to geographical coordinates (latitude, longitude) and then into the required system. The Lauf Method was shown in the National Mapping Council’s 1986, Special Publication 10, The Australian Geodetic Datum – Technical Manual, to be able to transform grid coordinates to AMG66 from Clarke 1858 to sub-metre accuracy avoiding the intermediate computations.
As the R502 series of maps and coordinate positions existed until the late 1960s project/survey work having Clarke 1858 grid coordinates in yards still likely exist. Given such historical project data (#), could GPS be used to now locate such historical points? Furthermore, as common coordinates generally only exist for geodetic positions how many points would be needed to achieve an acceptable accuracy for the transformed coordinates? (#: other reasons may also require GPS grid coordinates for historical R502 grid positions.)
Grid Transformation Evaluation
An R502 map sheet was considered to be an internally consistent unit of area for analysis. Such a 1.5 x 1 degree R502 map sheet area encompassed some 15 000 square kilometres. Lauf himself stated that four common stations are sufficient for surveys up to 10 000 kilometres square. For still larger areas, five common stations could be used, or else, the whole area should be broken down into two or more smaller areas. In a particular case, covering about 10 000 kilometres square with more than 100 common stations, it was found that there was no significant improvement in the results when a complex polynomial of degree four or higher was used. In the abovementioned example from the National Mapping Council’s 1986, Special Publication 10, The Australian Geodetic Datum – Technical Manual, the transformation to AMG66 from Clarke 1858 coordinates to sub-metre accuracy, the four common control points used encompassed an area of some 9 000 square kilometres.
The SH54-15 Broken Hill R502 reprint of 1981 was selected as the test area. The map sheet contained a spread of common geodetic control and in addition to the 10 000 yard TM grid had been overprinted (in cyan) with the UTM 10 000 metre grid. A point could be quickly visually inspected for coordinates in both systems. Twelve common points which had grid coordinates in the TM (Clarke 1858), AMG66 and GDA94 (GPS) systems were chosen. With 12 points it would enable 5, 4 and 3 points to be used as the common control with the remainder for accuracy checks on the transform.
The table and R502 map sheet below show the twelve point’s grid coordinates and their relative locations.
|
ID |
POINT NAME |
CLARKE 1858-E&N |
AGD66-E&N |
GDA94 (GPS)-E&N |
|||
|
1 |
GAP |
476057.573 |
1150713.828 |
569307.040 |
6558423.634 |
569428.708 |
6558602.154 |
|
2 |
MOORKAIE |
467100.193 |
1062560.232 |
561120.729 |
6477848.117 |
561242.520 |
6478026.837 |
|
3 |
ROBE |
434020.851 |
1083977.734 |
530885.091 |
6497423.679 |
531006.864 |
6497602.317 |
|
4 |
RAMPARTS |
541596.785 |
1112438.291 |
629212.456 |
6523439.104 |
629334.269 |
6523617.843 |
|
5 |
ACACIA |
500860.517 |
1112429.311 |
591978.472 |
6523430.403 |
592100.104 |
6523609.095 |
|
6 |
DERING |
463789.963 |
1104509.422 |
558094.666 |
6516190.791 |
558216.352 |
6516369.401 |
|
7 |
BLUFF |
482656.516 |
1131355.178 |
575339.150 |
6540729.223 |
575460.810 |
6540907.829 |
|
8 |
GAIRDNER |
473753.586 |
1089398.600 |
567201.850 |
6502379.210 |
567323.564 |
6502557.851 |
|
9 |
STEPHEN |
451370.213 |
1064856.192 |
546743.290 |
6479946.431 |
546865.064 |
6480125.106 |
|
10 |
SCROPES |
548652.630 |
1066664.602 |
635662.258 |
6481600.750 |
635784.078 |
6481779.489 |
|
11 |
FELSPAR |
416655.987 |
1021565.584 |
515014.374 |
6440377.047 |
515136.358 |
6440555.748 |
|
12 |
NTH BARRIER |
446654.173 |
1153819.113 |
542431.009 |
6561261.550 |
542552.664 |
6561440.165 |
Note that the position of Felspar, depicted in the map above, is correct in relation to the other points.
The points used for common control were selected from the twelve available on the basis of them generally encompassing the area containing some of the remaining points and leaving others as outliers to check if there was any increased inaccuracy away from the controlled region. Thus the quadrilateral formed by points Gap, Moorkaie, Robe and Ramparts enclosed an area of some 8 000 square kilometres. The point, Acacia, at about the centre of the previous quadrilateral was added as the fifth control point. Points, Gap, Robe and Ramparts, formed the common control for the three point option.
The analysis transformed TM (Clarke 1858) grid coordinates in yards to both AMG66 and GDA94. The transformations to GDA94 yielded little difference between the transformed coordinates on AMG66, thus the discussion below may be considered as applicable to both grid type transformations. Please refer to Annex B for the EXCEL worksheets and their use.
|
Common control points used – 4 |
Common control points used – 5 |
||||||||||||
|
ID |
Common |
Difference E |
Difference N |
ID |
Common |
Difference E |
Difference N |
ID |
Common |
Difference E |
Difference N |
||
|
1 |
GAP |
1 |
GAP |
1 |
GAP |
||||||||
|
3 |
ROBE |
2 |
MOORKAIE |
2 |
MOORKAIE |
||||||||
|
4 |
RAMPARTS |
3 |
ROBE |
3 |
ROBE |
||||||||
|
4 |
RAMPARTS |
4 |
RAMPARTS |
||||||||||
|
5 |
ACACIA |
||||||||||||
|
Required |
Required |
Required |
|||||||||||
|
6 |
DERING |
-0.07 |
-0.11 |
|
6 |
DERING |
0.1 |
-0.1 |
|
6 |
DERING |
-0.1 |
-0.2 |
|
7 |
BLUFF |
0.14 |
0.08 |
|
7 |
BLUFF |
0.2 |
0.1 |
|
7 |
BLUFF |
0.1 |
0.1 |
|
8 |
GAIRDNER |
-0.20 |
-0.05 |
|
8 |
GAIRDNER |
0.0 |
0.0 |
|
8 |
GAIRDNER |
-0.1 |
0.0 |
|
9 |
STEPHEN |
-0.04 |
0.13 |
|
9 |
STEPHEN |
0.1 |
-0.1 |
|
9 |
STEPHEN |
0.0 |
-0.1 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
10 |
SCROPES |
0.05 |
0.07 |
|
10 |
SCROPES |
-0.2 |
0.7 |
|
10 |
SCROPES |
-1.2 |
-0.1 |
|
11 |
FELSPAR |
0.13 |
0.40 |
|
11 |
FELSPAR |
-0.5 |
-1.0 |
|
11 |
FELSPAR |
-2.7 |
-2.1 |
|
12 |
NORTH BARRIER |
-0.19 |
0.14 |
12 |
NORTH BARRIER |
-0.1 |
-0.1 |
12 |
NORTH BARRIER |
-0.5 |
0.3 |
||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Mean RMS Value Pts 6,7,8&9 |
|
|
|
|
|
|
|
|
|
|
||
|
|
|
0.2 |
|
|
|
0.2 |
|
|
|
|
|
0.1 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Mean RMS Value Pts 10,11&12 |
|
|
|
|
|
|
|
|
|
|
||
|
|
|
0.3 |
|
|
|
0.8 |
|
|
|
|
|
2.1 |
|
During the three transformations, 7 points (IDs 6-12) were always considered as the required points. Of these points, IDs 6-9 were generally within or close to the area enclosed by the common control points, while points, IDs 10-12 were in the order of 30 kilometres outside the controlled region. The table above summarises the three transformations showing the fit of the transformed coordinates to all seven required points. To see any trend(s) Root Mean Square (RMS) values were calculated and averaged for the groups of points deemed to be enclosed or exterior to the controlled region. It was then clear that for points 6-9 that were generally within or close to the area enclosed by the common control points, that with five common control points the transformation was marginally more accurate. For the group of points exterior to the controlled region, however, with four and five common control points the transformation accuracy fell away significantly.
Closing remarks
This analysis would suggest that for achieving the most reliable results from Lauf’s Method that :
|
- |
four points be selected for use as common control such that they generally encompass the area containing the points whose coordinates are required to be transformed;
|
|
- |
If only three common control points are available they should also generally encompass the area containing the points whose coordinates are required to be transformed, otherwise their spread should be as wide as possible across the region of interest to give the best possible transformation accuracy;
|
|
- |
five common control points will give the most accurate transformation provided that none of the points whose coordinates are required to be transformed fall significantly outside the perimeter of the controlled area;
|
|
- |
in all cases it is important to try to get the best spread of the common control points as three well spread control points is better than five points in a cluster;
|
|
- |
In all cases if more than the minimum number of common control points are available then one or more excess points may be used as quasi unknowns to check transformation accuracy. |
Lauf in his paper commented on one of his examples where the unit of measurement in one grid was unknown and the second used metres, that from the values given in the First Order Divided Differences Table, the real value gave an indication of the unit of measurement used in that first grid. Above, the TM (Clarke 1858) grid coordinates in yards were transformed to the MGA94 metre grid. Noting that 1 yard is 0.9144 metres the real First Order Divided Difference was 0.91403 (to five places of decimal) in the above. However, the MGA94 grid applied a fixed scale factor of 0.9996 (exactly) to all distances meaning the factor 0.9144 was reduced to 0.91403. In old grids where today unknowingly feet, miles or even chains may have been used, such a scale conversion indicator might prove useful.
The above analysis used the coordinates of geodetic survey control for the common control and accuracy points; this being the best case test scenario. Coordinates may also be scaled from appropriate maps to use as the common control being aware, however, that scaling errors will carry through the transformation and impact the final result. In such cases the final GPS grid coordinates could be expected to be of the same order of accuracy as the scaled input coordinates irrespective of the accuracy of the GPS device used.
This work has been a desktop analysis only and results have not been checked on the ground.
Acknowledgement
The contribution of former AUSLIG surveyors Doug White and Dan Jaksa is acknowledged.
Compiled by Paul Wise, 2023.
Sources
Bomford, Guy (1971), Geodesy, 3rd Edn., Oxford at the Clarendon Press.
De Jong, Sybren Hendrik (1968), The Skew Mercator Projection in Rectangular Plane Coordinate Systems, Ohio State University Thesis, accessed at :
http://rave.ohiolink.edu/etdc/view?acc_num=osu1486647101628433
Hristov, Vladimir Kirilov (1943), Die Gauss-Krüger 'schen Koordinaten Auf Dem Ellipsoid (The Gauss-Krüger Coordinates on the Ellipsoid), Leipzig / Berlin : BG Teubner.
Lauf, Gordon Bertram and Young, F (1961), Conformal Transformations from One Map Projection to Another, using Divided Difference Interpolation, Bulletin Geodesique (International Union of Geodesy and Geophysics – IUGG, September 1961), No.61, pp.191-212.
Lauf, Gordon Bertram (1983), Transformation of Coordinates from One Map Projection to Another, Geodesy and Map Projections, pp.119-124, Commissioned by the Department of Surveying, Royal Melbourne Institute of Technology.
National Mapping Council (1986), The Australian Geodetic Datum – Technical Manual.
Wise, Paul Joseph (2017), Gordon Bertram Lauf : A brief profile.
Basic Operations with Complex Numbers
Addition : add together the real parts (without i) and add up the imaginary parts (with i).
(a+bi) + (c+di) =
(a+c) + (b+d)i
Subtraction : subtract the real parts and subtract the imaginary parts (with i).
(a+bi) - (c+di) =
(a-c) + (b-d)i
Multiplication : to multiply two complex numbers use the distributive law, noting i^2= -1 (i squared equals -1).
(a+bi)
* (c+di) = (ac-bd) + (ad+bc)i
Division : the division of two complex numbers is achieved by multiplying the numerator and denominator by the denominator's complex conjugate.
(a+bi)/(c+di)
= (a+bi)*(c-di)
/ (c+di)*(c-di) where (c+di)(c-di) becomes (c^2+d^2)
(^2 means squared) giving
(a+bi)/(c+di) = [(ac+bd) / (c^2+d^2)] + [(bc−ad) / (c^2+d^2)] i
Grid To Grid Transform by Lauf – EXCEL Workbook
[Please read the XNATMAP Disclaimer before use]
The EXCEL Workbook (please read the XNATMAP Disclaimer before use) contains three separate sheets. Depending on whether three, four or five common points are available there is a calculation for each option.
Each sheet contains red and blue cells. The green calculation cells are hidden to avoid overwriting.
The red cells are for user input and the blue cells are the results.
At the top, the red cells are for the grid coordinates of the common points. Below that there are several rows where the grid coordinates of the required points may be entered with the transformed coordinates appearing (in blue) alongside.