[This Guidance Note is available in its entirety below, but users are encouraged to download Guidance Note 7, in PDF format or in a MS Word 97-2000 & 6/95 - RTF format (zipped) for improved presentation.]

Guidance Note Number 7

POSC literature pertaining to

Coordinate conversions and coordinate transformations.

1. Map Projections and their Coordinate Conversion Formulas

1.1 Introduction

1.2 Identification of Projected Coordinate Reference Systems

3. Map Projection Parameters
4. Map Projection (Coordinate Conversion) Formulas
1. Lambert Conic Conformal
2. Lambert Conic Near-Conformal
3. Krovak Oblique Conic Conformal
4. Mercator
5. Cassini-Soldner
6. Transverse Mercator
7. Oblique Mercator
8. Stereographic
9. New Zealand Map Grid
10. Tunisia Mining Grid
11. American Polyconic
12. Lambert Azimuthal Equal Area
13. Albers Equal Area
14. Equidistant Cylindrical (Plate Carree)
15. Bonne

2. Formulas for Coordinate Operations other than Map Projections

1. Introduction
2. Coordinate Conversions other than Map Projections
1. Geographic/Geocentric Conversions
3. Coordinate Transformations
1. Offsets
2. Geocentric Translations
3. Molodenski Transformation
4. Helmert Transformation
5. Molodenski-Badedas Transformation
6. Interpolation Methods
4. Coordinate Operation Methods that can be Operations or Transformations
1. Polynomial Transformations
1. General case
2. Polynomial transformations with complex numbers
3. Polynomial transformation for Spain
2. Miscellaneous Linear Coordinate Transformations
1. Affine General Parametric Transformation
2. Affine General Geometric Transformation
3. Affine Orthogonal Geometric Transformation
4. Similarity Transformation

Part 1. Map Projections and their Coordinate Conversion formulas

1.1 Introduction

Setting aside the large number of map projection methods which may be employed for atlas maps, equally small scale illustrative exploration maps, and wall maps of the world or continental areas, EPSG provides reference parameter values for orthomorphic or conformal map projections which are used for medium or large scale topographic or exploration mapping. Here accurate positions are important and sometimes users may wish to scale accurate positions, distances or areas from the maps.

Small scale maps normally assume a spherical earth and the inaccuracies inherent in this assumption are of no consequence at the usual scale of these maps. For medium and large scale sheet maps, or maps and map projection coordinates held digitally to a high accuracy, it is essential that due regard is paid to the actual shape of the earth. Such coordinate reference systems are therefore invariably based on an ellipsoid and its derived map projections. The EPSG data set and this supporting conversion documentation considers only map projections for the ellipsoid.

Though not exhaustive the following list of named map projection methods are those which are most frequently encountered for medium and large scale mapping, some of them much less frequently than others since they are designed to serve only one particular country. They are grouped according to their possession of similar properties, which will be explained later. Except where indicated all are conformal.

Mercator Cylindrical

with one standard parallel

with two standard parallels

Cassini-Soldner (N.B. not conformal) Transverse Cylindrical

Transverse Mercator Group

Transverse Mercator (including south oriented version) Transverse Cylindrical

Universal Transverse Mercator

Gauss-Kruger

Gauss-Boaga

Oblique Mercator Group

Hotine Oblique Mercator Oblique Cylindrical

Oblique Mercator

Laborde Oblique Mercator

New Zealand Map Grid (May also be considered as type of Oblique Stereographic)

Space Oblique Mercator (Not strictly conformal)

Lambert Conical Conformal

with one standard parallel Conical

with two standard parallels

one standard parallel with truncated formulas (n.b. not conformal)

Stereographic Azimuthal

Polar

Oblique and equatorial

1.2. Identification of Map Projection method

If a map or coordinate list is provided for which an EPSG listed coordinate reference system is clearly identifiable, then its name or identifier together with EPSG dataset version number will address the required parameters including the coordinate conversion parameter values. If the coordinate reference system is not listed it will be necessary to create a new coordinate reference system with its own coordinate conversion (parameter) values.

It may often happen that one is presented with a coordinate list or map for which the author or compiler has regrettably failed to provide any indication of parameter values or properties:- no projection name, no grid definition and no statement of ellipsoid or datum. On the map there may be no grid or graticule, or indeed neither.

In order to adequately relate the digital or displayed map data to other data it is necessary to establish the properties of the data or given map from what may be gleaned from their appearance and other information. Geographical coordinates without qualifying information do not allow identification of the coordinate reference system other than that it is a geographic one. Projected or map grid coordinates may, by virtue of the actual and relative magnitudes of the Easting and Northing and knowledge of where in the world they relate to, provide clues as to the map projection. For example eastings between say 150000m and 850000m, allied with 6 or 7 figure northings correlated with latitude may indicate a UTM.

If the map bears neither grid nor graticule it will be useless unless one can identify a number of the point features shown for which one already has coordinate data. One may then be able to superimpose and fit a rectangular grid at appropriate scale from which other coordinate data may be read. If the map carries a grid then the numerical labelling of the grid lines, the assumption that it will be conformal or orthomorphic, and prior knowledge of approximately where in the world it covers may give some indication of the type of projection, but this may not be totally definitive. If the map bears a graticule the nature of the graticule lines will give some indication of the type of projection used in its compilation. For example straight meridians and concentric parallels would suggest a conical projection or, less frequently, a polar azimuthal. If the former, and assuming that it will be orthomorphic, then it will either be with one standard parallel or two and these will have been selected in relation to the latitudinal extent of the area, very possibly those in general use for that state's mapping. If the parallels are equally spaced it will be a simple equidistant conical projection. However for large scale mapping purposes the requirement that it is conformal will dictate that the parallels will not be equally spaced and it is more than likely that it will be some form of Lambert projection with either one or two standard parallels. Unfortunately there is no easy way of detecting which, nor the values of the standard parallels. The country it comes from and its national mapping system, if known, may suggest what these are. The EPSG data set will assist but is not exhaustive.

If both parallels of latitude and meridians are curved the projection has numerous possibilities but a form of Transverse Mercator may well be the most likely. One may attempt to identify the projection by computing the grid values of some of the graticule intersections for several possible projections in turn, plotting these to a rounded value for the estimated scale of the map e.g. 1:50000 or 1:100000, and attempting to fit the overlaid grid plot on the graticule. Repeating this for a number of potential projections for the area may be successful in obtaining a reasonable fit. But bear in mind that paper stretch may slightly distort scale from the nominal scale of the map, and the scale factor used in the graticule to grid conversions is another variable which may take only slightly different values e.g. a Gauss-Kruger takes a central meridian scale factor of unity while a UTM (like Gauss-Kruger, a Transverse Mercator) takes 0.9996.

Digital cartographic techniques make it relatively easy to plot grid and graticule for different projections with different parameters onto transparencies for "trial and error" overlays. The process can be time consuming so it is preferable to make maximum use of the clues which one may infer from the appearance of the map as initially presented, - its origins, its national area, and the conventional projections used for that area.

1.3. Map Projection parameters

A map projection grid is related to the geographical graticule of an ellipsoid through the definition of a coordinate conversion method and a set of parameters appropriate to that method. Differing conversion methods may require different parameters. Any one coordinate conversion method may take several different sets of associated parameter values, each set related to a particular map projection zone applying to a particular country or area of the world. Before setting out the formulas involving these parameters, which enable the coordinate conversions for the projection methods listed above, it is as well to understand the nature of these parameters.

The plane of the map and the ellipsoid surface may be assumed to have one particular point in common. This point is referred to as the natural origin. It is the point from which the values of both the geographical coordinates on the ellipsoid and the grid coordinates on the projection are deemed to increment or decrement for computational purposes. Alternatively it may be considered as the point which in the absence of application of false coordinates has grid coordinates of (0,0). For example, for projected coordinate reference systems using the Cassini-Soldner or Transverse Mercator methods, the natural origin is at the intersection of a chosen parallel and the chosen central meridian (see Figure 2 at end of section). The chosen parallel will frequently but not necessarily be the equator. For the stereographic projection the origin is at the centre of the projection where the plane of the map is imagined to be tangential to the ellipsoid.

Since the natural origin may be at or near the centre of the projection and under normal coordinate circumstances would thus give rise to negative coordinates over parts of the map, this origin is usually given false coordinates which are large enough to avoid this inconvenience. Hence each natural origin will normally have False Easting, FE and False Northing, FN values. For example, the false easting for the origins of all Universal Transverse Mercator zones is 500000m. As the UTM origin lies on the equator, areas north of the equator do not need and are not given a false northing but for mapping southern hemisphere areas the equator origin is given a false northing of 10,000,000m, thus ensuring that no point in the southern hemisphere will take a negative northing coordinate. Figure 4 illustrates the UTM arrangements.

These arrangements suggest that if there are false easting and false northing for the real or natural origin, there is also a Grid Origin which has coordinates (0,0). In general this point is of no consequence though its geographic position may be computed if needed. Sometimes however, rather than base the easting and northing coordinate reference system on the natural origin by giving it FE and FN values, it may be convenient to select a False Origin at a specific meridian/parallel intersection and attribute the false coordinates (0,0) or, more usually, EF and NF to this. The related easting and northing of the natural origin may then be computed if required.

The natural origin will always lie on a meridian of longitude. Longitudes are most commonly expressed relative to the Prime Meridian of Greenwich but some countries, particularly in former times, have preferred to relate their longitudes to a prime meridian through their national astronomic observatory, usually sited in or near their capital city, e.g. Paris for France, Bogota for Colombia. The meridian of the projection zone origin is known as the Longitude of Origin. For certain projection types it is often termed the Central Meridian or abbreviated as CM and provides the direction of the northing axis of the projected coordinate reference system.

Because of the steadily increasing distortion in the scale of the map with increasing distance from the origin, central meridian or other line on which the scale is the nominal scale of the projection, it is usual to limit the extent of a projection to within a few degrees of latitude or longitude of this point or line. Thus, for example, a UTM or other Transverse Mercator projection zone will normally extend only 2 or 3 degrees from the central meridian. For areas beyond this another zone of the projection, with a new origin and central meridian, needs to be used or created. The UTM system has a specified 60 numbered zones, each 6 degrees wide, covering the ellipsoid between the 84 degree North and 80 degree South latitude parallels. Other Transverse Mercator projection zones may be constructed with different central meridians, and different origins chosen to suit the countries or states for which they are used. A number of these are included in the EPSG dataset. Similarly a Lambert Conic Conformal zone distorts most rapidly in the north-south direction and may, as in Texas, be divided into latitudinal bands.

In order to further limit the scale distortion within the coverage of the zone or projection area, some projections introduce a scale factor at the origin (on the central meridian for Transverse Mercator projections), which has the effect of reducing the nominal scale of the map here and making it have the nominal scale some distance away. For example in the case of the UTM and some other Transverse Mercator projections a scale factor of slightly less than unity is introduced on the central meridian thus making it unity on two meridians either side of the central one, and reducing its departure from unity beyond these. The scale factor is a required parameter whether or not it is unity and is usually symbolised as k0.

Thus for projections in the Transverse Mercator group in section 1.1 above, the parameters which are required to completely and unambiguously define the projected coordinate system are:

Latitude of the natural origin

Longitude of the natural origin

Scale factor at the natural origin

False Easting

False Northing

Since the UTM zones obey set rules, it is sufficient to state only the UTM zone number (or central meridian). The remaining parameters from the above list are defined by the rules.

It has been noted that the Transverse Mercator projection is employed for the topographical mapping of longitudinal bands of territories, limiting the amount of scale distortion by limiting the extent of the projection either side of the central meridian. Sometimes the shape, general trend and extent of some countries makes it preferable to apply a single zone of the same kind of projection but with its central line aligned with the trend of the territory concerned rather than with a meridian. So, instead of a meridian forming this true scale central line for one of the various forms of Transverse Mercator, or the equator forming the line for the Mercator, a line with a particular azimuth traversing the territory is chosen, and the same principles of construction are applied to derive what is now an Oblique Mercator. This projection is sometimes referred to as the Hotine Oblique Mercator after the British geodesist who set out its formulas for application to Malaysian Borneo (East Malaysia) and also West Malaysia. Laborde had previously developed the projection system for Madagascar, and Switzerland uses a similar system derived by Rosenmund.

More recently (1974) Lee has derived formulas for a minimum scale factor projection for New Zealand known as the New Zealand Map Grid. The line of minimum scale follows the general alignment of the two main islands. This resembles an Oblique Mercator projection in its effect, but is not strictly an Oblique Mercator. The additional mathematical complexity of the projection enables its derivation via an Oblique Stereographic projection, which is sometimes the way it is classified. Because of its unique formulation inclusion of the New Zealand Map Grid within international mapping software was sporadic; as a consequence New Zealand has reverted to the frequently-encountered Transverse Mercator for its most recent mapping.

