[Coordinate System Table of Contents]

3.5 Geographic Coordinate System Transformations

The geographic coordinate system transformations are another class of coordinate transformations. See a diagram on coordinate transformations to see how this class fits into the total schema of coordinate transformations.

3.5.1 Geographic Coordinate System Transformation Methods

It has been noted earlier that it is frequently required to convert coordinates derived in one geographic coordinate system to values expressed in another. For example, land and marine seismic surveys are nowadays most conveniently positioned by GPS satellite in the WGS 84 geographic coordinate system, whereas the national geodetic system in use for the country concerned will probably be a much earlier coordinate system. It may therefore be necessary to convert the observed WGS 84 data to the national geodetic system in order to avoid discrepancies. This form of Epicentre coordinate transformation has to be between source and target geographical coordinate systems with the coordinates of both expressed in geographical terms. It is not executed by a transformation between projected coordinates.

The transformations may be most readily achieved by first expressing the observed or source geographical coordinates in terms of three dimensional XYZ cartesian values (geocentric coordinates) instead of the normal angular expressions of latitude and longitude. However in order to derive the three true cartesian coordinates of a point on the earth's surface one must recognize that as well as having a latitude and longitude the point will also have an elevation above the ellipsoid. While XYZ cartesian coordinates may be derived from a mere latitude and longitude by assuming that the point lies on the ellipsoid's surface, in fact very few earth points actually do. Therefore the height of the point must be taken into account and the height required is the height above the ellipsoid. It is this height which is delivered by GPS in relation to the WGS84 ellipsoid. Heights above other ellipsoids are not generally immediately available. Instead the height that is usually available is the height above the national vertical datum, normally Mean Sea Level at a particular coastal point measured over a particular period. Levels will normally have been derived by conventional surveying methods and values will relate to the geoid surface. Hence in order to derive heights above the ellipsoid it is necessary to know the height of the geoid above the ellipsoid or vice versa and for most areas of the world this is not known with great precision. There exist various mathematical models of the geoid which have been derived for individual countries or parts of the world or the entire world and as satellite and terrestrial gravity data accumulates they are being steadily refined. The best available geoidal data should be used to convert surveyed heights (surveyed other than by GPS) to ellipsoidal heights for use in Epicentre's geodetic coordinate transformation.

In the event that the true ellipsoidal height of a point cannot be obtained the assumption that it is zero will often allow a geodetic transformation without introducing significant error in the horizontal coordinates.

In the early days of satellite surveying when relationships between datums were not well defined and the data itself was not very precise, it was usual to apply merely a three parameter dX, dY, dZ shift to the XYZ coordinate set in one datum to derive those in the second datum. This assumed, generally erroneously, that the axial directions of the two ellipsoids involved were parallel. For localized work in a particular country or territory, the consequent errors introduced by this assumption were small and generally less than the observation accuracy of the data. Nevertheless, as knowledge and data has accumulated and surveying methods have become more accurate, it has become evident that a three parameter transformation is neither appropriate for world wide use, nor indeed for widespread national use if one is seeking the maximum possible accuracy from the satellite surveying and a single set of transformation parameters.

For petroleum exploration purposes it may well be that a three parameter transformation within a particular area of interest is quite adequate but it should not be assumed that the same transformation parameters are appropriate for use with different data in an adjoining area of interest or for another purpose.

The simplest transformation to implement involves applying shifts to the three geocentric coordinates. Molodenski developed a transformation which applies the geocentric shifts directly to geographical coordinates. Both of these methods assume that the axes of the source and target systems are parallel to each other. As indicated above this assumption may not be true and consequently these transformation methods result in only moderate accuracy, especially if applied over large areas.

Improved accuracy can be obtained by applying a Helmert 7-parameter transformation to geocentric coordinates. However there are two opposing sign conventions for the three rotation parameters. In Epicentre these are considered to be two different transformation methods, - either a Position Vector Transformation or a Coordinate Frame Transformation. (The Position Vector Transformation is also called the Bursa-Wolf transformation). It is crucial that the signs of the rotation parameters are consistent with the rotation convention being applied.

It is also possible to interpolate geographical or grid shifts for points on the basis of known shifts for a number of control points in a specific defined area. One such application is the coordinate transformation introduced to enable the conversion of coordinates expressed in the North American Datum of 1927 (including the Clarke 1866 ellipsoid) to coordinates expressed in the new North American Datum of 1983 which takes the GRS 1980 ellipsoid. Because the North American survey control network was built by conventional terrestrial survey observations and suffers from the inevitable instrumental and adjustment shortcomings of the time, the old network, based on the non geocentric Clarke 1866 ellipsoid and a single datum point at Meades Ranch in Kansas, is not wholly consistent when compared with data which can be more readily and accurately secured nowadays with the advantages of satellite technology, modern instrumentation, and electronic computational techniques. Hence to convert between coordinates on the old system to values in the new datum, it is not appropriate to merely apply the type of orthogonal transformation represented by the Molodenski or Helmert transformations described above. Different points in different parts of the North American continent need to undergo different positional or coordinate shifts according to their position within the continental network. This is known in Epicentre as the NADCON transformation method. The conversion is achieved by Bilinear gridded interpolation for the new NAD 83 Latitudes and Longitudes using US Coast & Geodetic Survey NADCON control point grids. Bilinear interpolation is also used in Canada for the same purpose using the National Transformation (NT) application and grids. As with the US, longitude differences are applied to longitudes which are positive west. The NAD27 and NAD83 geographical coordinate systems documented by EPSG use the positive east longitude convention. The Canadian gridded file format has recently been adopted by Australia.