The parameters required to define an Oblique Mercator projection are:

Latitude of projection centre (the origin point on the initial line)

Longitude of projection centre

Azimuth of initial line at the projection centre

Scale factor at projection centre

Angle from rectified to skewed grid at the natural origin

and then either

False Easting (easting at the projection natural origin)

False Northing (northing at the projection natural origin)

or

The angle from rectified to skewed grid is normally applied at the natural origin of the projection, that is where the initial line of the projection intersects the aposphere. In some circumstances, for instance in the Alaskan panhandle State Plane zone, this angle is taken to be identical to the azimuth of the initial line at the projection centre. This results in grid and true north coinciding at the projection centre rather than at the natural origin as is more usual.

It is possible to define the azimuth of the initial line through the latitude and longitude of two widely spaced points along that line. This approach is not currently followed by POSC/EPSG.

For Conical map projections, which for the normal aspect may be considered as the projection of the ellipsoid onto an enveloping cone in contact with the ellipsoid along a parallel of latitude, the parallel of contact is known as a standard parallel and the scale is regarded as true along this parallel. Sometimes the cone is imagined to cut the ellipsoid with coincidence of the two surfaces along two standard parallels. All other parallels will be concentric with the chosen standard parallel or parallels but for the Lambert Conical Conformal will have varying separations to preserve the conformal property. All meridians will radiate with equal angular separations from the centre of the parallel circles but will be compressed from the 360 longitude degrees of the ellipsoid to a sector whose angular extent depends on the chosen standard parallel, - or both standard parallels if there are two. Of course the normal longitudinal extent of the projection will depend on the extent of the territory to be projected and will never approach 360 degrees.

As in the case of the Transverse Mercator above it is sometimes desirable to limit the maximum positive scale distortion for the one standard parallel case by distributing it more evenly over the extent of the mapped area. This may be achieved by introducing a scale factor on the standard parallel of slightly less than unity thus making it unity on two parallels either side of it. This achieves the same effect as choosing two specific standard parallels in the first place, on which the nominal scale will be preserved. The projection is then a Lambert Conical Conformal projection with two standard parallels. Although, strictly speaking, the scale on a standard parallel is always the nominal scale of the map and the scale factor on the one or two standard parallels should be unity, it is sometimes convenient to consider a Lambert Conical Conformal projection with one standard parallel yet which has a scale factor on the standard parallel of less than unity. This provision is allowed for by POSC/EPSG, where the single standard parallel is referred to as the latitude of the natural origin. For an ellipsoidal projection the natural origin will fall slightly poleward of the mean of the latitudes of the two standard parallels.

A longitude of origin or central meridian will again be chosen to bisect the area of the map or, more usually, the total national map area for the country or state concerned. Where this cuts the one standard parallel will be the natural origin of the projected coordinate system and, as for the Transverse Mercator, it will be given a False Easting and False Northing to ensure that there are no negative coordinates within the projected area (see Figure 5). Where two standard parallels are specified a false origin may be chosen at the intersection of a specific parallel with the central meridian outside the mapped area. This point will be given easting at false origin and northing at false origin to ensure that no negative coordinates will result. Figure 6 illustrates these arrangements.

It is clear that any number of Lambert projection zones may be formed according to which standard parallel or standard parallels are chosen and this is clearly exemplified by those which are used for many of the United States State Plane coordinate zones. They are normally chosen either, for one standard parallel, to approximately bisect the latitudinal extent of the country or area or, for two standard parallels, to embrace most of the latitudinal extent of the area. In the latter case the aim is to minimise the maximum scale distortion which will affect the mapped area and various formulas have been developed by different mathematicians to select the appropriate standard parallels to achieve this. Kavraisky was one mathematician who derived a recipe for choosing the standard parallels to achieve minimal scale distortion. But however the selection of the standard parallels is made the same projection formulas apply. Thus the parameters needed to specify a projection in the Lambert projected coordinate system group will be as follows:

For a Lambert Conical Conformal with one standard parallel (1SP),

Latitude of the natural origin (the Standard Parallel)

Longitude of the natural origin (the Central Meridian)

Scale factor at the natural origin (on the Standard Parallel)

False Easting

False Northing

For a Lambert Conical Conformal with two standard parallels (2SP),

Latitude of the first standard parallel

Latitude of the second standard parallel

Latitude of the false origin

Longitude of the false origin (the Central Meridian)

Easting at the false origin

Northing at the false origin

where the order of the standard parallels is not material if using the formulas which follow.

The limiting case of the Lambert Conic Conformal having the apex of the cone at infinity produces a cylindrical projection, the Mercator. Here, for the single standard parallel case the latitude of natural origin is the equator. For the two standard parallel case the two parallels have equal latitude in the north and south hemispheres. In both one and two standard parallel cases, grid coordinates are for the natural origin at the intersection of the equator and the central meridian (see figure 1). Thus the parameters needed to specify a projection in the Mercator projected coordinate system group will be:

For a Mercator with one standard parallel (1SP),

Latitude of the natural origin (the Equator)

Longitude of the natural origin (the Central Meridian)

Scale factor at the natural origin (on the Equator)

False Easting

False Northing

For a Mercator with two standard parallels (2SP),

Latitude of the first standard parallel

Longitude of the natural origin (the Central Meridian)

False Easting (grid coordinate at the intersection of the CM with the equator)

False Northing

In the formulas that follow the absolute value of the first standard parallel must be used.

For Azimuthal map projections, which are only infrequently used for ellipsoidal topographic mapping purposes, the natural origin will be at the centre of the projection where the map plane is imagined to be tangential to the ellipsoid and which will lie at the centre of the area to be projected. The central meridian will pass through the natural origin. This point will be given a False Easting and False Northing.

The parameters needed to specify the Stereographic map projection are:

Latitude of natural origin

Longitude of natural origin (the central meridian for the oblique case)

Scale Factor at the natural origin

False Easting at the natural origin

False Northing at the natural origin

1.4. Map Projection formulas

Only formulas for computation on the ellipsoid are considered. Projection formulas for the spherical earth are simpler but the spherical figure is inadequate to represent positional data with great accuracy at large map scales for the real earth. Projections of the sphere are only suitable for illustrative maps at scale of 1:1 million or less where precise positional definition is not critical.

The formulas which follow are largely adapted from "Map Projections - A Working Manual" by J.P.Snyder, published by the U.S. Geological Survey as Professional Paper No.1395. As well as providing an extensive overview of most map projections in current general use, and the formulas for their construction for both the spherical and ellipsoidal earth, this excellent publication provides computational hints and details of the accuracies attainable by the formulas. It is strongly recommended that all those who have to deal with map projections for medium and large scale mapping should follow its guidance.

There are a number of different formulas available in the literature for map projections other than those quoted by Snyder. Some are closed formulas; others, for ease of calculation, may depend on series expansions and their precision will generally depend on the number of terms used for computation. Generally those formulas which follow in this chapter will provide results which are accurate to within a decimetre, which is normally adequate for exploration mapping purposes. Coordinate expression and computations for engineering operations are usually consistently performed in grid terms.

The importance of one further variable should be noted. This is the unit of linear measurement used in the definition of projected coordinate systems. For metric map projections the unit of measurement is restricted to this unit. For non-metric map projections the metric ellipsoid semi-major axis needs to be converted to the projected coordinate system linear unit before use in the formulas below. The relevant ellipsoid is obtained through the datum part of the projected coordinate reference system.

In the formulas for map projections which follow, the basic ellipsoidal parameters are represented by symbols and derived as follows:

a is the ellipsoidal semi-major axis

b is the ellipsoidal semi-minor axis

f is the flattening of the ellipsoid where 1/f = a/(a - b)

e is the eccentricity of the ellipsoid where e2 = 2f - f2

e' is the second eccentricity where e'2 = e2 /(1 -e2)

r is the radius of curvature of the meridian at latitude j ,

where r = a(1 - e2)/(1 - e2sin2j )3/2

n is the radius of curvature on the prime vertical (i.e. perpendicular to the meridian) at

latitude j , where n = a /(1 - e2sin2j )1/2

j is the latitude of the point to be converted, positive if North and negative if South of the equator

l is the longitude of the point to be converted, positive if East and negative if West of the prime meridian

j 0 is the latitude of the natural origin

l 0 is the longitude of the natural origin (with respect to the prime meridian)

j F is the latitude of the false origin

l F is the longitude of the false origin (with respect to the prime meridian)

j 1 is the latitude of first standard parallel

j 2 is the latitude of second standard parallel

k0 is the scale factor at the natural origin

E is the Easting measured from the grid origin

N is the Northing measured from the grid origin

FE is the false easting, the Eastings value assigned to the natural origin

FN is the false northing, the Northings value assigned to the natural origin

EF is the Eastings value assigned to the false origin

NF is the Northings value assigned to the false origin

EC is the Eastings value assigned to the projection centre

NC is the Northings value assigned to the projection centre

Different formulas are require for forward and reverse map projection conversions: the forward formula cannot be used for the reverse conversion. However both forward and reverse formulas are explicitly given in the sections below as parts of a single conversion method. As such, map projection methods are described by EPSG as being reversible. Forward and reverse formulas for each conversion method use the projection parameters appropriate to that method with parameter values unchanged.

1.4.1. Lambert Conic Conformal

For territories with limited latitudinal extent but wide longitudinal width it may sometimes be preferred to use a single projection rather than several bands or zones of a Transverse Mercator. The Lambert Conic Conformal may often be adopted in these circumstances. But if the latitudinal extent is also large there may still be a need to use two or more zones if the scale distortion at the extremities of the one zone becomes too large to be tolerable.

Conical projections with one standard parallel are normally considered to maintain the nominal map scale along the parallel of latitude which is the line of contact between the imagined cone and the ellipsoid. For a one standard parallel Lambert the natural origin of the projected coordinate system is the intersection of the standard parallel with the longitude of origin (central meridian). See Figure 5 at end of section 1.3. To maintain the conformal property the spacing of the parallels is variable and increases with increasing distance from the standard parallel, while the meridians are all straight lines radiating from a point on the prolongation of the ellipsoid's minor axis.

Sometimes however, although a one standard parallel Lambert is normally considered to have unity scale factor on the standard parallel, a scale factor of slightly less than unity is introduced on this parallel. This is a regular feature of the mapping of some former French territories and has the effect of making the scale factor unity on two other parallels either side of the standard parallel. The projection thus, strictly speaking, becomes a Lambert Conic Conformal projection with two standard parallels. From the one standard parallel and its scale factor it is possible to derive the equivalent two standard parallels and then treat the projection as a two standard parallel Lambert conical conformal, but this procedure is seldom adopted. Since the two parallels obtained in this way will generally not have integer values of degrees or degrees minutes and seconds it is instead usually preferred to select two specific parallels on which the scale factor is to be unity, as for several State Plane Coordinate systems in the United States.

The choice of the two standard parallels will usually be made according to the latitudinal extent of the area which it is wished to map, the parallels usually being chosen so that they each lie a proportion inboard of the north and south margins of the mapped area. Various schemes and formulas have been developed to make this selection such that the maximum scale distortion within the mapped area is minimised, e.g. Kavraisky in 1934, but whatever two standard parallels are adopted the formulas are the same.

1.4.1.1 Lambert Conic Conformal (2SP) (EPSG coordinate operation method code 9802)

To derive the projected Easting and Northing coordinates of a point with geographical coordinates (j ,l ) the formulas for the Lambert Conic Conformal two standard parallel case (EPSG coordinate operation method code 9802) are:

Easting, E = EF + r sin q

Northing, N = NF + rF – r cos q

where m = cosj /(1 – e2sin2j )1/2 for m1, j 1, and m2, j 2 where j 1 and j 2 are the latitudes of the standard parallels

t = tan(p /4 – j /2)/[(1 – e sinj )/(1 + e sinj )]e/2 for t1, t2, tF and t using j 1,j 2,j F and j respectively

n = (ln m1 – ln m2)/(ln t1 – ln t2)

F = m1/(nt1n)

r = a F tn for rF and r, where rF is the radius of the parallel of latitude of the false origin

q = n(l – l 0)

The reverse formulas to derive the latitude and longitude of a point from its Easting and Northing values are:

j = p /2 – 2arctan{t'[(1 – esinj )/(1 + esinj )]e/2}

l = q '/n +l 0

where

r' = ± {(E – EF)2 + [rF – (N – NF)]2}1/2 , taking the sign of n

t' = (r'/(aF))1/n

q ' = arctan [(E – EF)/(rF – (N – NF))]

and n, F, and rF are derived as for the forward calculation.

Note that the formula for j requires iteration. First calculate t' and then a trial value for j using

j = π/2-2arctan t'. Then use the full equation for j substituting the trial value into the right hand side of the equation. Thus derive a new value for j . Iterate the process until j does not change significantly. The solution should quickly converge, in 3 or 4 iterations.

Example:

For Projected Coordinate Reference System NAD27 / Texas South Central

Parameters:

Ellipsoid Clarke 1866 a = 6378206.400 metres = 20925832.16 US survey feet 1/f = 294.97870

then e = 0.08227185, e2 = 0.00676866

First Standard Parallel j 1 28o23'00"N = 0.49538262 rad

Second Standard Parallel j 2 30o17'00"N = 0.52854388 rad

Latitude False Origin j F 27o50'00"N = 0.48578331 rad

Longitude False Origin l F 99o00'00"W = -1.72787596 rad

Easting at false origin EF 2000000.00 US survey feet

Northing at false origin NF 0.00 US survey feet

Forward calculation for:

Latitude j 28o30'00.00"N = 0.49741884 rad

Longitude l 96o00'00.00"W = -1.67551608 rad

first gives :

m1 = 0.88046050 m2 = 0.86428642

t = 0.59686306 tF = 0.60475101

t1 = 0.59823957 t2 = 0.57602212

n = 0.48991263 F = 2.31154807

r = 37565039.86 rF = 37807441.20

q = 0.02565177

Then Easting E = 2963503.91 US survey feet

Northing N = 254759.80 US survey feet

Reverse calculation for same easting and northing first gives:

q ' = 0.025651765 r' = 37565039.86

t' = 0.59686306

Then Latitude j = 28o30'00.000"N

Longitude l = 96o00'00.000"W

1.4.1.2 Lambert Conic Conformal (1SP) (EPSG coordinate operation method code 9801)

The formulas for the two standard parallel can be used for the Lambert Conic Conformal single standard parallel case (EPSG coordinate operation method code 9801) with minor modifications. Then

E = FE + r sinq

N = FN + r0 - r cosq , using the natural origin rather than the false origin.

where

n = sin j 0

r = aFtn k0 for r0, and r

t is found for t0, j 0 and t, j and m, F, and q are found as for the two standard parallel case.

The reverse formulas for j and l are as for the two standard parallel case above, with n, F and r0 as before and

q ' = arctan{(E – FE)/[r0 – (N – FN)]}

r' = ± {(E – FE)2 + [r0 – (N – FN)]2}1/2

t' = (r'/(ak0F))1/n

Example:

For Projected Coordinate Reference System JAD69 / Jamaica National Grid

Parameters:

Ellipsoid Clarke 1866 a = 6378206.400 m 1/f = 294.97870

then e = 0.08227185, e2 = 0.00676866

Latitude Natural Origin j o 18o 00'00" N = 0.31415927 rad

Longitude Natural Origin l o 77o 00'00" W = -1.34390352 rad

Scale factor at origin ko 1.000000

False Eastings FE 250000.00 m

False Northings FN 150000.00 m

Forward calculation for:

Latitude j = 17o 55' 55.80" N = 0.31297535 rad

Longitude l = 76o 56' 37.26" W = -1.34292061 rad

first gives

m0 = 0.95136402 t0 = 0.72806411

F = 3.39591092 n = 0.30901699

r = 19643955.26 r0 = 19636447.86

q = 0.00030374 t = 0.728965259

Then Easting E = 255966.58 m

Northing N = 142493.51 m

Reverse calculation for the same easting and northing first gives

q ' = 0.000303736

t' = 0.728965259

m0 = 0.95136402

r' = 19643955.26

Then Latitude j = 17o 55' 55.800" N

Longitude l = 76o 56' 37.260" W

1.4.1.3 Lambert Conic Conformal (West Orientated) (EPSG coordinate operation method code 9826)

In older mapping of Denmark and Greenland the Lambert Conic Conformal is used with axes positive north and west. To derive the projected Westing and Northing coordinates of a point with geographical coordinates (j , l ) the formulas are as for the standard Lambert Conic Conformal (1SP) case above (EPSG coordinate operation method code 9801) except for:

W = FE – r.sin(q )

In this formula the term FE retains its definition, i.e. in the Lambert Conic Conformal (West Orientated) method it increases the Westing value at the natural origin. In this method it is effectively false westing (FW).

The reverse formulas to derive the latitude and longitude of a point from its Westing and Northing values are as for the standard Lambert Conic Conformal (1SP) case except for:

q ' = arctan[(FE – W)/{r0 – (N – FN)}]

r' = +/-[(FE – W)² + {r0 – (N – FN)}²]^0.5

1.4.1.4 Lambert Conic Conformal (2 SP Belgium) (EPSG coordinate operation method code 9803)

In 1972, in order to retain approximately the same grid coordinates after a change of geodetic datum, a modified form of the two standard parallel case was introduced in Belgium. The formulas for the regular Lambert Conic Conformal (2SP) case given above are used except for:

Easting, E = EF + r sin (q – a)

Northing, N = NF + rF - r cos (q – a)

and for the reverse formulas

l = [(q ' + a)/n] + l 0

where a = 29.2985 seconds.

In 2000 this modification was replaced through use of the regular Lambert Conic Conformal (2 SP) map projection with appropriately modified parameter values.

Example:

For Projected Coordinate Reference System Belge l972 / Belge Lambert 72

Parameters:

Ellipsoid International 1924 a = 6378388 metres

1/f = 297

then e = 0.08199189, e2 = 0.006722670

First Standard Parallel j 1 49o50'00"N = 0.86975574 rad

Second Standard Parallel j 2 51o10'00"N = 0.89302680 rad

Latitude False Origin j F 90o00'00"N = 1.57079633 rad

Longitude False Origin l F 4o21'24.983"E = 0.07604294 rad

Easting at false origin EF 150000.01 metres

Northing at false origin NF 5400088.44 metres

Forward calculation for:

Latitude j = 50o40'46.4610"N = 0.88452540 rad

Longitude l = 5o48'26.533"E = 0.10135773 rad

first gives :

m1 = 0.64628304 m2 = 0.62834001

t = 0.59686306 tF = 0.00000000

t1 = 0.36750382 t2 = 0.35433583

n = 0.77164219 F = 1.81329763

r = 37565039.86 rF = 0.00

a= 0.00014204 q = 0.01953396

Then Easting E = 251763.20 metres

Northing N = 153034.13 metres

Reverse calculation for same easting and northing first gives:

q ' = 0.01939192 r' = 548041.03

t' = 0.35913403

Then Latitude j = 50o40'46.4610"N

Longitude l = 5o48'26.533"E

1.4.2. Lambert Conic Near-Conformal (EPSG coordinate operation method code 9817)

The Lambert Conformal Conic with one standard parallel formulas, as published by the Army Map Service, are still in use in several countries. The AMS uses series expansion formulas for ease of computation, as was normal before the electronic computer made such approximate methods unnecessary. Where the expansion series have been carried to enough terms the results are the same as the above formulas to the centimetre level. However in some countries the expansion formulas were truncated to the third order and the map projection is not fully conformal. The full formulas are used in Libya but from 1915 for France, Morocco, Algeria, Tunisia and Syria the truncated formulas were used. In 1943 in Algeria and Tunisia, from 1948 in France, from 1953 in Morocco and from 1973 in Syria the truncated formulas were replaced with the full formulas.

To compute the Lambert Conic Near-Conformal the following formulas are used:

E = FE + r sinq

N = FN + M + r sinq tan (q / 2) using the natural origin rather than the false origin.

Compute constants for the ellipse:

n = (a – b)/(a+b)

A' = a [ 1– n + 5 (n² – n³ ) / 4 + 81 ( n⁴ – n⁵ ) / 64]*p /180

B' = 3 a [ n – n² + 7 ( n³ – n⁴ ) / 8 + 55 n⁵ / 64] / 2

C' = 15 a [ n² – n³ + 3 ( n⁴ – n⁵ ) / 4 ] / 16

D' = 35 a [ n³ – n⁴ + 11 n⁵ / 16 ] / 48

E' = 315 a [ n⁴ – n⁵ ] / 512

Then compute the meridional arc from the equator to the parallel.

so = A' j o – B' sin 2j o + C' sin 4j o – D' sin 6j o + E' sin 8j o, where j o in the first term is in degrees

s = A' j – B' sin 2j + C' sin 4j – D' sin 6j + E' sin 8j , where j in the first term is in degrees

m = s – so

A = 1 / (6 r o n o)

M = ko ( m + Am³ + Bm⁴+ Cm⁵+ Dm⁶ ) This is the term that is truncated to the third order.

M_S = M per second of arc = M / ((j – j o) * 3600)

q = (l – l 0) sin j o

ro = ko n o / tan j o

r = ro – M

The reverse formulas for j and l from E and N with ro and M_S as above:

j = M’/ (M_S *3600) + j o where j o and j are in degrees

l = l 0 + q ‘ / sin j o where l o and l are in radians

where

X = E – FE Y = N – FN

q ' = arctan [X / (ro – Y)]

r' = X / sin q '

M' = r' – ro

Example:

For Projected Coordinate Reference System: Deir ez Zor / Levant Zone

Parameters:

Ellipsoid Clarke 1880 (IGN) a = 6378249.2 m 1/f = 293.46602

then b = 6356515.000 n = 0.001706682563

Latitude Natural Origin j o 34o 39'00" N = 0.604756586 rad

Longitude Natural Origin l o 37o 21'00" E = 0.651880476 rad

Scale factor at origin ko 0.99962560

False Eastings FE 300000.00 m

False Northings FN 300000.00 m

Forward calculation for:

Latitude j 37o 31' 17.625" N = 0.654874806 rad

Longitude l 34o 08' 11.291" E = 0.595793792 rad

first gives

A = 4.1067494 * 10^-15 A' =111131.8633

B' = 16300.64407 C' = 17.38751 D' = 0.02308 E' = 0.000033

so = 3835482.233 s = 4154101.458 m = 318619.225

M = 318632.72 M_S = 30.82262319

q = -0.03188875 ro = 9235264.405 r = 8916631.685

Then Easting E = 15707.96 m (c.f. E = 15708.00 using full formulas)

Northing N = 623165.96 m (c.f. N = 623167.20 using full formulas)

Reverse calculation for the same easting and northing first gives

q ' = -0.03188875

r' = 8916631.685

M' = 318632.72

Latitude j = 0.654874806 rad = 37o 31' 17.625" N

Longitude l = 0.595793792 rad = 34o 08' 11.291" E

1.4.3. Krovak Oblique Conformal Conic (EPSG coordinate operation method code 9819)

The normal case of the Lambert Conformal conic is for the axis of the cone to be coincident with the minor axis of the ellipsoid, that is the axis of the cone is normal to the ellipsoid at a geographic pole. For the Oblique Conformal Conic the axis of the cone is normal to the ellipsoid at a defined location and its extension cuts the minor axis at a defined angle. This map projection is used in the Czech Republic and Slovakia under the name ‘Krovak’ projection. The map projection method is similar in principle to the Oblique Mercator (see section 1.4.5). The geographic coordinates on the ellipsoid are first reduced to conformal coordinates on the conformal (Gaussian) sphere. These spherical coordinates are then projected onto the oblique cone and converted to grid coordinates. The pseudo standard parallel is defined on the conformal sphere after its rotation, to obtain the oblique aspect of the projection. It is then the parallel on this sphere at which the map projection is true to scale; on the ellipsoid it maps as a complex curve. A scale factor may be applied to the map projection to increase the useful area of coverage.

The defining parameters for the Krovak oblique conformal conic map projection are:

j _c = latitude of centre of the projection

l _c = longitude of centre of the projection

j ₁ = latitude of pseudo standard parallel

k_c = scale factor on the pseudo standard parallel

E_c = False Easting of the centre of the projection at the apex of the cone

N_c = False Northing of the centre of the projection at the apex of the cone

From these the following constants for the projection may be calculated :

B = ((1 + e² cos⁴(j _c)) / (1 – e² ))^0.5

A = a (1 – e² )^0.5 / ( 1 – e² sin² (j _c))

g ₀ = asin(sin (j _c) / B)

t₀ = tan(p / 4 + g ₀/ 2) . ((1 + e sin (j _c)) / (1 – e sin (j _c)))^e.B/2 / tan(p / 4 + j _c/ 2) ^B

n = sin j ₁

r0 = k_c A / tan (j ₁)

To derive the projected Easting and Northing coordinates of a point with geographical coordinates (j ,l ) the formulas for the oblique conic conformal are:

Easting: E = Ec + r cos q

Northing: N = Nc + r sinq

where

U = 2 (atan ( k tan^B(j / 2 + p / 4 ) ((1 + e sin (j )) / (1 – e sin (j )))^e.B/2 ) – p / 4)

V = B (l _c – _l )

S = asin ( cos (a _c) sin ( U ) + sin (a _c) cos (U) cos (V) )

D = asin ( cos ( U ) sin ( V ) / cos ( S ) )

q = n D

r = r0 tan(p / 4 + g ₀/ 2) / tan ⁿ ( S/2 + p / 4 )

Note that the terms Easting and Northing here refer to the two map grid ordinates. Their actual geographic direction depends upon the azimuth of the centre line.

The reverse formulas to derive the latitude and longitude of a point from its Easting and Northing values are:

j _j = 2*(atan(k^{-1/ B} tan ^{1/ B} ( U’/2 + p / 4 ) ((1 + e sin ( j _j-1) / ( 1 – e sin ( j _j-1)) ^e/2 ) – p / 4)

where j = 1,2 and the latitude is found by iteration.

l = l c– V’ / B

where

r' = [(E – EF)2 + (N – NF)2]1/2

q ' = arctan [(E – Ec)/(N – Nc)]

D' = q ' / sin ( j ₁)

S' = 2*(atan((r0 / r )^1/n tan(p / 4 + j ₀/ 2)) - p / 4)

U' = asin ( cos (a _c) sin (S') – sin (a _c) cos (S') cos (D') )

V' = asin ( cos (S') sin (D') / cos (U'))