Alternatively a polynomial expression with listed coefficients for both latitude and longitude may be used as the transformation method. A transformation, applicable to offshore Norway to effect the transformation between coordinates expressed in the imperfect European Datum 1950 (ED50) and the newer ED87 uses this approach. Statens Kartwerk, the Norwegian Survey authority, publishes a document which lists 15 coefficients for each of separate latitude and longitude polynomial conversion formulas, involving up to fourth order expressions of latitude and longitude expressed in degrees.

If and when other countries decide to convert to using Geocentric Datums to suit satellite positional data acquisition methods, coordinate transformations similar to NADCON or the Norwegian example may well be employed to facilitate the conversion between the old terrestrial survey derived coordinates and the new geocentric datum values.

Note that it is very important to ensure that the signs of the parameters used in the transformations are correct in respect of the transformation being executed. Preferably one should always express transformations in terms of "From".........."To"........... thus avoiding the confusion which may result from interpreting a dash as a minus sign or vice versa.

3.5.2 Geographic Coordinate System Transformation Formulas

3.5.2.1 Geographic/Geocentric Conversions

Latitude, j, and Longitude, l, in terms of Geographic Coordinate System may be expressed in terms of a geocentric (earth-centered) cartesian coordinate system (X, Y, Z) with the Z axis corresponding with the Polar axis positive northwards, the X axis through the intersection of the Greenwich meridian and equator, and the Y axis through the intersection of the equator with longitude 90°E.

Geocentric coordinate systems are conventionally taken to be defined with the X axis through the intersection of the Greenwich meridian and equator. This requires that the equivalent geographic coordinate system is based on the Greenwich meridian. In application of the formulas below, geographic coordinate systems based on a non-Greenwich prime meridian should first be transformed to their Greenwich equivalent.

If the earth's spheroidal semi major axis is a, semi minor axis b, and inverse flattening 1/f, then

where is n the prime vertical radius of curvature at latitude j and is equal to

j and l are respectively the latitude and longitude (related to the prime meridian) of the point

h is height above the ellipsoid, (see note below), and

e is the eccentricity of the ellipsoid where

(Note: that h is the height above the ellipsoid.. This is the height value which is delivered by GPS satellite observations but is not the gravity-related height value which is normally used for national mapping and levelling operations. The gravity-related height (H) is usually the height above mean sea level or an alternative level reference for the country. If one starts with a gravity-related height H, it will be necessary to convert it to an ellipsoid height h before using the above transformation formulas. h = H + N, where N is the geoid height above the ellipsoid at the point and is sometimes negative. The geoid is a gravity surface approximating mean sea level. For the WGS84 ellipsoid the value of N, representing the height of the geoid relative to the ellipsoid, can vary between values of -100m in the Sri Lanka area to +60m in the North Atlantic. Geoid heights of points above the nationally used ellipsoid may not be readily available.)

For the reverse conversion, Cartesian coordinates in the geocentric coordinate system may then be converted to the geographical coordinates in terms of geographic coordinate system by:

by iteration

where l_B is relative to Greenwich. If the geographic system has a non-Greenwich prime meridian, the Greenwich value of the local prime meridian should be applied to longitude.

See section 3.5.2.4 below for an example.

3.5.2.2 Offsets

Several transformation methods utilising offsets in coordinate values are recognised. These include longitude rotations, geographical coordinate offsets and vertical offsets.

Mathematically, if the origin of a one-dimensional coordinate system is shifted along the positive axis and placed at a point with ordinate A (where A>0), then the transformation formula is:

  X_new = X_old - A

However it is common practice in coordinate system transformations to apply the shift as an addition, with the sign of the shift parameter value having been suitably reversed to compensate for the practice. Since 1999 this practice has been adopted by EPSG. Hence transformations allow calculation of coordinates in the target system by adding a correction parameter to the coordinate values of the point in the source system:

  X_t = X_s + A

where X_s and X_t are the values of the coordinates in the source and target coordinate systems and A is the value of the transformation parameter to transform source coordinate system coordinate to target coordinate system coordinate.

3.5.2.3 Geocentric Translations

If we assume that the minor axes of the ellipsoids are parallel and that the prime meridian is Greenwich, then shifts dX, dY, dZ in the sense from source geocentric coordinate system to target geocentric coordinate system may then be applied as

Xt	=	Xs + dX
Yt	=	Ys + dY
Zt	=	Zs + dZ

Example:

This example combines the geographic/geocentric transformation of section 2.2.1 above with a geocentric transformation. Consider a North Sea point with coordinates derived by GPS satellite in the WGS84 geographical coordinate system, with coordinates of:

latitude js	=	53 deg 48 min 33.82 sec N
longitude ls	=	2 deg 07 min 46.38 sec E
ellipsoidal height hs	=	73.0 m

whose coordinates are required in terms of the ED50 geographical coordinate system which takes the International 1924 ellipsoid. The three parameter geocentric translations transformation from WGS84 to ED50 for this North Sea area is given as dX = +84.87m, dY = +96.49m, dZ = +116.95m. The WGS84 geographical coordinates first convert to the following geocentric values using the formulas from section 3.5.2.1:

Xs =	3 771 793.97 m
Ys =	140 253.34 m
Zs =	5 124 304.35 m

Applying the quoted geocentric translations to these, we obtain new geocentric values now related to ED50:

Xt	=	3 771 793.97 + 84.87	=	3 771 878.84 m
Yt	=	140 253.34 + 96.49	=	140 349.83 m
Zt	=	5 124 304.35 + 116.95	=	5 124 421.30 m

Using the reverse conversion given in section 3.5.2.1 above, these convert to ED50 values on the International 1924 ellipsoid as:

latitude jt	=	53 deg 48 min 36.565 sec N
longitude lt	=	2 deg 07 min 51.477 sec E
ellipsoidal height ht	=	28.02 m

Note that the derived height is referred to the International 1924 ellipsoidal surface and will need a further correction for the height of the geoid at this point in order to relate it to Mean Sea Level.

3.5.2.4 Abridged Molodenski Transformation

As an alternative to the above computation of the new latitude, longitude and height above ellipsoid in discrete steps, the changes in these coordinates may be derived directly by formulas derived by Molodenski. Abridged versions of these formulas, which are quite satisfactory for three parameter transformations, are as follows:

where the dX, dY and dZ terms are as before, and r and n are the meridian and prime vertical radii of curvature at the given latitude j on the first ellipsoid (see section 3.4.4), da is the difference in the semi-major axes of the target and source ellipsoids (a_t - a_s) and df is the difference in the flattening of the two ellipsoids (f_t - f_s).

The formulas for dj and dl indicate changes in j and l in arc-seconds.

3.5.2.5 Helmert Transformation

Transformation of coordinates from one geographic coordinate system into another (also colloquially known as a “datum transformation”) is usually carried out as an implicit concatenation of three transformations:

[geographical to geocentric >> geocentric to geocentric >> geocentric to geographic]

The middle part of the concatenated transformation, from geocentric to geocentric, is usually described as a simplified 7-parameter Helmert transformation, expressed in matrix form with 7 parameters, in what is known as the "Bursa-Wolf" formula:

The parameters are commonly referred to defining the transformation "from source coordinate system to target coordinate system", whereby (X_S, Y_S, Z_S) are the coordinates of the point in the source geocentric coordinate system and (X_T, Y_T, Z_T) are the coordinates of the point in the target geocentric coordinate system. But that does not define the parameters uniquely; neither is the definition of the parameters implied in the formula, as is often believed. However, the following definition, which is consistent with the “Position Vector Transformation” convention, is common E&P survey practice:

(dX, dY, dZ): Translation vector, to be added to the point's position vector in the source coordinate system in order to transform from source system to target system; also: the coordinates of the origin of source coordinate system in the target coordinate system

(R_X, R_Y, R_Z): Rotations to be applied to the point's vector. The sign convention is such that a positive rotation about an axis is defined as a clockwise rotation of the position vector when viewed from the origin of the Cartesian coordinate system in the positive direction of that axis; e.g. a positive rotation about the Z-axis only from source system to target system will result in a larger longitude value for the point in the target system. Although rotation angles may be quoted in any angular unit of measure, the formula as given here requires the angles to be provided in radians.

M: The scale correction to be made to the position vector in the source coordinate system in order to obtain the correct scale in the target coordinate system. M = (1 + dS*10^-6), whereby dS is the scale correction expressed in parts per million.

Example:

This example combines the geographic/geocentric transformation of section 3.5.2.1 above with a position vector transformation.

Transformation from WGS 72 to WGS 84 (EPSG transformation code 1238). Transformation parameters:

dX =	0.000 m
dY =	0.000 m
dZ =	+4.5 m
RX =	0.000 seca
RY =	0.000 seca
RZ =	+0.554 seca
dS =	+0.219 ppm

Input point coordinate system: WGS 72 (geographic 3D)

Latitude jS	=	55 deg 00 min 00 sec N
Longitude lS	=	4 deg 00 min 00 sec E
Ellipsoid height hS	=	0 m

Using the geographic to geocentric transformation method given in section 3.5.2.1, this transforms to Cartesian geocentric coordinates:

XS =	3 657 660.66 m
YS =	255 768.55 m
ZS =	5 201 382.11 m

Application of the 7 parameter Position Vector Transformation (code 1238) results in

XS =	3 657 666.78 m
YS =	255 778.43 m
ZS =	5 201 387.75 m

Using the reverse formulas for the geographic/geocentric transformation method given in section 3.5.2.1 this converts into:

Latitude jT	=	55 deg 00 min 00.090 sec N
Longitude lT	=	4 deg 00 min 00.554 sec E
Ellipsoid height hT	=	+3.22 m

Although being common practice particularly in the European E&P industry, the Position Vector Transformation sign convention is not universally accepted. A variation on this formula is also used, particularly in the USA E&P industry. That formula is based on the same definition of translation and scale parameters, but a different definition of the rotation parameters. The associated convention is known as the "Coordinate Frame Rotation" convention.

The formula is:

and the parameters are defined as:

(dX, dY, dZ): Translation vector, to be added to the point's position vector in the source coordinate system in order to transform from source system to target system; also: the coordinates of the origin of source coordinate system in the target coordinate system

(R_X, R_Y, R_Z): Rotations to be applied to the coordinate frame. The sign convention is such that a positive rotation of the frame about an axis is defined as a clockwise rotation of the coordinate frame when viewed from the origin of the Cartesian coordinate system in the positive direction of that axis, that is a positive rotation about the Z-axis only from source coordinate system to target coordinate system will result in a smaller longitude value for the point in the target coordinate system. Although rotation angles may be quoted in any angular unit of measure, the formula as given here requires the angles to be provided in radians.

M: The scale factor to be applied to the position vector in the source coordinate system in order to obtain the correct scale of the target coordinate system: M = (1 + dS*10^-6), whereby dS is the scale correction expressed in parts per million.

In the absence of rotations the two formulas are identical; the difference is solely in the rotations. The name of the second method reflects this.

Note that the same rotation that is defined as positive in the first method is consequently negative in the second and vice versa. It is therefore crucial that the convention underlying the definition of the rotation parameters is clearly understood and is communicated when exchanging datum transformation parameters, so that the parameters may be associated with the correct coordinate transformation method (algorithm).

The same example as for the Position Vector Transformation can be calculated, however the following transformation parameters have to be applied to achieve the same input and output in terms of coordinate values:

dX =	0.000 m
dY =	0.000 m
dZ =	+4.5 m
RX =	0.000 seca
RY =	0.000 seca
RZ =	-0.554 seca
dS =	+0.219 ppm

Please note that only the rotation has changed sign as compared to the Position Vector Transformation.

3.5.3 Polynomial Transformations

Note: In the sections that follow, the general mathematical symbols X and Y representing the axes of a coordinate system must not be confused with the specific axis abbreviations or axis order in particular coordinate systems.

3.5.3.1 General Case

Polynomial transformations between two coordinate systems are typically applied in cases where one or both of the coordinate systems exhibits lack of homogeneity in orientation and scale. The small distortions are then approximated by polynomial functions in latitude and longitude or in easting and northing. Depending on the degree of variability in the distortions, approximation may be carried out using 2^nd, 3_{rd, or higher degree polynomials. In the case of transformations between two projected coordinate systems, the additional distortions resulting from the application of two map projections and a datum transformation can be included in a single polynomial approximation function.}

Polynomial approximation functions themselves are subject to variations, as different approximation characteristics may be achieved by different polynomial functions. The simplest of all polynomials is the general polynomial function. In order to avoid problems of numerical instability this type of polynomial should be used after reducing the input parameters, usually coordinate offsets U and V relative to a central evaluation point, to 'manageable' numbers, between -10 and +10 at most.

U = X_S - X_S0 in defined units (which may not be those of the coordinate system)
V = Y_S - Y_S0

Then

(X_T - X_T0) = (X_S - X_S0) + dX
(Y_T - Y_T0) = (Y_S - Y_S0) + dY

or

X_T = X_S - X_S0 + X_T0 + dX
Y_T = Y_S - Y_S0 + Y_T0+ dY

where

X_T, Y_T are coordinates in the target coordinate system,
X_S, Y_S are coordinates in the source coordinate system,
X_S0, Y_S0 are coordinates of the evaluation point in the source coordinate system,
X_T0, Y_T0 are coordinates of the evaluation point in the target coordinate system.

and where

dY =	B0	+B1U+B2V+B3U2+B4U.V+B5V2	(to 2nd order)
		+B6U3+B7U2.V+B8U.V2+B9V3	(3rd order terms)
		+B10U4+B11U3.V+B12U2.V2+B13U.V3+B14V4	(4th order terms)

Other (arguably better) approximating polynomials are described in mathematical text books such as “Theory and applications of numerical analysis”, by G.M. Phillips and P.J. Taylor (Academic Press, 1973).

Polynomial Reversibility

Approximation polynomials are in a strict mathematical sense not reversible, i.e. in principle the same polynomial coefficients cannot be used to execute the reverse transformation. However, in practical application there are exceptions when applied to the approximation of small differences between the geometry of a set of points in two different coordinate systems. An example is the polynomial approximation between ED50 and ED87 geographical coordinate systems in the North Sea. The typical difference in coordinate values between ED50 and ED87 is in the order of 2 metres (»10^-6 degrees), with a maximum of about 3 metres. This holds for the application area of 52ºN to 62ºN and 3ºW to 9ºE. The polynomial functions are evaluated about a central point, 55ºN, 0ºE. The offsets with respect to this central or evaluation point, in degrees, are used as input parameters to the polynomial functions. The output of the polynomial functions are corrections to the input coordinates of about 10^-6 degrees. This difference of several orders of magnitude between input and output values is the property that makes a polynomial function reversible in practice (although not in a formal mathematical sense).

Suppose one would want to check the transformation from ED50 to ED87 at a point in which both coordinate sets are known. If the proper input coordinate offsets are used, the result will be a correction vector, which, if added to the ED50 coordinates, would yield the ED87 coordinates. The correction itself is in the order of 1 to 2 metres

Let us evaluate what happens if the user accidentally uses the ED87 coordinate offsets as input instead of the ED50 offsets. The input values now contain an error in the order of 1 to 2 metres, or approximately 10^-6 degrees. Consequently the output, the corrections, will be 10^-6 *10^-6 degrees or 1 micron in error. This is a perhaps somewhat complex way of saying that the corrections from ED50 to ED87 at two neighboring points are approximately equal.

Now what happens if the transformation from ED50 to ED87 is applied with the correct input values, derived from the ED50 coordinates of the point, but with the signs of the coefficients reversed? The result would be a correction of the same magnitude, but with changed sign.

The final step is to evaluate what would happen if the polynomial transformation were applied with reversed sign coefficients and with the ED87 coordinates of the point as input. The result would be a correction vector with reversed sign, but which would be wrong by a factor of 10^-6 or 1 micron which is not significant.

The last described scenario is the reverse transformation, from ED87 to ED50.

The error made by the polynomial approximation formulas in calculating the reverse correction is of the same order of magnitude as the ratio of output versus input:

output error     output value
---------------  »   --------------
input error       output error

As long as the input values (the coordinate offsets from the evaluation point) are orders of magnitude larger than the output (the corrections), and provided the coefficients are used with changed signs, the polynomial transformation may be considered to be reversible.

EPSG recognises two classes of general polynomial functions, distinguished by whether or not the coefficients may be used in both forward and reverse transformations, i.e. are reversible. Note that in the reversible polynomial class, the signs of the coefficients but not of the evaluation point coordinates are reversed.

Because the coordinate differences are small and the coordinates in the source and target systems are very similar in magnitude, in reversible polynomial transformations it is possible to adopt the same coordinate values for the evaluation point in both source and target systems. As the evaluation point coordinates are the same in source and target coordinate systems, the points are at two physically different locations. The difference in position between these is accounted for within the zero order coefficient values of the polynomial.

Then

U = X_S - X₀ in defined units (which may not be those of the coordinate system)
V = Y_S - Y₀

Then

X_T = X_S + dX
Y_T = Y_S + dY

Example: 4^th order reversible polynomial

For geodetic transformation ED50 to ED87 (1).

Offset unit: degree
Ordinate 1 of evaluation point X₀ = 55º 00' 00.000"N = +55 deg.
Ordinate 2 of evaluation point Y₀ =  0º 00' 00.000"E =  +0deg

Parameters:

A0 = -5.56098E-06	A1 = -1.55391E-06	…	A14 = -4.01383E-09
B0 = +1.48944E-05	B1 = +2.68191E-05	…	B14 = +7.62236E-09

Forward calculation for:

Latitude jED50	=	XS	=	52º 30' 30" N	=	+52.508333333 deg
Longitude lED50	=	YS	=	2º E	=	+2.0 deg

U = X_S - X₀ = j_ED50 - X₀ = 52.508333333 - 55.0 = -2.4916666667 deg
V = Y_S - Y₀ = l_ED50 - Y₀ = 2.0 - 0.0 = 2.0 deg

dX	=	A0 + A1U +...+ A14V4
	=	-5.56098E-06 +(-1.55391E-06 * -2.491666667) + ... + (-4.01383E-09 * 2.04)
	=	-3.12958E-06 degreesTD>

dY	=	B0 + B1U +...+ B14V4
	=	+1.48944E-05 +(2.68191E-05 * -2.491666667) + ... + (7.62236E-09 * 2.04)
	=	+9.80126E-06 degrees

Then

Latitude jED87	=	XT	=	XS + dX
			=	52.508333333 - 3.12958E-06 deg
			=	52º 30' 29.9887" N
Longitude lED87	=	YT	=	YS + dY
			=	2º 00' 00.0353" E

Reverse Calculation:

Transformation ED50 to ED87 (1)

The transformation method for the ED50 to ED87 (1) transformation, 4th-order reversible polynomial, is reversible. The same formulas may be applied for the reverse calculation, but coefficients A₀ through A₁₄ and B₀ through B₁₄ are applied with reversal of their signs. Sign reversal is not applied to the coordinates of the evaluation point. Thus

Ordinate 1 of evaluation point X₀ = 55º 00' 00.000"N = +55 deg.
Ordinate 2 of evaluation point Y₀ =  0º 00' 00.000"E =  +0deg

Parameters:

A0 = +5.56098E-06	A1 = +1.55391E-06	…	A14 = +4.01383E-09
B0 = -1.48944E-05	B1 = -2.68191E-05	…	B14 = -7.62236E-09

Reverse calculation for

Latitude jED87	=	XS	=	52º 30' 29.9887" N	=	+52.5083301944 deg
Longitude lED50	=	YS	=	2º 00' 00.0353" E	=	+2.0000098055 deg

U = 52.5083301944 - 55.0 = -2.4916698056 deg
V = 2.0000098055 - 0.0 = 2.0000098055 deg

dX	=	A0 + A1U +...+ A14V4
	=	+5.56098E-06 +(1.55391E-06 * -2.491666667) + ... + (4.01383E-09 * 2.04)
	=	+3.12957E-06 degreesTD>

dY	=	B0 + B1U +...+ B14V4
	=	-1.48944E-05 +(-2.68191E-05 * -2.491666667) + ... + (-7.62236E-09 * 2.04)
	=	-9.80124E-06 degrees

Then

Latitude jED50	=	XT	=	XS + dX
			=	52.5083301944 + 3.12957E-06 deg
			=	52º 30' 30.000" N
Longitude lED50	=	YT	=	YS + dY
			=	2º 00' 00.000" E

3.5.3.2 Polynomial Transformation with Complex Numbers

The relationship between two projected coordinate systems may be approximated more elegantly by a single polynomial regression formula written in terms of complex numbers. The advantage is that the dependence between the 'A' and 'B' coefficients (for U and V) is taken into account in the formula, resulting in fewer coefficients for the same order polynomial. A third-order polynomial in complex numbers is used in Belgium. A fourth-order polynomial in complex numbers is used in The Netherlands for transforming coordinates referenced to the Amersfoort / RD system to and from ED50 / UTM.

(dX + i.dY) =	(A1+i.A2)(U+iV) + (A3+i.A4)(U+iV)2	(to 2nd order)
	+(A5+i.A6)(U+iV)3	(3rd order terms)
	+(A7+i.A8)(U+iV)4	(4th order terms)

where U = (X_S - X_X0).10^-5
and V = (Y_S - Y_X0).10^-5

The fourth order polynomial can alternatively be expressed in matrix form as

Then as for the general polynomial case in section 3.5.3.1 above,

X_T = X_S - X_S0 + X_T0 + dX
Y_T = Y_S - Y_S0 + Y_T0+ dY

and where, as above,

X_T, Y_T are coordinates in the target coordinate system,
X_S, Y_S are coordinates in the source coordinate system,
X_S0, Y_S0 are coordinates of the evaluation point in the source coordinate system,
X_T0, Y_T0 are coordinates of the evaluation point in the target coordinate system.

Note that the zero order coefficients of the general polynomial, A₀ and B₀, have apparently disappeared. In reality they are absorbed by the different coordinates of the source and of the target evaluation point, which in this case, are numerically very different because of the use of two different projected coordinate systems for source and target.

The transformation parameter values (the coefficients) are not reversible. For the reverse transformation a different set of parameter values are required, used within the same formulas as the forward direction.

Example:
Projected coordinate transformation: Amersfoort / RD New to ED50 / UTM zone 31N (1)

For input point:

Easting,   X_{AMERSFOORT/RD} = X_S = 200,000.00 m
Northing, Y_{AMERSFOORT/RD} = Y_S = 500,000.00 m

U = (X_S-X_S0).10^-5 = (200,000 - 155,000).10^-5 = 0.45
V = (Y_S-Y_S0).10^-5 = (500,000 - 463,000).10^-5 = 0.37

dX = -1,240.050
dY = 1,468.748

Then

Easting	EED50/TM31N	= XT	= XS - XS0 + XT0 + dX
			= 200,000 - 155,000 + 663,395.607 + (-1240.050)
			= 707,155.557 m
Northing	NED50/TM31N	= YT	= YS - YS0 + YT0 + dY
			= 500,000 - 463,000 + 5,781,194.380 + 1,468.748
			= 5,819,663.128 m

3.5.3.3 Polynomial Transformation for Spain

The original geographic coordinate system for the Spanish mainland was based on Madrid 1870 datum, Struve 1860 ellipsoid, with longitudes related to the Madrid meridian. Three second-order polynomial expressions have been empirically derived by El Servicio Geográfico del Ejército to convert geographical coordinates based on this system to equivalent values based on the European Datum of 1950 (ED50). The polynomial coefficients derived can be used to convert from Madrid 1870 (Madrid) geographic coordinate system to the ED50 system. Three pairs of expressions have been derived: each pair is used to calculate the shift in latitude and longitude respectively for

1. a mean for all Spain,
2. a better fit for the north of Spain,
3. a better fit for the south of Spain.

The polynomial expressions transformations are:

dj seca = A₀ + (A₁*j_s) + (A₂*l_s) + (A₃*H_s)
dl seca = B₀₀ + B₀ + (B₁*j_s) + (B₂*l_s) + (B₃*H_s)

where latitude j_s and longitude l_s are in decimal degrees referred to the Madrid 1870 (Madrid) geographic coordinate system and H_s is the height in metres. B₀₀ is the longitude (in arc-seconds) of the Madrid meridian measured from the Greenwich meridian. It is the value to be applied to a longitude relative to the Madrid meridian to transform it to a longitude relative to the Greenwich meridian.

The results of these expressions are applied through the formulas:

	jED50	= jM1870(M) + dj
and	lED50	= lM1870(M) + dl

Example:

Input point, coordinate system: Madrid 1870 (Madrid) (geographic 2D)

Latitude js	=	42 deg 38 min 52.77 sec N
	=	+42.647992 dega
Longitude ls	=	3 deg 39 min 34.57 sec E of Madrid
	=	+3.659603 dega from the Madrid meridian.
Gravity-related height Hs	=	0 m

For the north zone transformation:

	B00 = -13276.58
A0 = 11.3287790	B0 = 2.5079425
A1 = -0.1674	B1 = 0.8352
A2 = -0.03852	B2 = -0.00864
A3 = 0.0000379	B3 = -0.0000038

The resulting calculations:

	dj	=	+4.05 seca
Then latitude	jED50	=	42 deg 38 min 52.77 sec N + 4.05 seca
		=	42 deg 38 min 56.82 sec N
	dl	=	-13270.54 seca = -3 deg 41 min 10.54 sec
Then longitude	lED50	=	3 deg 39 min 34.57 sec E - 3 deg 41 min 10.54 sec
		=	0 deg 01 min 35.97 sec W of Greenwich

3.5.4 Miscellaneous Linear Transformations

An affine 2D transformation is used for transforming a coordinate system possibly with non-orthogonal axes and possibly different units along the two axes to an isometric coordinate system (i.e. a system of which the axes are orthogonal and have equal scale units). The transformation therefore involves a change of origin, differential change of axis orientation and a differential scale change.

EPSG distinguishes four transformation methods to implement this class of transformation:

1. the parametric representation
2. the geometric representation, which is described in
   1. a general case,
   2. a simplified case in which the axes are constrained to be orthogonal, and
   3. a further simplified case known as the similarity transformation.

3.5.4.1 The affine parametric transformation

Mathematical and survey literature usually provides the parametric representation of the affine transformation. The parametric algorithm is commonly used for rectification of digitised maps. It is often embedded in CAD software and Geographical Information Systems where it is frequently referred to as "rubber sheeting." Although the application of this algorithm falls outside the scope of the EPSG geodesy data set it is presented here for reasons of clarity.

The formula in matrix form is as follows:

or, using algebraic coefficients:

where

X_T, Y_T are the coordinates of a point P in the target coordinate sytem;
X_S, Y_S are the coordinates of P in the source coordinate sytem;

This form of describing an affine transformation is analogous to the general polynomial transformation formulas (section 3.5.3.1). Although it is somewhat artificial, an affine transformation can be considered to be a first order general polynomial transformation but without the reduction to source and target evaluation points.

Reversibility

The parameter values for an affine transformation cannot be used for the reverse transformation. However, the reverse transformation is another affine transformation using the same formulas but with different parameter values. The reverse parameter values, indicated by a prime ('), can be calculated from those of the forward transformation as follows:

D = A₁B₂ - A₂B₁
A₀' = (A₂B₀ - B₂A₀)/D
B₀' = (B₁A₀ - A₁B₀)/D
A₁' = +B₂/D
A₂' = -A₂/D
B₁' = -B₁/D
B₂' = +A₁/D

Geometric representation of the affine transformation. Please note that to prevent cluttering of the figure, the scale parameters of the X_S and Y_S axes have been omitted.

From the diagram above, it can be seen that

Scale of the two coordinate systems:

The scaling of both source and target coordinate systems adds some complexity to this formula.

The transformation will often be applied to transform a local coordinate system to a projected 2D coordinate system. The local coordinate system, e.g. a seismic bin grid, may have different units of measure on its two axes. These have scale ratios of dS_X and dS_Y respective to the axes of the projected coordinate system.

The projected 2D coordinate system is nominally defined to be in well-known units, e.g. metres. However, the distortion characteristics of the map projection only preserve true scale along certain defined lines or curves, hence the projected coordinate system's unit of measure is strictly speaking only valid along those lines or curves. Everywhere else its scale is distorted by the map projection. For conformal map projections the distortion at any point can be expressed by the point scale factor 'k' for that point. Please note that this point scale factor 'k' should NOT be confused with the scale factor at the natural origin of the projection, denominated by 'k₀'. For non-conformal map projections the scale distortion at a point is bearing-dependent and will not be described in this document. It has developed as working practice to choose the origin of the source (local) coordinate system as the point in which to calculate this point scale factor 'k', although for local coordinate systems with a large coverage area a point in the middle of the coverage area may be a better choice.

Adding the scaling between each pair of axes and dropping the suffix for point P, after rearranging the terms we have the geometric representation of the general affine transformation:

or, in matrix form:

where:

XT0, YT0	=	the coordinates of the origin point of the source coordinate system, expressed in the target coordinate system.
dSX, dSY	=	the length of one unit of the source axis, expressed in units of the target axis, for the first and second source and target axis pairs respectively;
k	=	point scale factor of the target coordinate system at a chosen reference point;
qX, qY	=	the angles about which the source coordinate system axes XS and YS must be rotated to coincide with the target coordinate system axes XT and YT respectively (counter-clockwise being positive).

Comparing the algebraic representation with the parameters of the parameteric form in section 3.5.4.1 above it can be seen that the parametric and geometric forms of the affine transformation are related as follows:

A₀ = X_T0
A₁ = k dS_X cosq_X
A₂ = k dS_Y sinq_Y
B₀ = Y_T0
B₁ = -k dS_X sinq_X
B₂ = k dS_Y cosq_Y

Reversibility:

The parameters for an affine transformation cannot be used for the reverse transformation. However, the reverse transformation is another affine transformation using the same formulas but with different parameters. The reverse parameter values can be calculated from the formulae provided above and applying those to same algorithm. Alternatively the reverse transformation can be described by a different formula, as shown below, using the same parameters as the forward transformation:

where D = cos(q_X-q_Y), or, algebraically:

3.5.4.3 The affine orthogonal geometric transformation

If the source coordinate system happens to have orthogonal axes, that is both axes are rotated through the same angle to bring them into the direction of the orthogonal target coordinate system axes, i.e. q_X = q_Y = q, then the general geometric representation can be simplified to the orthogonal geometric affine transformation. In matrix form this is:

or algebraically:

where

XT0, YT0	=	the coordinates of the origin point of the source coordinate system, expressed in the target coordinate system.
dSX, dSY	=	the length of one unit of the source axis, expressed in units of the target axis, for the X axes and Y axes respectively;
k	=	point scale factor of the target coordinate system at a chosen reference point;
q	=	the angle through which the source coordinate system axes must be rotated to coincide with the target coordinate system axes (counter-clockwise is positive). Alternatively, the bearing (clockwise positive) of the source coordinate system YS-axis measured relative to target coordinate system north.

For the reverse case the general geometric affine formulas given above can be similarly simplified by replacing q_X and q_Y with q.

See below for example.

3.5.4.4 The similarity transformation

If the source coordinate system happens to have axes of the same scale, that is both axes are scaled by the same factor to bring them into the scale of the target coordinate system axes (i.e. dS_X = dS_Y = dS) then the orthogonal geometric affine transformation can be simplified further to a similarity transformation:

Similarity transformation

From the above diagram, the similarity transformation in algebraic form is:

Dropping the suffix for point P and rearranging the terms:

or, in matrix form:

where

XT0, YT0	=	the coordinates of the origin point of the source coordinate system, expressed in the target coordinate system.
1+dS	=	the length of one unit in the source coordinate system expressed in units of the target coordinate systm;
q	=	the angle about which the axes of the source coordinate system axes must be rotated to coincide with the axes of the target coordinate system axes (counter-clockwise is positive). Alternatively, the bearing of the source coordinate system YS-axis measured relative to target coordinate system north.

The similarity transformation can also be described as a special case of the parametric affine transformation where coefficients A₁ = B₂ and A₂ = -B₁, and the expression, A₁B₂-A₂B₁ = (1+dS).

Reversibility

In contrast with the affine transformation, the similarity transformation parameters are reversible, but only on the condition that the scale difference between the two coordinate systems is small (order of several parts per million). Then dS is the deviation from unity of the ratio of the units of measure of the two coordinate systems. In these cases the reverse transformation would require a scale correction of 1/(1+dS) » (1-dS). This enables usage of the same scale and rotation parameters, but with reversed sign, for the reverse transformation. (The rotation angle of +q becomes -q, which is valid for all q). However for the reverse transformation the translation parameters, X_T0 and Y_T0, take entirely different values.

When to use the similarity transformation

Similarity transformations can be used when source and target systems

• each have orthogonal axes
• each have the same scale along both axes
  and
• both have the same units of measure.

for example between engineering plant grids and projected coordinate systems.

A transformation between two systems that

• each have orthogonal axes,
  but
• significantly different scales or units of measure,

should be defined as an orthogonal affine transformation and not as a similarity transformation. This may be the case with seismic bingrids with square bins.

Transformations between two coordinate systems where in either system either the scale along axes differ or the axes are not orthogonal should be defined as a general affine transformation in either the parametric or geometric form.

Examples:

Similarity transformation method

Source coordinate system:	Astra Minas Grid (local coordinate system)
Target coordinate system:	Campo Inchauspe / Argentina 2 (projected 2D system)

Note that for the Astra Minas Grid, the coordinate axes are:
X (positive axis oriented north)
Y (positive axis oriented west) and coordinates are quoted in that order.

whereas for Campo Inchauspe / Argentina 2 the axes are:
X (positive axis oriented north)
Y (positive axis oriented east) and coordinates are quoted in that order.

Thus the Astra Minas grid X and Y axes map to the Campo Inchauspe / Argentina 2 Y and X axes respectively. With respect to the symbols in the formulas,
X_S = Astra Minas X
Y_S = Astra Minas Y
X_T = Campo Inchauspe / Argentina 2 Y
Y_T = Campo Inchauspe / Argentina 2 X

Parameters of the similarity transformation:
X_T0 = 2610200.48 m
Y_T0 = 4905282.73 m
q = 271º 05' 30" = 271.0916667 dega
k = 0 thus (1 + k) = 1.0

Forward calculation for Astra Minas point:

X(north) = 10000 m
Y(west) = 50000 m

X_S = Astra Minas X = 10000
Y_S = Astra Minas Y = 50000

Gauss-Kruger zone 2 Easting (Y)	=	XT = XT0 + XS.dS.cosq + YS.dS.sinq
	=	2610200.48 + (50000*1.0*cos(271.0916667deg)) + (10000*1.0*sin(271.0916667deg))
	=	2601154.90 m
Gauss-Kruger zone 2 Northing (X)	=	YT = YT0 - XS.dS.sinq + YS.dS.cosq
	=	4905282.73 - (50000*1.0*sin(271.0916667deg)) + (10000*1.0*cos(271.0916667deg))
	=	4955464.17 m

Affine transformation (orthogonal geometric) method

Source coordinate system: imaginary 3D seismic acquisition bin grid. The two axes are orthogonal, but the unit on the I-axis is 25 metres, while the unit on the J-axis is 12.5 metres. The target projected coordinate system is WGS 84 / UTM Zone 31N and the origin of the bin grid (centre of bin 0,0) is defined at E = 456781.0, N = 5836723.0. The projected coordinate system point scale factor at the bin grid origin is 0.99984.

The map grid bearing of the I and J axes are 110º and 20º respectively. Thus the angle through which both the positive I and J axes need to be rotated to coincide with the positive Easting axis and Northing axis respectively is +20 degrees

Hence:

XT0	=	456 781.0 m
YT0	=	5 836 723.0 m
dSX	=	25
dSY	=	12.5
k	=	0.99984
q	=	+20 dega

Forward calculation for center of bin with coordinates: I = 230, J = 247:

XT	=	Easting	=	XT0 + XS.k.dSX.cosq + YS.k.dSY.sinq
			=	464 855.62 m
YT	=	Northing	=	YT0 - XS.k.dSX.sinq + YS.k.dSY.cosq
			=	5 837 055.90 m

Reverse calculation for this point:

XS	=	[(XT - XT0).cosqY - (YT - YT0).sinqY]/[k.dSX.cos(qX-qY)]
	=	230 bins
YS	=	[(XT - XT0).sinqX + (YT - YT0).cosqY]/[k.dSX.cos(qX-qY)]
	=	162 bins

---

[Coordinate System Table of Contents]

---

Last modified: 13 Sept 2000
© Copyright 1998-2000 POSC. All rights reserved.