Example:

For Projected Coordinate Reference System: S-JTSK (Ferro) / Krovak

N.B. Krovak projection uses Ferro as the prime meridian. This has a longitude with reference to Greenwich of 17 deg 40 min West. To apply the formulas the defining longitudes must be corrected to the Greenwich meridian.

Parameters:

Ellipsoid Bessel 1841 a = 6377397.155m 1/f = 299.15281

then e = 0.081696831 e2 = 0.006674372

Latitude of projection centre (j _c) 49o 30'00" N = 0.604756586 rad

Longitude of Origin 42o 30'00" East of Ferro

Longitude of Ferro is 17o 40'00" West of Greenwich

Longitude of Origin (l _c) 24o 50'00" East of Greenwich = 0.433423431 rad

Latitude of pseudo standard parallel (j ₁) 78o 30'00" N

Azimuth of centre line (a _c) 30o 17' 17.303"

Scale factor on pseudo standard parallel (ko) 0.99990

Easting at projection centre (E_c) 0.00 m

Northing at projection centre (N_c) 0.00 m

Projection constants:

B = 1.000597498

A = 6380703.61

g ₀ = 0.863239103

t₀ = 1.003419164

n = 0.979924705

r0 = 1298039.005

Forward calculation for:

Latitude j 50o 12' 32.4416" N = 0.876312566 rad

Longitude l 16o 50' 59.1790" E = 0.294083999 rad

Gives

U = 0.875596949

V = 0.139422687

S = 1.386275049

D = 0.506554623

q = 0.496385389

r = 1194731.014

Then "Easting" E = 1050538.643 m

"Northing" N = 568990.997 m

where "Easting" increases southwards and "Northing" increases westwards.

Reverse calculation for the same "Easting" and "Northing" gives

r' = 1194731.014

q ' = 0.496385389

D' = 0.506554623

S' = 1.386275049

U' = 0.875596949

V' = 0.139422687

j 1 = 0.876310601

j 2 = 0.876312560

j 3 = 0.876312566

Latitude j = 0.876312566 rad = 50o 12' 32.4416" N

Longitude l = 0.595793792 rad = 16o 50' 59.1790" E

1.4.4. Mercator (EPSG coordinate operation method codes 9804 and 9805)

The Mercator map projection is a special limiting case of the Lambert Conic Conformal map projection with the equator as the single standard parallel. All other parallels of latitude are straight lines and the meridians are also straight lines at right angles to the equator, equally spaced. It is the basis for the transverse and oblique forms of the projection. It is little used for land mapping purposes but is in almost universal use for navigation charts. As well as being conformal, it has the particular property that straight lines drawn on it are lines of constant bearing. Thus navigators may derive their course from the angle the straight course line makes with the meridians.

In the few cases in which the Mercator projection is used for terrestrial applications or land mapping, such as in Indonesia prior to the introduction of the Universal Transverse Mercator, a scale factor may be applied to the projection. This has the same effect as choosing two standard parallels on which the true scale is maintained at equal north and south latitudes either side of the equator.

The formulas to derive projected Easting and Northing coordinates are:

For the two standard parallel case, k0 is first calculated from

k0 = cosj 1/(1 – e2sin2j 1)1/2

where j 1 is the absolute value of the first standard parallel (i.e. positive).

Then, for both one and two standard parallel cases,

E = FE + ak0(l – l 0)

N = FN + ak0 ln{tan(p /4 + j /2)[(1 – esinj )/(1 + esinj )]e/2 }

where symbols are as listed above and logarithms are natural.

The reverse formulas to derive latitude and longitude from E and N values are:

j = c + (e2/2 + 5e4/24 + e6/12 + 13e8/360) sin(2c )

+ (7e4/48 + 29e6/240 + 811e8/11520) sin(4c )

+ (7e6/120 + 81e8/1120) sin(6c ) + (4279e8/161280) sin(8c )

where c = p /2 – 2 arctan t

t = B (FN-N)/(ak0⁾

B = base of the natural logarithm, 2.7182818...

and for the 2 SP Case, k0 is calculated as for the forward transformation above.

l = [(E – FE)/ak0] + l 0

Examples

1. Mercator (with two standard parallels) (EPSG coordinate operation method code 9805)

For Projected Coordinate Reference System Pulkovo 1942 / Mercator Caspian Sea

Parameters:

Ellipsoid Krassowski 1940 a = 6378245.00m 1/f = 298.300

then e = 0.08181333, e2 = 0.00669342

Latitude first SP j 0 42o00'00"N = 0.73303829 rad

Longitude Natural Origin l 0 51o00'00"E = 0.89011792 rad

False Eastings FE 0.00 m

False Northings (at equator) FN 0.00 m

then natural origin at latitude of 0oN has scale factor k0= 0.74426089

Forward calculation for:

Latitude j 53o00'00.00"N = 0.9250245 rad

Longitude l 53o00'00.00"E = 0.9250245 rad

gives Easting E = 165704.29 m

Northing N = 5171848.07 m

Reverse calculation for same easting and northing first gives :

t = 0.33639129

c = 0.92179596

Then Latitude j = 53o00'00.000"N

Longitude l = 53o00'00.000"E

2. Mercator (with the equator as standard parallel) (EPSG coordinate operation method code 9804)

For Projected Coordinate Reference System Makassar / NEIEZ

Parameters:

Ellipsoid Bessel 1841 a = 6377397.155 m 1/f = 299.15281

then e = 0.08169683

Latitude Natural Origin j o 00o00'00"N = 0.0000000 rad

Longitude Natural Origin l o 110o00'00"E = 1.91986218 rad

Scale factor ko 0.997

False Eastings FE 3900000.00 m

False Northings FE 900000.00 m

Forward calculation for:

Latitude j = 3o00'00.00"S = -0.05235988 rad

Longitude l = 120o00'00.00"E = 2.09439510 rad

gives

Easting E = 5009726.58 m

Northing N = 569150.82 m

Reverse calculation for same easting and northing first gives :

t = 1.0534121

c = -0.0520110

Then Latitude j = 3o00'00.000"S

Longitude l = 120o00'00.000"E

1.4.5. Cassini-Soldner (EPSG coordinate operation method code 9806)

The Cassini-Soldner projection is the ellipsoidal version of the Cassini projection for the sphere. It is not conformal but as it is relatively simple to construct it was extensively used in the last century and is still useful for mapping areas with limited longitudinal extent. It has now largely been replaced by the conformal Transverse Mercator which it resembles. Like this, it has a straight central meridian along which the scale is true, all other meridians and parallels are curved, and the scale distortion increases rapidly with increasing distance from the central meridian.

The formulas to derive projected Easting and Northing coordinates are:

Easting, E = FE + n [A – TA3/6 – (8 – T + 8C)TA5/120]

Northing, N = FN + M – M0 + n tanj [A2/2 + (5 – T + 6C)A4/24]

where A = (l – l 0)cosj

T = tan2j

C = e2 cos2j /(1 - e2)

n = a /(1 – e2sin2j )1/2

and M, the distance along the meridian from equator to latitude j , is given by

M = a[1 – e2/4 – 3e4/64 – 5e6/256 –....)j – (3e2/8 + 3e4/32 + 45e6/1024 +....)sin2j

+ (15e4/256 + 45e6/1024 +.....)sin4j – (35e6/3072 + ....)sin6j + .....]

with j in radians.

M0 is the value of M calculated for the latitude of the chosen origin. This may not necessarily be chosen as the equator.

To compute latitude and longitude from Easting and Northing the reverse formulas are:

j = j 1 – (n 1tanj 1/r 1)[D2/2 – (1 + 3T1)D4/24]

l = l 0 + [D – T1D3/3 + (1 + 3T1)T1D5/15]/cosj 1

where

n 1 = a /(1 – e2sin2j 1)1/2

r 1 = a(1 – e2)/(1 – e2sin2j 1)3/2

j 1 is the latitude of the point on the central meridian which has the same Northing as the point whose coordinates are sought, and is found from:

j 1 = m 1 + (3e1/2 – 27e13/32 +.....)sin2m 1 + (21e12/16 – 55e14/32 + ....)sin4m 1

+ (151e13/96 +.....)sin6m 1 + (1097e14/512 – ....)sin8m 1 + ......

where

e1 = [1 – (1 – e2)1/2]/[1 + (1 – e2)1/2]

m 1 = M1/[a(1 – e2/4 – 3e4/64 – 5e6/256 – ....)]

M1 = M0 + (N – FN)

T1 = tan2j 1

D = (E – FE)/n 1

For Projected Coordinate Reference System Trinidad 1903 / Trinidad Grid

Parameters:

Ellipsoid Clarke 1858 a = 20926348 ft = 31706587.88 links

b = 20855233 ft

then 1/f = 294.97870, e2 = 0.00676866

Latitude of natural origin j 0 10o26'30"N = 0.182241463 rad

Longitude of natural origin l 0 61o20'00"W = -1.07046861 rad

False Eastings FE 430000.00 links

False Northings FN 325000.00 links

Forward calculation for:

Latitude j 10o00'00.00" N = 0.17453293 rad

Longitude l 62o00'00.00"W = – 1.08210414 rad

A = -0.01145876 C = 0.00662550

T = 0.03109120 M = 5496860.24

n = 31709831.92 M0 = 5739691.12

Then Easting E = 66644.94 links

Northing N = 82536.22 links

Reverse calculation for same easting and northing first gives :

e1 = 0.00170207 D = -0.01145875

T1 = 0.03109544 M1 = 5497227.34

n 1 = 31709832.34 m 1 = 0.17367306

j 1 = 0.17454458 r 1 = 31501122.40

Then Latitude j = 10o00'00.000"N

Longitude l = 62o00'00.000"W

1.4.6. Transverse Mercator

1.4.6.1 General Case (EPSG coordinate operation method code 9807)

The Transverse Mercator projection in its various forms is the most widely used projected coordinate system for world topographical and offshore mapping. All versions have the same basic characteristics and formulas. The differences which distinguish the different forms of the projection which are applied in different countries arise from variations in the choice of values for the the coordinate conversion parameters, namely the latitude of the natural origin, the longitude of the natural origin (central meridian), the scale factor at the natural origin (on the central meridian), and the values of False Easting and False Northing, which embody the units of measurement, given to the origin. Additionally there are variations in the width of the longitudinal zones for the projections used in different territories.

The following table indicates the variations in the coordinate conversion parameters which distinguish the different forms of the Transverse Mercator projection and are used in the EPSG Transverse Mercator map projection operations:

The most familiar and commonly used Transverse Mercator in the oil industry is the Universal Transverse Mercator (UTM) whose natural origin lies on the equator. However, some territories use a Transverse Mercator with a natural origin at a latitude of natural origin closer to that territory.

In EPSG the coordinate conversion method is considered to be the same for all forms of the Transverse Mercator projection. The formulas to derive the projected Easting and Northing coordinates for the normal case (EPSG coordinate operation method code 9807) are in the form of a series as follows:

Northing, N = FN + k0{M – M0 + n tanj [A2/2 + (5 – T + 9C + 4C2)A4/24 +

(61 – 58T + T2 + 600C – 330e'2)A6/720]}

where T = tan2j

C = e2 cos2j /(1 – e2)

A = (l – l 0)cosj , with l and l 0 in radians

n = a /(1 – e2sin2j )1/2

M = a[(1 – e2/4 – 3e4/64 – 5e6/256 –....)j – (3e2/8 + 3e4/32 + 45e6/1024+....)sin2j

+ (15e4/256 + 45e6/1024 +.....)sin4j – (35e6/3072 + ....)sin6j + .....]

with j in radians and M0 for j 0, the latitude of the origin, derived in the same way.

The reverse formulas to convert Easting and Northing projected coordinates to latitude and longitude are:

j = j 1 – (n 1tanj 1/r 1)[D2/2 – (5 + 3T1 + 10C1 – 4C12 – 9e'2)D4/24

+ (61 + 90T1 + 298C1 + 45T12 – 252e'2 – 3C12)D6/720]

l = l 0 + [D – (1 + 2T1 + C1)D3/6 + (5 – 2C1 + 28T1 – 3C12 + 8e'2 + 24T12)D5/120] / cosj 1

where

n 1 = a /(1 – e2sin2j 1)1/2

r 1 = a(1 – e2)/(1 – e2sin2j 1)3/2

j 1 may be found as for the Cassini projection from:

j 1 = m 1 + (3e1/2 – 27e13/32 +.....)sin2m 1 + (21e12/16 – 55e14/32 + ....)sin4m 1

+ (151e13/96 +.....)sin6m 1 + (1097e14/512 – ....)sin8m 1 + ......

and where

e1 = [1 – (1 – e2)1/2]/[1 + (1 – e2)1/2]

m 1 = M1/[a(1 – e2/4 – 3e4/64 – 5e6/256 – ....)]

M1 = M0 + (N – FN)/k0

T1 = tan2j 1

C1 = e'2cos2j 1

e'2 = e2 /(1 – e2)

D = (E – FE)/(n 1k0)

For areas south of the equator the value of latitude j will be negative and the formulas above, to compute the E and N, will automatically result in the correct values. Note that the false northings of the origin, if the equator, will need to be large to avoid negative northings and for the UTM projection is in fact 10,000,000m. Alternatively, as in the case of Argentina's Transverse Mercator (Gauss-Kruger) zones, the origin is at the south pole with a northings of zero. However each zone central meridian takes a false easting of 500000m prefixed by an identifying zone number. This ensures that instead of points in different zones having the same eastings, every point in the country, irrespective of its projection zone, will have a unique set of projected system coordinates. Strict application of the above formulas, with south latitudes negative, will result in the derivation of the correct Eastings and Northings.

Similarly, in applying the reverse formulas to determine a latitude south of the equator, a negative sign for j results from a negative j 1 which in turn results from a negative M1.

For Projected Coordinate Reference System OSGB 1936 / British National Grid

Parameters:

Ellipsoid Airy 1830 a = 6377563.396 m 1/f = 299.32496

then e'2 = 0.00671534 e2 = 0.00667054

Latitude Natural Origin j o 49o00'00"N = 0.85521133 rad

Longitude Natural Origin l o 2o00'00"W = -0.03490659 rad

Scale factor ko 0.9996013

False Eastings FE 400000.00 m

False Northings FN – 100000.00 m

Forward calculation for:

Latitude j 50o30'00.00"N = 0.88139127 rad

Longitude l 00o30'00.00"E = 0.00872665 rad

A = 0.02775415 C = 0.00271699

T = 1.47160434 M = 5596050.46

M0 = 5429228.60 n = 6390266.03

Then Easting E = 577274.99 m

Northing N = 69740.50 m

Reverse calculations for same easting and northing first gives :

e1 = 0.00167322 m 1 = 0.87939562

M1 = 5599036.80 n 1 = 6390275.88

j 1 = 0.88185987 D = 0.02775243

r 1 = 6372980.21 C1 = 0.00271391

T1 = 1.47441726

Then Latitude j = 50o30'00.000"N

Longitude l = 00o30'00.000"E

1.4.6.2 Transverse Mercator Zoned Grid System (EPSG coordinate operation method code 9824)

When the growth in distortion away from the projection origin is of concern, a projected coordinate reference system cannot be used far from its origin. A means of creating a grid system over a large area but also limiting distortion is to have several grid zones with most defining parameters being made common. Coordinates throughout the system are repeated in each zone. To make coordinates unambiguous the easting is prefixed by the relevant zone number. This procedure was adopted by German mapping in the 1930’s through the Gauss-Kruger systems and later by American military mapping through the Universal Transverse Mercator (or UTM) grid system. (Note: subsequent civilian adoption of the systems usually ignores the zone prefix to easting. Where this is the case the formulas below do not apply: use the standard TM formula separately for each zone).

The parameter Longitude of Natural Origin (l 0) is changed from being a defining parameter to a derived parameter, replaced by two other defining parameters, the Initial Longitude (the western limit of zone 1) (l I) and the Zone Width (W). Each of the remaining four Transverse Mercator defining parameters – Latitude of Natural Origin, Scale factor at natural origin, False easting and False Northing – have the same parameter values in every zone.

The standard Transverse Mercator formulas above are modified as follows:

Zone number, Z, = int((Long + LongI + W) / W) with Long, LongI and W in degrees.

where (LongI) is the Initial Longitude of the zoned grid system

and W is the width of each zone of the zoned grid system.

If Long < 0, Long = (Long + 360) degrees.

Then,

Long0 = [Z * W] – [LongI + (W/2)]

For the forward calculation,

Easting, E = Z*10^6 + FE + k0.nu[A + (1 – T + C)A^3/6 + (5 – 18T + T^2 + 72C – 58e'^2)A^5/120]

and in the reverse calculation for longitude,

D = (E – [FE + Z*10^6])/(nu1.k0)

1.4.6.3 Transverse Mercator (South Orientated) (EPSG coordinate operation method code 9808)

For the mapping of southern Africa a south oriented Transverse Mercator map projection method is used. Here the coordinate axes are called Westings and Southings and increment to the West and South from the origin respectively. See Figure 3 for a diagrammatic illustration. The standard Transverse Mercator above formulas need to be modified to cope with this arrangement with

Westing, W = FE – k0 n [A + (1 – T + C)A3/6 + (5 - 18T + T2 + 72C - 58e'2)A5/120]

Southing, S = FN – k0{M – M0 + n tanj [A2/2 + (5 – T + 9C + 4C2)A4/24 +

(61 - 58T + T2 + 600C - 330e'2)A6/720]}

In these formulas the terms FE and FN retain their definition, i.e. in the Transverse Mercator (South Orientated) method they increase the Westing and Southing value at the natural origin. In this method they are effectively false westing (FW) and false southing (FS) respectively.

For the reverse formulas, those for the standard Transverse Mercator above apply, with the exception that:

M1 = M0 – (S – FN)/k0

and D = – (W – FE)/(n 1k0), with n 1 = n for j 1

1.4.7. Oblique Mercator and Hotine Oblique Mercator

The initial line central to the map area of given azimuth a _c passes through a defined centre of the projection (j _c , l _c ) . The point where the projection of this line cuts the equator on the aposphere is the origin of the (u , v) co-ordinate system. The u axis is parallel to the centre line and the v axis is perpendicular to (90° clockwise from) this line.

In applying the formulas for the (Hotine) Oblique Mercator the first set of co-ordinates computed are referred to the (u, v) co-ordinate axes defined with respect to the azimuth of the centre line. These co-ordinates are then ‘rectified’ to the usual Easting and Northing by applying an orthogonal conversion. Hence the alternative name as the Rectified Skew Orthomorphic. In the special case of the projection covering the Alaskan panhandle the azimuth of the line at the natural origin is taken to be identical to the azimuth of the initial line at the projection centre. This results in grid and true north coinciding at the projection centre rather than at the natural origin as is more usual.

To ensure that all co-ordinates in the map area have positive grid values, false co-ordinates are applied. These may be given values (E_c , N_c) if applied at the projection centre [EPSG Oblique Mercator method] or be applied as false easting (FE) and false northing (FN) at the natural origin [EPSG Hotine Oblique Mercator method].

The formulas can be used for the following cases:

Alaska State Plane Zone 1

Hungary EOV

Madagascar Laborde Grid

East and West Malaysia Rectified Skew Orthomorphic grids

Swiss Cylindrical projection

The Swiss and Hungarian systems are a special case where the azimuth of the line through the projection centre is 90 degrees. These therefore give similar but not exactly the same results as a conventional Transverse Mercator projection.

Specific references for the formulas originally used in the individual cases of these projections are:

The defining parameters for the [Hotine] Oblique Mercator projection are:

j _c = latitude of centre of the projection

l _c = longitude of centre of the projection

a _c = azimuth (true) of the centre line passing through the centre of the projection

g _c = rectified bearing of the centre line

k_c = scale factor at the centre of the projection

and either

for the Oblique Mercator

E_c = False Easting at the centre of projection

N_c = False Northing at the centre of projection

or for the Hotine Oblique Mercator:

FE = False Easting at the natural origin

FN = False Northing at the natural origin

From these defining parameters the following constants for the map projection may be calculated :

B = [1 + e² cos^4j _c / (1 – e² )]^0.5

A = a B k_c (1 – e² )^0.5 / (1 – e² sin² j _c)

t_o = tan(p / 4 – j _c / 2) / [(1 – e sin j _c) / (1 + e sin j _c)]^e/2

D = B (1 – e² )^0.5 / [cosj _c (1 – e² sin² j _c)^0.5]

To avoid problems with computation of F, if D < 1 make D² = 1

F = D + (D² – 1)^0.5 . SIGN(j _c)

H = F t_o^B

G = (F – 1/F) / 2

g _o = asin[sin (a _c) / D]

l _o = l _c – [asin(G tang _o)] / B

Then compute the (u_c , v_c) co-ordinates for the centre point (j _c , l _c).

v_c = 0

In general

u_c = (A / B) atan[(D² – 1)^0.5 / cos (a _c) ]. SIGN(j _c)

but for the special cases where a _c = 90 degrees (e.g. Hungary, Switzerland) then

u_c = A (l _c – l _o )

Forward case: To compute (E,N) from a given (j ,l ), for both the Hotine Oblique Mercator method and the Oblique Mercator method:

t = tan(p / 4 – j / 2) / [(1 – e sin j ) / (1 + e sin j )]^e/2

Q = H / t^B

S = (Q – 1/Q) / 2

T = (Q + 1/Q) / 2

V = sin(B (l – l _o))

U = (– V cos(g _o) + S sin(g _o)) / T

v = A ln[(1 – U) / (1 + U)] / 2 B

Then

either (a) for the Hotine Oblique Mercator (where the FE and FN values have been specified with respect to the natural origin of the (u , v) axes):

u = A atan{(S cosg _o + V sing _o) / cos[B (l – l _o)]} / B

The rectified skew co-ordinates are then derived from:

E = v cos(g _c) + u sin(g _c) + FE

N = u cos(g _c) – v sin(g _c) + FN

or (b) for the Oblique Mercator (where the false easting and northing values (E_c , N_c) have been specified with respect to the centre of the projection (j _c , l _c) then :

u = (A atan{(S cosg _o + V sing _o) / cos[B (l – l _o)]} / B) – (u_c . SIGN(l – l _c))

The rectified skew co-ordinates are then derived from:

E = v cos(g _c) + u sin(g _c) + E_c

N = u cos(g _c) – v sin(g _c) + N_c

Reverse case: To compute (j ,l ) from a given (E,N) :

For the Hotine Oblique Mercator:

v’ = (E – FE) cos(g _c) – (N – FN) sin(g _c)

u’ = (N – FN) cos(g _c) + (E - FE) sin(g _c)

or for the Oblique Mercator:

v’ = (E – E_c) cos(g _c) – (N – N_c) sin(g _c)

u’ = (N – N_c) cos(g _c) + (E – E_c) sin(g _c) + u_c

then for both cases:

Q’ = e^{- (B v ‘/ A)} where e is the base of natural logarithms.

S’ = (Q’ – 1 / Q’) / 2

T’ = (Q’ + 1 / Q’) / 2

V’ = sin (B u’ / A)

U’ = (V’ cos(g _o) + S’ sin(g _o)) / T’

t’ = {H / [(1 + U’) / (1 – U’)]^0.5}^{1 / B}

c = p / 2 – 2 atan(t’)

• = c + sin(2c ).( e² / 2 + 5 e⁴ / 24 + e⁶ / 12 + 13 e⁸ / 360)

l = l ₀ – atan [(S’ cosg _c – V’ sing _c) / cos(B u’ / A)] / B

For Projected Coordinate Reference System Timbalai 1948 / R.S.O. Borneo (m) using the Oblique Mercator method.

Parameters:

Ellipsoid: Everest 1830 (1967 Definition)

a = 6377298.556 metres 1/f = 300.8017

then e = 0.081472981 e2 = 0.006637847

Latitude of Projection Centre j c 4o00'00"N = 0.069813170 rad

Longitude of Projection Centre l c 115o00'00"E = 2.007128640 rad

Azimuth of central line ac 53o18'56.9537" = 0.930536611 rad

Rectified to skew gc 53o07'48.3685" = 0.927295218 rad

Scale factor ko 0.99984

Easting at Projection Centre Ec 590476.87 m

Northings at Projection Centre Nc 442857.65 m

Forward calculation for:

Latitude j 5o 23' 14.1129"N = 0.094025313 rad

Longitude l 115o 48' 19.8196"E = 2.021187362 rad

t_o = 0.932946976 g_o = 0.927295218

D = 1.002425787 l_o = 1.914373469

D² = 1.004857458

u_c = 738096.09 v_c = 0.00

S = 0.093990763 T = 1.004407419

V = 0.106961709 U = 0.010967247

v = -69702.787 u = 163238.163

Then Easting E = 679245.73 m

Northing N = 596562.78 m

Reverse calculations for same Easting and Northing first gives :

v’ = -69702.787 u’ = 901334.257

Q’ = 1.011028053

S’ = 0.010967907 T’ = 1.000060146

V’ = 0.141349378 U’ = 0.093578324

t’ = 0.910700729 c = 0.093404829

Then Latitude j = 5o23'14.113"N

Longitude l = 115o48'19.820"E

If the same projection is defined using the Hotine Oblique Mercator method then:

False Easting Ec 0.00 m

False Northing Nc 0.00 m

u = 901334.257

and all other values are as for the Oblique Mercator method.

1.4.8. Stereographic

The Stereographic projection may be imagined to be a projection of the earth's surface onto a plane in contact with the earth at a single tangent point from a projection point at the opposite end of the diameter through that tangent point.

This projection is best known in its polar form and is frequently used for mapping polar areas where it complements the Universal Transverse Mercator used for lower latitudes. Its spherical form has also been widely used by the US Geological Survey for planetary mapping and the mapping at small scale of continental hydrocarbon provinces. In its transverse or oblique ellipsoidal forms it is useful for mapping limited areas centred on the point where the plane of the projection is regarded as tangential to the ellipsoid., e.g. the Netherlands. The tangent point is the origin of the projected coordinate system and the meridian through it is regarded as the central meridian. In order to reduce the scale error at the extremities of the projection area it is usual to introduce a scale factor of less than unity at the origin such that a unitary scale factor applies on a near circle centred at the origin and some distance from it.

The coordinate conversion from geographical to projected coordinates is executed via the distance and azimuth of the point from the centre point or origin. For a sphere the formulas are relatively simple. For the ellipsoid the parameters defining the conformal sphere at the tangent point as origin are first derived. The conformal latitudes and longitudes are substituted for the geodetic latitudes and longitudes of the spherical formulas for the origin and the point.

1.4.8.1 Oblique and Equatorial Stereographic cases (EPSG coordinate operation method code 9809)

Given the geodetic origin of the projection at the tangent point (j 0, l 0), the parameters defining the conformal sphere are:

R= Ö ( r 0 . n 0)

n= Ö {1+[(e2 cos4j 0) / (1 – e2)]}

c= (n + sinj 0) (1 – sinc 0) / [(n – sinj 0) (1 + sin(c 0)]

where: sinc 0 = (w1 – 1)/(w1 + 1)

w1 = [S1.(S2)e]n

S1 = (1 + sinj 0)/(1 – sinj 0)

S2 = (1 – e sinj 0)/(1 + e sinj 0)

The conformal latitude and longitude of the origin (c 0,L 0) are then computed from :

c 0 = sin-1[(w2 – 1)/( w2 + 1)]

where S1 and S2 are as above and w2 = c [S1.(S2)e]n = c.w1

L 0 = l 0

For any point with geodetic coordinates (j ,l ) the equivalent conformal latitude and longitude

( c , L ) are then computed from

L = n( l – L 0 ) + L 0

and

c = sin-1[(w – 1)/(w + 1)]

where w = c [Sa.(Sb)e]n

Sa = (1 + sinj ) / (1 – sinj )

Sb = (1 – e.sinj ) / (1 + e.sinj )

Then

E = FE + 2 R k0 cosc sin(L – L 0 ) / B

and

N = FN + 2 R k0 [sinc cosc 0 – cosc sinc 0cos(L – L 0 )] / B

where B = [1+sinc sinc 0 + cosc cosc 0cos(L – L 0 )]

The reverse formulas to compute the geodetic coordinates from the grid coordinates involves computing the conformal values, then the isometric latitude and finally the geodetic values.

The parameters of the conformal sphere and conformal latitude and longitude at the origin are computed as above. Then for any point with Stereographic grid coordinates (E,N) :

c = c 0 + 2 tan-1{[(N – FN) – (E – FE) tan (j/2)] / (2 Rk0)}

L = j + 2 i + L 0

where g = 2 Rk0 tan (p /4 – c 0/ 2 )

h = 4 Rk0 tan c 0 + g

i = tan-1 {(E – FE) / [h + (N – FN)]}

j = tan-1 {(E – FE) / [g – (N – FN)]} - i

Geodetic longitude l = (L – L 0 ) / n + L 0

Isometric latitude y = 0.5 ln {(1 + sinc ) / [ c (1 – sinc )]} / n

First approximation j 1 = 2 tan-1 eY - p / 2 where e=base of natural logarithms.

y i = isometric latitude at j i

where y i= ln{[tan(j i/2+p / 4] [(1 – e sinj i)/(1 + e sinj i)]e/2}

Then iterate j i+1 = j i – ( y i – y ) cos j i ( 1 – e2 sin2j i) / (1 – e2)

until the change in j is sufficiently small.

If the projection is the equatorial case, j 0 and c 0 will be zero degrees and the formulas are simplified as a result, but the above formulas remain valid.

For the polar version, j 0 and c 0 will be 90 degrees and the formulas become indeterminate. See below for formulas for the polar case.

For stereographic projections centred on points in the southern hemisphere, the signs of E, N, l 0 and l must be reversed to be used in the equations and j will be negative anyway as a southerly latitude.

For Projected Coordinate Reference System RD / Netherlands New

Parameters:

Ellipsoid Bessel 1841 a = 6377397.155 m 1/f = 299.15281

then e = 0.08169683

Latitude Natural Origin j o 52o09'22.178"N = 0.910296727 rad

Longitude Natural Origin l o 5o23'15.500"E = 0.094032038 rad

Scale factor k0 0.9999079

False Eastings FE 155000.00 m

False Northings FN 463000.00 m

Forward calculation for:

Latitude j 530N = 0.925024504 rad

Longitude l 60E = 0.104719755 rad

first gives the conformal sphere constants:

r 0 = 6374588.71 n 0 = 6390710.613

R = 6382644.571 n = 1.000475857 c = 1.007576465

where S1 = 8.509582274 S2 = 0.878790173 w1 = 8.428769183

sin c 0 = 0.787883237

w2 = 8.492629457 c 0 = 0.909684757 L 0 = l 0 = 0.094032038 rad

For the point (j ,l ) c = 0.924394997 L = 0.104724841 rad

hence B = 1.999870665

and E = 196105.283 m N = 557057.739 m

Reverse calculation for the same Easting and Northing (196105.28 E, 557057.74 N) first gives:

g = 4379954.188 h = 37197327.960 i = 0.001102255 j = 0.008488122

then L = 0.10472467 whence l = 0.104719584 rad = 60E

Also c = 0.924394767 and y = 1.089495123

Then j 1 = 0.921804948 y 1 = 1.084170164

j 2 = 0.925031162 y 2 = 1.089506925

j 3 = 0.925024504 y 3 = 1.089495505

j 4 = 0.925024504

Then Latitude j = 53o00'00.000"N

Longitude l = 6o00'00.000"E

An alternative approach is given by Snyder, where, instead of defining a single conformal sphere at the origin point, the conformal latitude at each point on the ellipsoid is computed. The conformal longitude is then always equivalent to the geodetic longitude. This approach is a valid alternative to the above, but gives slightly different results away from the origin point. It is therefore considered by EPSG to be a different coordinate operation method to that described above.

1.4.8.2 Polar Stereographic (EPSG coordinate operation method code 9810).

For the forward conversion from latitude and longitude,

E = FE + r sin (l – l o)

N = FN – r cos (l – l o)

where

r = 2 a ko t /{[(1+e)(1+e) (1 – e)(1-e)]0.5}

t = tan (p /4 – j /2) / [(1 – esinj ) / (1 + e sinj )]e/2

For the reverse conversion,

j = c + (e2/2 + 5e4/24 + e6/12 + 13e8/360) sin(2c )

+ (7e4/48 + 29e6/240 + 811e8/11520) sin(4c )

+ (7e6/120 + 81e8/1120) sin(6c ) + (4279e8/161280) sin(8c )

l = l o + arctan [(E – FE) / (FN – N)]

where c = p /2 - 2 arctan t

t = r {[(1+e)(1+e) (1– )1-e]0.5} / 2 a ko

r = [(E-FE)2 + (N – FN)2]0.5

1.4.9 New Zealand Map Grid (EPSG coordinate operation method code 9811)

This projection system typifies the most recent development in the design and formulation of map projections where, by more complex mathematics yielding formulas readily handled by modern computers, it is possible to maintain the conformal property and minimise scale distortion over the total extent of a country area without being tied to a particular zone distributed about a set line such as a central meridian or standard parallel. Thus both North and South Islands of New Zealand, previously treated not very satisfactorily in two zones of a Transverse Mercator projection can now be projected as one zone of what resembles most closely a curved version Oblique Mercator but which, instead of being based on a minimum scale factor straight central line, has a central line which is a complex curve roughly following the trend of both North and South Islands. The projected coordinate reference system achieves this by a form of double projection where a conformal projection of the ellipsoid is first made to say an oblique Stereographic projection and then the Cauchy-Riemann equations are invoked in order to further project the rectangular coordinates on this to a modification in which lines of constant scale can be made to follow other than the normal great or small circles of Central meridians or standard parallels. The mathematical treatment of the New Zealand Map Grid is covered by a publication by New Zealand Department of Lands and Survey Technical Circular 1973/32 by I.F.Stirling.

1.4.10 Tunisia Mining Grid (EPSG coordinate operation method code 9816)

l_Paris (grads) = 7.83445 + [(E – 270) * 0.012185], where E is in kilometres.

The reverse formulas are:

E (km) = 270 + [(l_Paris – 7.83445) / 0.012185] where l_Paris is in grads.

N (km) = 360 + [(j – 36.5964) / B]

where j is in grads and B = 0.010015 if j > 36.5964, else B = 0.01002.

l_Greenwich (degrees) = [(l_Paris + 2.5969213) * 0.9] where l_Paris is in grads.

j _g (grads) = (j _d / 0.9) where j _d is in decimal degrees.

l_Paris (grads) = [(l_Greenwich / 0.9) – 2.5969213)] where l_Greenwich is in decimal degrees.

Example:

For grid location 302598,

Latitude j = 36.5964 + [(598 – 360) * A]. As N > 360, A = 0.010015.

Longitude l = 7.83445 + [(E – 270) * 0.012185, where E = 302.

• = 8.22437 grads east of Paris = 9.73916 degrees east of Greenwich.

1.4.11 American Polyconic (EPSG coordinate operation method code 9818)

This projection was popular before the advent of modern computers due to its ease of mechanical construction, particularly in the United States. It is neither conformal nor equal area, and is distortion-free only along the longitude of origin. A modified form of the polyconic projection was adopted in 1909 for the International Map of the World series of 1/1,000,000 scale topographic maps. A general study of the polyconic family of projections was published in 1919 (and reprinted in 1934) by Oscar Adams of the US Geological Survey. The mathematical treatment of the American Polyconic is covered on pages 124-131 of USGS professional paper 1395 by J.P.Snyder.

1.4.12 Lambert Azimuthal Equal Area (EPSG coordinate operation method code 9820)

To derive the projected coordinates of a point, geodetic latitude (Lat) is converted to authalic latitude (ß). The formulae to convert geodetic latitude and longitude (Lat,Lon) to Northing (Y) and Easting (X) are:

Northing, Y = FN + (B / D) . {(cos ß0 . sin ß) – [sin ß0 . cos ß . Cos(Lon – Lon0)]}

Easting, X = FE + {(B . D) . [cos ß . sin(Lon – Lon0)]}

where

B = Rq . (2 / {1 + sin ß0 . sin ß + [cos ß0 . cos ß . Cos(Lon – Lon0)]})^(1/2)

D = a . [cos Lon0 / (1 – e^2 sin^2 Lat0)^(1/2)] / (Rq . cos ß0)

Rq = a . (qp / 2)^(1/2)

ß = arcsin (q / qp)

ß0 = arcsin (q0 / qp)

q = (1 – e^2) . ([sin Lat / (1 – e^2 sin^2 Lat)] – {[1/(2e)] . ln [(1 – e sin Lat) / (1 + e sin Lat)]})

q0 = (1 – e^2) . ([sin Lat0 / (1 – e^2 sin^2 Lat0)] – {[1/(2e)] . ln [(1 – e sin Lat0) / (1 + e sin Lat0)]})

qp = (1 – e^2) . ([sin Latp / (1 – e^2 sin^2 Latp)] – {[1/(2e)] . ln [(1 – e sin Latp) / (1 + e sin Latp)]})

and where

f is the flattening of the ellipsoid

e is the eccentricity of the ellipsoid where e^2 = 2f – f^2

Lat is the latitude of the point to be converted, positive if North and negative if South of the equator

Lon is the longitude of the point to be converted, positive if East and negative if West of the prime meridian (Greenwich)

Lat0 is the latitude of the natural origin

Lon0 is the longitude of the natural origin (with respect to the prime meridian Greenwich)

X is the Easting measured from the grid origin

Y is the Northing measured from the grid origin

FE is the false easting, the Eastings value assigned to the natural origin

FN is the false northing, the Northings value assigned to the natural origin

The reverse formulas to derive the geodetic latitude and longitude of a point from its Northing and Easting values are:

Lat = ß’ + [(e^2/3 + 31e^4/180 + 517e^6/5040) . sin 2ß’] + [(23e^4/360 + 251e^6/3780) . sin 4ß’] + [(761e^6/45360) . sin 6ß’]

Lon = Lon0 + arctan {(X-FE) . sin C / [D. rho . cos ß0 . cos C – D^2. (Y-FN) . sin ß0 . sin C]}

where

ß’ = arcsin{(cosC . sin ß0) + [(D . (Y-FN) . sinC . cos ß0) / rho]}

C = 2 . arcsin(? / 2 . Rq)

rho = {[(X-FE)/D]^2 + [D . (Y –FN)]^2}^(1/2)

and D, Rq, and ß0 are as in the forward equations.

See USGS Professional Paper 1395, "Map Projections - A Working Manual" by John P. Snyder.

1.4.12.2 Lambert Azimuthal Equal Area (Spherical) (EPSG coordinate operation method code 9821)

The US National Atlas uses the spherical form of the map projection method, so exceptionally EPSG documents this. See USGS Professional Paper 1395, "Map Projections - A Working Manual" by John P. Snyder for formulas and example.

1.4.13 Albers Equal Area (EPSG coordinate operation method code 9822)

To derive the projected coordinates of a point, geodetic latitude (Lat) is converted to authalic latitude (ß). The formulae to convert geodetic latitude and longitude (Lat,Lon) to Northing (Y) and Easting (X) are:

Northing, Y = Nf + rho0 – (rho . cos theta)

Easting, X = Ef + (rho . sin theta)

where

theta = n . (Lon - Lon0)

rho = [a . (C – n . q)^(1/2)] / n

rho0 = [a . (C – n . q0)^(1/2)] / n

C = m12 + (n . q1)

n = (m1^2 – m2^2) / (q2 – q1)

m1 = cos Lat1 / (1 – e^2 sin^2Lat1)^(1/2)

m2 = cos Lat1 / (1 – e^2 sin^2Lat2)^(1/2)

q = (1 – e^2) . ([sin Lat / (1 – e^2 sin^2Lat)] – {[1/(2e)] . ln [(1 – e sinLat) / (1 + e sinLat)]})

q0 = (1 – e^2) . ([sin Lat0 / (1 – e^2 sin^2Lat0)] – {[1/(2e)] . ln [(1 – e sinLat0) / (1 + e sinLat0)]})

q1 = (1 – e^2) . ([sin Lat1 / (1 – e^2 sin^2Lat1)] – {[1/(2e)] . ln [(1 – e sinLat1) / (1 + e sinLat1)]})

q2 = (1 – e^2) . ([sin Lat2 / (1 – e^2 sin^2Lat2)] – {[1/(2e)] . ln [(1 – e sinLat2) / (1 + e sinLat2)]})

and where

f is the flattening of the ellipsoid

e is the eccentricity of the ellipsoid where e^2 = 2f – f^2

Lat is the latitude of the point to be converted, positive if North and negative if South of the equator

Lon is the longitude of the point to be converted, positive if East and negative if West of the prime meridian (Greenwich)

Lat1 and Lat2 are the latitudes of the standard parallels

Lat0 is the latitude of the projection (not grid) origin

Lon0 is the longitude of the projection origin (with respect to the prime meridian Greenwich)

X is the Easting measured from the grid origin

Y is the Northing measured from the grid origin

Ef is the false easting, the Eastings value assigned to the projection origin

Nf is the false northing, the Northings value assigned to the projection origin

The reverse formulas to derive the geodetic latitude and longitude of a point from its Northing and Easting values are:

Lat = ß’ + [(e^2/3 + 31e^4/180 + 517e^6/5040) . sin 2ß’] + [(23e^4/360 + 251e^6/3780) . sin 4ß’] + [(761e^6/45360) . sin 6ß’]

l = l0 + (theta / n)

where

ß’ = arcsin (q’ / {1 – [(1-e^2)/(2 . e)] . ln [(1-e) . (1+e)]})

q’ = [C – (rho^2 . n^2 / a^2)] / n

theta = arctan[(X-Ef) / [rho0 – (Y-Nf)]}

rho = {(X- Ef)^2 + [rho0 – (Y-Nf)]^2}^(1/2)

and C, n and rho0 are as in the forward equations.

See USGS Professional Paper 1395, "Map Projections - A Working Manual" by John P. Snyder.

1.4.14 Equidistant Cylindrical (EPSG coordinate operation method code 9823)

This method has one of the simplest formulas available. If the latitude of natural origin (j _O) is at the equator the method is also known as Plate Carrée. It is not used for rigorous topographic mapping because its distortion characteristics are unsuitable. Formulas are included to distinguish this map projection method from an approach sometimes mistakenly called by the same name and used for simple computer display of geographic coordinates – see Pseudo Plate Carrée below.

For the forward calculation of the Equidistant Cylindrical method:

X = R . (l - l_O) . cos(j _O)

Y = R . j

where R = [a² * (1 – e²)]^0.5

and j and l are expressed in radians.

For the reverse calculation:

j = Y / R

l = l_O + (X / R cos(j _O))

where R is as for the forward method.

1.4.14.1 Pseudo Plate Carrée (EPSG coordinate operation method code 9825)

This is not a map projection in the true sense as the coordinate system units are decimal degrees and therefore of variable linear scale. It is used only for depiction of graticule (latitude/longitude) coordinates on a computer display. The origin is at latitude (j ) = 0, longitude (l) = 0. See above for the formulas for the proper Plate Carrée map projection method.

For the forward calculation:

X = l

Y = j

where j and l are expressed in radians.

For the reverse calculation:

j = Y

l = X

1.4.14 Bonne (EPSG coordinate operation method code 9827)

The Bonne projection was frequently adopted for 18^th and 19^th century mapping, but being equal area rather than conformal its use for topographic mapping was replaced during the 20^th century by conformal map projection methods.

The formulas to convert geodetic latitude and longitude (j , l) to Easting and Northing are:

E = (r . sin T) + FE

N = (a . m0 / sin j 0 – r . cos T) + FN

where

m = cosj / (1 – e2sin2j )1/2

with j in radians and m0 for j 0, the latitude of the origin, derived in the same way.

M = a[(1 – e2/4 – 3e4/64 – 5e6/256 –....)j – (3e2/8 + 3e4/32 + 45e6/1024+....)sin2j

+ (15e4/256 + 45e6/1024 +.....)sin4j – (35e6/3072 + ....)sin6j + .....]

with j in radians and M0 for j 0, the latitude of the origin, derived in the same way.

r = a . m0 / sin j 0 + M0 – M

T = a . m (l – l 0) / r with l and l 0 in radians

For the reverse calculation:

X = E – FE

Y = N – FN

r = ± [X² + (a . m0 / sin j 0 – Y)²]^1/2 taking the sign of j 0

M = a . m0 / sin j 0 + M0 – r

μ = M / [a (1 – e²/4 – 3e⁴/64 – 5e⁶/256 – …)]

e₁ = [1 – (1 – e²)^1/2] / [1 + (1 – e²)^1/2]

j = m + (3e1/2 – 27e13/32 +.....)sin2m + (21e12/16 – 55e14/32 + ....)sin4m

+ (151e13/96 +.....)sin6m + (1097e14/512 – ....)sin8m + ......

m = cos j / (1 – e² sin² j )^1/2

If j 0 is not negative

l = l 0 + r {arctan[X / (a . m0 / sin j 0 – Y)]} / a . m

but if j 0 is negative

l = l 0 + r {arctan[– X / (Y – a . m0 / sin j 0 )]} / a . m

In either case, if j = ±90°, m = 0 and the equation for l is indeterminate, so use l = l 0.

1.4.15.1 Bonne (South Orientated) (EPSG coordinate operation method code 9828)

In Portugal a special case of the method with coordinate system axes positive south and west has been used for older mapping. The formulas are as for the general case above except:

W = FE – (r . sin T)

S = FN – (a . m0 / sin j 0 – r . cos T)

In these formulas the terms FE and FN retain their definition, i.e. in the Bonne (South Orientated) method they increase the Westing and Southing value at the natural origin. In this method they are effectively false westing (FW) and false southing (FS) respectively.

For the reverse formulas, those for the standard Bonne method above apply, with the exception that:

Part 2: Formulas for Coordinate Operations other than Map Projections

2.1. Introduction

It has been noted earlier that it is frequently required to convert coordinates derived in one geographic coordinate reference 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 reference system, whereas the national geodetic reference system in use for the country concerned will probably be a much earlier coordinate reference system. It may therefore be necessary to transform the observed WGS 84 data to the national geodetic reference system in order to avoid discrepancies. This form of coordinate transformation has to be between source and target geographical coordinate reference systems with the coordinates of both expressed in geographical terms. It is not executed by a transformation between projected coordinate reference systems.

The coordinate 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 recognise 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 gravity-related 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 gravity-related heights (surveyed other than by GPS) to ellipsoidal heights for use in a 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 coordinate transformation without introducing significant error in the horizontal coordinates.

In the early days of satellite surveying when relationships between geographic coordinate reference systems 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 geographic coordinate reference system to derive those in the second geographic coordinate reference system. This assumed, generally erroneously, that the axial directions of the two ellipsoids involved were parallel. For localised 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 parameter values.

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

The simplest coordinate transformation to implement involves applying shifts to the three geocentric coordinates. Molodenski developed a transformation method which applies the geocentric shifts directly to geographical coordinates. Both of these methods assume that the axes of the source and target coordinate reference 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 EPSG these are considered to be two different transformation methods, either a Position Vector transformation or a Coordinate Frame transformation. (The Position Vector transformation is elsewhere called the Bursa-Wolf transformation). It is crucial that the signs and interpretation of the transformation parameters are consistent with the convention being applied. However these methods may suffer from high correlation between the translations and rotations. An alternative approach that avoids this correlation is the Molodenski-Badekas transformation.

It is also possible to interpolate geographical or grid coordinate differences 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 (based on the Clarke 1966 ellipsoid) to coordinates expressed in the newer 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 EPSG as the Bi-linear Gridded Interpolation transformation technique. The transformation is achieved by bilinear 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 also been adopted by Australia and New Zealand. In Great Britain, bi-linear interpolation of gridded easting and northing differences is used.

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 transformation formulas, involving up to fourth order expressions of latitude and longitude expressed in degrees.

If and when other countries decide to convert to using coordinate reference systems based on 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 parameter values 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.

2.2 Coordinate Conversions other than Map Projections

2.2.1 Geographic/Geocentric Conversions

Latitude, j , and Longitude, l , in terms of a geographic coordinate reference system may be expressed in terms of a geocentric (earth centred) 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 prime meridian and equator, and the Y axis through the intersection of the equator with longitude 90° E.

Geocentric coordinate reference 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 reference 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 ellipsoidal semi major axis is a, semi minor axis b, and inverse flattening 1/f, then

X = (n + h) cos j cos l

Y = (n + h) cos j sin l

Z = ((1 – e2) n + h) sin j

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

n = a /(1 – e2sin2j )1/2,

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 e2 = (a2 – b2)/a2 = 2f – f2

(Note that h is the height above the ellipsoid.. This is the height value that 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 WGS 84 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 reference system may then be converted to the geographical coordinates in terms of geographic coordinate reference system by:

j = arctan (Z + e2n sin j ) / (X2 + Y2)1/2 by iteration

l = arctan Y/X

h = X sec l sec j – n

where l is relative to the Greenwich prime meridian.

See section 2.3.2 below for an example.

2.3 Coordinate Transformation Formulas

2.3.1 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 reference system coordinate to target coordinate reference system coordinate.

2.3.1.1 Vertical Offset (EPSG coordinate operation method code 9616)

If the coordinate reference system axes are positive in opposite directions, for instance in the transformation of heights in the source vertical CRS to depths in the target vertical CRS, or in differing units, the above formula is modified to:

X_t = [(X_s * U_s) + (A * U_A)] * (m / U_t)

where

X_t = value in the target vertical coordinate reference system.

X_s = value in the source vertical coordinate reference system.

A is the value of the origin of the target system in the source system.

m is unit direction multiplier (m=1 if both systems are height or both are depth; m = –1 if one system is height and the other system is depth; the value of m is implied through the vertical coordinate reference system type attribute).

U_s U_t and U_A are unit conversion ratios to metres for the source and target systems and the offset value respectively.

2.3.2 Geocentric Translations (EPSG coordinate operation method code 9603)

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 reference system to target geocentric coordinate reference system may then be applied as

Xt = Xs + dX

Yt = Ys + dY

Zt = Zs+ dZ

This example combines the geographic/geocentric conversion of section 2.2.1 above with the geocentric translation method.

Consider a North Sea point with coordinates derived by GPS satellite in the WGS84 geographical coordinate reference system, with coordinates of:

latitude j _s = 53° 48' 33.82" N,

longitude l _s = 2° 07' 46.38" E,

and ellipsoidal height h_s = 73.0m,

whose coordinates are required in terms of the ED50 geographical coordinate reference system which takes the International 1924 ellipsoid. The three parameter geocentric translations method's parameter values from WGS84 to ED50 for this North Sea area are given as dX = +84.87m, dY = +96.49m, dZ = +116.95m.

The WGS 84 geographical coordinates first convert to the following geocentric values using the formulas from section 2.2.1:

Xs = 3771 793.97 m

Ys = 140 253.34 m

Zs = 5124 304.35 m

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

Xt = 3771 793.97 + 84.87 = 3771 878.84 m

Yt = 140 253.34 + 96.49 = 140 349.83 m

Zt = 5124 304.35 + 116.95 = 5124 421.30 m

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

latitude j _t = 53° 48' 36.565" N,

longitude l _t = 2° 07' 51.477" E,

and ellipsoidal height h_t = 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.

2.3.3 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 (EPSG coordinate operation method code 9604). Abridged versions of these formulas (EPSG coordinate operation method code 9605), which are quite satisfactory for three parameter transformations, are as follows:

dj " = (– dX.sinj .cosl – dY.sinj .sinl + dZ.cosj + [a .df + f .da ].sin2j ) / (r . sin1")

dl " = (– dX.sinl + dY.cosl ) / (n.cosj . sin 1")

dh = dX.cosj .cosl + dY.cosj .sinl + dZ.sinj + (a .df + f .da).sin^2j – da

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 1.4), da is the difference in the semi-major axes of the target and source ellipsoids (at – as) and df is the difference in the flattening of the two ellipsoids (ft – fs).

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

2.3.4 Helmert transformation

and the parameters are defined as:

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

(R_X, R_Y, R_Z) : Rotations to be applied to the coordinate reference 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 reference frame when viewed from the origin of the Cartesian coordinate reference system in the positive direction of that axis, that is a positive rotation about the Z-axis only from source coordinate reference system to target coordinate reference system will result in a smaller longitude value for the point in the target coordinate reference 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 reference system in order to obtain the correct scale of the target coordinate reference system. M = (1+dS*10^-6), where 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:

Transformation parameters Coordinate Frame Rotation convention:

dX = 0.000 m

dY = 0.000 m

dZ = +4.5 m

R_X = 0.000 sec

R_Y = 0.000 sec

R_Z = – 0.554 sec

dS = +0.219 ppm

Please note that only the rotation has changed sign as compared to the Position Vector Transformation. The Position Vector convention is used by IAG and recommended by ISO 19111.

2.3.5 Molodenski-Badekas transformation

2.4 Coordinate Operation Methods that can be Operations or Transformations

In theory, certain coordinate operation methods do not readily fit the ISO 19111 classification of being either a coordinate conversion (no change of datum involved) or a coordinate transformation. These methods change coordinates directly from one coordinate reference system to another and may be applied with or without change of datum, depending upon whether the source and target coordinate reference systems are based on the same or different datums. In practice, most usage of these methods does in fact include a change of datum. EPSG follows the general mathematical usage of these methods and classifies them as transformations.

2.4.1 Polynomial Transformations

2.4.1.1 General Case

Polynomial transformations between two coordinate reference systems are typically applied in cases where one or both of the coordinate reference 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 coordinate operations between two projected coordinate reference 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.

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_S0 , Y_S0 are coordinates of the evaluation point in the source coordinate reference system,

X_T0 , Y_T0 are coordinates of the evaluation point in the target coordinate reference system.

and where

dX = A₀ + A₁.U + A₂.V + A₃.U² + A₄.U.V + A₅.V² (to 2^nd order)

dY = B₀ + B₁.U + B₂.V +B₃.U² +B₄.U.V +B₅.V²

+ B₆.U³ +B₇.U².V +B₈.U.V² +B₉.V³

+ B₁₀.U⁴ + B₁₁.U³.V + B₁₂.U².V² + B₁₃.U.V³ + B₁₄.V⁴

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).

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 reference systems. An example is the polynomial approximation between ED50 and ED87 geographical coordinate reference 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 coordinate 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 neighbouring 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 coordinate 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:

¾ ¾ ¾ ¾ ¾ » ¾ ¾ ¾ ¾ ¾ (» 10^-6)

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 reference 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

and

X_T = X_S + dX

Y_T = Y_S + dY

where dX and dY are derived from the polynomial expressions for the general polynomial.

Example: 4^th order reversible polynomial

For geodetic coordinate transformation ED50 to ED87 (1)

Offset unit: degree

Ordinate 1 of evaluation point X₀ = 55° 00' 00.000"N = +55 degrees

Ordinate 2 of evaluation point Y₀ = 0° 00' 00.000"E = +0 degrees

Parameters:

A₀ = – 5.56098E-06 A₁ = – 1.55391E-06 ... A₁₄ = – 4.01383E-09

B₀ = +1.48944E-05 B₂ = +2.68191E-05 ... B₁₄ = +7.62236E-09

Forward calculation for:

Latitude j _ED50 = X_s = 52° 30' 30" N = +52.508333333 degrees

Longitude l _ED50 = Y_s = 2° E = +2.0 degrees

U = X_S – X_{0 =} j _ED50 – X₀ = 52.508333333 – 55.0 = -2.491666667 degrees

V = Y_S – Y₀ = l _ED50 – Y₀ = 2.0 – 0.0 = 2.0 degrees

dX = A₀ + A₁.U + ... + A₁₄.V⁴

= – 5.56098E-06 + (– 1.55391E-06 * – 2.491666667) + ... + (– 4.01383E-09 * 2.0^4)

= – 3.12958E-06 degrees

dY = B₀ + B₁.U + ... + B₁₄.V⁴

= +1.48944E-05 + (2.68191E-05 * – 2.491666667) + ... + (7.62236E-09 * 2.0^4)

= +9.80126E-06 degrees

Then Latitude j _ED87 = X_T = X_S + dX

Longitude l _ED87 = Y_T = Y_S + dY

= 2° 00' 00.0353" E

Reverse calculation for coordinate transformation ED50 to ED87 (1).

The transformation method for the ED50 to ED87 (1) coordinate 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 degrees

Ordinate 2 of evaluation point Y₀ = 0° 00' 00.000"E = +0 degrees

A₀ = +5.56098E-06 A₁ = +1.55391E-06 ... A₁₄ = +4.01383E-09

B₀ = – 1.48944E-05 B₁ = – 2.68191E-05 ... B₁₄ = – 7.62236E-09

Reverse calculation for:

Latitude j _ED87 = X_S = 52° 30' 29.9887"N = +52.5083301944 degrees

Longitude l _ED87 = Y_S = 2° 00' 00.0353" E = +2.0000098055 degrees

U = 52.5083301944 – 55.0 = – 2.4916698056 degrees

V = 2.0000098055 – 0.0 = 2.0000098055 degrees

dX = A₀ + A₁.U + ... + A₁₄.V⁴

= +5.56098E-06 + (1.55391E-06 * -2.491666667) + ...

= +3.12957E-06 degrees

dY = B₀ + B₁.U + ... + B₁₄.V⁴

= – 1.48944E-05 + (– 2.68191E-05 * – 2.491666667) + ...

= – 9.80124E-06 degrees

Then Latitude j _ED50 = X_T = X_S + dX

Longitude l _ED50 = Y_T = Y_S + dY

= 2° 00' 00.000" E

2.4.1.2 Polynomial transformation with complex numbers

The relationship between two projected coordinate reference 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.

The fourth order polynomial (EPSG coordinate operation method code 9632) can alternatively be expressed in matrix form as

Then as for the general polynomial case in section 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_S0 , Y_S0 are coordinates of the evaluation point in the source coordinate reference 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 reference 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:

Coordinate transformation: Amersfoort / RD New to ED50 / UTM zone 31N (1):

Formula Parameter

EPSG coordinate transformation parameter name symbol value Unit

------------------------------------------------------------------- ----------- ----------------- --------

ordinate 1 of the evaluation point in the source CRS X_S0 155,000.000 metre

ordinate 2 of the evaluation point in the source CRS Y_S0 463,000.000 metre

ordinate 1 of the evaluation point in the target CRS X_T0 663,395.607 metre

ordinate 2 of the evaluation point in the target CRS Y_T0 5,781,194.380 metre

A1 A₁ – 51.681 coefficient

A2 A₂ +3,290.525 coefficient

A3 A₃ +20.172 coefficient

A4 A₄ +1.133 coefficient

A5 A₅ +2.075 coefficient

A6 A₆ +0.251 coefficient

A7 A₇ +0.075 coefficient

A8 A₈ – 0.012 coefficient

For input point:

Easting, X_{AMERSFOORT/RD} = X_S = 200,000.00 metres

Northing, Y_{AMERSFOORT/RD} = Y_S = 500,000.00 metres

U = (X_S – X_S0) .10⁵ = (200,000 – 155,000) .10^-5 = 0.45

V = (Y_S – Y_S0) .10⁵ = (500,000 – 463,000) .10^-5 = 0.37

dX = – 1,240.050

dY = 1,468.748

Then

Easting, E_ED50/UTM31N = X_T = X_S – X_S0 + X_T0 + dX

Northing, N_ED50/UTM31N = Y_T = Y_S – Y_S0 + Y_T0 + dY

2.4.1.3 Polynomial transformation for Spain (EPSG coordinate operation method code 9617)

The original geographic coordinate reference 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 transform 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 transform coordinates from the Madrid 1870 (Madrid) geographic coordinate reference 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 (i) a mean for all Spain, (ii) a better fit for the north of Spain, (iii) a better fit for the south of Spain.

The polynomial expressions are:

dj (arc sec) = A₀ + (A₁*j _s) + (A₂*l _s) + (A₃*H_s)

dl (arc sec) = 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 reference system and H_s is gravity-related height in metres. B₀₀ is the longitude (in 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:

j _ED50 = j _M1870(M) + dj

and l _ED50 = l _M1870(M) + dl .

Example:

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

Latitude j _s = 42° 38' 52.77" N

= +42.647992 degrees

Longitude l _s = 3° 39' 34.57" E of Madrid

= +3.659603 degrees from the Madrid meridian.

Gravity-related height H_s = 0 m

For the north zone transformation:

A₀ = 11.328779 B₀₀ = -13276.58

A₁ = -0.1674 B₀ = 2.5079425

A₂ = -0.03852 B₁ = 0.8352

A₃ = 0.0000379 B₂ = -0.00864

B₃ = -0.0000038

Then latitude j _ED50 = 42° 38' 52.77" N + 4.05"

dl = -13270.54 seconds = -3° deg 41' min 10.54" sec

Then longitude l _ED50 = 3° 39' 34.57" E - 3° 41' 10.54"

2.4.2 Miscellaneous Linear Coordinate Operations

2.4.2.1 The Affine Parametric Transformation (EPSG coordinate operation method code 9624)

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:

X_T = A₀ + A₁. X_S + A₂.Y_S

Y_T = B₀ + B₁. X_S + B₂.Y_S

where

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

The parameter values for an affine transformation cannot be used for the reverse operation. However, the reverse operation 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 operation as follows:

D = A_{1 .} B₂ – A_{2 .} B₁

A₀' = (A_{2 .} B₀ – B_{2 .} A₀) / D

B₀' = (B_{1 .} A₀ – A_{1 .} B₀) / D

A₁' = +B₂ / D

A₂' = – A₂ / D

B₁' = – B₁ / D

B₂' = +A₁ / D

2.4.2.2 The Affine General Geometric Transformation

(EPSG coordinate operation method code 9623)

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

The operation will often be applied to transform an engineering coordinate reference system to a projected coordinate reference system. The engineering coordinate reference 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 reference system.

The projected coordinate reference 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 reference 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 (engineering) coordinate reference system as the point in which to calculate this point scale factor ‘k’, although for engineering coordinate reference 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: