SURV3510 Geodesy

Contents

Chapter 1 Introduction to Geodesy

Chapter 2 Gravity of the Earth

Chapter 3 Coordinate System

Chapter 4 Geodetic Reference System

Chapter 5‐6 Geodetic Datums and Geodetic datum of Australia

Chapter 7 Map projection

Chapter 8 Curves on the Surface of the Ellipsoid

Chapter 9 Reduction of Measured Distances to the Ellipsoid

Chapter 10 Geodetic Techniques and Applications

 1 Introduction to Geodesy
The practice of geodesy has changed dramatically in the last few decades, mainly through the
impact of satellite methods. For example, positioning for most geodetic applications is now
almost exclusively performed with GPS (Global Positioning System) and the finer structure of
the earth’s gravity field is being determined over the oceans by satellite altimetry and over both
the land and ocean by recent satellite gravity field missions (CHAMP (Challenging
MiniSatellite Payload for Geophysical Research and Applications) launched on July 15, 2000;
GRACE (Gravity Recovery and Climate Experiment) launched on March 17, 2002; and GOCE
(Gravity Field and Steady-State Ocean Circulation Explorer) scheduled for launch in 2005).
Geodesy today is truly an interdisciplinary science, overlapping with fields such as geophysics,
oceanography, meteorology, and glaciology to name a few.

1.1 What is Geodesy?

Geodesy is a science, the oldest earth or geoscience in fact. The origin of the terminology can
be traced back to the times of the Greeks, in the time of Aristotle (384-322 BC). The Greek
word “” (geodesia) means dividing the earth and came from  (Earth) and  (I
divide). The ancient Greeks were able to use astronomy in order to determine the size and shape
of the Earth (they knew it was round!) and gained a basic understanding of the Earth's motion in
relation to other celestial bodies. In around 600 BC, Pythagoras was the first to consider a
spherical earth while some years later, Aristotle continued this notion and also first discussed
the concept of gravity. One of the most famous geodesists at the end of this era of discovery
was Ptolemy who published many geodetic and astronomical findings along with the first
extensive version of a world map.

1.1.1 Definitions of Geodesy

To understand the essential characteristics of the discipline of geodesy, it is useful to look at the
historical definitions used during its development…
1. The task of geodesy is the determination of the potential function (Bruns, 1878).
2. Geodesy is the science of measuring and portraying (mapping) the Earth’s surface
(Helmert, 1880).
3. Geodesy is a branch of science which investigates methods to accurately measure
elements of the Earth’s surface and to determine from them geographic positions of
points on this surface and which studies the figure of the Earth from a theoretical point
of view and by evaluating results of measurements (Zakatov, 1957).
4. Geodesy is both theoretical and practical. Its theoretical function is to determine the
size and shape of the Earth and, in conjunction with other Earth sciences, to study the
structure of the Earth’s crust and of the immediately underlying layers. Its practical
function is to perform the measurements and computations that will give the
coordinates of selected control points on the Earth’s surface, i.e., to fix their positions
on the Earth’s surface. (Heiskanen and Vening-Meinesz, 1958).
5. Au sens étymologique du mot, la géodésie est la science qui a pour objet la mesure des
dimensions de la Terre. Déterminer, d’une part, la forme et les dimensions précises de
la planète; réaliser, d’autre part, principalement au moyen de triangulations, la
mensuration des territoires terrestres pour permettre d’en dresser des cartes exactes et
fournir des données géométriques précises pour les diverses enterprises de l’ingénieur,
sont en effet les buts principaux, scientifiques et practiques de l’activité des
géodésiens (Dupuy et Dufour, 1969).
6. Geodesy is the discipline that deals with the measurement and representation of the
earth, including its gravity field, in a three-dimensional time varying space (Vaníček
and Krakiwsky, 1986).

1-1
 7. Geodesy is considered as a discipline which deals mainly with the mapping of the
Earth and the monitoring of variations at its surface. From the very beginning those
tasks were connected with the gravity vector g (absolute value and direction) (Groten,
1979).
8. The problem of geodesy is to determine the figure and the external gravity field of the
Earth and of the other heavenly bodies as a function of time; as well as to determine
the mean Earth ellipsoid from parameters observed on and exterior to the Earth’s
surface (Torge, 1980).
9. The study of geodesy is the study of the shape and gravity field of the Earth. The
methods used to carry out these investigations include terrestrial triangulation,
levelling, and gravity surveys and the methods based on the technologies developed
for space research, such as the tracking of artificial satellites and the Moon. …. With
the broad definition of geodesy to include the study of crustal motion, the spatial and
temporal variations in the gravity field and the tidal and rotational deformations of the
Earth, the geodetic observations play an important role in the study of the structure
and evolution of the Earth. …. Geodesy is high frequency geology and low frequency
seismology (Lambeck, 1989).

As can be appreciated from the above definitions, the original focus of geodesy has expanded
from the classical definition of F.R. Helmert (1880), from being the ‘science of the
measurement and mapping of the earth’s surface’, to one that now recognises that the Earth is
dynamic and that the discipline has expanded to include applications in oceanography,
planetary and atmospheric sciences, geophysics, glaciology, geodynamics, etc.

1.1.2 Modern Geodesy

The main focus of all modern geodesy is directed towards the determination of :
(i) the size and shape of the Earth and other celestial bodies
(ii) the Earth’s (and planet’s) gravity field(s)
(iii) the temporal variations of (i) and (ii)
(iv) the orientation of the Earth in space
In this extended definition, geodesy is considered to be part of the geosciences and engineering
sciences, including navigation and geomatics. As a rough generalisation, it can be said that
there are three principal types of activity in geodesy:
 terrestrial or classical geodesy based on measurements of the geometry and gravity
of points on the surface of the Earth;
 space geodesy based on observations of and from artificial satellites
or extraterrestrial radio sources (quasars);
 theoretical geodesy, which is concerned with the analysis and
interpretations of measurements and fundamental theory.
Today, most practical geodesy is performed using space geodetic techniques. The primary
techniques, in order of highest to lowest accuracy/precision, are:
 Very Long Baseline Interferometry (VLBI),
 Satellite Laser Ranging (SLR),
 Lunar Laser Ranging (LLR),
 Global Positioning System (GPS),
 DORIS (Doppler Orbitography by Radiopositioning Integrated on Satellite),
 GLONASS (Russian version of GPS),
 PRARE (Precise Range and Range Rate Experiment).

1-2
The boundaries of our measurement accuracies and precisions are always been pushed to new
limits. Relative positioning accuracies of a few parts in 1011, which is equivalent to a few
millimetres in 100,000 km (the Earth’s radius is about 6,371 km), have been achieved with
VLBI and SLR methods. Relative accuracies of few parts in 1010 (1 cm in a 100,000 km line)
are routinely achievable with GPS and DORIS.

1.1.3 Sub-Disciplines of Geodesy

(i) Geometrical Geodesy – relates to the representation of measurements made on the
Earth’s surface in a consistent reference system/mathematical model.
(ii) Physical Geodesy – involves the study of the Earth’s gravity field and its temporal
variations.
(iii) Mathematical Geodesy – deals with geodetic boundary value problems incorporating a
geometric (figure of the Earth) and a physical (gravity field) formulation of the
problem.
(iv) Dynamic/Satellite Geodesy – involves study of the orbital motion of planets and
artificial satellites, and space geodetic techniques.
(v) Geodetic Astronomy – based on the rules of spherical trigonometry, this discipline
deals with the determination of the orientation of the local gravity vector
(astronomical latitude  and longitude ) and the astronomical azimuth A of a
terrestrial mark fixed in relation to natural celestial bodies, particularly fixed stars and
the Sun. With the advent of space geodesy, the amount of geodetic astronomy now
practised is minimal.

There are no sharp or distinct separations of the above fields of geodesy, they are used more as
a generic classification. An alternative classification can be made based on spatial scales, in
which geodesy may be sub-divided into the following areas:
(i) Global Geodesy - which includes the determination of the size and shape of the Earth,
its orientation in space, its external gravity field and the temporal variations of the
dynamic Earth.
(ii) Regional/National Geodesy - involves regional/national geodetic surveys in which the
fundamentals for the determination of the Earth’s surface and gravity field within the
region are established within or related to the global reference system. This is realised
by coordinate and/or gravity values of a sufficiently large number of control points,
arranged in geodetic and gravimetric networks. In this case, curvature and the gravity
field of the Earth must be considered. This area involves geodetic surveys.
(iii) Local Geodesy - or ‘plane’ surveying is usually done on a horizontal plane and the
variations of the Earth’s gravity field are ignored.

There is a close relation between global geodesy, geodetic surveying, and plane surveying. The
connections between these spatial groupings follow the classical hierarchical system. Geodetic
surveys are linked to reference frames (networks) established by global geodesy; these surveys
adopt the global parameters for the figure of the Earth and the gravity field. Plane surveys, in
turn, are generally referenced to control points established by geodetic surveys. Plane surveys
are used extensively in the development of the national and state map series, digital cadastral
information systems, spatial databases for GIS applications, and in engineering projects. The
measurement and data evaluation methods used in national geodetic surveys are often similar to
those used in global geodetic work. These days, even local surveys make much use of GPS as a
measurement tool.
The example of how this hierarchical system works for Australia is the following (also see
http://www.ga.gov.au/earth-monitoring/geodesy.html):
 ITRF (International Terrestrial Reference Frame) – global reference frame.

1-3
 AFN (Australian Fiducial Network) – based on a set of 10 fundamental stations on
mainland Australia and offshore territories, fixed within the ITRF92.
 ANN (Australian National Network) – a ~ 500 km network of stations across
Australia, related to the AFN.
 TFN (Tasmanian Fiducial Network) – a network of control within Tasmania (approx.
at 100 km spacing), based on the control from the AFN/ANN.
 Local State networks and subsequent densification surveys.

ITRF2000 Network

1-4
 AFN ANN

Figure 1.1 – Global – Regional – National - Local Networks

1.2 Goals of Geodesy

1. Establishment and maintenance of precise global, regional and local 3-D geodetic
networks, recognising the time dependent nature of the networks. Points can be
positioned individually or as part of a network of points, and defined in a specified
reference frame (‘absolute’) or in relation to other points (relative).
2. Determination of the Earth’s gravity field and functions of this field (temporal
variations included). This goal is important in order to transform the geodetic
observations made in physical space into a geometric space.
3. Measurement and modelling of geodynamic phenomena, e.g. polar motion, Earth tides,
crustal motion, postglacial rebound, length-of-day changes (lod), ocean and atmospheric
loading, etc. Such geodynamic phenomena can be secular, periodic or abrupt and can
occur over a wide range of spatial scales (from global to local). Today our geodetic
techniques are capable of detecting such changes at the few mm level and hence we
must consider these phenomena as time dependent variables.
A short summary of possible applications of these goals for the different areas of geodesy is:
 Global Geodesy
- general shape of the Earth’s figure and gravity field
- dimensions of the mean Earth ellipsoid
- establishment of a global terrestrial reference frame
- determination of the geoid as a reference surface on land and at sea
- connection between different geodetic datums
- connection of national datums within a global geodetic datum.

1-5
 National/Regional Geodesy
- establishment of geodetic control for national networks
- maintenance of 3-D homogeneous/heterogeneous networks
- analysis and improvement of existing terrestrial networks
- establishment of geodetic connections between land areas
- densification of existing networks
 Applied & Plane Geodesy
- detailed plane surveying (land registers, GIS, town planning, etc)
- specialised networks (e.g., for engineering tasks)
- terrestrial control for photogrammetry and remote sensing
 Geodynamics
- control points for crustal motion studies
- polar motion, earth rotation, lod, etc studies
- solid earth tides, isostatic rebound studies
 Navigation & Marine Geodesy
- precise navigation of land, sea, and air vehicles
- precise positioning for marine mapping, exploration, hydrography, oceanography,
marine geology, geophysics
- connection of tide gauges (unification of height systems)
 Related Geophysical Fields
- position and velocity determination for different platforms – air, sea, land
- ice motion, ice sheet mass balance studies
- sea level change in the global oceans using satellite altimeters

1.3 International Geodesy

1.3.1 The International Union of Geodesy & Geophysics (IUGG)

The IUGG was founded in 1919 and is the umbrella organisation for international cooperation
between the seven member associations. The associations covered by the IUGG include the
International Association of Geodesy (IAG), the IASPEI (Seismology), the IAVCEI
(Volcanology), the IAPSO (Physical Oceanography), the IAMAP (Meteorology), the IAGA
(Geomagnetism and Aeronomy) and the IAHS (Hydrological Sciences).

1.3.2 The International Association of Geodesy (IAG)

[see http://www.iag-aig.org]
The IAG is an organisation devoted to the study of scientific problems of geodesy, and the
promotion of international cooperation in this field. Since 1999 the IAG has been undergoing a
major restructure resulting in a more focussed and organised approach to global collaboration in
geodetic studies.
"The Mission of the Association is the advancement of geodesy, an Earth
science that includes the study of the planets and their satellites. The IAG
implements its mission by advancing geodetic theory through research and
teaching, by collecting, analysing, and modelling observational data, by
stimulating technological development and by providing a consistent
representation of the figure, rotation, and gravity field of the Earth and planets
and their temporal variations".

1-6
The IAG is overseen by the Bureau, which consists of the President (Gerhard Beutler
(Switzerland)), the Vice President (Michael Sideris (Canada)) and the Secretary General (C. C.
Tscherning (Denmark)). There is also an executive committee that is comprised of the Bureau,
Past President, Commission Presidents, Three Service Representatives, President of the
Communication and Outreach Branch and Two members at large (geographical coverage).
Membership of the Association consists of Countries, Candidate Members, Individual Members
and Fellows. Academic links to IAG in Australia are through the Australian Academy of
Sciences, and the National Committee for Solid Earth Sciences. For the latest Australian IAG
report see a reference paper or the web site at http://iag.dgfi.badw.de/fileadmin/IAG-
docs/NationalReports2011/Australian_National_Report_2008-2011.pdf.

The IAG participates in the IUGG General Assemblies, which are held every 4 years. Between
each of the General Assembly meetings, the Associations usually hold a General Scientific
Meeting with scientific presentations on several selected topics. The IAG used to publish two
quarterly journals, Bulletin Géodésique and Manuscripta Geodaetica. In 1996, these two
publications were merged to produce a single monthly journal, the Journal of Geodesy (see
http://link.springer.de/link/service/journals/00190/index.htm for free electronic access to
articles). Every four years the IAG also produces ‘The Geodesist’s Handbook’, which contains
useful information about the IAG structure, summaries of all fundamental constants, etc.
The physical structure of the IAG is organised into the following components:
 Commissions
 Services
 Inter-commission Committees
 IAG Projects
 Communication and Outreach Branch

1-7
There are 4 Commissions within the IAG that represent the major interests of the association.
These Commissions are:
1) Reference Frames
2) Gravity Field
3) Earth Rotation and Geodynamics
4) Positioning and Applications
Each Commission has a President and Steering Committee who define the appropriate sub-
structure of the Commission, which may consist of Sub-commissions, Study/Working Groups,
Commission Projects and Inter-Commission Projects.
The Services of the IAG are (generally) independent organisations which provide specific
information for various sections of the Geodetic community and include the:
 IERS (International Earth Rotation and Reference Systems Service)
 IGS (International GPS Service)
 ILRS (International Laser Ranging Service)
 IVS (International VLBI Service for Geodesy and Astrometry)
 IGFS (International Gravity Field Service)
 IDS (International DORIS Service)
 BGI (International Gravimetric Bureau)
 IGES (International Geoid Service )
 ICET (International Centre for Earth Tides) (Belgium)
 PSMSL (Permanent Service for Mean Sea Level)
 BIPM (Bureau International de Poids and Measure - time section)
 IBS (IAG Bibliographic Service).

Inter-commission Committees were created during the restructure to handle important and
permanent tasks involving all Commissions (e.g. geodetic theory, geodetic conventions, space
geodesy) and therefore requiring input from and communication to each of the Commissions.
IAG Projects are also a new feature of the IAG structure. This section was designed to allow
IAG Projects of a broad scope, long duration (> 10 years) and of highest interest and
importance for the entire field of geodesy to be designed and implemented. Currently there is
only one project, the Integrated Global Geodetic Observing System (IGGOS).
The Communication and Outreach Branch is the third of the new elements of the IAG structure
and was designed to “provide the Association with communication, educational/public in-
formation and outreach links to the membership, to other scientific Associations and to the
world as a whole.”

1.4 Australian Geodesy

In Australia, geodetic activities in support of the National geodetic infrastructure are
coordinated through the Federal and State government agencies by the Intergovernmental
Committee on Surveying and Mapping (ICSM) – see http://www.icsm.gov.au/ . This committee
has established a working level Geodesy Group, with representatives from Federal and State
agencies, which plans and implements approved national projects (e.g.,
http://www.icsm.gov.au/geodesy/index.html). In practice, most of Australia’s regional geodetic
program is managed by Geoscience Australia (federal agency). Research in geodesy in
Australia is mostly undertaken through Universities, principally funded by grants from the
Australian Research Council (ARC).

1-8
1.4.1 Geoscience Australia (GA)
[see http://www.ga.gov.au/earth-monitoring/geodesy.html ]
Geoscience Australia is part of the Industry, Tourism and Resources government portfolio and
is the national agency for geoscience research and geospatial information.
The National Mapping Division and the Earth Monitoring Group represent the geodetic arm of
GA and are responsible for all of the geodetic surveys and space geodetic activities within
Australia. The main focus of Geodesy in GA is in regard to Australia’s contribution to the
determination of changes in the shape, size, position, and orientation of the Earth in space and
assessing the significance of these changes in the Australian context. The function of the
geodetic program is to make accurate measurements of these parameters for:
1. the development and maintenance of the fundamental national positioning
(geospatial) infrastructure in Australia, including positioning infrastructure for
navigation (for example, GPS)
2. monitoring of global change, including climate change (e.g. sea-level rise) and
tectonic movements of, and within, the Australian continent.

Examples of Geodetic Programs Undertaken by GA:

 Positioning
- Data/solutions for definition of the International Terrestrial Reference Frame (ITRF)
- Geocentric Datum of Australia (GDA)
- Australian Fiducial Network (AFN)
- Australian National Network (ANN)
- Australian Regional GPS Network (ARGN)
- Local Survey Connections at Fundamental Geodetic Sites
- National Geodetic Data Base
- Australian Baseline Sea Level Monitoring Array (ABSLMA)
- Inertial Positioning

The ARGN
 Advanced Space Technology
- VLBI
- SLR/LLR
- GPS
- Geodetic Positioning and Monitoring of Tide Gauge Datums
- Altimeter Calibrations

1-9
  Determination of the Gravity Field
- Gravity Networks
- Geoid Computations
 General Theory and Methodology
- Local Scale Parameter Method
- Fundamental Constants - Refractive Index of Air
- Geoid Related Studies
- GPS
- High Precision Engineering Surveys
- Linear Estimation
 Geodynamics
- Crustal Movements
- Satellite Altimetry
 Contributions to the International GPS Service (IGS)

1.4.2 Australian Universities

There are several Australian universities carrying out high-level geodetic research. These
projects are usually overseen by a member of the academic staff and worked on by a group of
staff and post-graduate students. Many of the university geodesy groups work together in
collaboration on some projects and/or share equipment. Each university has its area of expertise
or specialisation within the geodesy discipline and overall the Australian research quality is of a
very high, international standard. The major university geodesy groups in Australia include (in
no particular order):
 UNSW (NSW) - http://www.gmat.unsw.edu.au/
 Curtin (WA) - http://geodesy.curtin.edu.au/
 RMIT University (Vic) - http://www.rmit.edu.au/SPACE
 University of Tasmania (Tas) - http://www.utas.edu.au/spatial/
 ANU(RSES) (ACT) - http://rses.anu.edu.au/
 University of Newcastle - http://www.newcastle.edu.au/school/engineering/

Note: For a good overview of the different aspects of geodesy check out the web tutorial
developed by National Resources Canada at
http://www.geod.nrcan.gc.ca/edu/geod/index_e.php.

1-10
2 Gravity Field of the Earth

In this chapter we will introduce the theory of the Earth's gravity field, its determination and its
relationship to geometrical investigations.
Remember that gravity is not just a good idea… it’s the law!

2.1 Newton's Laws of Gravitation

Since in Geodesy we are dealing with stationary or slow moving objects, we need to apply the laws of
Newton’s gravitational theory:

−k m1 m2
F= (2.1)
s2

where
k is the gravitational constant (~ 6.67259 x 10-11)
m1 and m2 are masses of attracting bodies separated by distance s (a vector)
F is the gravitational force
This attraction is symmetric in opposite directions and propagates with the speed of light.
If a unit mass is assigned to the attracted mass (i.e. m2 = 1) and the attracting mass is m (i.e. m1 = m) then
(2.1) becomes

−k m
F= (2.2)
s2
The negative sign indicates a convention in the direction of the vector force, i.e. the force points from the
attracted to the attracting mass (m2 → m1).

Now if the mass m has coordinates (Xm, Ym, Zm) and the unit mass has cords (X, Y, Z) in a 3D Cartesian
system, then our force can be written as:
r r r
F = FX i + FY j + FZ k (2.3)
where
FX, FY, FZ are components of F parallel to axes X, Y, Z
r r r
i , j , k are unit vectors along the axes

Therefore, we can write:

X − Xm
FX = − km
s3
Y −Y
FY = − km 3 m
s (2.4)
Z − Zm
FZ = − km
s3
2-1
or
⎡ X − Xm ⎤
⎢ s ⎥
⎢ ⎥
km ⎢ Y − Ym ⎥
F =− 2 (2.5)
s ⎢ s ⎥
⎢ ⎥
⎢ Z − Zm ⎥
⎣⎢ s ⎦⎥
where
s = [(X - Xm)2 + (Y - Ym)2 + (Z - Zm)2]½

2.1.1 Gravitational Potential

The amount of work required to move a unit mass from infinity to a distance s from the attracting mass is
called gravitational potential (V).
km
V= (2.6)
s
or
s
ds
V = km∫ (2.7)
∞
s2
Substituting (2.6) into (2.7) and taking partial derivatives with respect to X, Y, Z gives:

dV d ⎛1⎞ ⎧ d d ds ⎫
= km ⎜ ⎟ ⎨ = ⋅ ⎬
dX dX ⎝ s ⎠ ⎩ dX ds dX ⎭
−km ds
= 2 ⋅
s dX
−km X − X m
= 2 ⋅ = FX
s s
(2.8)

dV − km Y − Ym
= 2 ⋅ = FY
dY s s

dV − km Z − Z m
= 2 ⋅ = FZ
dZ s s

Hence F and V are related by gradient, i.e. F is a gradient vector of a scalar function V:
⎛ dV dV dV ⎞
F = ( FX , FY , FZ ) = grad V = ∇V = ⎜ , , ⎟ (2.9)
⎝ dX dY dZ ⎠
The replacement of vector F by scalar V is helpful when considering several attracting masses acting on a
point mass at the same time.
We can say that the potential of all attracting masses is:
kmi
VT = ∑ Vi = ∑ (2.10)
si

2-2
Considering a solid body of mass M (i.e. a system of infinitesimal point masses), the discrete summation
is replaced by a volume integral over the body (ρ = density):

dm ρ dv
V( X ,Y , Z ) = k ∫∫∫ = k ∫∫∫
M
s earth
s (2.11)
↑ applies for a body at rest

2.2 Gravity Field of the Earth

2.2.1 Gravity Potential & Gravity Force

In order to account correctly for the earth's potential, the earth's rotation must be taken into account.
The earth's rotation introduces an additional force that acts on all parts, called centrifugal force.
In a 3D Cartesian system with the origin at the centre of mass of the body and the z-axis in the plane of
(coincident with) the rotation axis, the centrifugal force is given by:
q = ω 2 ( X 2 + Y 2 )1 / 2 = ω 2 p (2.12)
−

X
p f

p
Y

ω = angular velocity of earth

p = perpendicular distance from Pt to Z axis

Figure 2.1 - Forces Acting on a Point on the Earth.

• Centrifugal force (f) acts parallel to the equatorial plane (or, in the same direction as p) and
• Centrifugal force (f) is subject to variations over time in both direction and magnitude.

If q is regarded as a gradient vector of a scalar potential, the centrifugal potential is:

ω2 1
Φ= p2 = ω 2 ( X 2 + Y 2 ) (2.13)
2 2

2-3
The potential of gravity on the earth (W) is the sum of the gravitational potential (V) and centrifugal
potential (Φ):
W =V +Φ
ρ dV 1 (2.14)
= k ∫∫∫ + ω2 ( X 2 +Y2 )
V
s 2
The resultant of F and f is called the gravity force.
• Gravity force is stronger at the poles than at the equator because of the decrease towards the pole of
the centrifugal force.

{Note: W = direction of the plumbline}

Z

W
V
Y

Figure 2.2. The Potential of Gravity (W) as a Resultant of the Gravitational Potential (V) and the
Centrifugal Potential (Φ).

2-4
The gradient vector of gravity (W) is the gravity force vector g:
⎛ dW dW dW ⎞ ⎛ dW r dW r dW r ⎞
g = ∇W = ⎜ , , ⎟=⎜ i+ j+ k⎟ (2.15)
⎝ dX dY dZ ⎠ ⎝ dX dY dZ ⎠

The components of g are:

dW X − Xm
gX = = −k ∫∫∫ 3
ρ dV + ω X2
dX V
s
dW Y − Ym
gY = = − k ∫∫∫ 3
ρ dV + ωY2 (2.16)
dY V
s
dW Z − Zm
gZ = = − k ∫∫∫ 3
ρ dV
dZ V
s

• The gravity force g is a total force acting on a unit mass of coordinates (X, Y, Z) and is the resultant
of the gravitational and centrifugal forces of F and q.
• The magnitude of g is called gravity and has dimensions of acceleration (eg. ms-2).
• Gravity is traditionally measured in units of gals, where
1 gal = 1 cm sec-2
1 mgal = 1-3 gal = 10-5 m sec-2
The most common term used is mgals.
• Gravity meters now have a resolution of order 1-2 μgal (i.e. 10-6 gal).
• The gravity vector is often called gravity acceleration. The direction of g at a point is the direction
of the plumbline or direction of the vertical.

• Approximate values of gravity are:

gequator = 978 gals
gpoles = 983 gals.

2.2.2 Equipotential Surfaces

• The simplest way to use the gravity potential (W) to characterise the gravity field is to use
equipotential surfaces and lines of force (i.e. plumblines).

o Definition:

An equipotential surface is the surface, on which gravity potential W(X,Y,Z) is a

constant, i.e.

W(X,Y,Z) = constant

• In reality there are an infinite number of potential surfaces.

• In geodesy we use a few of these potential surfaces for fundamental reference definition
o at W(X, Y, Z) = W0 = constant, we have approximate sea level.

2-5
• The equipotential surfaces can be determined by evaluating the relation of the gravity potential on
the Earth (Eq. 2.14), if the density (ρ) and centrifugal force (i.e. angular momentum - ω) are
known.
o However, this is not very convenient, mainly because we have imprecise knowledge of ρ.
• Properties of equipotential surfaces (or geops) are:
- They never cross each other.
- They are never parallel to each other (i.e. dH≠constant).
- The distance between geopotential cannot be 0.
- They are closed surfaces within or exterior to a volume (i.e. no discontinuities).
- The local radii of curvature vary smoothly from point to point.
- The surfaces are complex everywhere.

Moving along a geop gives no change in potential. The total differential of the gravity potential is:

dW dW dW
dW = dX + dY + dZ
dX dY dZ
⎡ dX ⎤
⎡ dW dW dW ⎤ ⎢ ⎥
=⎢ dY (2.17)
⎣ dX dY dZ ⎥⎦ ⎢ ⎥
⎢⎣ dZ ⎥⎦
= ∇W ⋅ dx
= g ⋅ dx
where
g is the gravity force vector,
dx is a position vector, and
dW is the change in potential between two points separated by (dX, dY, dZ).

• If two points occupy the same equipotential surface then dW = 0 or g⋅dx = 0.

o If the scalar product of two vectors is 0, then vectors g and dx are orthogonal to each
other.
o This implies that the gravity force vector g at a point is normal to the equipotential surface
passing through the point.
Because of the curvature of geops, the lines that are normal to them also have slight curvature and
possibly some torsion (twist). These lines are the plumblines or lines of force.

2-6
 plumblin

geops

Eart

d
H

plumblines

geops

Mountain

Figure 2.3 – Equipotential Surfaces (geops) and Plumblines

The spacing of the equipotential surfaces is directly related to the magnitude of gravity
• i.e. the closer the surfaces are, the stronger the gravity field.
Change in g is nearly a difference between potential surfaces. Over infinitesimally small line elements
along a plumbline, |dx| ~ dH where dx is the distance along the plumbline and dH is the linear distance
between the start/end of the line element. If H is considered positive upwards and g positive downwards:
dW = g ⋅ dx
= g ⋅ dH ⋅ cos( g ⋅ dx)
14243 (2.18)
~180°

= − g ⋅ dH
Eq. (2.18) relates the change in potential to a change in orthometric height (measured by linear staff
increments). This equation is fundamental to the theory of height determination. Eq. (2.18) can be written
as:
dW
g=− (2.19)
dH

From this equation (and Fig. 2.4 below), gravity cannot be constant on the same equipotential surface.

2-7
 W + dW

dH1 dH2
W

g1 P2
P1 g2

dW dW
g1 = − ≠ g2 = −
dH1 dH 2

Fig.2.4 – Relationship between change in potential and change in height

2.2.3 Geoid, Geoid Height & Deflection of the Vertical

The geoid is the gravity equipotential surface that best approximates Mean Sea Level.
• The mean sea level is a real, global surface.
• The geoid ≡ MSL to a first order of ~1 – 2 m.
• It is important to remember that the actual ocean surface is not an equipotential surface due to the
presence of ocean currents, density differences between ocean basins (i.e. temperature and salinity),
tides, waves, etc.
• The geoid is a fundamental reference surface and is used as a datum for height values.
• The potential of the geoid is denoted by W0 (but Wgeoid ≠ 0).
• The distance measured along the plumbline from the geoid to a point is called the orthometric
height (H).

At any point, the position of a particular ellipsoid with respect to the geoid is defined by the distance
between the two surfaces, measured along the normal to the ellipsoid, and/or the relative slope of the two
surfaces (see Fig.2.5).

2-8
 θ Note: Values for either N or θ can be used to define the geoid surface.

geoid
N

g
ref. ellipsoid

Fig.2.5 - Fundamental Surfaces of Geodesy

N is called either the geoid (or geoidal) height, or the geoid-ellipsoid separation, or geoidal
undulation.

• If the reference ellipsoid is geocentric, then N is 'absolute', if the ellipsoid is regional (or non-
geocentric), then the N values are relative.
• N is taken to be positive above the ellipsoid with a range of values from ± 200 m.

The deflection of the vertical (θ) is the angle between the ellipsoid normal and the local vertical (gravity
vector or plumbline) at a point.
• The deflection of the vertical defines the relative slope of the two surfaces (geoid/ellipsoid).
• Typical values of θ are order
o 1" – 20" in most cases
o 60" –100" in mountainous areas, i.e. areas with high density (gravity) contrast.

Since the deflection of the vertical is a space angle we normally separate it into two orthogonal
components in the meridional (longitudinal) direction (ξ) and prime vertical (latitudinal) direction (η)

• θ = ξ 2 +η2
• The deflection angle along a line of ellipsoidal azimuth (α) is θα = ξ cos α + η sin α .
• The difference in N between adjacent stations on an ellipsoid (say A and B) can be calculated as:

ΔN ≅ −θα AB ⋅ s AB { for small angles} (2.20)

This is called astro-geodetic levelling.
o The first Australian geoid was computed this way
o The accuracy depends on precision of astronomical observations.

To get profiles of N, we can integrate the changes in ΔN along the route:

2-9
 B
N A − N B = − ∫ θ ds (2.21)
A

This technique can be used for local area surveys but today astro-geodetic techniques would not generally
be used. Instead either GPS levelling or models derived from satellite geopotential fields are used.

2.2.4 Normal Gravity Field

As a first approximation, the earth is considered an ellipsoid of revolution.
As it is difficult to directly compute the actual gravity field of earth, we try to determine the gravity field
of an ellipsoid, because the ellipsoid is easy to handle mathematically. This is called the “normal
gravity field”.
• The deviations of the actual gravity field from the ellipsoid "normal" field are small so that they can
be considered linear.
• This splitting of the earth's gravity field into a "normal" and a remaining, small "disturbing" field
considerably simplifies the problem of its determination.
• The problem could hardly be solved without introduction to the normal gravity field.

We therefore assume that the normal figure of the earth is a level ellipsoid, i.e. an ellipsoid of revolution
which is an equipotential surface of a normal gravity field.
Denoting the potential of the normal gravity field by:
U = U ( x, y , z ) (2.22)
The level ellipsoid is a surface with U = constant
• It exactly corresponds to the geoid, defined as a surface where W = constant.
• By defining the given ellipsoid is an equipotential surface of the normal gravity field, and by
prescribing the total mass M, we completely and uniquely determine the normal potential U as:
↓ centrifugal potential
U =V +Φ (2.23)
↑ gravitational potential
• The normal potential function U(x, y, z) is completely determined by a system of 4 constants:
- the shape of the ellipsoid (a, f)
- the total mass (M)
- the angular velocity (ω)
Note that other systems of four independent parameters may also be used as defining constants (e.g. U =
f(a, J2, kM, ω)).

2-10
The normal gravitational potential V, like the gravitation potential of the earth, can be given by a
spherical harmonic series. This harmonic expansion for V plus the centrifugal potential Φ leads to the
following expression for the potential of normal gravity:

kM ⎡ ∞
⎛a⎞
n
⎤ ω2 2 2
U (r ) = ⎢1 + ∑ ⎜ ⎟ Cn ,o Pn ,o (cosθ ) ⎥ + r sin θ (2.24)
r ⎢⎣ n = 2 ⎝ r ⎠ ⎥⎦ 2

where
r = distance from geocentre to point
a = semi-major axis of ellipsoid
n = degree of the spherical harmonic expansion (even number since symmetric mass
distribution)
Cn,o = harmonic coefficients (zonal harmonics)
Pn,o(cosθ) = Legendre polynomials
θ = co-latitude (i.e. 90 - ϕ)
ω2
r 2 sin 2 θ = the centrifugal potential expressed in spherical coordinates
2

The normal gravity vector (γ) is given by:

γ = gradU = ∇U (2.25)
• Differentiating equation (2.24) with respect to r, yields the magnitude of normal gravity.
• Normal gravity is usually denoted by γ and is a function of both distance from the centre of mass of
the earth and attitude.
• Normal gravity is rotationally symmetric and hence does not depend on longitude.

Normal gravity (γ0) on the level ellipsoid (i.e. h = 0)

It is given by the formula of Somigliana (1929):
a ⋅ γ e cos 2 ϕ + b ⋅ γ p sin 2 ϕ
γ0 = (2.26)
a 2 cos 2 ϕ + b 2 sin 2 ϕ
where
γe and γp are the normal gravity at the equator and poles respectively;
a and b are the semi-axes of the ellipsoid.
γe and γp can be obtained by closed formulae from the defining parameters of the level ellipsoid.

By defining the gravity flattening (β) as:

γ p −γe
β= (2.28)
γe

2-11
an expression for normal gravity on the ellipsoid is:
γ 0 = γ e (1 + β sin 2 ϕ + higher order terms ) (2.32)

An international gravity formula is:

γ 0 = γ e (1 + B sin 2 ϕ + C sin 4 ϕ + D sin 6 ϕ + E sin 8 ϕ (2.33)
This is related to an ellipsoid GRS80, for which:
a = 6378137 m
b = 6356752.3
f = (a-b)/a = 1/298.257 722 101
The constants are:

γ e = 9.780 326 7715 ms-2

B = 0.005279 044
C = 0.000 023 2718
D = 0.000 000 1262
E = 0.000 000 0007

Normal gravity (γ0) above the level ellipsoid (i.e. h ≠ 0)

The normal gravity at height h is:
⎛ 2 ⎞
γ h = γ 0 ⎜1 −
⎝ a
(1 + f + m − 2 f sin 2 ϕ ) h ⎟
⎠
(2.35)

2.2.5 Disturbing Potential and Gravity Anomaly

1. Disturbing Potential
The disturbing potential (T) is the difference between the actual gravity potential (W) of the earth and
the normal gravity potential (U) as
T=W–U (2.36)
T is a small but irregular potential.
T is a function of s.
ΔT = 0 in the exterior space.

2-12
 P

Geoid (W = W0)

N
gP = measured gravity

Q
Ellipsoid (U = W0)

γQ = gravity on ellipsoid

Fig. 2.6 –Geoid and reference ellipsoid

2. Gravity Anomaly
We assume that both surfaces have the same potential (as per definition), i.e.:

WP = W ( x, y, z ) P = W0
U Q = U ( x, y , z )Q = U 0
(2.37)
with
U 0 = W0 ( because equal mass of geoid and ellipsoid )

The difference in magnitude between the observed gravity on the geoid (gP) and the normal gravity on the
ellipsoid (γQ) is called the gravity anomaly, Δg, at P:
Δg = g P − γ Q (2.38)
The difference in directions between the observed gravity on the geoid (gP) and the normal gravity on the
ellipsoid (γQ) is called the deflection of the vertical, N.

The magnitude of gravity can be measured using gravimeters. Gravity networks are surveyed using both
absolute and relative gravimeters.
In order to make possible a comparison between gravity values, they need to be reduced to a common
surface (e.g. the geoid).
One of the most common reductions is a free-air reduction based on the vertical gradient of gravity
(assumed linear with height) between the observation point and the geoid.

2-13
The reduction is:
∂g ∂γ
δ g FA = − ⋅H ≈ − ⋅H
∂H ∂h
where
δ g FA is the free air reduction ( correction )
∂g (2.39)
is the vertical gravity gradient
∂H
∂γ
is the normal vertical gravity gradient
∂h
H is the orth om etric height

∂γ
For most purposes, can be taken as +0.3086 mgal/m. The free-air gravity anomaly then becomes:
∂h
Δg FA = g p − (γ Q − 0.3086 ⋅ H ) (2.40)

The free-air reduction corresponds to a parallel displacement of the topographic masses by the height of
the computation point and provides good approximate boundary values on the geoid.

3. Geoid
The geoid undulation N can be related to the disturbing potential by means of Brun’s Formula:
T
N=
γ

With reference to Fig.2.6, we can also show that:

⎛ ∂U ⎞
U P = UQ + ⎜ ⎟ N = UQ − γ N
⎝ ∂N ⎠
∂U
where is the partial differentiation in the (2.41)
∂N
direction of the ellipsoid normal ,
and is equal to γ
Also

WP = U P + TP
= U Q − γ N + Tp
because WP = U Q = W0 , we get :
(2.42)
T
T =γ N or N =
γ
• Eq. (2.42) is a very important equation in physical geodesy.
• Eq. (2.42) relates the disturbing potential (T) to the geoidal undulation (N) and the normal gravity
(γ).

2-14
• Eq. (2.42) is called Brun's Formula.

4. Spherical Approximation of T
From Eq. (2.42), to obtain the geoid height N or the deflection of the vertical of a point on the geoid, we
need to know the disturbing potential (T) at that point.
The disturbing potential is resolved using the relationship between Δg and T.

Using a spherical approximation, the gravity anomaly can be then written as:
⎛ ∂T ⎞ ⎛ 2 ⎞
Δg = − ⎜ ⎟ − ⎜ ⎟T (2.42)
⎝ ∂r ⎠ ⎝ r ⎠

where r = 6371 km is the mean radius of the Earth.

• Eq. (2.42) is called the fundamental equation of physical geodesy.
• Its importance lies in relating a quantity that can be measured (Δg) to the unknown disturbing
potential (T).
• It is now a differential equation which can be solved.
• T cannot be observed, but can be computed using Eq.(2.42), which is the geodetic boundary value
problem.

2.2.6 Geodetic Boundary Value Problems (GBVP)

The GBVP consists of the determination of
• the physical surface of the earth or the geoid, and
• the external gravity field from the gravity potential (W) and the gravity acceleration (g) on the
earth's surface.

The normal gravity potential (U) can be calculated in terms of spherical harmonics (e.g. Eq. 2.24).
The problem then is to find T given values of a linear combination of T and its derivative normal to the
surface.
This type of problem is termed the third boundary value problem of potential theory.
If T can be determined, then any other of its functionals (i.e. N, ξ, η, Δg) can be calculated and vice-
versa.
There are many solutions of the GBVP for the geoid such as:
- The classical solution
- Spherical harmonic solutions
- Molodenskii's solution
- Statistical solutions
- Combined solutions
Spherical harmonic solutions of the geoid are very commonly used (e.g. GEM, JGM, OSU) with the
model coefficients being easily obtained from satellite orbit analyses, satellite altimetry and gravimetric
terrestrial data.

2-15
 3 Coordinate Systems

Appropriate, well-defined and reproducible reference coordinate systems are essential

for such things as the description of satellite motion, the modelling of our
measurements, the representation of position and interpretation of results on the Earth’s
surface. The increasing accuracy of many satellite observation techniques requires a
corresponding increase in the accuracy of the reference systems due to the time-varying
nature of the geodynamic processes influencing our measurements.
There are a number of different coordinate systems which are typically used in geodetic
and mapping-type applications. The use of different systems creates a need for accurate
transformation techniques so that coordinates may be converted from one system to
another, thereby allowing information to be easily exchanged between applications.

3.1 Basic Coordinate Systems and Transformations

3.1.1 Cartesian Coordinate Systems

In a Cartesian coordinate system, with axes (x, y, z), the position of a point P is
determined by its position vector:

z
⎡ xp ⎤ P
⎢ ⎥
X P = ⎢ yp ⎥ Xp
zp
⎢zp ⎥
⎣ ⎦
y
xp
yp
x

Figure 3.1 – Cartesian Coordinate System

3.1.2 Cylindrical Coordinate Systems

In a cylindrical coordinate system, with axes ( r , λ , z ), the position of a point P is

determined by its position vector:

⎡ rp ⎤ z
⎢ ⎥ P
X P = ⎢λ p ⎥
⎢ zp ⎥ Xp
⎣ ⎦
o zp
rp
r λp

Figure 3.2 – Cylindrical Coordinate System

3-1
A coordinate surface is a surface along which one of the three coordinates is constant.
As can be seen in Figure 3.3, if we hold r = constant and let λ and z vary, the locus of
P(r, λ, z) is then a right circular cylinder of radius r and axis along OZ. The locus of λ =
constant is a vertical plane containing the z-axis and making an angle λ with the xz
plane. Similarly, if z = constant, we create a horizontal plane perpendicular to the z-axis.

Figure 3.3 – Cylindrical Coordinate Surfaces

3.1.3 Spherical Coordinate Systems

In a spherical coordinate system, with axes ( r , φ , λ ), the position of a point P is

determined by its position vector:

⎡ rp ⎤
⎢ ⎥
X P = ⎢φ p ⎥ P
⎢λ p ⎥
⎣ ⎦ rp

o φp
λp

Figure 3.4 – Spherical Coordinate System

In this type of system, every point in the whole space can be given spherical coordinates
restricted to the ranges:
r ≥ 0, -π/2 ≤ φ ≤ π/2, 0 ≤ λ ≤ 2π

3-2
Because of the analogy between the surface of a sphere and the Earth’s surface, the
vertical axis of this system is sometime called the polar or rotation axis, while φ is
commonly referred to as latitude and λ as longitude. It is also common to discuss a
spherical coordinate system in terms of meridians, parallels, and hemispheres. When
looking at coordinate surfaces in a spherical system (see Figure 3.5), we can see that
holding r = constant and letting λ and φ vary creates a sphere of radius r. If φ = constant,
the surface created is a cone, the height of which is determined by the range of r values.
The locus of λ = constant is again a vertical plane containing the vertical-axis and
making an angle λ with the horizontal axis.

Figure 3.5 - Spherical Coordinate Surfaces

3.1.4 Transformations Between Basic Coordinate Systems

In practice it is common to obtain coordinates that are stated in terms of one coordinate
system and for some reason have a need to transform these coordinates into another
system more useful to your application. Thus understanding and utilising coordinate
transformations is a basic and essential skill of geodetic studies.
3.1.4.1 Cartesian to Cartesian
A transformation from one Cartesian coordinate system (x,y,z) to another (x',y',z') -
when both systems have the same origin and z-axis (this is purely for convenience) but
different x and y axes - simply requires a rotation of angle γ around the z axis. This
transformation can be realised through the standard matrix operation:

X p' = Rz (γ ) Xp
frame 2 rotation frame1 (3.1)
coords matrix coords

3-3
with the rotation matrix about the z axis being given by:

⎡ cos γ sin γ 0⎤
Rz (γ ) = ⎢⎢ − sin γ cos γ 0 ⎥⎥ (3.2)
⎢⎣ 0 0 1 ⎥⎦

x 'p = x p cos γ + y p sin γ

ie. y 'p = − x p sin γ + y p cos γ
z 'p = z p

Derivation:

y'
ypsinγ

γ
yp y

xp' γ
xpcosγ
xp

yp P
γ

yp' xpsinγ

ypcosγ
x x'

3-4
The equivalent rotation matrices (required if z axes not aligned) about the x axis (Rx)
and y axis (Ry) can similarly be derived:

⎡1 0 0 ⎤
Rx (α ) = ⎢⎢ 0 cos α sin α ⎥⎥
⎢⎣ 0 − sin α cos α ⎥⎦
(3.3)
⎡cos β 0 − sin β ⎤
Ry ( β ) = ⎢⎢ 0 1 0 ⎥⎥
⎢⎣ sin β 0 cos β ⎥⎦

This representation of the Cartesian to Cartesian transformation is valid for a right-

handed coordinate system. Also, in these equations a counter-clockwise rotation (when
viewed towards the origin from the axis) is taken as positive.
When dealing with Cartesian systems, any required coordinate transformation can be
realised through combination of the above rotations. The complete transformation from
one system (X) to another system (X') can therefore be written as:

X p' = Rx (α ) Ry ( β ) Rz (γ ) X p = R X p (3.4)

where the combined matrix, R, for rotating the coordinates of XP by angles (α, β, γ) is:

⎡ cos β cos γ cos β sin γ − sin β ⎤

⎢
R = ⎢sin α sin β cos γ − cos α sin γ sin α sin β sin γ + cos α cos γ sin α cos β ⎥⎥
⎢⎣cos α sin β cos γ + sin α sin γ cos α sin β sin γ − sin α cos γ cos α cos β ⎥⎦

Therefore, the forwards and backwards transformation between the position vectors XP
and X P' in two arbitrarily rotated Cartesian coordinate systems can be written as:

X p' = R X p and X p = RT X p' (3.5)

The mathematical properties of rotation matrices apply to all of the above equations,
with the following rules being important:
(a) Rotation does not change the length of the position vector, i.e., scale is invariant
to rotation.
(b) Matrix multiplication is not commutative, Ri ( μ ) R j (υ ) ≠ Ri (υ ) R j ( μ )
(c) Matrix multiplication is associative,
Ri ( μ ) { R j (θ ) Rk (υ )} = { Ri ( μ ) R j (θ )} Rk (υ )
(d) Rotation around the same axes are additive, Ri ( μ ) Ri (υ ) = Ri ( μ + υ )
(e) The inverse and transpose matrices are equal since they are orthogonal
Ri−1 ( μ ) = RiT ( μ ) = Ri (− μ )
and also
( Ri ( μ ) R j (υ )) −1 = R −j 1 (υ ) Ri−1 ( μ )

3-5
In geodesy, the rotation angles between systems are most often small and hence we can
use the linearised form of the rotation matrix (R) in our computations. For small rotation
angles (θ), we can make the assumptions that cosθ ~ 1 and sinθ ~ θ (where the angle θ
is in radians), if we also neglect higher order terms (i.e., products of rotation angles) we
can therefore write the rotation matrix for very small angles as:

⎡1 γ −β ⎤
R = ⎢⎢ −γ 1 α ⎥⎥ (3.6)
⎢⎣ β −α 1 ⎥⎦

Another useful transformation allows us to change the orientation of the coordinate axes
in our Cartesian system (i.e., we can change from a right-handed coordinate system to a
left-handed coordinate system or vice versa). This can be achieved using the following
reflection matrices:

⎡ −1 0 0 ⎤ ⎡1 0 0 ⎤ ⎡1 0 0 ⎤
S x = ⎢⎢ 0 1 0 ⎥⎥ S y = ⎢⎢ 0 −1 0 ⎥⎥ S z = ⎢⎢0 1 0 ⎥⎥ (3.7)
⎣⎢ 0 0 1 ⎦⎥ ⎣⎢ 0 0 1 ⎦⎥ ⎣⎢0 0 −1⎦⎥

3.1.4.2 Cylindrical to Cartesian

The transformation of coordinates between these two systems is quite straight forward.
In this case, the z axes and origin of the systems are already co-incident and the x and y
axes are simply defined at λ = 0° and 90° respectively (see Figure 3.6).

z
P

Xp zp

o y
rp
xp λp
x
yp

Figure 3.6 - Transformation from Cylindrical to Cartesian Coordinates

The equations for the transformation from cylindrical (polar) to Cartesian (rectangular)
and vice versa are therefore:

⎡ x2 + y 2 ⎤
⎡ x ⎤ ⎡ r cos λ ⎤ ⎡r ⎤ ⎢ ⎥
⎢ y ⎥ = ⎢ r sin λ ⎥ and ⎢λ ⎥ = ⎢ tan −1 y ⎥ (3.8)
⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ x ⎥
⎢⎣ z ⎥⎦ ⎢⎣ z ⎥⎦ ⎢⎣ z ⎥⎦ ⎢ ⎥
⎢⎣ z ⎥⎦

3-6
3.1.4.3 Spherical to Cartesian
In the transformation between spherical and Cartesian coordinates we align the z-axis
with the rotational/polar axis of the spherical system and again define the x and y axes
from the origin of the spherical system at λ = 0° and 90° respectively (see Figure 3.7)

z
P
rp
zp
φp y
o
xp λp
x
yp

Figure 3.7 - Transformation from Spherical to Cartesian Coordinates

The equations for the transformation from spherical to Cartesian (rectangular) and vice
versa are therefore:

⎡ x ⎤ ⎡ r cos φ cos λ ⎤
⎢ y ⎥ = ⎢ r cos φ sin λ ⎥
⎢ ⎥ ⎢ ⎥
⎢⎣ z ⎥⎦ ⎢⎣ r sin φ ⎥⎦

and (3.9)

r = x2 + y 2 + z 2
⎛ y⎞
λ = tan −1 ⎜ ⎟
⎝ ⎠ x
⎛ z ⎞ ⎛z⎞
φ = tan −1 ⎜ ⎟ or sin −1 ⎜ ⎟
⎜ 2 2 ⎟
x +y ⎠ ⎝r⎠
⎝

3-7
3.2 Ellipsoidal Coordinate Systems
The earth’s surface may be closely approximated by a rotational, oblate ellipsoid (i.e.
flattened poles). As a result, many geodetic calculations are done using geometrically-
defined ellipsoidal systems. It is therefore important to have a full understanding of the
definition and use of these systems.

3.2.1 Geometric Parameters of the Ellipse

The rotational ellipsoid is created by rotating the meridian ellipse about its minor axis.
The shape of the ellipsoid is generally defined by two parameters, the semi-major axis
(a) and semi-minor axis (b).

P
meridian ellipse
Z

O
P' O
F1 F2 Y

reference cartesian axes

Figure 3.8 – Geometry of the Ellipse

where: F1 and F2 are the focii of the ellipse (i.e. sum of the distances from F1 and F2 to
any point on the ellipse is constant).
OP' = a
OP = b
OF2 = E = OF1
F1P' = a(1 – e)
reference X , Y and Z Cartesian axes are defined with respect to the origin
and axes of the ellipse (see inset)

and the equation of the ellipsoid is written as

x2 + y2 z 2
+ 2 =1 (3.10)
a2 b

3-8
Generally, the parameter b is replaced by one of a number of smaller quantities that are
more suitable for mathematical calculations, i.e. they lend themselves to series
expansion. These quantities are:
(i) the geometric flattening (f)
a −b
f = or b = a (1 − f )
a
(ii) the linear eccentricity (E)

E = a 2 − b2
(iii) the first eccentricity (e)

E a 2 − b2
e= =
a a2
(iv) the second eccentricity (e')

E a 2 − b2 e2
( e ')
2
e' = = or =
b b2 1 − e2
Other useful relations between these quantities are
b 1
= 1 − f = 1 − e2 =
1 + ( e ')
a 2

and
e2 = 2 f − f 2

Generally you will find that the shape of an ellipsoid will be defined using the a and f
parameters, although theoretically f must be derived from other fundamental constants.
Typical values for the above parameters when referring to an earth based reference
system are (e.g. Geodetic Reference System (GRS 80) values):
a = 6 378 137 m
b = 6 356 752.3141 m
E = 521 854.0097 m
f-1 = 298.257 222 101
e2 = 0.006 694 380 022 90
(e')2 = 0.006 739 496 775 48

3-9
3.2.2 Coordinate Definition on the Ellipsoid

A system of Cartesian coordinates ( X , Y , Z ) which provide one means of defining

points in the ellipsoidal system is located with its origin at the centre of the ellipsoid and
the Z axis coincident with the minor axis of the ellipsoid. The X axis is usually
directed to the ellipsoidal zero meridian and the Y axis then completes the right-handed
system.
The second means of defining coordinates in an ellipsoidal system is to use the
geographic (or geodetic or curvilinear) coordinates of ellipsoidal latitude (ϕ), ellipsoidal
longitude (λ) and ellipsoidal height (h). These coordinates are ellipsoid dependent and
will vary with changes in the ellipsoidal definition.
Figure 3.9 illustrates the relationship between the rectangular and curvilinear ellipsoidal
coordinate systems. The transformation between these coordinate systems will be
discussed further later.

Z P

b P'

ϕ
O
a Y
λ

Figure 3.9 – Ellipsoidal Coordinate Systems

It is also sometimes useful to define coordinates in terms of the radius of the circle of
latitude, p , where:

Z P where:
h'p h
P' p = X P2 + YP2
p
X P = p cos λ
YP = p sin λ
Z p' Z P = Z P' + hP'

Y
Figure 3.10 – Definition of Ellipsoidal Coordinates

3-10
3.2.3 Types of Latitude

earth surface

Z P

sphere (radius a)
h'p P'' h

P'
ellipsoid
p
a

r
b a
dp
Z p' tangent at P'
dz
ϕp
βp ϕp 90-ϕp

O
Y (X )
a

ϕ p = geodetic (ellipsoidal) latitude

ϕ p = geocentric latitude geodetic coords of P (on surface)
β p = reduced (spherical) latitude (ϕ , λ , h )
a = semi-major axis
r = geocentric radius geodetic coords of P' (on ellipsoid)
h = ellipsoidal height (ϕ , λ )

Figure 3.11 – Definition of Latitude

Using the geocentric latitude ϕ p and geocentric radius r, we can define the quantities

p = r cos ϕ p
Z p' = r sin ϕ p (3.11)

r= (p 2
+ Z p'
2
)
and the equation of the ellipse is
' 2
p2 Z p
+ =1 (3.12)
a 2 b2

3-11
 b
We can show that, since the ratio of ellipsoidal coordinates to circular coordinates is ,
a
we have
p = a cos β p
⎧ b⎫ (3.13)
Z p' = b sin β p ⎨=a ⋅ sinβ p ⋅ ⎬
⎩ a⎭

The slope of the ellipsoidal tangent at P' is

∂Z p'
= − cot ϕ p (3.14)
∂p
Also, from circular coordinates

∂p ∂Z p'
= − a sin β p and = b cos β p
∂β p ∂β p
(3.15)
∂Z p' ∂Z p' ∂β p b cos β p ⎛b⎞
∴ = ⋅ = = − ⎜ ⎟ cot β p
∂p ∂β p ∂p − a sin β p ⎝a⎠

Thus, equating Error! Reference source not found. and

Error! Reference source not found. gives

⎛b⎞
tan β p = ⎜ ⎟ tan ϕ p or tan β p = (1 − f ) tan ϕ p
⎝a⎠ (3.16)
( reduced β ) p ( geodetic ϕ ) p

From equations (3.13) and (3.16) and from Figure 3.11 we can see that that relationship
between geocentric and geodetic latitude can be expressed as:
Z p' b sin β p
(geocentric) tan ϕ p = =
p a cos β p
⎛b⎞
= ⎜ ⎟ tan β p (3.17)
⎝a⎠
2
⎛b⎞
= ⎜ ⎟ tan ϕ p or
⎝a⎠
(1- e ) tan ϕ ( geodetic )
2
p

Using series expansions for the differences between the various latitude parameters and
relations above we can get approximate equations such as:
e2 e4 e6
ϕ −ϕ = sin 2ϕ + sin 2 ϕ sin 2ϕ + sin 4 ϕ sin 2ϕ + K
2 2 2 (3.18)
= 2 (ϕ − β )

3-12
Numerical Example

if f −1 ≈ 300
then e2 = 2 f − f 2
2
2 ⎛ 1 ⎞
≈ −⎜ ⎟
300 ⎝ 300 ⎠
1
≈
150

e2
∴ ϕ −ϕ ≈ sin 2ϕ
2
1
= sin 2ϕ
300

the maximum difference between latitude types occurs at

sin 2ϕ = 1 i.e. at ϕ = 45°

∴ max ϕ − ϕ ≈ 11'30"
∴ max ϕ − β ≈ 5'45"

3.2.4 Radii of Curvature on the Ellipsoid

Curvature is the change in the direction of the normal (or tangent) to a curve (dθ) over a
distance along the curve (ds).

1 −∂θ
ds curvature = ≈ (3.19)
dθ R ∂s
R

where R = radius of the curve

dθ

Figure 3.12 – Curvature

3-13
For an ellipsoid, the principal radii of curvature are:
• in the meridian plane (i.e. a north-south section) denoted as ρ
• in the prime vertical plane (i.e. an east-west section) denoted as ν

meridian
P

prime vertical

Figure 3.13 – The Meridian and Prime Vertical Planes at Point P on an Ellipsoid
(sourced from http://www.mentorsoftwareinc.com/CC/gistips/TIPS0899.HTM)

The radius of curvature of a plane curve whose equation, in rectangular coordinates, is

z = f(p) is given by
3
⎡ ⎛ ∂z ⎞ 2 ⎤ 2

⎢1 + ⎜ ⎟ ⎥
⎢ ⎝ ∂p ⎠ ⎥⎦
ρ=⎣ (3.20)
⎛ ∂2 z ⎞
⎜ 2⎟
⎝ ∂p ⎠

In the case of a meridional section, the radius of curvature is obtained using expressions
2
⎛b⎞ ⎡ ⎛ ∂z ⎞ ⎤
− ⎜ ⎟ ⎢ z − p ⎜ ⎟⎥
⎝ ∂p ⎠ ⎦
2
∂ z ⎝a⎠ ⎣
=
∂p 2 z2
(3.21)
2
∂z ⎛b⎞ p
= −⎜ ⎟ = − cot ϕ
∂p ⎝a⎠ z

Substitution of equations (3.21) into (3.20) yields an expression for ρ.

The radius of curvature in the prime vertical at a point P is psecϕ and substitution of
this into the above equations gives an expression for ν.

a (1 − e 2 ) a
ρ= ν=
(1 − e sin 2 ϕ ) (1 − e sin 2 ϕ )
3 1
2 2 2 2
(3.22)

( radius to centre of curvature ) ( radius to rotation axis )

3-14
Using typical ellipsoidal parameters of a = 6378136.3 m and f-1 = 298.257 m at
latitude ϕ = 45° gives the following radii of curvature values:
ρ = 6,367,381.109 m and ν = 6,388,837.597 m (difference ~ 21 km)

From the above equations and example we can see that ν is greater than ρ except at the
poles (i.e. ϕ = 90°) when ν and ρ are equal:

ν = ρ = a (1 − e 2 )
−1
2

At the equator (i.e. ϕ = 0°) we have:

ν =a and ρ = a (1 − e 2 )
The mean curvature at any point on the ellipsoid is given as

1⎛ 1 1⎞
⎜ + ⎟
2 ⎜⎝ ρ ν ⎟⎠

Any plane which contains the ellipsoidal normal at a point P intersects the ellipsoid in a
normal section . The radius of curvature of an arbitrary normal section with an azimuth
of α' is computed according to Euler’s formula by:

1 cos 2 α sin 2 α
= + (curvature)
Rα ' ρ ν
or (3.23)

ρν
Rα ' = (radius)
ν cos α + ρ sin 2 α
2

3.2.5 Radii of Spherical Approximation

Various spherical approximations may be made to simplify calculations on the ellipsoid.
The determination of the radius of the chosen spheroids is dependent on the type of
approximation used, the most common are:
• Mean Radius (Gauss Radius)

ρ ⋅ν = a (1 − e 2 ) 2 (1 − e 2 sin 2 ϕ )
1
RM =
• Radius of a sphere having the mean radius of the biaxial ellipsoid
2a + b
RS =
3
• Radius of a sphere having the same area as the ellipsoid
Area
RA = (Area = area of ellipsoid)
4a
• Radius of a sphere having the same volume as the ellipsoid

(
RV = a 1 − e 2 ) 1
6

3-15
3.2.6 Arc Lengths and Areas on the Surface of an Ellipsoid
The arc element of a meridian arc can be written as
S meridian = ρ ∂ϕ (3.24)

In order to find the length of arc between two latitudes ϕ1 and ϕ2, we integrate the above
equation
∂ϕ
S meridian = ∫ ρ ∂ϕ = a (1 − e 2 ) ∫
ϕ2 ϕ2
(3.25)
(1 − e2 sin 2 ϕ )
ϕ1 ϕ1 3
2

This representation is complex since the form of the integral is one of an elliptic integral
of the second kind. The computations may be achieved through numerical integration.
Alternatively, the denominator can be expanded by binomial expansion and
subsequently integrated term by term. The series method leads to the formula for the
meridian arc length (correct to 0.5 mm at latitude 45°) as:
S meridian = a ( Ao (ϕ 2 − ϕ1 ) − A2 ( sin 2ϕ 2 − sin 2ϕ1 ) + A4 ( sin 4ϕ 2 − sin 4ϕ1 ) − A6 ( sin 6ϕ 2 − sin 6ϕ1 ) + K) (3.26)
where a is semi-major axis of the ellipsoid
ϕ1 and ϕ2 are the geodetic latitudes (in radians)
A0, … A6 are constants given as:

e2 3e 4 5e 6 3⎛ e 4 15e6 ⎞
A0 = 1 − − − A2 = ⎜ e 2 + + ⎟
4 64 256 8⎝ 4 128 ⎠
15 ⎛ 4 3e 6 ⎞ 35
A4 = ⎜e + ⎟ A6 = e6
256 ⎝ 4 ⎠ 3072

The arc length of a parallel on the ellipse is

S parallel = ν cos ϕ dλ (3.27)

and the length of the arc of a circle of latitude between geodetic longitudes λ1 and λ2 is
given by
λ2
S parallel = ν cos ϕ ∫ ∂λ
λ1 (3.28)
= ν cos ϕ ( λ2 − λ1 )
Thus, for an area in an ellipsoid between latitudes ϕ1 and ϕ2 and longitudes λ1 and λ2 we
can write:
ϕ2 λ2
Area = ∫ ∫λ ρν cos ϕ ∂ϕ ∂λ (3.29)
ϕ1 1

3-16
 4 Geodetic Reference Systems

Reference coordinate systems used in geodesy are typically global and geocentric (i.e.,
the origin of the coordinate system is situated at the Earth’s centre of mass, called the
geocentre). One reason why the coordinate system is geocentric is that the artificial
satellite’s motion refers to the centre of mass of the Earth. In contrast, terrestrial
measurements are usually made over a small region on the Earth and are hence best
described in local reference coordinate systems. The relationship between both systems
must be known with sufficient accuracy for use in geodetic activities. Since the relative
position and orientation of reference frames changes with time, appropriate modelling
must be done. The establishment of precise transformation formulae between various
types of reference systems is one of the most important tasks of geodesy.
There is a distinction made between the terms ‘reference frame’ and ‘reference system’.
A reference system embodies the conceptual idea of a reference system definition,
including the fundamental theory and standards, whereas a reference frame relates to the
practical realisation of a reference system through measurements and the adoption of a
particular set of station coordinates. An example of this is the current Conventional
Terrestrial Reference System standards (see http://www.iers.org/iers/products/itrs/) and
the latest International Terrestrial Reference Frame (ITRF2005), defined by a particular
set of station coordinates and velocities given at a fixed epoch. Note that any set of
ITRF coordinates will have an associated time or epoch at which the reference frame is
fixed.
The positions of points on the Earth’s surface are usually given either in Cartesian,
ellipsoidal (geodetic) or natural (geographical) coordinates. Spherical and cylindrical
coordinates are also sometimes used, particularly with astronomy and space techniques.
Over small regions, local topocentric coordinates are typically used. The relationships
between each of these coordinate systems are very important in geodetic study and need
to be understood.

4.1 Geodetic Surfaces

Before we discuss the various reference systems used in Geodesy it is important to
understand the various geodetic surfaces and their relationships in regard to determining
a three-dimensional position for a point of interest.
The first point to note is that, for geodetic positioning, the physical shape of the real
earth is closely approximated by the mathematical surface of the rotational ellipsoid.
The ellipsoidal surface provides a smooth and convenient mathematical reference
surface for horizontal coordinates in geodetic networks.
On the other hand, the ellipsoid is much less suitable as a reference surface for vertical
coordinates (heights). Instead the geoid is used, where the geoid is defined as 'that level
surface of the Earth’s gravity field that best fits mean sea level on a global basis'. The
geoid is an equipotential surface (i.e. the gravity vector is perpendicular to the geoid at
all points on a potential surface) that is a function of the lateral density variations within
the Earth. As such, there are an infinite number of potential surfaces, with the geoid
being only one choice. However, the geoid is an ‘imaginary’ surface and its physical
realisation is conventionally defined in relation to MSL, as the undisturbed ocean
surface closely approximates (within about 1-2m) a potential surface.
The relationship between the geoid and ellipsoid is shown in Figure 4.1. The vertical
separation between the geoid and a particular reference ellipsoid is called the geoid
undulation or geoid height (N). The numerical value of the geoid height depends on the
particular ellipsoid chosen. For a global reference ellipsoid, typical N values are up to

1
+/-100 m. The geometrical relation between N, the height above the ellipsoid (h) and
the orthometric height (H) (mostly determined from spirit levelling) at a given point, is:
h=H+N (4.1)

geoid

ellipsoid

deflection of the vertical

θ

zenith (direction of gravity)

H ground surface
h

N
ellipsoid

ellipsoidal normal

geoid

Figure 4.1 – Relationships Between Geodetic Surfaces

It can be seen that N must be known when observations from satellite geodesy (h) and
from terrestrial geodesy (H) are used in a combined adjustment. h is derived primarily
from geometric considerations whereas H is defined in the Earth’s gravity field. The
details of the determination of height was discussed in Chapter 2.
The angle θ between the directions of the ellipsoidal normal and of the plumbline
(direction of gravity or zenith) at point P is called the deflection of the vertical (see
Figure 4.1). Usually θ is divided into two components, defined as:
ξ=Φ-ϕ meridian component (N-S)
η = (Λ - λ) cos ϕ prime vertical component (E-W)
where (Φ, Λ) are astronomical latitude and longitude values (i.e. relative to zenith) and
(ϕ, λ) are geodetic latitude and longitude (i.e. relative to ellipsoid).

Figure 4.1

Sign convention
for deflection of
vertical
components.

2
3.1 Overview of Geodetic Reference Systems
(i) Global Ellipsoidal – when an ellipsoid is centred on the earth's geocentre
and is designed to 'fit' the earth’s shape it may be described as a global
ellipsoidal (or geodetic) reference system. In this system we usually define
the position of terrestrial points in terms of curvilinear coordinates (ϕ, λ, h),
however a Cartesian representation (X,Y,Z) is also often used. There are
many different global ellipsoids which are used in geodesy, each with
slightly different values for the a and f parameters required to define the size
and shape of the reference surface. Therefore, when stating coordinates in a
particular global ellipsoidal datum, it is also necessary to identify the name
of the datum and/or the parameters which define it.
(ii) Global Cartesian – global Cartesian reference systems use rectangular
coords (X, Y, Z) to define position and are designed such that the:
• origin is at the geocentre (earth’s centre of mass)
• Z-axis directed towards the conventional definition of the north pole
(NP), or more exactly, towards the conventional terrestrial pole (CTP)
as defined by the International Earth Rotation Service (IERS).
• X-axis passes through a point of zero longitude (the Greenwich
meridian) and latitude.
• Y-axis forms the orthogonal triad for right-handed coordinate system.
(iii) Regional Ellipsoidal – An ellipsoidal reference system that is not geocentric
is described as a regional system. These systems are usually defined to
provide a 'best fit' to a region of interest on the earth’s surface (e.g. for an
individual country). This type of ellipsoidal datum is defined by fixing eight
parameters – two to specify the dimensions of the ellipsoid (a, f); three to
specify the location of its centre (fixing ϕ, λ, h of point at centre of
network); three to specify the orientation of the ellipsoid (fixing an azimuth,
deflection of vertical components at the origin point).
(iv) Regional Cartesian (geodetic) - a Cartesian reference system that is non-
geocentric and non-topocentric (see below), represents a regional coordinate
system (u, v, w). A regional Cartesian datum will be defined with respect to
the regional reference ellipsoid from which it was derived, i.e.:
• origin is at the centre of the reference ellipsoid used for defining the
datum in question.
• w-axis coincides with the semi major axis b of the ellipsoid (whether
the w-axis is also parallel to the CTP-axis is a matter of choice when
defining the datum)
• u-axis lies in the geodetic meridian λ = 0° and equatorial plane ϕ = 0°
of the ellipsoid
• v-axis completes the right-handed system.
(v) Local Cartesian – this type of frame uses coordinates (x', y', z') and may be
defined as follows
• origin at any point P (may or may not be on the earth’s surface)
• the x', y', z' axes are parallel to any geocentric terrestrial frame x, y, z.
(vi) Local Astronomical - relative positioning of points can use local
astronomical frames (LA). These systems are Cartesian frames with their
origin at an observation point, P on the earth’s surface. This local coordinate

3
 system is also called a local topocentric frame. Local astronomical frames
have the following definition:
• origin at any point P on the surface of the earth
• z' axis is directed to the astronomical zenith (i.e. local vertical or
‘plumbline’)
• x' axis is directed to the astronomical north meridian
• y' axis is directed to the east (i.e. perpendicular to x')
This definition indicates that a local astronomical system is therefore a left-
handed system and is not ellipsoid dependent.
(vii) Local Geodetic - relative positioning of points may also use local geodetic
(LG) frames as a reference. In these frames, positions are usually stated in
terms of east, north, up (or height (h)) coordinates relative to a local origin at
P. Local geodetic reference frames are defined as follows:
• origin at any point P (ϕ, λ, h) referred to a given ellipsoid
• u (or h) axis is normal through P – i.e. contains the ellipsoidal normal
through the point P. Convention is +ve in the outward or up direction.
• e-axis is normal to u and to the geodetic meridian plane, +ve east in
direction of increasing λ.
• n-axis is perpendicular to e and u axes, +ve north in the direction of
increasing φ.
The (e, n, u) system is therefore also a left-handed, local topocentric frame
however, LG frames are ellipsoid dependent and the orientation of the
frame changes depending on the geodetic position of P.

4.2 Inertial and Terrestrial Reference Coordinate Systems

To define relative positions of points on the Earth’s surface and to establish the motion
of a planet as a whole in space, it is necessary to introduce two fundamental coordinate
reference systems. These reference systems are:
- A space-fixed, inertial reference system for the description of satellite
motion, termed Conventional Inertial System (CIS) or Conventional
Celestial Reference System (CCRS).
- An earth-fixed, terrestrial reference system for the positions of the
observation stations, termed Conventional Terrestrial Reference System
(CTS or CTRS) and also called Earth Centred Earth Fixed System
(ECEF).
Either reference system can be defined by a set of three-dimensional coordinates
associated with a polyhedron of points or a set of curvilinear/angular coordinates.

4.2.1 Inertial frame

Newton’s laws of motion are only valid in an inertial reference system, i.e., a coordinate
system at rest or in a state of uniform rectilinear motion without any acceleration. The
theory of motion for artificial satellites is developed with respect to such a system. The
equatorial system, at a given epoch to yields a rather good approximation to an inertial
reference system.
The coordinates of points in an inertial frame correspond (usually) to stellar sources
whose positions are assumed to be fixed in space. At the present time, these positions
are realised through a catalogue of positions and proper motions for a given number of

4
fundamental stars (FK5 catalogue) in addition to a system of defined astronomical
constants – listed by the International Earth Rotation Service (IERS).
The origin of the inertial system is supposed to coincide with the geocentre (M). The
positive Z axis is oriented towards the north pole and the positive X axis to the First
Point of Aries, γ. The Y axis completes a right-handed system. The geocentre actually
undergoes small accelerations because of the annual motion of the earth around the sun
and so the reference system is sometimes termed quasi-inertial.
The realisation of the celestial reference system (ICRS) is the IERS Celestial Reference
Frame (ICRF) (see http://www.iers.org/iers/products/icrs/). The ICRF is based on the
coordinates of extragalactic objects determined by VLBI and is a realisation of a system
of directions which are consistent with those of the FK5. The origin of right ascension
of the frame (γ) is implicitly defined by adoption of the right ascensions of 23 radio
sources in catalogues obtained by various geodetic agencies, using a fixed right
ascension value for 3C273B at its FK5 value. The IERS origin for right ascension and
the FK5 are consistent within the uncertainties of the latter (~0.05"). The ICRF polar
axis points in the direction of the mean pole at J2000.0 (standard Julian Date epoch at
Yr 2000.0) as defined by the International Astronomical Union (IAU) conventional
models for precession and nutation. The FK5 system related to J2000, has been used as
the CIS since January 1, 1988. Improvements in the accuracy of the CIS are expected
with the astrometric satellite missions VSOP and HIPPARCOS. The angular
momentum axis serves as the reference Z axis and its free position in space is the
Celestial Ephemeris Pole (CEP).

4.2.2 Terrestrial frame

A suitable earth-fixed reference system must be connected in a well-defined way with
the earth’s crust. Such a CTRS can be realised through a set of Cartesian coordinates of
fundamental stations within a global network. Currently it is established through the
conventional direction to the mean orientation of the polar axis (or mean position of
earth’s rotational axis), called the Conventional Terrestrial Pole (CTP) and a zero
longitude (λ = 0°) on the equator.
Global earth deformations occur at the level of cm's/yr and there is a need to establish a
terrestrial network with an internal consistency of the same order or better.
Deformations of the terrestrial frame are produced by tectonic and geodynamic
processes, as well as components of the motion of the Earth’s rotation axis, etc. These
deformations are largely unpredictable and it is necessary to constantly monitor the
relations between the two fundamental reference frames.
The conventional terrestrial reference frame is defined by a CTS having the following
criteria:
(a) geocentric, the centre of mass being defined for the whole Earth,
including the oceans and atmosphere.
(b) its scale is that of a local Earth frame, in the meaning of a
relativistic theory of gravitation.
(c) its orientation is given by the Bureau International de l’Heure
(BIH) orientation at 1984.0.
(d) its time evolution in orientation will create no residual global
rotation with regards to the crust.
The geocentre is not fixed relative to the crust since the seasonal movements of the
oceans and atmosphere will result in origin shifts of the order of a few centimetres. The
direction of the z-axis was previously defined by the average position of the earth’s
rotation axis between the years 1900-1905. This position was known as the
conventional international origin (CIO) pole. In recent years, a new definition for a

5
reference pole has been used – the CTP – which is defined by IERS/BIH using a variety
of space techniques (SLR, LLR, VLBI, GPS). The x-axis passes through the point of
zero longitude and the earth’s mean equatorial plane, again defined by IERS/BIH. The
orientations are defined by adopting IERS Earth orientation parameters at a particular
reference epoch. The y-axis forms the completed right-handed system.
When one wants to realise such a CTRS through a reference frame (i.e., a network of
station coordinates), a selected set of station coordinates is chosen or prescribed by
IERS at a chosen epoch to. The reference frame defined by IERS from various solutions
is known as the International Terrestrial Reference Frame (ITRF). Currently we use
ITRF2000 station coordinates. Coordinates are usually specified by Cartesian
coordinates (X, Y, Z) but if geodetic coordinates are required, the use of the GRS80
ellipsoid is recommended.
Within the ITRF, the (instantaneous) position of a point on the surface of the solid Earth
should be expressed as
X (t ) = X o + V (t − to ) + ∑ Δ X i (t ) (4.2)
i

where Δ X i are corrections for various time changing effects (solid Earth tides, ocean
loading, atmosphere loading, etc), X o and V o are position and velocity components of
the point at time to and X (t ) is the position at time t. V o is usually given to first order
from a horizontal plate velocity model (NUVEL-1).

4.2.2.1 The ITRF2000 and ITRF2005

The following paragraph was obtained directly from the IERS technical report on the
ITRF2000 - see http://www.iers.org/iers/publications/tn/tn31/.
"The ITRF2000 solution is intended to be a standard solution for a wide
application community (geodesy, geophysics, navigation, etc.). It is achieved
by combining simultaneously station positions and velocities using full
variance-covariance information provided, in SINEX format, by the IERS
analysis centers. The individual solutions incorporated in the combination
are free from any external constraints, yielding actual space geodesy
estimates in terms of station positions and velocities. The ITRF2000 origin is
established by a weighted average of 5 over 7 submitted SLR solutions. Its
scale is obtained by a weighted average of 5 SLR and 3 VLBI solutions. The
ITRF2000 orientation is aligned to ITRF97 at epoch 1997.0 and its
orientation time evolution is aligned to the geophysical model NNR-NUVEL-
1A. The orientation and its time evolution were implemented using a
consistent geodetic way, anchored over 50 sites of high geodetic quality. The
ITRF2000 contains about 800 stations located on about 500 sites. Figure 1
shows the coverage of these sites, underlying the collocated space geodesy
techniques."
The actual information regarding the realisation of the ITRF2000 that geodesists require
for implementing the reference frame include station positions and velocities (at a given
epoch), transformation parameters between ITRF2000 and previous ITRF solutions as
well as the post-fit residuals and statistics from the current solution.
The most recent reference frame update, the ITRF2005, has just been introduced.
Transformation parameters between the ITRF2000 and ITRF2005 can be found at
http://itrf.ensg.ign.fr/ITRF_solutions/2005/tp_05-00.php.
For more information on the IERS and ITRF refer to http://www.iers.org/ and
http://itrf.ensg.ign.fr/.

6
 4.2.3 Transition Between Inertial and Terrestrial Systems
The transition from the space-fixed equatorial system (CIS) to the conventional
terrestrial system (CTS) is realised through a sequence of rotations that account for
- precession RP
- nutation RN
- earth rotation (GAST sidereal time) RS }
- polar motion RM } usually combined RR
For an arbitrary vector x, the transformation is given by

x[CTS ] = R (4.3)
M S N P
R R R x[CIS ]
These rotations are commonly referred to as Earth Orientation Parameters or EOP's and
are available from http://www.iers.org/iers/products/eop/ . The values are usually given
as
- Universal time (UT1) – is the time of the earth clock, which performs one
revolution in about 24h. It is practically proportional to sidereal time.
- Coordinates of the pole – x and y are the coordinates of the Celestial Ephemeris
Pole (CEP) relative to the IRP, the IERS Reference Pole.
- Celestial pole offsets – offsets are described in the IAU Precession and Nutation
models. The observed differences with respect to the conventional celestial pole
position defined by the models are monitored and reported by the IERS.
The CIS, defined at the standard epoch J2000.0, is transformed into the instantaneous or
true system at observation epoch by applying the corrections due to precession and
nutation. The z-axis of the true CIS represents the free position of the angular
momentum axis and thus points to the Celestial Ephemeris Pole (CEP). Rotating this
system around the z axis of the CIS system by the matrix RS does not change the CEP.
The CEP is rotated into the CIS (z-axis aligned) by RM, completing the transformation.
Unlike precession and nutation, the earth rotation parameters (RM, RS) cannot be
described through theory but must be determined through actual observations.

mean axis of rotation

Figure 4.2 – True Instantaneous and Mean Conventional Terrestrial System (space fixed)

7
 4.3 Ellipsoidal Reference Systems
A specifically oriented reference ellipsoid constitutes a geodetic datum. Over the years,
large numbers of geodetic (horizontal) datums have been created by various agencies
for surveying and mapping purposes over the entire earth or within particular regions. In
the past, the term geodetic datum has primarily been a horizontal datum and separate
vertical datums are established for orthometric heights. On each datum, a network of
control stations is established to provide access to accurate coordinates.

4.3.1 Global Geodetic Datums

As mentioned in section 4.1, a global ellipsoidal datum is one in which the centre of the
reference system is coincident with the earth's geocentre, the semi-major axis is
approximately equal to the equatorial radius of the earth and the 0° longitude is aligned
with the Greenwich meridian. In this system we usually define the position of terrestrial
points in terms of curvilinear coordinates (ϕ, λ, h) – although these may be converted to
rectangular coordinates.
There are many different global ellipsoids which are used in geodesy, each with slightly
different values for the a and f parameters required to define the size and shape of the
reference surface. Therefore, when stating curvilinear coordinates in a particular global
ellipsoidal datum, it is also necessary to identify the name of the datum and/or the
parameters which define it.

4.3.2 Regional Geodetic Datums

Eight parameters are required to define a regional geodetic datum: two to specify the
dimensions of the ellipsoid (a and f); three to specify the location of its centre
[( X o ,Yo , Z o ) or (φo , λo , ho )] ; and three to specify the orientation of the ellipsoid (R1,
R2, R3). The definition of the datum was usually effected by adopting a reference
ellipsoid of a particular shape; fixing the latitude, longitude and height of an initial point
or datum origin located near the centre of the geodetic network being established; fixing
the azimuth of a line from the origin point; and fixing the deflection of the vertical at the
initial point.
The ellipsoids of the regional datums are usually non-geocentric so that the ellipsoid can
conform as closely as possible to the geoid over the region (the area of interest for
mapping) rather than the whole globe – see Figure 4.3. A desirable condition for the
best-fitting ellipsoid is that the magnitudes of the geoid undulations be minimised, i.e.
n

∑N i
2
= minimum.
i =1

8
 REGIONAL GEODETIC
DATUM (ANS)

AUSTRALIA

Figure 4.3 – A Regional Ellipsoid (e.g. the ANS).

[Image obtained from http://www.ga.gov.au/geodesy/datums/aboutdatums.jsp]

The result is that not only do the centres of different local reference ellipsoids not
coincide, but the ellipsoids may have different a and f values, and may be slightly
rotated with respect to each other. Table 4-1 lists a few parameters for various
regional/global datums.

Table 4-1 - Parameters of Some Adopted Reference Ellipsoids

Coordinate System (Datum) Reference Ellipsoid a (m) l/f

GDA GRS-80 6,378,137 298.257222101

WGS-84 WGS-84 6,378,137 298.257223563
AGD ANS (or SA-69) 6,378,160 298.25
ED-79 International 6,378,388 297
GEM-8 GEM-8 6,378,145 298.255
GEM-9 (or GEM-10) GEM-9 (or GEM-10) 6,378,140 298.255
GEM-10B GEM-10B 6,378,138 298.257
GEM-T1 GEM-T1 6,378,137 298.257
NAD-27 Clarke 1866 6,378,206.4 294.9786982
NAD-83 GRS-80 6,378,137 298.257222101
NWL-9D = NSWC-9Z2 WGS-66 6,378,145 298.25
SA-69 SA-69 (or AN) 6,378,160 298.25
WGS-72 WGS-72 6,378,135 298.26
Note: AGD = Australian geodetic datum; ANS = Australian national spheriod; ED = European datum; GEM =
Goddard earth model; GRS = geodetic reference system; NAD = North American datum; NSWC = Naval surface
warfare center; NWL = Naval Weapons Laboratory; SA = South American; WGS = World geodetic system.

Of the datums listed in Table 4-1, those most commonly encountered in Australia are
the GDA (GRS80), AGD and WGS84. The GDA and WGS84 datums are geocentric
whilst the AGD is the regional datum for Australia.

9
 4.3.2.1 The AGD
For Australia, the geodetic datum (horizontal) was previously defined by the Australian
Geodetic Datum (AGD). The origin of this network is at Johnston Geodetic Station
where
φo = −25o56′54.5515′′
λo = 133o12′30.0771′′
ho = 571.2 m.
The minor axis of the AGD ellipsoid is termed the Australian National Spheroid (ANS)
and has parameters of:
a = 6,378,160 m
f-1 = 298.25
The N (geoid-ellipsoid separation) value at Johnston was taken as +4.9 m (for AGD84
and 0 for AGD66) with deflection of vertical components of
ξ o = −4.2 and ηo = −4.5 .
′′ ′′
Coordinates of points in this geodetic datum where given in either ellipsoidal or grid
values:
AGD84 φ, λ, h coords
AMG84 grid coords
In Tasmania, the State network operated in the AGD66 system, whereas most mainland
states worked in the AGD84 system.

4.3.2.2 The GDA

As of 2000, the Australian Geodetic Datum (AGD) has been replaced by the Geocentric
Datum of Australia (GDA) as the official Australian coordinate system. GDA is a 3-D,
geocentric reference system, defined within the ITRF and compatible with the WGS84
system (see Table 4-2). For full details of the GDA and its implementation see
http://www.ga.gov.au/geodesy/datums/gda.jsp.

Table 4-2 - GDA Specifications (from http://www.ga.gov.au/geodesy/datums/gda.jsp)

TERMINOLOGY
Datum: Geocentric Datum of Australia (GDA)
Geographical coordinate set:
Geocentric Datum of Australia 1994 (GDA94)
(latitude and longitude)
Map Grid of Australia 1994 (MGA94)
Grid coordinates : (Universal Transverse Mercator, using the
GRS80 ellipsoid)
DEFINITION
Reference Frame : ITRF92
Epoch : 1994.0
Ellipsoid : GRS80
Semi-major axis (a) : 6,378,137.0 metres
Inverse flattening (1/f) : 298.257222101

10
The GDA is defined by the sites of the Australian Fiducial Network (AFN), see Figure
1.1. The adopted coordinates of the AFN were computed in terms of the International
Terrestrial Reference Frame (ITRF) and from these AFN coordinates the geodetic
network that covers the whole continent was recomputed.
The main reason for adopting a 3-D geocentric datum was to accommodate the
increased use of GPS by both scientists and the general public. As GPS is based on a
geocentric datum, the introduction of GDA eliminates the need for complex coordinate
conversions.

4.4 Local (Topocentric) Coordinate Systems

4.4.1 Local Astronomical Coordinate Systems

When we make terrestrial geodetic observations (e.g. with a theodolite) all of our
readings (with the exception of slant ranges) are related to the local gravity vector g. It
is therefore easiest to describe our coordinates from these observations in terms of a
local reference coordinate system that is tied to the direction of the plumb line at the
observation point.

Az
Az

Figure 4.4 – Local Astronomical (Topocentric) Coordinate System

As shown in Figure 4.4, we define the local system with point P (the observation point)
as the origin and the orientation of the unit normal vector n (or the z'-axis) is usually
determined from astronomical observations of the astronomical (geographical) latitude
Φ and the astronomical longitude Λ, and can be written in terms of direction cosines as

⎡ cos Φ cos Λ ⎤
n = ⎢⎢ cos Φ sin Λ ⎥⎥ (4.4)
⎢⎣ sin Φ ⎥⎦
The axes of this system are therefore defined as:
- z' axis directed to the astronomical zenith (local vertical)
- x' axis directed to the astronomical north meridian
- y' axis directed to the east or zonal direction
Note that this coordinate system is a left-handed coordinate system.

11
The location of an (observed) point Pi within the local astronomical coordinate system is
usually derived from terrestrial observations of astronomical azimuth Az, slant range s,
and zenith angle z (see Figure 4.4). The Cartesian coordinates in the local topocentric
system for Pi can therefore be written as
⎡ xP' i ⎤ ⎡cos Az sin z ⎤
⎢ ' ⎥
X Pi = ⎢ yPi ⎥ = s ⎢⎢ sin Az sin z ⎥⎥
'
(4.5)
⎢ z P' ⎥ ⎢⎣ cos z ⎥⎦
⎣ i⎦
The relationship between the local astronomical and the global astronomical Cartesian
coordinate systems is also illustrated in Figure 4.4.

4.4.2 Local Geodetic Coordinate Systems

The most important difference between the two types of topocentric coordinate systems
is that in a local astronomical system observations are referred to the local vertical (i.e.
gravity), whilst in a local geodetic coordinate system observations are referred to the
local ellipsoidal normal (see Figure 4.5). As this coordinate system is referred to an
ellipsoidal parameter, the definition of the local system is therefore dependent on the
selection of the reference ellipsoid (i.e. when quoting coordinates, also need to specify
the reference ellipsoid).

Figure 4.5 – Local Geodetic (Topocentric) Coordinate System

In a local geodetic (ellipsoidal) system we therefore define our axes (see Figure 4.5):
- origin at the observation point P
- u axis in the direction of the ellipsoidal vertical
- n axis perpendicular to the u axis and directed toward the
ellipsoidal north (i.e. the geodetic meridian)
- e axis perpendicular to u and n and directed to the east
Again, note that this system is a left-handed system.

12
We can define the location of a second point (Pi) in a local ellipsoidal system via the
following quantities (see Figure 4.5)
- slant range s
- ellipsoidal azimuth α
- ellipsoidal zenith angle ζ
Given this information (s, α, ζ), we can determine the local geodetic coordinates of Pi
using the following relations:

⎡n⎤ ⎡cos α sin ζ ⎤

⎢ e ⎥ = s ⎢ sin α sin ζ ⎥ (4.6)
⎢ ⎥ ⎢ ⎥
⎢⎣u ⎥⎦ ⎢⎣ cos ζ ⎥⎦
It is also possible to determine the inverse relations (observations in terms of
coordinates) for a local ellipsoidal system:

⎛e⎞
α = tan −1 ⎜ ⎟
⎝ ⎠
n
⎛u⎞
ζ = cos −1 ⎜ ⎟ (4.7)
⎝s⎠
s = n2 + e2 + u 2
The relationship between the local geodetic and the global geodetic coordinate systems
is also illustrated in Figure 4.5 and will be discussed further later.

4.5 Transformations Between Datums

4.5.1 Curvilinear Geodetic <> Cartesian Geodetic (ϕ,λ,h <> X,Y,Z)

The point P' on the ellipsoid (see Fig 3.7) is obtained by projecting the surface point P
along the ellipsoidal normal (Helmert’s projection). From Figure 4.6 we can see that the
relationship between the curvilinear and Cartesian coordinates for P' can be expressed
as:

⎡ X P' ⎤ ⎡ cos ϕ cos λ ⎤

⎢ ⎥ ⎢ ⎥
P ' = ⎢ YP' ⎥ = ν ⎢ cos ϕ sin λ ⎥ (4.8)
⎢ ⎥
⎣⎢(1 − e ) sin ϕ ⎦⎥
⎢Z ' ⎥ 2
⎣ P⎦

where ν is the prime vertical radius of curvature at point P'.

When considering a surface point, P, at a distance h along the ellipsoidal normal either
above or below the ellipsoid, the transformation from curvilinear to Cartesian
coordinates becomes:
⎡ ⎤
⎡ X P ⎤ ⎢ (ν + h ) cos ϕ cos λ ⎥
⎢ Y ⎥ = ⎢ ν + h cos ϕ sin λ ⎥
⎢ P⎥ ⎢( )
⎥
(4.9)

(( ) )
⎢⎣ Z P ⎥⎦ ⎢ ν 1 − e 2 + h sin ϕ ⎥
⎣ ⎦

13
 Z Z P
p
p'
h
P'
Y p' ϕ
p' ϕ
λ
X Y
ν

p'
cos ϕ =
ν
∴ p ' = ν cos ϕ
⎛ dz ⎛ b ⎞ p'⎞
2
⎛b⎞
2

⎟⎟ and ⎜ ⎟ = (1 − e )
2
X P ' = p 'cos λ = ν cos ϕ cos λ from eqn 2.21 ⎜ = −⎜ ⎟
⎜ dp ⎝a⎠ z ⎠ ⎝a⎠
⎝
YP ' = p 'sin λ = v cos ϕ sin λ
⇒ Z P ' = (1 − e 2 ) tan ϕ ⋅ p ' = ν (1 − e 2 ) sin ϕ

Figure 4.6 – Proof of Curvilinear to Cartesian Transformation

The inverse problem of determining ϕ, λ, h from X,Y,Z is not as straight forward and
requires iteration to solve for ϕ and h. Fortunately, the system of equations converges
quickly since h << v. From Figure 4.6, we can see:
p
h= −v
cos ϕ
⎡z⎛ 2⎛ v ⎞⎞ ⎤
−1

ϕ = arctan ⎢ ⎜1 − e ⎜ ⎟⎟ ⎥ (4.10)
⎢⎣ p ⎝ ⎝ v + h ⎠ ⎠ ⎥⎦

where: p = X 2 + Y 2
In this transformation, the solution for λ is exact and given as:
⎛Y ⎞
λ = arctan⎜ ⎟ (4.11)
⎝X⎠
Another approach is to use the following rigorous, non-iterative method developed by
Bowring (1985 - see Survey Review, 28: 202-206).

14
 Y ⎫
tan λ = ⎪
X
⎪
Z (1 − f ) + e 2 a sin 3 μ ⎪
tan ϕ = ⎪
(1 − f ) ⎡⎣ p − e2 a cos3 μ ⎤⎦ ⎪
a2 ⎪
h = p cos ϕ + z sin ϕ − ⎪
v ⎬ ( X , Y , Z ) → (ϕ , λ , h )
⎪
⎪
where e2 = 2 f − f 2 ⎪
⎪
⎡ ⎤ ⎪
z⎢ e2 a ⎥
tan μ = (1 − f ) + ⎪
p⎢ 2 ⎥
( X + Y + Z ) ⎥⎦
1
2 2 2 ⎪
⎢⎣ ⎭⎪
Recently Pollard (Journal of Geodesy (2002) 76: 36-40) has suggested a vector method
which in some cases provides better results (especially height) than Bowring's however,
for most applications the results (especially for latitude) are of a similar quality.
It is important to note that using these transformations, we can convert from (X, Y, Z) to
any ellipsoid system (ϕ, λ, h) given the parameters of the ellipsoid (a, f). Hence we can
transform from (ϕ, λ, h)ellipsoid1 to (ϕ, λ, h)ellipsoid2 by using the steps

(ϕ , λ , h )ellipsoid 1 → ( X , Y , Z ) → (ϕ , λ , h )ellipsoid 2
However, this assumes that the origins and orientations of ellipsoids and XYZ reference
frames coincide. If any origin shifts are involved we need to use different formulae but
it can still be broken down to simple steps as above. There are also other formulae
which directly transform coordinates between ellipsoid datums and we will deal with
these later in this chapter.

4.5.1.1 Transforming coordinate changes

Other useful equations allow transformation between changes in global Cartesian and
curvilinear coordinates, since it is usually easier to visualise changes in latitude,
longitude and height rather than changes in global Cartesian coordinates (e.g. useful for
error analysis or for 3-D eccentricity computations that refer eccentric station
observations back to the main control station). The following equations can be obtained
by differentiating equation (4.9) and it can be shown (after a fair bit of algebra) that:

⎡ dX ⎤ ⎡ − ( ρ + h ) cos λ sin ϕ − ( v + h ) cos ϕ sin λ cos ϕ cos λ ⎤

⎢ dY ⎥ = ⎢ − ρ + h sin λ sin ϕ
⎢ ⎥ ⎢ ( ) ( v + h ) cos ϕ cos λ cos ϕ sin λ ⎥⎥ ⋅ [ dϕ d λ dh ]
T

⎢⎣ dZ ⎥⎦ ⎢⎣ ( ρ + h ) cos ϕ 0 sin ϕ ⎥⎦

⎡ dϕ ⎤
= J ⎢⎢ d λ ⎥⎥
⎢⎣ dh ⎥⎦

with precisions carried through by variance propagation QXYZ = JQϕλh J T .

15
Rearranging the above equation we get:

⎡ dϕ ⎤ ⎡ dX ⎤
⎢ d λ ⎥ = J −1 ⎢ dY ⎥
⎢ ⎥ ⎢ ⎥
⎢⎣ dh ⎥⎦ ⎢⎣ dZ ⎥⎦

⎡ − sin ϕ cos λ − sin ϕ sin λ cos ϕ ⎤

⎢ ρ +h ρ +h ρ + h⎥
⎢ ⎥
⎢ − sin λ cos λ ⎥
where J −1 = ⎢ 0 ⎥
⎢(
ν + h ) cos ϕ (ν + h ) cos ϕ ⎥
⎢ cos ϕ cos λ cos ϕ sin λ sin ϕ ⎥
⎢ ⎥
⎣ ⎦
We can also write a similar equation (easily proved using ellipsoidal geometry) for
transforming between global curvilinear and local geodetic coordinate changes:

⎡ Δn ⎤ ⎡( ρ + h ) 0 0 ⎤ ⎡ Δφ ⎤
⎢ Δe ⎥ = ⎢ 0 ⎥
0 ⎥ ⎢⎢ Δλ ⎥⎥
⎢ ⎥ ⎢ (ν + h ) cos φ
⎢⎣ Δu ⎥⎦ ⎣⎢ 0 0 1 ⎦⎥ ⎢⎣ Δh ⎥⎦

4.5.2 Local <> Global

4.5.2.1 Astronomical
Observed coordinate differences can be transformed from the local topocentric
astronomical system into the global Cartesian system using the transformation

ΔX Global = A ΔX local (4.12)

where A = Rz (180° − Λ) Ry (90° − Φ) S y

The matrix Sy changes the orientation of the y' axis and converts the left-handed into a
right-handed coordinate system. Doing the matrix multiplication for A gives

⎡ − sin Φ cos Λ − sin Λ cos Φ cos Λ ⎤

A = ⎢⎢ − sin Φ sin Λ cos Λ cos Φ sin Λ ⎥⎥ (4.13)
⎢⎣ cos Φ 0 sin Φ ⎥⎦ evaluated at po int P

The inverse transformation is:

ΔX local = A−1ΔX Global = AT ΔX Global (4.14)

The above formulae are most useful when wanting to combine results from local
terrestrial observations and from satellite techniques.
4.5.2.2 Ellipsoidal
For the transformation of coordinate differences from the local topocentric to the global
ellipsoidal system we obtain the relation:

16
 ΔX global = A ΔX local
or (4.15)
ΔX global = Rz (180o − λ ) Ry ( 90o − ϕ ) S y ΔX n ,e,u

where

⎡ − sin ϕ cos λ − sin λ cos ϕ cos λ ⎤

A = ⎢⎢ − sin ϕ sin λ cos λ cos ϕ sin λ ⎥⎥ (4.16)
⎢⎣ cos ϕ 0 sin ϕ ⎥⎦ evaluated at point P

Using the above expressions it is also possible to express observations in the local
system (from point P looking to point Pi) in terms of the global Cartesian coordinate
differences:
⎛ − sin λΔX + cos λΔY ⎞
α = tan −1 ⎜ ⎟
⎝ − sin ϕ cos λΔX − sin ϕ sin λΔY + cos ϕΔZ ⎠
⎛ cos ϕ cos λΔX + cos ϕ sin λΔY + sin ϕΔZ ⎞
ζ = cos −1 ⎜ ⎟ (4.17)
⎝ s ⎠
s = ΔX 2 + ΔY 2 + ΔZ 2

Proof: by expanding terms for ΔX n ,e ,u i.e. ΔX n ,e,u = AT ΔX global

⎡ Δn ⎤ ⎡ − sin ϕ cos λΔX − sin ϕ sin λΔY + cos ϕΔZ ⎤

⎢ Δe ⎥ = ⎢ − sin λΔX + cos λΔY ⎥ (4.18)
⎢ ⎥ ⎢ ⎥
⎢⎣ Δu ⎥⎦ ⎢⎣cos ϕ cos λΔX + cos ϕ sin λΔY + sin ϕΔZ ⎥⎦
and substituting equation (4.18) into equation (4.7) gives the relations above (equation
4.17).

Figure 4.7 – Local and Global Geodetic Reference Systems

17
 4.5.3 Local Geodetic <> Local Astronomical (n,e,u <> x',y',z')

⎡n⎤ ⎡ x '⎤
⎢ e ⎥ = R ⎢ y '⎥
⎢ ⎥ [ ]⎢ ⎥ (4.19)
⎢⎣u ⎥⎦ ⎢⎣ z ' ⎥⎦

where R = Rx ( 90° − ϕ ) Rz ( 90° + λ )

⎡1 0 0 ⎤ ⎡cos θ 0 − sin θ ⎤ ⎡ cos θ sin θ 0⎤

Rx (θ ) = ⎢⎢0 cos θ sin θ ⎥⎥ Ry (θ ) = ⎢⎢ 0 1 0 ⎥⎥ Rz (θ ) = ⎢⎢ − sin θ cos θ 0 ⎥⎥
⎢⎣0 − sin θ cos θ ⎥⎦ ⎢⎣ sin θ 0 cos θ ⎥⎦ ⎢⎣ 0 0 1 ⎥⎦

⎡ − sin ϕ cos λ − sin ϕ sin λ cos ϕ ⎤

∴ R = ⎢⎢ − sin λ cos λ 0 ⎥⎥
⎢⎣ cos ϕ cos λ cos ϕ sin λ sin ϕ ⎥⎦

Referring back to equation (4.16) we can see that in fact the rotation matrix R is equal to
the transpose of the A matrix used in previous transformations. Therefore we have:
ΔX LG = AT ΔX LA (4.20)

Conversely if (n,e,u) are known and (x',y',z') are required we can use
⎡ x' ⎤ ⎡n⎤
⎢ '⎥ ⎢ ⎥
⎢ y ⎥ = R (ϕ , λ ) ⎢ e ⎥
−1
(4.21)
⎢ z' ⎥ ⎢⎣u ⎥⎦
⎣ ⎦

⎡ − sin ϕ cos λ − sin λ cos ϕ cos λ ⎤

where R = ⎢⎢ − sin ϕ sin λ
−1
cos λ cos ϕ sin λ ⎥⎥
⎢⎣ cos ϕ 0 sin ϕ ⎥⎦

Using equation (4.21) in equation (4.7) we can get alternate expressions for the
ellipsoidal quantities α, ζ and s in terms of Δx', Δy', Δz' and the geodetic position of P.

⎛ − sin λΔx '+ cos λΔy ' ⎞

α = tan −1 ⎜ ⎟
⎝ − sin ϕ cos λΔx '− sin ϕ sin λΔy '+ cos ϕΔz ' ⎠
⎛ ⎞
⎜ cos ϕ cos λΔx '+ cos ϕ sin λΔy '+ sin ϕΔz ' ⎟
ζ = cos −1
(4.22)
⎜⎜ ⎟⎟
( )
1
2 2 2 2
⎝ Δx ' + Δy ' + Δz ' ⎠
s = ( Δx '2 + Δy '2 + Δz '2 )
1
2

18
 4.5.4 Cartesian <> Cartesian
The most commonly used transformation between any two Cartesian systems (e.g. from
(u,v,w) to (X,Y,Z)) is a seven parameter (similarity) transformation (see equation 4.23
below). This transformation includes 3 translations, 3 rotations and 1 scale factor to
account for the difference in location of the geometrical centre of each reference system,
the alignment of the axes and a scale change between the systems.

⎡ X ⎤ ⎡ ΔX ⎤ ⎡1 γ −β ⎤ ⎡ u ⎤
⎢ Y ⎥ = ⎢ ΔY ⎥ + 1 + δ s ⎢ −γ α ⎥⎥ ⎢⎢ v ⎥⎥
⎢ ⎥ ⎢ ⎥ ( )⎢ 1 (4.23)
⎢⎣ Z ⎥⎦ ⎢⎣ ΔZ ⎥⎦ ⎢⎣ β −α 1 ⎥⎦ ⎢⎣ w ⎥⎦

where ΔX, ΔY, ΔZ are the coordinates of the origin of the (u,v,w) frame in the (X,Y,Z)
frame, α, β, γ are differential rotations (small) around the (u,v,w) axes respectively to
establish parallelism with the (X,Y,Z) frame and δs is a differential scale change
between the two frames.
Typical examples of the use of this technique in the Australian context are
transformations from the non-geocentric Australian Geodetic Datum (AGD66 or
AGD84) to the geocentric WGS84 or GDA datums.
In order to employ a similarity transformation, information is required on the 7
parameters defined above. These parameters are generally well defined between
commonly used datums and are widely available. An example of these parameters are
those which are used convert from AGD84 (and AGD66) to GDA94 (expected to have
an accuracy of about 1 metre – see the GDA technical manual
http://www.icsm.gov.au/icsm/gda/gdatm/gdav2.2.pdf ) and WGS84:

Table 4.3 - National Table 4.3 - Tasmanian

parameters AGD84 to GDA parameters AGD66 to GDA

ΔX -117.763 m ΔX -120.271 m
ΔY -51.510 m ΔY -64.543 m
ΔZ 139.061 m ΔZ 161.632 m
α -0.292" α -0.217"
β -0.443" β 0.067"
γ -0.277" γ 0.129"
δs -0.191 ppm δs 2.499 ppm

Table 4.4 – Parameters for AGD66/AGD84 to WGS84

Parameter AGD66 → WGS84 AGD84 → WGS84
ΔX -123.314 m -116.00 m
ΔY -47.223 m -50.47 m
ΔZ 136.594 m 141.69 m
α -0.264" -0.23"
β -0.332" -0.39"
γ -0.270" -0.344"
δs -1.384 ppm 0.0983 ppm

19
If you want to transform coordinates in the other direction i.e. from WGS84 (or GDA)
to AGD84 (or AGD66), then the SIGNS of all of the above transformation parameters
are reversed and equation (4.23) is used as before.

4.5.6 Curvilinear Geodetic <> Curvilinear Geodetic

For a transformation between two different ellipsoidal reference systems with a
difference in size and shape of ellipsoids, we can use several approaches. In this
transformation it is necessary to account for the difference in the location of the
geometrical centre of each reference ellipsoid, the difference in size and shape of the
two ellipsoids, and the difference in orientation and scale.
Consider the ellipsoids with sizes and shapes defined by (a1, f1) and (a2, f2) and the
differences in location of each ellipsoid centre defined by a set of translations (ΔX, ΔY,
ΔZ), rotation angles (α, β, γ) and scale change (δs).
If we have coordinates (ϕ1, λ1, h1) on datum 1 and wish to find the coordinates (ϕ2, λ2,
h2) on datum 2, we can use several methods to find a solution.

4.5.6.1 Similarity Transformation

The first approach, a direct approach, involves the use of formula given in this chapter
and is best described in diagrammatic form:

Curvilinear Coords (ϕ1, λ1, h1) Curvilinear Coords (ϕ2, λ2, h2)

(ϕ, λ, h) to (X, Y, Z) transformation (X, Y, Z) to (ϕ, λ, h) transformation

Cartesian Coords (X1, Y1, Z1) Cartesian Coords (X2, Y2, Z2)

7-parameter transformation

Datum 1 Datum 2

Figure 4.8 – 7-parameter Transformation from Ellipsoid 1 to Ellipsoid 2

This approach relies on a knowledge of the transformation parameters between the two
Cartesian reference frames.

20
Numerical Example
Given E, N, H coordinates in AMG84, transform to ϕ, λ, h in WGS84.
AMG84 (k0 = 0.9996, λ0 = 147°, E0 = 500,000 m, N0 = 10,000,000 m)
E: 526 559.1976 m N: 5 249 816.8013 m H: 8.087 m (AHD)
Step 1
Convert E & N to ϕ & λ on ANS (a = 6378160, f-1 = 298.25)
(e.g. using Redfearn's method)
Gives geodetic curvilinear coords on AGD84:
ϕANS: -42° 54′ 13.1379″ λANS: 147° 19′ 31.2137″
Step 2
Now we have (ϕ, λ), we need to determine geoid height NANS on the same
ellipsoid (ANS) to define the full set of curvilinear coords. (Note: Can use
program WINTER to get N value – or if need NGRS80 can use
http://www.ga.gov.au/geodesy/ausgeoid/nvalcomp.jsp)
hANS = H + NANS = 8.087 + 18.509 = 26.596 m
Step 3
We now have (ϕ, λ, h) on the ANS, next step is to calculate (X, Y, Z)AGD84 (e.g.
using Bowring's formula or iterating direct formula)
X: -39388784.777 m Y: 2526195.908 m Z: -4319700.532 m
Step 4
Next, we want to use 7-parameter transformation to go from (X,Y,Z)AGD84 to
(X,Y,Z)WGS84 - using parameters given above in Table 4.4.
X: -3938913.545 m Y: 2526143.935 m Z: -4319549.00 m
Step 5
To convert these Cartesian coordinates to (ϕ, λ, h)WGS84, we need a & f values
for the WGS84 ellipsoid
a = 6378137 f-1 = 298.257223563
The solution is therefore
ϕWGS84: -42° 54 07′.7495″ λWGS84: 147° 19′ 36.2064″ hWGS84: 5.010 m

Points to note:
- Curvilinear coords require knowledge of ellipsoid parameters (a & f)
- Orthometric height (H) is invariant to datum definition.
- Geoid ht (N) and ellipsoidal ht (h) are ellipsoid dependent, and hence values
must be on the same ellipsoid (i.e. hANS ≠ NWGS84 + H)
- When changing datum, we can sometimes make use of the fact that the
difference in ellipsoidal height is equivalent to the geoid height change
between datums (ignoring rotations and scale change), i.e.:
hWGS84 = NWGS84 + H and hANS = NANS + H
and therefore since H is constant:
Δh(WGS84 - ANS) = ΔN(WGS84 – ANS).

A common problem that arises with this method of transformation is that one of the
geoid heights (N) may not be readily available for converting our H value to an h value
(e.g may not have access to regional datum values). We will now go through the above
example again, this time without NANS i.e. skip step 2.

21
In step 3, use an initial ellipsoidal height of hANS = 0 and determine Δh between the two
ellipsoids. Then determine NANS from Δh and NWGS84 and re-do steps 3-5 with the
correct hANS value.

4.5.6.2 Soler and Hothem Equations

(ref. Journal of Geodesy (1998) 72: 482-490)
Transformation of curvilinear geodetic coordinates between datums (say, from Datum 1
to Datum 2) can be done by a general differential transformation as

⎡ϕ ⎤ ⎡ϕ ⎤ ⎡ dϕ ⎤
⎢ λ ⎥ = ⎢ λ ⎥ + ⎢ dλ ⎥ (4.24)
⎢ ⎥ ⎢ ⎥ ⎢ ⎥
⎢⎣ h ⎥⎦ D ⎢⎣ h ⎥⎦ D ⎢⎣ dh ⎥⎦
2 1

or, in terms of Cartesian coordinates:

⎡(ν + h ) cos ϕ ∂λ ⎤ ⎧⎧∂u ⎫ ⎧∂u ⎫ ⎫

⎢ ⎥ ⎪⎪ ⎪ ⎪ ⎪ ⎪
⎢ ( ρ + h ) ∂ϕ ⎥ = [ R ] ⎨⎨∂v ⎬ + ⎨∂v ⎬ ⎬ (4.25)
⎢⎣ ∂h ⎥⎦ ⎪⎪∂w⎪ ⎪∂w⎪ ⎪
⎩⎩ ⎭7 − param ⎩ ⎭∂a ,∂f ⎭
where

⎡ ∂u ⎤ ⎡ ΔX ⎤ ⎡ ∂s ∂ω −∂ψ ⎤ ⎡ u ⎤
⎢ ∂v ⎥ = ⎢⎢ ΔY ⎥⎥ + ⎢⎢ −∂ω ∂s ∂ε ⎥⎥ ⎢⎢ v ⎥⎥ (4.26)
⎢ ⎥
⎢⎣ ∂w⎥⎦ 7 − param ⎢⎣ ΔZ ⎥⎦ ⎢⎣ ∂ψ −∂ε ∂s ⎥⎦ ⎢⎣ w ⎥⎦

and

⎡ ∂u ⎤
⎢ ∂v ⎥ ⎡∂a ⎤
⎢ ⎥ = − [ D ] ⎢ ∂f ⎥ (4.27)
⎢⎣ ∂w⎥⎦ ∂a ,∂f ⎣ ⎦

where

⎡ ⎤
⎢ cos ϕ cos λ a (1 − f ) sin 2 ϕ cos ϕ cos λ ⎥
⎢ ⎥
( ) (1 − e2 sin 2 ϕ ) 2
1 3
⎢ 1 − e 2
sin 2
ϕ 2
⎥
⎢ ⎥
cos ϕ sin λ a (1 − f ) sin ϕ cos ϕ sin λ ⎥
2

[ D ] = ⎢⎢ ⎥
⎢ (1 − e sin ϕ ) (1 − e2 sin 2 ϕ ) 2
1 3
2 2 2
⎥
⎢ ⎥
⎢ (1 − e ) sin ϕ
2

⎢ ( ρ sin ϕ − 2ν ) (1 − f ) sin ϕ ⎥⎥
2

⎢⎣ (1 − e sin ϕ )
1
2 2 2
⎥⎦ Datum1

and

∂f = f D 2 − f D1
∂a = aD 2 − aD1 ⇒ ∂a = ∂a − aD1∂s
where ∂s = sD 2 − sD1

22
 4.5.6.3 Molodenskii formulae
For transforming between ellipsoid datums, the (more straight-forward) alternative to
the Soler and Hothem equations is to use the Molodenskii formulae:

⎛ ∂a (ν e sin ϕ cos ϕ ) ⎡ ρ ⎛ a ⎞ + ν ⎛ b ⎞⎤ sin ϕ cos ϕ ⎞

2

⎜ −ΔX sin ϕ cos λ − ΔY sin ϕ sin λ + ΔZ cos ϕ + + ∂f

⎢⎣ ⎜⎝ b ⎟⎠ ⎜⎝ a ⎠⎦⎟⎥ ⎟
∂ϕ =
⎝ a ⎠
( ρ + h)
( −ΔX sin λ + ΔY cos λ )
∂λ = (4.28)
(ν + h ) cos ϕ
∂h = ΔX cos ϕ cos λ + ΔY cos ϕ sin λ + ΔZ sin ϕ − ∂a ⎜
⎛ a ⎞+δ ⎛ b ⎞ν sin 2 ϕ
⎟ f ⎜ ⎟
⎝ν ⎠ ⎝a⎠

where dϕ and dλ are in radians dh is in metres

ϕ, λ, ν, ρ, a, b, h are values of the point/ellipsoid (Datum 1)
ΔX, ΔY, ΔZ are the shifts between the ellipsoid centres of D1 → D2.

Note that you should be able to validate equations (4.28) by expanding the matrix
equations (4.24) to (4.27).
To simplify calculation, a series of Abridged Molodenskii equations were developed:

∂ϕ =
( −ΔX sin ϕ cos λ − ΔY sin ϕ sin λ + ΔZ cos ϕ + ( a∂f + f ∂a ) sin 2ϕ )
ρ
( −ΔX sin λ + ΔY cos λ )
∂λ =
ν cos ϕ
∂h = ΔX cos ϕ cos λ + ΔY cos ϕ sin λ + ΔZ sin ϕ + ( a∂f + f ∂a ) sin 2 ϕ − ∂a

The Molodenskii formulae only are applicable if there is no scale change and no
rotation between Datum 1 and Datum 2, otherwise need to use Soler and Hothem
equations or 7-parameter transformation technique. Using these Molodenskii formulae
will generally produce results of lower accuracy than the other methods (depending on
the accuracy of the parameters in the 7-parameter transformation).

Numerical Example.
Using the data from the previous example we have curvilinear coordinates for a point in
the AGD84 datum (a = 6378160 m, f-1 = 298.25).
ϕ -42° 54′ 13.1379″ λ 147° 19′ 31.2137″ h 26.596 m
Now, using the Soler and Hothem equations and the Molodenskii Formulae, we want to
solve for the curvilinear cords of the point in the WGS84 datum (a = 6378137 m, f-1 =
298.257223563)
Using Soler & Hothem equations:

Using Molodenskii formulae:

23
From these solutions it can be seen that scale uncertainty has pronounced effect on h,
but minimal on ϕ and λ.

4.5.7 ITRF?? <> ITRF??

Over the last two decades, there have been different versions of the ITRF published.
Although these datums are all called 'ITRF' they should in fact be considered as
different reference frames. Therefore, transforming between ITRF frames is a
specialised form of Cartesian to Cartesian transformation that is commonly encountered
in geodetic work. Information on the required transformation and parameters is
available at ftp://lareg.ensg.ign.fr/pub/itrf/ITRF.TP sections of which are reproduced
below. See http://itrf.ensg.ign.fr/ITRF_solutions/2005/tp_05-00.php for ITRF200/2005
transformation parameters.
----------------------------------------------------------------------------------------------------------
Transformation Parameters and their Rates from ITRF2000 to Previous Frames

SOLUTION T1 T2 T3 D R1 R2 R3 EPOCH Ref.

cm cm cm ppb .001" .001" .001" IERS Tech.Note#
RATES T1 T2 T3 D R1 R2 R3
cm/y cm/y cm/y ppb/y .001"/y .001"/y .001"/y
ITRF97 0.67 0.61 -1.85 1.55 0.00 0.00 0.00 1997.0 27
Rates 0.00 -0.06 -0.14 0.01 0.00 0.00 0.02
ITRF96 0.67 0.61 -1.85 1.55 0.00 0.00 0.00 1997.0 24
Rates 0.00 -0.06 -0.14 0.01 0.00 0.00 0.02
ITRF94 0.67 0.61 -1.85 1.55 0.00 0.00 0.00 1997.0 20
Rates 0.00 -0.06 -0.14 0.01 0.00 0.00 0.02
ITRF93 1.27 0.65 -2.09 1.95 -0.39 0.80 -1.14 1988.0 18
Rates -0.29 -0.02 -0.06 0.01 -0.11 -0.19 0.07
ITRF92 1.47 1.35 -1.39 0.75 0.00 0.00 -0.18 1988.0 15
Rates 0.00 -0.06 -0.14 0.01 0.00 0.00 0.02
ITRF91 2.67 2.75 -1.99 2.15 0.00 0.00 -0.18 1988.0 12
Rates 0.00 -0.06 -0.14 0.01 0.00 0.00 0.02
ITRF90 2.47 2.35 -3.59 2.45 0.00 0.00 -0.18 1988.0 9
Rates 0.00 -0.06 -0.14 0.01 0.00 0.00 0.02
ITRF89 2.97 4.75 -7.39 5.85 0.00 0.00 -0.18 1988.0 6
Rates 0.00 -0.06 -0.14 0.01 0.00 0.00 0.02
ITRF88 2.47 1.15 -9.79 8.95 0.10 0.00 -0.18 1988.0 IERS An.Rep.1988
Rates 0.00 -0.06 -0.14 0.01 0.00 0.00 0.02

Note : These transformation parameters should be used with the standard model given
below and are valid at the indicated epoch.

⎡ XS ⎤ ⎡ X ⎤ ⎡ T 1 ⎤ ⎡ D - R3 R 2 ⎤ ⎡ X ⎤
⎢ YS ⎥ = ⎢ Y ⎥ + ⎢T 2 ⎥ + ⎢ R3 D - R1⎥ ⎢ Y ⎥
⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥⎢ ⎥
⎢⎣ ZS ⎥⎦ ⎢⎣ Z ⎥⎦ ⎢⎣T 3 ⎥⎦ ⎢⎣- R 2 R1 D ⎥⎦ ⎢⎣ Z ⎥⎦
Where X,Y,Z are the coordinates in ITRF2000 and XS,YS,ZS are the coordinates in the
desired frame.
Before using this transformation however, it may be necessary to adjust the parameters
to obtain their value at the observation epoch (t). This adjustment is done using the
following equation

P (t ) = P( EPOCH ) + P& * (t - EPOCH )

24
where EPOCH is the epoch indicated in the above table and P& is the rate of that
parameter.
The procedure for a 7-parameter transformation and 14-parameter transformation (i.e. 7
parameters and their rates/velocities) between the ITRF and GDA94 is outlined in the
Geoscience Australia publication “International Terrestrial Reference Frame (ITRF) to
GDA94 Coordinate Transformation” (Dawson & Steed, 2004), which is available from
the GA website at http://www.ga.gov.au/geodesy/datums/calcs.jsp.

4.5.8 Parallelism Conditions for Astronomical and Ellipsoidal

Systems (Laplace Equation)
Let us now examine the relations among the four major geodetic coordinate systems:
astronomical → CTRS (CT) and Local Astronomical (LA)
ellipsoidal → Geocentric (G) and Local Geodetic (LG)
A unit vector (or coordinate differences) can be rotated from the LA to CT system as
follows (from equations in previous sections):
ΔX CT = RΔX G and ΔX G = A ΔX LG
From Figure 4.9 we can see from inspection that
ΔX LG = RZ ( ΔAz ) RY ( −ξ ) RX (η ) ΔX LA
where
ΔAz = difference in Azimuth = angle between xLG & xLA axes
ξ = meridion component of deflection of vertical
η = prime vertical component of deflection of vertical

Figure 4.9 – Physical Relationship Between the Four Major Geodetic Coordinate Systems

25
Thus, combining the above three equations we get:
ΔX CT = R A RZ ( ΔAz ) Ry ( −ξ ) Rx (η ) ΔX LA
≡ fn (α , β , γ , ϕ , λ , ΔΑz , ξ ,η )
To maintain consistency among the four systems, this equation must be equivalent to
equation (4.12)
ΔX CTS = A ΔX LA
≡ fn(Φ, Λ )
The rotational values, deflection components and Λ − λ , Φ − ϕ would normally be
very small quantities. Retaining only first order terms in power series expansion of the
trig. functions, we get the results

⎡ ΔAz ⎤ ⎡ ( Λ − λ ) sin ϕ ⎤ ⎡α ⎤
⎢ ξ ⎥=⎢ ⎥ o ⎢ ⎥
Φ −ϕ ⎥ − RY (ϕ − 180 ) RZ ( λ − 180 ) ⎢ β ⎥
o
⎢ ⎥ ⎢
⎢⎣ −η ⎥⎦ ⎢⎣ − ( Λ − λ ) cos ϕ ⎥⎦ ⎢⎣γ ⎥⎦

If the G system is exactly aligned with the CT system (i.e. α , β , γ = 0o ), then the 2nd
term in the above equation disappears and we have a system of three very simple, but
very important equations.

ΔAz = Az − α = ( Λ − λ ) sin ϕ = η tan ϕ → Laplace Equation

ξ = Φ −ϕ
→ Deflection of Vertical Components
η = (Λ − λ )cos ϕ
These three equations constitute the topocentric conditions for parallelism for the CT
and G system. For parallelism all three equations must be satisfied for all points on the
earth’s surface.
An alternate derivation to the above steps can be done via spherical trig using the
following diagram

CTP
Δλ=Λ−λ
90-ϕ
Pi
α ζ

u η
(ellipsoidal
normal)
−ξ
Az z
θ

z'
(local vertical)

Figure 4.10 – Spherical Trigonometry of Parallelism Equations

26
By applying spherical trig to the various triangles above, one can derive the following
relations (see H&M, 1967, p186-190):
(1) Az − α = ( Λ − λ ) sin ϕ + (ξ sin α − η cos α ) cot ζ

(2) ξ = Φ −ϕ
(3) η = (Λ − λ )cos ϕ
(4) ζ = z + ξ cos α + η sin α or z = ζ - (ξ cos α + η sin α )
where
z, Az , Φ, Λ are (observed) astronomical observations, and
ζ , α , ϕ , λ are (computed) ellipsoidal values
α is defined as the geodetic azimuth between two normal sections. The
second term in the ΔAz equation accounts for a correction to the
observed astronomical azimuth (Az) to relate it to the same ellipsoidal
normal as the geodetic azimuth α.

If these parallelism conditions are satisfied, the transformations ϕ ↔ Φ and λ ↔ Λ

are given by relations (2) and (3) and the transformations α ↔ Az and ζ ↔ z by
relations (1) and (4).
These equations are most commonly employed for reducing observed quantities (Az
and z) to corresponding ellipsoidal quantities (α and ζ).

27
 1.6.1 Introduction
z A Geodetic Datum is a framework that enables
us to define coordinate systems.
5. The geodetic datums and » A reference ellipsoid is required to which the points
Geocentric Datum of Australia are referred.
» Coordinates may be used as
5.1 Introduction – 3D Cartesian coordinates X,Y and Z (no requests for
ellipsoid definition)
5.2 The AGD – Geographical coordinates φ, λ and h
5.3 The GDA – Map projection coordinates E, N and H

5.4 The Australian Geoid - AUSGeoid98 » Examples:

– AGD has ANS associated with it, and
– GDA has GRS80 associated with it.
1 2

Introduction (Cont.) Introduction (Cont.)

z An ellipsoid as a permanent basis for a z 8 parameters are required for positioning

geodetic coordinate system will be a ellipsoid
included as part of datum establishment. » The semi-major axis (a) and flattening (f)
z The shape and size of the ellipsoid must define the shape and size (of the Earth).
be selected, and its position in the Earth » 3 parameters define its position.
has to be located and oriented. » 3 parameters orient it with respect to 3
Cartesian axes.

3 4

Introduction (Cont.) Introduction (Cont.)

z Geocentric Datum
z Regional/Local Datum » It became apparent that GPS can measure the positions
» Background more accurately than conventional surveying methods.
– Lower precision of geoid (few metres) at the time. » GPS coordinates are related to the centre of mass of the
– Most geodetic surveying was done using earth (ie. geocentre).
conventional terrestrial methods » Better geoid models
» The ellipsoid is fitted to the regional geoid as » The ellipsoid is fitted to the geoid best in a global scale.
closely as possible to reduce effects of N, ξ » Eg. WGS84 and GDA94.
and η on observations. z However, the Geocentric Datum may not locally fit
» Eg. AGD66, AGD84 the ellipsoid well.
5 6

1
 1.6.2 Local Ellipsoids and Datums Local Ellipsoids and Datums (Cont.)

z As stated, a regional reference ellipsoid is z As the geoid is an undulating surface, the

chosen so as to fit the geoid, as closely best-fitting ellipsoid in one region is not
as possible, over the area to be mapped. necessarily the best fit in another.
z This approach enables subsequent z Consequently, different reference
geodetic data, which are collected on the ellipsoid s have been defined in different
physical surface of the earth, to be regions.
reduced to the surface of this ellipsoid » The Australian National ellipsoid (ANS) is
without the introduction of significant used in Australia
horizontal scale error. » The Airy 1830 ellipsoid is used in Great
Britain.
7 8

Local Ellipsoids and Datums (Cont.) Local Ellipsoids and Datums (Cont.)

z Regional ellipsoids have generally been z Once the best fitting regional ellipsoid has been
defined, it is adopted as the reference surface
established using astronomical for geodetic positions in that region.
observations to define the deflection of z The available geodetic observations are then
the vertical to be zero. adjusted, in a least squares sense, to form the
z The orientation and scale of the ellipsoid regional datum based on that reference
ellipsoid.
is defined using further geodetic » In Australia the AGD66 and AGD84 datums are
observations. based on the ANS.
» In Great Britain, the OSGB36 and OS(SN)80 datums
are based on the Airy 1830 ellipsoid.
9 10

Summary 5.2 The AGD

z For Australia, the geodetic datum (horizontal)
z A regional or geocentric datum comprises was previously defined by the Australian
geographical coordinates that are Geodetic Datum (AGD).
referenced to the surface of a particular » The origin of this network is at Johnston Geodetic
Station where
reference ellipsoid.
z The minor axis of the AGD ellipsoid is termed
z Therefore, a single ground point can have the Australian National Spheroid (ANS) and
different geographical coordinates by has parameters of:
» a = 6,378,160 m
virtue of the datum and ellipsoid used.
» f = 1/298.25

11 12

2
 The Origin Point of AGD The AGD
• S 25° 56' 54.5515",
• E 133° 12' 0.0771" z Coordinates of points in this geodetic
h 571.2 m (ellipsoid height) datum where given in either ellipsoidal or
z The N (geoid-ellipsoid grid values:
separation) value at
» AGD84 φ, λ, h coords
Johnston was taken as:
» +4.9 m (for AGD84) » AMG84 grid coords
» 0 m (for AGD66)
z With deflection of vertical
components of
Johnston Geodetic Station
13
ξ o = −4.2′′ and ηo = −4.5′′ Stone marker 14

Establishing A Regional/Local Datum (1) Establishing the Origin Point

z Observations are carried out at the point to determine φA
z Objectives and λA.
» Find the most appropriate parameters a and f for the z Geodetic latitude (φG) and longitude (λG) are set equal to
region of interest φA and λA.
» Fit the geoid as closely as possible to this ellipsoid z Measure astronomical azimuth, AA and constrains AA= AG
z The Origin Point of a local datum is a station » This will define the ellipsoid’s Z axis to be parallel to the rotational
where axis of the Earth at the time of observation).
» The coordinates are fixed z Above 3 conditions satisfy the Laplace’s Equation
» The ellipsoid and geoid coincide AG= AA – (λA - λG) sin φ
» Astronomical and geodetic coordinates are equal.
z Observed longitude gives the position of X, and then Y.
» For AGD66. The origin point is Johnston Geodetic
Station in the Northern Territory z But, this does not make that the ellipsoid fits the geoid over
a continent.
15 16

(2) Determining the Deflection of the Vertical (2) Determining the Deflection of the Vertical
over the Continent over the Continent (Cont.)

z Laplace and other astronomical observations z Using a least squares adjustment, we can
are measured over the continent. position and orient the ellipsoid through
z Then, at each station, we obtain the minimising ξ and η over the entire continent.
components of deflection of the vertical: z A similar processing can be used to minimise
ξ = φA - φG the geoid spheroid separation, N.
η = (λA - λG) cos φ z This is the main method to establish AGD66.
ξ and η can be obtained at all stations. z But, the orientation of the spheroid may be a bit
z The values of ξ and η are relative to the datum skew!
chosen, but the geoid does not changed.
17 18

3
 (3) Determining Orientation Parameters of
Laplace Station the Geodetic Datum

z A Laplace station is a survey station at which z Now we use Laplace stations to constrain
astronomical observations have been made the azimuth in the adjustment of the
and used in Laplace's equation to determine geodetic network.
the deviation of the vertical (the difference » This improves the precision of the network.
between the local vertical and the normal to the » This also defines the Z axis, so that it is
spheriod of reference). parallel to the Earth’s rotational axis.
z Within triangulation networks, Laplace stations z The Z axis is actually parallel to the mean
help to control the errors in scale and azimuth rotational axis of the earth related to a
that would otherwise accumulate. particular epoch!!
z Australia has ~275 Laplace stations for AGD66.
19 20

(3) Determining Orientation Parameters of

the Geodetic Datum (Cont.)

z The earth’s rotational axis is in motion, gradually

changing the pole positions in space.
» The motion is due to the phenomena of precession and
nutation caused by the mutual attractions among the
earth, sun, moon and planets.
» The polar motion corrections have to be applied to
astronomical observations at Laplace stations and
reduced them to the Conventional International Origin
(ie. mean pole of the earth’s rotational axis).
z The procedure defines two of three orientation
21
Distribution of 275 Laplace Stations used in AGD66
22 parameters of the geodetic datum.

(4) Determining the Orientation Parameter of (5) Determining a and f Parameters of the
the Geodetic Datum Geodetic Datum

z The third orientation parameter is the

z The size and geometrical shape of the ellipsoid
position of the X axis of the ellipsoid with are defined by its semi-major axis (a) and
respect to the Earth. flattening (f=1-b/a).
» The X axis passes through Greenwich
meridian. z The a and f can be determined using either
astro-geodetic or satellite means.
» This can be done by astronomical
determination of λ at the origin point, but the z Note: we do not determine other 3 parameters
precision would be lower. defining the position with respect to the
» So, the longitude was usually measured at an geocentre.
observatory with high quality instruments (eg. » Local or regional datums are usually not geocentric
Mount Stromlo observatory in Australia). datums!
23 24

4
 The AGD66 The AGD66

z The least squares adjustment of the z In AGD66, a local ellipsoid, Australian

Australian geodetic network performed in National Spheroid (ANS), was used to
March 1966 used the Australian Geodetic define the adopted size and shape of the
Datum. earth, and the position of the origin point
z This adjustment produced a set of (ie. Johnston Geodetic Station).
coordinates which, in the form of latitudes
and longitudes, is known as the z ANS has:
Australian Geodetic Datum 1966 » a = 6,378,160 m and f = 1/298.25
coordinate set (AGD66).
25 26

The AGD84 5.3 The GDA

z In 1982 a new national adjustment, referred to as the z In 2000, the Australian Geodetic Datum (AGD)
Geodetic Model of Australia 1982 (GMA82), was
performed using all data previously included in the has been replaced by the Geocentric Datum of
1966 adjustment as well as additional, modern Australia (GDA) as the official Australian
terrestrial and space-based observations. coordinate system.
z This new adjustment also used the gazetted Australian
Geodetic Datum.
z GDA is a 3-D, geocentric reference system,
z The coordinate set resulting from this adjustment was defined within the ITRF and compatible with the
accepted by the National Mapping Council in 1984. WGS84 system.
z This is known as the Australian Geodetic Datum 1984 z For full details of the GDA and its
(AGD84).
implementation see
http://www.ga.gov.au/geodesy/datums/gda.jsp
27 28

GDA Specifications Horizontal Datums in Australia

TERMINOLOGY
z The Australian Map Grid (MG84 and AMG66) is a UTM
Datum: Geocentric Datum of Australia (GDA) projection of AGD geographical coordinates from the ANS
Geographical coordinate set: spheroid.
Geocentric Datum of Australia 1994 (GDA94)
(latitude and longitude)
Map Grid of Australia 1994 (MGA94)
z The Map Grid of Australia (MGA) is a UTM projection of
Grid coordinates : (Universal Transverse Mercator, using the GDA geographical coordinates from the GRS80 spheroid.
GRS80 ellipsoid)
DEFINITION datum ellipsoid Map grid a (m) 1/f
Reference Frame : ITRF92
Epoch : 1994.0 AGD66 ANS AMG66 6378160 298.25
Ellipsoid : GRS80 AGD84 ANS AMG84 6378160 298.25
Semi-major axis (a) : 6,378,137.0 metres
Inverse flattening (1/f) : 298.257222101
GDA94 GRS80 MGA94 6378137 298.257222101
29 30

5
 What is the difference between
GDA and AGD?
z GDA is based on a global
spheroid, which best fits to the
geoid over the entire earth.
» Its origin coincide with the centre of
mass of the Earth.
z AGD uses a regional reference
spheroid, which best fits to the
geoid for Australia only.
» Its origin does not coincide with the
centre of mass of the Earth.
z GDA’s coordinates on the Earth’s
surface depart from those from
GDA ~ 200 m north east.
31 32

33 34

Background to GDA Background to GDA (Cont.)

z In 1992, eight GPS permanent stations across Australia z In 1988, the Inter-governmental Committee on Surveying
formed the Australian Fiducial Network (AFN). and Mapping (ICSM), which was established following the
» as part of the world-wide International GPS Service (IGS) campaign.
cessation of the National Mapping Council (NMC), initiated
the conversion from the regional AGD66 and AGD84
z During 1992-1994, 70 GPS sites with a 500 km spacing datums to a geocentric datum (ICSM 1990, 1991a, 1991b).
across Australia formed the Australian National Network
(ANN). z The ICSM resolved to adopt the Geocentric Datum of
Australia 1994 (GDA94) at its November 1994 meeting.
z Combined both GPS observations in a single regional GPS » The gradual implementation of this datum Australia-wide by 1st
solution January, 2000.
» in terms of ITRF92. » The GDA is intended to supersede the AGD66 and AGD84
» to an epoch of 1994.0. » The GDA will lead to a homogeneous national coordinate datum for
future mapping.
» position accuracy: AFN ±2 cm; ANN ±5cm.
z Positions of the AFN sites were used to define the GDA94.
35
z GDA94 was released on 6 September 1995. 36

6
 Background to GDA (Cont.) Geocentric Datum of Australia (Cont.)

A global sub-network of accurately

z z The AFN is a continent-wide GPS network
positioned ground stations, tied to the of eight accurately positioned geocentric
International Terrestrial Reference Frame stations, which is, in turn, tied geodetically
(ITRF), has been established using GPS, to the IGS.
and is called the International GPS » The AFN Accuracy: one part in one hundred
Service for Geodynamics (IGS). million or a few cm at stations.
z In 1992, the Australian Surveying and Land
Information Group (AUSLIG) took the
opportunity to simultaneously establish the
37 Australian Fiducial Network (AFN). 38

Australian Fiducial Network (AFN) Geocentric Datum of Australia (Cont.)

z The AFN consists of 8

permanent, z The Australian National Network (ANN)
continuously operating is a further sub-network of the AFN, and
Geodetic GPS uses GPS-derived positions at
receivers.
z These stations are
approximately 500-kilometre intervals
accurate to 0.02 - 0.04 over the continent.
ppm. » The ANN accuracy: one part in ten million.
z AFN (and hence
GDA94) station z The AFN and ANN coordinates form the
coordinates are based backbone of the new Geocentric Datum
on the ITRF92 at of Australia.
epoch 1994.0 .
39 40

Australian National Network (ANN) GDA94

z The ANN network
of 78 sites is z At mean time, state and Territory GPS
constrained to the networks were established and
AFN. connected to the ANN.
z AFN and ANN
form the z All these possible GPS observations
framework for the combined with observations of Australian
GDA. Geodetic Network are constrained to the
z It has an ANN and AFN.
estimated
accuracy of z This established the GDA94 throughout
0.1 ppm. Australia.
41 42

7
 Summary - GDA94 Global Geocentric Ellipsoid – WGS84

z The GRS80 is the ellipsoid associated with the z The most recent global geocentric ellipsoid that is
national geocentric datum GDA94. widely used is the World Geodetic System 1984
» f=1/298.257222101 (WGS84).
» a=6378137 m » The WGS84 ellipsoid is slightly different from the
z Its X, Y and Z axes coincide with axes of GRS80 ellipsoid:
ITRF92. – f=1/298.257223563 (WGS84), while
– f=1/298.257222101 (GRS80)
» A satellite datum of the ITRF provides for a
continental geocentric datum. » This slight difference in flattening is due to rounding
» Coordinated determined from this datum are related errors.
to the epoch 1994.0 » However, the effect on three-dimensional coordinates is
z The GDA94 has been used since 1 Jan 2000. less than a millimetre and can therefore be neglected.
43 44

Global Geocentric ellipsoid – WGS84 Terrestrial Reference Frame (ITRF)

z A global network of accurate coordinates (and their

z The GPS provides positions that are velocities) having the earth’s centre of mass as its
referenced to the WGS84 ellipsoid. origin.
z ITRF is maintained by the International Earth
z All GPS receivers compute and store Rotation Service.
coordinates in terms of WGS84, then z The ITRF is positioned relative to the geocentre
transform to other datums when using a variety of space-based geodetic
techniques, such as
information is displayed. » Satellite Laser Ranging (SLR),
z But in GDA94, the ellipsoid used is » Very Long Baseline Interferometry (VLBI),
» GPS, and
GRS80! » Doppler Orbitography and Radio positioning Integrated
by Satellite (DORIS).
45 46

ITRF and WGS84

z The ITRF92 is the reference frame of the

Geocentric Datum of Australia 1994.
z The ellipsoid used in GDA94 is GRS80.
z In fact, the WGS84 ellipsoid was modified in
early 1994 and is now coincident with the
ITRF92 at the 10cm level.
z Therefore, for all practical purposes, the
GDA is fully compatible with GPS in terms of
ellipsoid and datum.
47

8
7. Map Projection
7.1 Introduction
The fundamental coordinate reference system for surveying and mapping is a set of geodetic
coordinates related to a particular datum. It is then necessary to consider how to arrange the data so
that it can be placed on a flat surface.

There are two reasons doing this.

1. For presentation. Whether the data is to be shown on a paper map or on a computer screen,
it must of necessary be presented in a 2-D format.
2. For computation. Even a simple concept such as the distance between two points becomes
complex when expressed in ellipsoidal formulae, and wherever possible it is more desirable
to carry out computation in a simple 2-D Cartesian system.

A map projection, then, is defined as an ordered system of meridian and parallels on a flat surface,
or

A map projection is a mathematical algorithm to transform locations defined on the curved surface
of the earth into locations defined on the flat surface of a map.

It should be immediately realised that it is impossible to convert a sphere or an ellipsoid into a flat
plane without in some way distorting or cutting it.

It then follows that there is no single method for doing this, and hence the classification of types of
map projection:

• Conic: projection surface is a conic surface. Lamp at the center of the earth.
o Examples: Albers Equal Area, Lambert Conformal Conic. Good for East-West
land areas.

1
• Cylindrical: projection surface is a cylindrical surface. Lamp at the center of the earth.
o Examples: (Transverse Mercator). Good for North-South land areas.

Tangent Secant

Transverse Oblique

• Azimuthal: projection surface is a flat surface tangent to the earth. Lamp at the center of the
earth (gnomonic), at the other side of the earth (stereographic), or far from the earth
(orthographic).
o Examples: Lambert Azimuthal Equal Area. Good for global views.

2
7.2 Fundamental Concepts
7.2.1 Grids and graticules
Meridians and parallels appear on maps depended on the type of projection, but in general they
are in arrangement of straight or curved lines.

In Fig.1, the set of parallels and meridians is known as the graticule.

Fig. 1 An example of the graticule

• This graticule does not constitute the basis of a coordinate reference system.
• The graticule does not show on some maps.

In Fig.2, the set of straight lines related to a rectangular coordinate reference system is known as
the grid.

Fig.2 A grid superimposed on the graticule

This grid is given coordinates that may be labelled y and x, or easting (E) and northing (N), or
some other convention.

3
7.2.2 Scale factor
Features on the surface of a sphere or an ellipsoid will be distorted when projected onto a plane.
So, it is necessary to have a precise definition of the amount of distortion that has resulted.

The scale factor, k, is defined for this purpose.

distance on the projection

k= (7.1)
distance on the globe

• k is different at each point on the projection.

• k in many cases have different values in each direction.
• Eq. (7.1) in general can apply to only a short distance (in theory infinitesimal short).
• For longer lines, the relevant parameter is the integrated mean of the point scale factors along
the whole length of the line.
• It is important to understand that k is unrelated to the scale of map (a ration such as 1:50 000),
because k results purely from the act of projecting to a flat surface.
• The ideal value is k=1, meaning no distortion.
• The distortion here is not the same as an ‘error’ in the map.

In general, the distortion in the direction of the parallels will be different to the distortion in the
direction of the meridian.

Let kP represent the scale factor along a parallel, and kM represent the scale factor along a meridian.
A small square (1x1) is then projected as rectangle of dimensions (kP×kM).

We define two scale factors:

1. Point scale factor is the ratio of an infinitesimal distance in a particular direction on the map
to the corresponding distance on the globe.
2. Line scale factor is the ratio of the length of a projected (finite) line to the length of the
corresponding line on the globe.

7.2.3 Geometric form of the projection surface

Historically, projections were derived by first projecting from the ellipsoid to an intermediate
surface.
• The surface has a nature that features can be recovered without distortion.
• Such a surface is called as the projection surface.

There are three main types:

• Azimuthal or Zenithal, with a plane as projection surface.
• Cylindrical, with a cylinder as projection surface.
• Conical, with a cone as projection surface.

4
7.2.4 Preserved features
It is not possible to devises a projection without introducing distortion. The usual approach is to
preserve one of features, such as the shape, area and size on the surface of the ellipsoid.

For examples, some preserved features are:

• Equidistant projection: kM=1. The distances along all meridians should remain undistorted.
• Equal area projection: kMkP=1. The area of the projected unit square remains equal to one.
o e.g., Lambert’s Cylindrical Equal-area projection.
• Orthomorphic or conformal projection: kM=kP. The shape feature remains undistorted.
o The conformal projection is also preserving angles.
o This means that angles measured on the ground can be transferred to the projection
for use in computations.
o It is thus one of frequently used and almost all large scale mapping.
o e.g., Mercator projection.

In summary, most projection methods may be classified according to:

1. the shape of the projection surface, which is dictated primarily by the geographical area to be
mapped but also in part by the function of map, and
2. the features on the sphere (or ellipsoid) that are to be preserved on the projection.

7.2.5 Sphere and ellipsoid

As stated, the fundamental coordinate system is an ellipsoidal one related to a geodetic datum using
an ellipsoid. The relative positions of points in such a coordinate system are different to those on a
sphere.

However, in discussions followed, diagrams and equations, which are derived from a sphere with
radius (ie. R=1 or R= mean radius of the Earth), are sometimes used to give approximate
indications the scale factor and so on.

There are some reasons to do this.

1. The flattening of most ellipsoids is <1/300. Thus, the shape of a particular projection when
applied to the sphere is very similar to the shape when it is applied to the ellipsoid.
• The significant difference will appear on a large scale map, but the sphere is useful for
finding how the resulting map has been distorted.
2. The ellipsoid is not a conformable surface.

7.3 Cylindrical Projections

A cylindrical projection can be imagined in its simplest form:
• A cylinder has been wrapped around a globe at the equator.
• If the graticule of latitude and longitude are projected onto the cylinder and the cylinder
unwrapped, then a grid-like pattern of straight lines of latitude and longitude would
result.

5
 Fig. 3 Cylindrical projection surface (normal case)

After projection:
• The meridians of longitude would be equally spaced.
• The parallels of latitude would remain parallel but may not appear equally spaced
anymore.

In reality cylindrical map projections are not so simply constructed.

There are three aspects of the cylindrical projections:

• Tangent or secant to equator is termed regular or normal.

• Tangent or secant to a meridian is the transverse aspect.
• Tangent or secant to any point on the globe is called oblique.

Fig. 4 Transverse cylindrical projection surface

6
 Fig. 5 Oblique cylindrical projection surface

7.3.1 Normal cases

(1). Gall’s Projection

Gall's stereographic cylindrical projection results from projecting the earth's surface from the
equator onto a secant cylinder intersected by the globe at 45 degrees north and 45 degrees south.
This projection moderately distorts distance, shape, direction, and area.

7
Scale factors:

1 + cos 45o
kM =
1 + cos φ
cos 45o
kP =
cos φ

(2) Mercator’s Projection

The Mercator projection has straight meridians and parallels that intersect at right angles.

It is the conformal version. The projection looks like:

Scale Factors (cf. course notes on p.116-117):

dX d   π φ 
kM = = ln tan +  = secφ
dφ dφ   4 2 
2π
kP = = sec φ
2π cosφ

The general features of Mercator are:

1. kM=kP=secφ at the any point and in any direction.

2. In consequence, the pole is now of infinite size and at infinite distance from the equator, and
hence cannot be represented on the projection.

8
The meridians are parallel to each other and that the angles are preserved. This fact makes it an
ideal projection for navigation.

A line drawn between 2 points A and B has a constant azimuth wrt the meridians, which can
be read directly from the map.

o But this line is not the shortest line due to the variation of scale factor within the
projection. The shortest line between A and B should be the great circle, which
should be projected as a curved line.

Fig.6 A line of constant azimuth from A to B on a MT

3. An individual map sheet will usually represent a small part of the world (because
kM=kP=sec90°=infinite).

9
7.3.2 Transverse Cases

So far, the cylindrical projections discussed were formed by placing a cylinder in contact with the
equator. Because they are often used to plot the Earth as a whole, they are optimised for use in the
equatorial regions.

For those parts of the Earth that do not lie close to the equator, an alternative is used.

Transverse Mercator projections turns the cylinder onto its side and make it contact (or tangent to)
a particular meridian known as the “central meridian” (CM).

Note: For Surveyor’s projection and mapping purposes, we project from a sphere of radius R
(ie. the mean radius of the Earth) instead of the R=1 to obtain projection coordinates (in m) when
dealing with Transverse cases.

A set of rules can again be proposed to produce equal area, equidistant, or conformal projections.
The most important of these is the Transverse Mercator projection.

(1). Transverse Mercator projection

The central meridian is on 0 degree of the longitude of origin.

The important features of this projection are:

1. The projection is conformal.

2. The scale factor at each point is the same in any direction.
3. The meridians are no longer parallel to each other, and in fact are no longer straight lines.

The exception is the central meridian, which is the straight line.

At any general point, the meridian makes an angle with the central meridian, which is the
direction of grid north. This angle is termed the grid convergence, γ.

For the sphere, γ can be computed by

γ ≈ (λ − λcentre of zone ) sin φ (in second)

10
 Fig.7 Examples of varying values of convergence within a TM projection. Note that true north is
always towards the CM, and that the magnitude of convergence increases with distance from the
CM.

The Transverse Mercator projection is commonly used for a narrow band of no more than ±3° on
either side of the central meridian.

This is because the scale factor distortion in the E-W direction increases with distance from
the CM. In a surveying projection, it is required to minimise the scale factor distortion. This
is done by limiting the projection to a narrow N-S ‘strip’ symmetrical about the CM.

The strips is termed the ‘zones’ and <6°.

When using a 6° zone, Fig. here shows that the scale factor varies between 1 on the CM and
1.0014 on the edge of the projection (black curve in Fig.8).

To minimise the scale factor distortion, we use a central meridian scale factor. A typical
value of medium scale topographic mapping would be 0.9996.

This means that the scale factor across the projection ranges from 0.9996 on the CM to 1.0010
(=0.9996sec3°) on the edge (red curve in Fig.8)

11
 Fig.8 Scale factor for TM projections

When an ellipsoid rather than a sphere is used as the model of the Earth, the formulae for
converting ellipsoidal coordinates to grid coordinates become more complex.

The ellipsoidal version is called the Gauss-Krüger.

Ellipsoidal coordinates are first converted to spherical coordinates before being further
converted into Cartesian grid coordinates.

Transverse Mercator maps are often used to portray areas with larger north-south than east-west
extent. Distortion of scale, distance, direction and area increases away from the central meridian.

(2). Universal Transverse Mercator (UTM) Coordinate System

The TM zone concept may be extended to the whole Earth. It is the basis of a world-wide
projection system known as UTM.

The UTM projection divides the surface of the Earth into 6 degree zones, each mapped by the TM
projection with a central meridian in the centre of the zone.

• The zones are numbered from 1 starting at λ=180°E, and increase eastwards (60 zones in
total).

12
• Each zone is divided into 2, for northern and southern hemisphere. This is to ensure that
coordinates in the southern hemisphere can be made positive through introducing a false
northing.

• Within each zone, the origin λ lies on the meridian along the centre of the zone.

o UTM zone 1 has the central meridian 177°W

o UTM zone 2 has the central meridian 171°W

• The true origin of each zone is at the intersection of the CM with the equator. But this makes
some coordinates (E) negative.

To avoid negative coordinates, false origin coordinates are used when φ<0° (e.g. in Australia)
as:

• False easting = 500,000 m

• False northing = 10,000,000 m
• Scale factor at origin = 0.9996

• All parameters of the UTM system are listed in Table below:

13
 Table projection parameters and their values for global zoned TM systems

λ0):
Computing Zone Number (N) and CM (λ

N = integer part of [(180+λ)/6+1]

λ0 = 6(N-1)-177

14
 Grid Coordinates

Geodetic coordinates (φ, λ) are represented on a map, by
Map Projection 2 mathematically "projecting" them onto a surface.

The UTM adopted uniform scale factor, false origins and
zone size and numbering for the TM projection.
By Dr X Deng for SURV3510
The globe is divided into narrow 66 longitude zones, which
are projected onto a Transverse Mercator projection.

1 2

UTM Grid Coordinates UTM Grid Coordinates

Expressed as a distance
Expressed as a distance in
in meters to the east, meters to the north, referred
referred to as the to as the “northing“ (N).
"easting“ (E).
Measured relative to the
equator.

Min & max values are: » For locations north of the
» 160,000 and 834,000 m at equator the equator is
equator assigned N=0 m.
» For locations south of the
» 465,000 and 516,000 m at
equator the equator is
84°N. assigned N= 10,000,000 m.
3 4

UTM Grid Coordinates Latitude Zones

UTM zones extend from
Latitude is also divided into zones, but less regularly.
a latitude of 80° S to 84°
Zones are lettered from A at the South Pole to Z at the North.

The circle south of 80° is divided into two zones, A and B.
N. Thereafter zones are 8° wide.

In the polar regions the
Zone M is just south of the equator and N is north.
Universal Polar
Zone T, between 40 and 48 degrees north

Zone X, from 72°N to 84°N, is 12° wide.
Stereographic (UPS) grid

Zones Y and Z cover the north polar region north of 84°.
system is used.
I and O are not used to avoid the confusion with numbers.

5 6

1
 Grid Coordinates (Cont.) Grid North

Since meridians converge, the UTM grid

coincides exactly with true north only along the
central meridian of each zone.

Elsewhere, the grid makes an angle to the
meridians, especially at the edges of the zone.

Therefore:
» The true north is along the meridian directions.
» The grid north is along the north-south lines of the
grid.

7 60 longitude zones, 6º wide, zone 1 at 180ºW. 8

Grid Convergence Grid Convergence

The spherical angle between the grid and true
The grid convergence, γ can be
meridians at the point is called the “grid
convergence, γ”.
calculated using:

This angular quantity is to be added to an azimuth
to obtain a grid bearing: (
tan γ = tan λ − λcentre of zone sin φ )
» Grid Bearing = Azimuth + Grid Convergence.

In the southern hemisphere, grid convergence is
» positive for points east of the central meridian (grid north
After some approximation, it becomes.
is west of true north) and
» negative for points west of the central meridian (grid north
is east of true north).
( )
γ " ≈ λ − λcentre of zone " sin φ
9 10

Australia is covered by zones 49 - 56 Horizontal Datums in Australia

The Australian Map Grid (AMG84 and AMG66) is a UTM
projection of AGD geographical coordinates from the ANS
spheroid.

The Map Grid of Australia (MGA) is a UTM projection of
GDA geographical coordinates from the GRS80 spheroid.

datum spheroid Map grid a (m) 1/f

AGD66 ANS AMG66 6378160 298.25

AGD84 ANS AMG84 6378160 298.25

GDA94 GRS80 MGA94 6378137 298.257222101
11 12

2
 Australian National Grid (ANG) Australian National Grid (ANG)

Prior to 1966, "the rectangular grid
Parameters used in ANG are:

coordinate system", used in conjunction » Projection: Transverse Mercator
with the Clarke 1858 spheroid, was called » Zone Width: 5 degrees
the Australian National Grid (ANG). » Longitude of Origin: CM of each zone
» Latitude of Origin: Equator (0 degrees)
» False Easting: 400 000 yards
» False Northing: 4,915,813.467 yards
» Central Scale Factor: 1.0
» Ellipsoid: Clarke 1858
13 14

Australian Map Grid (AMG66) Integrated Survey Grid (ISG)

The ISG is a projection used in NSW only

In March 1966, the least squares adjustment (produced in 1976).
of the Australian geodetic network performed » It was introduced to minimise scale factor corrections
(mainly in cadastral surveys).
using the AGD. » Coordinates are derived from the AGD66
» This adjustment produced a set of coordinates of
Parameters used in ISG are:
(φ, λ), known as the AGD66 coordinate set. » Projection: Transverse Mercator
» Zone Width: 2 degrees

The grid coordinates derived from a UTM » Longitude of Origin: CM of each zone
projection of the AGD66 coordinates, using » Latitude of Origin: Equator (zero degrees)
» False Easting: 300 000 m
the ANS, is now known as the Australian » False Northing: 5 000 000 m
Map Grid 1966 coordinate set (AMG66). » Central Scale Factor: 0.99994
15 16 » Ellipsoid: ANS

Australian Map Grid (AMG84) Australian Map Grid (AMG84)

In 1982 a new national adjustment was
Coordinates on the AMG are derived from a
UTM projection of (φ, λ) on the AGD.
performed using all data previously included in

Parameters used for the AMG are the same as
the 1966 adjustment as well as additional, AMG66:
modern terrestrial and space-based » Projection: Transverse Mercator
observations. » Longitude of Origin: CM of each zone

This new adjustment used the gazetted AGD, » Latitude of Origin: Equator (0°)
» Zone width: 6°
and produced the coordinate set of the AGD84. » Central scale factor: 0.9996

The equivalent UTM grid coordinates, projected » False easting: 500,000 m
using the ANS, are known as the Australian » False northing (φ<0°): 10,000,000 m
Map Grid 1984 (AMG84). » Ellipsoid: ANS a = 6,378,160 m 1/f = 298.25
17 18

3
 Map Grid of Australia (MGA94) Map Grid of Australia (MGA94)

The UTM projection was used with the
Parameters of UTM are:
GRS80 ellipsoid and GDA94 (φ, λ) to » Projection: Transverse Mercator
produce grid coordinates, known as the » Longitude of Origin: CM of each zone
Map Grid of Australia 1994 coordinates » Latitude of Origin: Equator (0°)
(MGA94). » Zone width: 6°
» Central scale factor: 0.9996
» False easting: 500,000 m
» False northing (φ<0°): 10,000,000 m
» Ellipsoid: GRS80
19 20

Conversion between (φ
φ, λ) and (E,N)

Redfearn's formulae are used to convert

between (φ, λ), and AMG grid coordinates
(E,N).
» These formulae are shown with worked
examples in Chapter 5 in “GDA Technical
Manual”. Questions?

Note that Grid Calculation formulae can also

be found in Chapter 6 in “GDA Technical
Manual”.
21 22

4
 SURV3510 Topic: Introduction to the UTM and Map Grid of Australia

The University of Newcastle

Faculty of Engineering and Built Environment
School of Engineering
SURV3510 Geodesy
Introduction to the UTM and Map Grid of Australia

Notes and References

Reference: Wikipedia: see Mercator Projection, Universal Transverse Mercator.

A map projection is a way of creating plane easting and northing coordinates from ‘curvilinear’
latitudes and longitudes. Different map projections have different properties! Notice: surveyors don’t
really have to be creating a paper topographic map to want a map projection.

The Mercator map projection system mathematically projects geodetic coordinates (latitude /
longitude) onto a cylinder which makes contact along the equator, to create grid coordinates (e.g.
easting / northing). The Mercator projection can cover much of the earth with a single ‘map’, but the
distortions can get large.

The Transverse Mercator system projects geodetic

coordinates onto a concentric cylinder which makes
contact along one (central) meridian. Note: in surveying,
the cylinder is wrapped around an ellipsoid (not a
sphere). This projection covers a ‘strip’ of the earth,
from pole to pole, but then it uses different strips to cover
the earth: see diagram.

Shape distortion exists, and increases with distance from the central meridian. To minimise distortion,
the cylinder is "rotated" to create different zones, bringing a different central meridian into contact
with the ellipsoid. The width of the bands (and therefore the number of zones) can be chosen.
Distortion increases with zone width, but on the other hand narrow zones increase the complexity by
having more zones and more crossings to deal with.

This projection is used around the world, and is called Universal TM, with zones 6° in width,
numbered to east from Zone 1 at 3° longitude. Some zones relevant to Australia are:

Central Meridian Zone Number

135° 53
141° 54
147° 55
153° 56

TM Characteristics (crucial)
• There is a scale factor at any point. This is not map scale; this is the scale between ellipsoid
and map, whereas map scale is between map grid and ‘paper’. This point scale factor is not
constant, and varies, mainly in the direction across the zone and the variation is greater with
greater zone widths.
• But the scale factor does not vary with direction at the point, meaning that angles on the
ellipsoid (or the ground?) are the same angles on the grid. This is important for surveyors
1
 SURV3510 Topic: Introduction to the UTM and Map Grid of Australia

doing calculations of eastings and northings with their angles. In a small area the UTM is
conformal.
• Note: meridians on the ellipsoid are not north-south on the grid (except the central
meridian), i.e. grid north is not parallel to true north, azimuths are not bearings (but their
differences are correct: angles on the ellipsoid are the same angles on the grid).
• Over a larger area, shapes are distorted, and long lines on the ground are not straight lines
on the projection. Along a long line, the point scale factor can vary, and the average scale
factor is a line scale factor.

The UTM is used for the Map Grid of Australia (MGA), as follows.
• The origin for each zone is the intersection of the equator and the contacting meridian (the
central meridian), but a false origin is often used to avoid negative coordinates. False origin
is at -500,000 m east and -10,000,000 m north.
• The scale factor is normally exactly 1.000 on the central meridian, and increases away from the
CM, but for the MGA it is set at less than 1.000, so that the average scale across the zone is
close to 1.000. The central scale factor k0 (scale factor on the central meridian) is set at 0.9996.

The UTM was used for the old AMG, which has
• false origin values of 500,000 east and 10,000,000 north
• zones of 6 degrees
• central scale factor: 0.9996

The TM was also for the old ISG, which has

• false origin values 300,000 and 5,000,000 m
• zones of 2 degrees with ¼° overlap
• central scale factor: 0.99994. Note that this is closer to 1.000 than is the MGA or AMG value,
and is made possible by the narrower zone.

Calculations to interchange coordinates between ellipsoid and grid are covered in a later Topic.

Tutorial Questions
Bold font questions will be taken up in the tutorial class.

1. Explain the two most important characteristic properties of the UTM map projection for
surveyors? Do other projections have these properties? (This can take a paragraph).

2. Why do we need a map projection? Why not just work with latitudes and longitudes?

3. Why do we have latitudes and longitudes anyway? Why not just have a system of
eastings and northings, without a basis in latitudes and longitudes?

4. What are the zone widths on the MGA? Is this a good value? 6°

5. Give the longitude and zone number of one central meridian on MGA. e.g. No.56: 153°

6. How do you obtain true MGA coordinates from MGA based on the false origin?
Subtract East 500,000 and North 10,000,000.

7. What is the central scale factor for the MGA? 0.9996

Why is it not 1.0000? (This last part can take a paragraph).

2
 8 Curves on the Surface of the Ellipsoid

8.1 The Normal Section and its Properties

All normals to the surface of a sphere contain the centre of the sphere and therefore, in
this case, a normal plane at point A containing point B also contains the normal at B.
However on an ellipsoid, the normal planes at A and B do not generally intersect (unless
A and B are on the same meridian or parallel) because the ellipsoidal normals at A and
B are skewed and are not co-planar.
We define the (forward) normal section on an ellipsoid as the curve of the intersection
of the vertical plane, containing the ellipsoidal normal at a point A and the projected
point (B1) on the ellipsoid of a point B at some elevation above the ellipsoid, with the
ellipsoidal surface (see Figure 8.1). If the target and theodolite stations were switched,
we would define another plane (containing the ellipsoidal normal at point B and the
projected point A1 on the ellipsoid) and therefore another (or the reverse) normal
section (between B1 and A1) where this plane intersects with the ellipsoid surface. That
is, the two normal planes containing A1, B1 and their respective ellipsoidal normals are
inclined to each other and cut the surface of the ellipsoid in different curves (see Figure
8.1 inset).

target station
theod. station

normal section at A

points on ellipsoid

B1
ellipsoid normals A1

N
M

Figure 8.1 – Normal Section Definition on the Ellipsoid

A theodolite levelled at point A with the vertical axis along the ellipsoidal normal NA1A
(ignoring the deflection of the vertical) directed towards a target at B defines a plane
that contains the normal at A and passes through B (the dashed line (---) in Figure 8.1).

8-1
The intersection of this plane with the ellipsoid represents the reduced line of sight from
A to B, and defines an azimuth αA1B1 (on the ellipsoid surface) called the geodetic
azimuth. However, to be able to do computations on the ellipsoid, we need to use the
normal section azimuth α'A1B1 (i.e. the azimuth of the solid line A1B1) (see section
4.2).
The existence of two different normal sections between two points on the ellipsoid (i.e.
from A1 to B1 and from B1 to A1) creates problems when directions are being used in
computations. For example, in Figure 8.2, the observed angles (θ1,θ2,θ3) do not belong
to a closed figure.
C

θ3 NS BC
NS CA
NS CB
NS AC
θ2

normal section BA B
θ1
A
normal section AB

Figure 8.2 - Observed Angles and Normal Sections

In the configuration of three points, one would have six normal sections. Thus the use
of normal sections introduces an added ambiguity in the definition of an ellipsoidal
triangle.
The above problems are overcome by making two corrections to the data:
(i) angular corrections for skew normals (i.e. αAB to α'AB)
(ii) normal section to geodesic correction, called the geodetic line
correction.
All of the corrections that we will deal with in this section arise from the peculiarities of
the geometry of the biaxial ellipsoid.

8.2 Azimuth Correction Due To The Height Of The Observed Point

Directions measured on the surface of the earth are related to the local vertical.
Corrections to these directions to take account of the deflection of the vertical can be
done via the equations given in section 3.6.7 (p. 3-27). Once the directions observed are
corrected to lie along the ellipsoidal normal, we have the problem caused by the
geometry seen in Figure 8.3, i.e. if points PI and Pj are not at the same height, a
correction for the effect of the skew normals needs to be made.

8-2
 Figure 8.3 – Effect of Skew Normals

With reference to Figure 8.3:

αij is the geodetic azimuth of line PiPj (shaded plane)
α'ij is the normal section azimuth – P'iP'j on ellipsoid
hj is the ellipsoidal height of Pj (= H+N )
sij is the ellipsoidal distance between Pi and Pj.

The skew normal correction (α'ij-αij) is given by:

hj
( e ')
2
∆α = α 'ij − α ij = sin 2α ij cos 2 ϕ m
2 ρm
where ρm = mean curvature in meridian between Pi and Pj

2(
=1 ρi + ρ j )
ϕm = mean latitude

2(
=1 ϕi + ϕ j )
e' = second eccentricity (see p. 2-9)

The correction is sometimes referred to as the height of target correction. For all except
high precision calculations, a working formula for ∆α can be used, where

∆α ≅ 0.108′′h j ( ) sin 2α ij cos 2 ϕ m

km

For example, if h j = 3000m, α ij = 45o and ϕ m = 0o , ∆α ~ 0.3′′.

The effect of skew normals is always small and most often insignificant. It should be
noted that the correction is independent of the height of the theodolite station (hi) and of
the distance sighted (sij). The height hj need only be accurate to the nearest 100 m or so

8-3
and errors of a few degrees in αij and ϕm give no significant errors in ∆α except in
critical cases. The general range of skew normal corrections in Australia is of the order
± 0.00 – 0.06″ (for heights < 1000 m).
Figure 8.4 gives an example of skew normal correction magnitudes for a number of
different cases.

Figure 8.4 – An Example of Skew Normal Correction for ϕi = 40°°, ϕj = 41°°

8.3 Forward and Reverse Normal Sections

It is necessary to determine if the difference between the forward and reverse normal
sections is significant for geodetic lines. The linear separation between the forward and
reverse normal sections can be summarised as:

Pi d P'j

Figure 8.5 – Linear Separation Between Normal Sections

(s )
3

d ~ ( e ')
2 ij
2
sin 2α ij ⋅ cos 2 ϕ m
12 R

Therefore, at ϕm = 45°, αij = 45°:

length of line (sij) 200 km 100 km 50 km

displacement of normal sections (d) 55 mm 7 mm 0.9 mm

8-4
The azimuth separation of the forward and reverse normal sections can be illustrated as

North

tangent P'j
α'ji

α'ij
∆α
Pi

Figure 8.6 – Azimuth Separation of Normal Sections

( e ') ( sij )
2 2
cos 2 ϕ m sin 2α 'ij
∆α ~
4ν m

Therefore, at ϕm = 0°, α'ij = 45°:

sij 200 km 100 km 50 km

∆α 0.399″ 0.085″ 0.021″

∆α is maximum at the ends of the line and zero in the middle, i.e. the normal sections
are parallel at the mid-point.
In general, for distances of 20-25 km or less, it is not necessary to consider the angle of
separation ( ∆α ). For longer distances, the difference becomes significant for first order
geodetic networks. The neglect of any correction will therefore cause systematic errors
in the computations. To make these computations simpler, we need to be able to define
a unique curve between Pi and P'j. Such a curve is called a geodesic.

8.4 Geodesic Curve and its Properties

The normal section has a disadvantage as a curve on the surface of the ellipsoid – it is
not unique between two points. The geodesic is a unique curve, which gives the shortest
distance on any surface between two points. If the surface is a plane, the geodesic is a
straight line. On a sphere, the geodesic is a great circle whereas on the ellipsoid, the
geodesic is a curve having double curvature and is thus not a plane curve.
Generally, a surface has, at a point, a curvature that is different in different directions.
A normal plane, containing the normal to the surface at Pi, is the osculating plane of the
normal section made by this plane. It is the curvature of this normal section that defines
the curvature of the surface at Pi in the direction of the normal plane. This relation is,
however, generally valid only at Pij; the curve for which this relation is valid
everywhere is the geodesic. The geodesic has no curvature in the tangent plane, it is
locally straight on the surface; but is the line of maximum curvature between two
points.

8-5
It can be shown (by Clairaut’s theorem) that for any point on a geodesic:
ν123
cos ϕ sin α = constant
p

However, a parallel of latitude satisfies Clairaut’s theorem but it is not a geodesic unless
at ϕ = 0°.

8.5 Geodetic line correction

Figure 8.7 summaries the geometry of the geodetic line correction.

Figure 8.7 – Difference Between Normal Sections and Geodesic Curve

The angle between the geodesic and the normal section at a point is given by:
1
∆α = α 'ij − α ij ≈ ∆α
3
The geodesic approximately trisects the angle between the normal sections at a point.
Since the geodesic is the (locally) shortest length of all curves drawn on the ellipsoid
and its osculating plane contains the ellipsoidal normal at every point along its path, we
can derive the correction term between the geodesic azimuth α ij ( ) and normal

( )
section azimuth α 'ij , required to be added to the forward normal section azimuth as

8-6
 ∆α = α ij − α 'ij

( e ') ( sij ) ( e ') ( sij )

2 2 2 3
sin 2α ij cos 2 ϕ m sin α ij sin 2ϕ m
=− +
12ν m2 48ν m3
2
 s ( km ) 
~ −0.028′′ sin 2α ij cos ϕ m  ij 2
 100 
 

where νm = radius of curvature in prime vertical

= 1 (ν i + ν j )
2
sij = length of geodesic ≡ length of normal section.

For most geodetic surveys, this correction is rarely needed; exceptions are for VLBI
lines. Some approximate magnitudes are given below and also summarised in Figure
8.8.
At ϕm = 45°, α = 45°
s 30 km 60 km 100 km 120 km
0.001″ 0.005″ 0.014″ 0.020″
∆α
↓
~ 12 mm or 1 in 107

Figure 8.8 - Normal Section To Geodesic Correction For ϕm = 45°°.

8.6 Distances on the geodesic

The geodesic distance is a curve contained within the envelope formed by the normal
sections. Its value can be developed as a difference from the normal section distance
using a power series expansion with negligible distortion.

8-7
The difference in length between the normal section (sij) and the geodesic (sg) is

− (e ) s g  s g 
4 4

s g − sij =   sin 2 2α ij cos 4 ϕ m + L

360  ν m 
For s = 100 km, the difference between the geodesic and normal section is of the order
10-9 m and for a distance of 1000 km, the correction term is ~ 4 x 10-5 m. It is ~ 4 mm at
2500 km. Hence this difference can be ignored in almost all circumstances and distances
reduced to the ellipsoid can be considered as either the normal section length or
geodesic length.

8.7 Summary of Reduction from Observed to Ellipsoidal Azimuth

Observed azimuth Azij - Astronomical Azimuth observed at Earth’s surface in
relation to the local vertical, i.e. direction of gravity.
Geodetic azimuth αij - Azimuth at Earth’s surface from Pi → Pj in the normal section
containing both points, i.e. observed azimuth corrected for the deflection of the vertical.
Normal section azimuth α'ij - Azimuth on the ellipsoid containing the normals
through both points, i.e. geodetic azimuth corrected for skew normals.

Geodesic azimuth α ij - Azimuth of geodesic curve (also called ellipsoidal azimuth)

between two points on the ellipsoid, i.e. normal section azimuth corrected for the
geodetic line correction.

The geodetic azimuth αij is obtained from the observed azimuth using the Laplace
equation (see p 3-27):

α ij = Azij − ηi tan ϕi − (ξi sin α ij − ηi cos α ij ) cot ζ ij

1
424
3 14444
4244444
3
constant at deflection of vertical
stn Pi components in azimuth α ij
at zenith angle ζ ij
C1
C2

Typical variations in the correction terms (C1 and C2 in above equation) are shown in
Figure 8.9. These corrections are due to the effects of the earth’s gravity field.

Figure 8.9 – The Effects of the Components of the Laplace Correction ∆α.

Reductions for the above angular quantities are applied to both directions and azimuths
in the same way. However, an observed horizontal angle is reduced to the ellipsoid by
applying the corrections to the two directions of the arms of the angle – i.e. one is left
with the differences of the corrections.

8-8
9 Reduction of Measured Distances to the Ellipsoid

9.1 Remarks on Distance Reductions

Before observed distances can be used in geodetic computations, it is necessary to apply
physical and geometrical reductions. Most terrestrial distance measurements are done
using EDM devices and hence the effect of the propagation medium (the atmosphere) must
be considered (i.e. physical) as well as geometric reductions to ellipsoidal distances.
For EDM distances, these physical corrections may be summarised as:
(i) First Velocity Correction (see Rüeger, Ch. 6). The distance (d') actually displayed
on the distance meter is:
c0 ∆t ' ∆t ' nREF
d'= or = −d'
nREF 2 2 c0
where c0 is the velocity of light in a vacuum
nREF is the reference refractive index of the instrument
∆t' is the measured ‘two way’ travel time of the signal from the EDM
to reflector.
However, the actual length of the wave path (d) is
c0 ∆t '
d=
n 2
where n is the actual refractive index (ideally the integrated mean value
along the wave path).
In order to obtain the (true) wave path length (d) from the read-out distance (d'),
we can derive the equation:

n 
d =  REF d '
 n 
Since we normally deal with differential corrections, we define the first velocity
correction k' as:
n −n
k ' = d − d ' =  REF  d ' ≈ ( nREF − n ) d '
 n 
Thus the corrected distance (d) is given as
d = d '+ k '
First velocity corrections for electro-optical short range EDMs are generally given
in a form like:

 ∆p 11.27e  −6
k ' = c − +  10 d
 ( 273.15 + t ) ( 273.15 + t ) 
≈ f ( t , p, e ) good to order 0.5 ppm
i.e. a function of temperature (t), atmospheric pressure (p) and partial water
vapour pressure (e).

9-1
 (ii) Second Velocity Correction. Usually, a mean first velocity correction is obtained
from meteorological observations at both ends of the line. However another small
correction needs to be made for the effect of differences in curvature between a
radius (R) containing the end points of the line and the actual wave path radius (r)
as this has an effect on the correct refractive index.

For the actual wave path, the correct refractive index would be (assuming linear
gradients of temperature and pressure):

n = 1 (n1 + n2 ) + ∆n
2

where ∆n =
( k − k ) ( d ')
2 2

12 Rα2
The second velocity correction k" is thus given by
k ′′ = − d ' ∆n

( d ')
3

or k ′′ = − ( k − k 2 )
12 Rα2
where k is the coefficient of refraction
1 curvature of ray path
k= r
1 curvature of ellipsoid
R
for light waves kL ~ 0.13 rL = 8R
for microwaves km ~ 0.25 rm = 4R
Rα is the mean radius of curvature of the ellipsoid
along the line.
The second velocity correction k ′′ is more important for microwaves than for light
waves – e.g. k" is of order –50 mm for d' ~ 50 km using a microwave instrument.
Therefore, the wave path length d1 is obtained as:
d1 = d '+ k '+ k "

9-2
After the first and second velocity corrections are applied, the arc distance (d1) must be
reduced geometrically for geodetic computations. If a 2-D adjustment is to be done, the
distance must be reduced to the ellipsoid (common datum). If a 3-D network adjustment is
being carried out, the wave path (d1) must be reduced to the wave path chord (d2) – this
provides the space distance between the two terminals of the line.
There are basically two different methods of reduction – one using known or measured
heights, the other measured zenith angles. The first method is employed in conjunction
with medium to long range distance measurements (~5-100 km) and the second method is
used for short distances (≤5-10 km) where the heights of target stations are usually not
available.

9.2 Reductions using Measured Heights

In what follows, we assume that other geometrical corrections for telescope- and
theodolite-mounted EDMs have been made.

9-3
The ellipsoid is always replaced by a sphere of radius R for the reduction of distances. This
leads to an intersection of the ellipsoidal normals, which were originally skew. Ellipsoidal
heights should always be used for high order geodetic surveys. In low order and short
range surveys, the approximation H ~ h is usually made although the effects of this
approximation are not always negligible. Note using GPS, we can determine h directly and
hence given H, we can determine the geoid-ellipsoid separation N and perhaps produce
more accurate local geoid maps.
The step-by-step 'measured heights' reduction procedure for distances can be summarised
as follows:
{d1 = wave path arc distance}
(i) d1 → d2 (arc-to-chord correction): accounting for curvature of the
electromagnetic wave.

( d1 )
3 ( d1 )
3 5

d 2 = d1 − k 2
2
+k 4
~ 1 mm for d1 ~ 39 km
24 R 640 R 4
 ρν 
where R ≈ Rα =  
 ν cos α + ρ sin α 
2
2

This equation is derived from the relation:

β d1  d1 
d 2 = 2r sin
2
= 2r sin
2r using β ~ r 

 d1  R
and by carrying out series expansion of sin   , remembering r =
 2r  k

Note: d2 is needed for 3-D network (e.g. GPS baseline) adjustment calculations.

{d2 = wave path chord distance}

(ii) d2 → d3 (chord-to-chord correction): chord distance (d2) in space must be

reduced to a chord distance (d3) on the ellipsoid, i.e., reducing all distances to a
common datum. Made up of a slope correction plus a sea-level correction.
Using cosine rule in ∆ P1 P2 C:

d 22 = ( R + h1 ) + ( R + h2 ) − 2 ( R + h1 )( R + h2 ) cos γ
2 2

γ γ d3
and then using cos γ = 1 − 2 sin =
2
and sin we get
2 2 2R
d 22 − (h2 − h1 )
2
d3 =
 h1  h2 
1 + 1 + 
 R  R

The above equation is rigorous and accounts for a slope correction and sea-level
correction, which is an alternative route to reductions using two separate
corrections.

9-4
 Note: The normal expression used for the reduction of horizontal distances (d2') at
hm '
station P1 or P2 is given as − d 2 . This expression is not always sufficient and it
R
is better to use the equation given above or more exact expressions for the sea-level
 
(
hm 2
d 2 − ∆h 2 ) 2  . This is because the slope correction does
1
correction  e.g . −
 R 
not lead to the chord at mean height hm, which would be the true horizontal
distance needed.

{d3 = ellipsoidal chord distance}

(iii) d3 → d4 (chord-to-arc correction): chord distance (d3) must be reduced to the

surface of the ellipsoid to give an arc distance (d4).
γ  d 
d 4 = Rα γ = 2 Rα   = 2 Rα arcsin 3  *exact*
2  2R 
and performing a series expansion of the arcsin term to give an approximation
formula:
d3 3d35
d{4 = d3 + 3 2 + +K
{ 24 R 640 R 4
arc chord 1442443
chord − to − arc corrn.

this formula is correct to an order of ~1 mm at d3 = 200 km

The distance d4 is the distance along the surface of the ellipsoid, denoted
previously as s (i.e. the geodesic or normal section length), and is used in 2-D
calculations of position on the ellipsoid.
An alternate, closed form solution for d4, without any approximations, is given by:

 d k  k2
Rα2 sin 2  1  − (h2 − h1 )
2

s = d 4 = 2 Rα arcsin  2 Rα  4
k (Rα + h1 )(Rα + h2 )
2

derived from
d3
d 4 = 2 Rα arcsin
2 Rα
where
d 22 − ( h2 − h1 )
2

d3 =
 h1   h2 
1 +  1 + 
 Rα   Rα 
and
d1 R
d 2 = 2r + sin and r=
2r k

{d4 = ellipsoidal arc distance}

9-5
9.2.1 Analysis of Errors
- If we assume h ~ H (i.e. geoid-ellipsoid separation term is ignored) then, if N1 = N2,
the resulting distance will be in error by 1 ppm per 6.4 m of N (1 ppm is equal to 0.1
 ∆H 
m in 100 km). If N1 ≠ N2, an additional error will occur of order   ⋅ ∆N .
 d 
- Over long lines (e.g. d1 ~ 20 km, ∆H = 1000 m, hm = 500 m) ∆h must be known to
order 0.1 m, (R + hm) to ~ 1.6 m, R to ~ 1.6 m and k to ~ 2.3 for the d4 calculation to
be better than 5 mm. Hence we must use Rα and heights from geodetic levelling
when desiring results of this precision.

9.3 Reduction using Measured Zenith Angles

9-6
From ∆ Q P1 P2 , we can show that

QP1 = d 2 cos ( z1 + θ1 + δ )
QP2 = d 2 sin ( z1 + θ1 + δ )

In ∆QCP2 ,

 d 2 sin ( z1 + θ1 + δ ) 
γ = arctan  
 R + h1 + d 2 cos ( z1 + θ1 + δ ) 
Since d4 = γ ⋅ R
 d 2 sin ( z1 + θ1 + δ ) 
= R arctan  
 R + h1 + d 2 cos ( z1 + θ1 + δ ) 
d1k
where δ =
2R
For distances ≤ 30 km, this equation for d4 can be simplified (to order 1 mm) as:

  dk 
 d1 sin  z1 + θ1 + 1  
d 4 = R arctan   2R 

 R + h + d cos  z + θ + d1k  
  1 1 
2 R  
1 1

where
z1 = measured zenith angle at P1 (in radians)
θ1 = deflection of the vertical at P1 in azimuth P1P2 (in radians)
k = coefficient of refraction of light = 0.13
R = radius of curvature of ellipsoid along line P1P2
h1 = ellipsoidal height at P1
d1 = wave path length (assumed d2 ≅ d1)

Note: (h2 – h1) from ellipsoidal zenith angles (ζ) (need defl. of vertical) = ellipsoidal ∆h
(H2 – H1) from measured zenith angles (z) = orthometric ∆H

9.3.1 Zenith Angle Measurements/Heights

The observed zenith angles are referred to the local direction of the plumbline and the
tangent to the curved ray at the observation point.
Considering Figure 5.3, z1 is the observed zenith at P1 to P2, however for height
determinations with respect to the ellipsoid we need a zenith angle with reference to the
ellipsoidal normal and the chord distance d2. If the angle of refraction is δ and the
component of the deflection of the vertical in the azimuth of the line P1P2 is θ, then the
required zenith angle ζ1' is given by:
ζ 1' = z1 + δ + θ1

9-7
The difference in ellipsoidal heights between stations P1 and P2 is computed from Jordan’s
formula for heights:

( d4 )
2
 h +h 
∆h12 = h2 − h1 = d 4 1 + 1 2  cot ζ 1' +
 2R  2 R sin 2 ζ 1'
where R is the mean radius of curvature (= ½ [R1 + R2]) and Ri is equivalent to Rα defined
previously on p. 5-4.
This equation is analogous to the one-way trig. heighting formula. δ depends on the
coefficient of refraction k, and hence on the meteorological conditions, particularly the
vertical gradient of temperature, along the line. δ is usually derived from observations at
the end points of the line leading to errors of a few seconds of arc or more for distances >
2-3 km. The applicability of using this formula over longer lines is very dubious. For lines
of up to 10 km and using reciprocal zenith angles (simultaneous), the effect of the
uncertainties in refraction can be reduced.
Using the equation
d4
γ= = ζ 1' + ζ 2' − 180o
R
and the application of the tangent law to the plane triangle CP1P2 gives the equation

 h +h d2   ζ ' −ζ ' 
∆h12 = d 4 1 + 1 2 + 4 2  ⋅ tan  2 1 
 2R 12 R   2 
Since ζ i' = zi + δ i + θ i only the differences in δ and θ at each end of the line appear in this
equation and so if δ 1 = δ 2 & θ1 = θ 2 (likely on short lines) the height difference is only
affected by random errors.
Usually, orthometric height differences related to the geoid are mostly required. This
height difference (∆H = H2 – H1) can be obtained using zenith angles referred to the local
vertical:

∆H12 = H 2 − H1 = d 2 cos z1 +
(1 − k ) (d sin z1 )
2
(+ hI − hT )
2
2R
i.e., the one-way trig levelling formula – this is approximate only and does not take
account of any change in the radius of curvature along the line (~ 12 mm over 10 km) nor
the approximations in deriving the equation (~2 m over 10 km line at 20° elevation).
The exact formula is

 γ 
cos  z1 − (1 − k ) 
∆H12 = d 2  2 
γ
cos
2
but this equation requires a determination of k.
Observations of zenith angles are usually made for two reasons
(i) to determine differences in height
(ii) to reduce slope distances to horizontal etc.
For both purposes, the observed zenith angle should be corrected for curvature of the line
of sight caused by atmospheric refraction.

9-8
If zenith angles are being used to do calculations on a plane, the zenith angles should be
corrected for curvature and refraction. If a true XYZ system is to be achieved, the zenith
angles should be corrected in relation to an adopted Z axis (usually the centroid of the
network) and the determined height differences adjusted for curvature.
For zenith angles, the refraction correction (with correct sign) is given by:
d1k
δ = ≈ +2.1′′ / km
2R
and the correction for curvature is:
γ −d3
= ≈ −16.2′′ / km.
2 2R
Thus the combined curvature and refraction correction to zenith angles is –14.1″/km. We
would apply this correction to zenith angles before making slope distance reductions /
height difference calculations using plane trigonometry calculations.

9-9
 10 Geodetic Techniques and Applications
This chapter will look at some of the major geodetic techniques and give an overview of
their use in geodetic applications.
GPS
GPS (The Global Positioning System) was initially designed as a military satellite
positioning system and is now the most commonly used (by lay-people, surveyors and
geodesists) of all the positioning systems (overtaking the compass!)
In this system, a series of ~30 satellites are orbiting the earth in known orbits. GPS
receivers are designed to pick up the signals from these satellites and therefore
determine a travel time for the signals and consequently calculate a distance between
the particular satellite and the point on the earth. To solve for the 3D coordinate of the
point on the earth the receiver needs to simultaneously acquire signals from at least 3
satellites (usually up to 12 satellites are available) but in reality we need four satellites
to allow us to also solve for a time offset in the signals.
For geodetic applications, GPS has opened up many doors especially in the area of
geophysics. Some examples of the information which can now be determined from GPS
observations (besides coordinates) include:
- Tropospheric water vapour monitoring
(see http://cepheus.cosmic.ucar.edu/indexGlobal.html)

10-1
- Ionospheric Mapping
(see http://iono.jpl.nasa.gov/latest_rti_global.html, http://iono.jpl.nasa.gov/gaim/rtgaim.html)

10-2
- Sea-level monitoring: via direct measurements of sea level (e.g. GPS buoys, tide gauge
installations) and indirect measurements (e.g. crustal uplift, ice mass balance).

INTERACTION OF SEA LEVEL, TIDE GAUGES AND GPS

• In order to estimate ‘absolute’ sea level and its variations, we
need to measure motion of the land in addition to relative sea level
changes wrt a local BM.

• We can achieve accuracies within the 10-20mm level in a global

reference frame. We require observations over a number of years
(~10) to resolve motion to the ±1-1.5 mm/yr level.

Tide Gauge
GPS Site

Tide Gauge
VLBI
Fiducial GPS
Site

Geoid Equipotential Surface

(Ideally considered MSL)

Mathematical
Ellipsoid reference surface

- Determination of earth rotation parameters: via solutions of a global network of

sites (e.g. IGS final solution)

Sample
Accuracy Latency Updates Archive locations
Interval

PM 0.05 mas
CDDIS(US-MD)
daily
PM <0.2 ~13 SOPAC(US-CA)
Final weekly (1200
rate mas/day days IGN(FR)
UTC)
IGS CB(US-CA)
LOD 0.02 ms

10-3
- Plate motion/Crustal deformation: via permanent networks

- Glacial motion and ice mass balance

10-4
- Rift mechanics and iceberg calving

- Earthquake & Volcanic monitoring: via dense, permanent networks in

earthquake/volcanic zones

10-5
VLBI
(the following is directly from http://www.ga.gov.au/geodesy/sgc/vlbi/vlbitech.jsp)
What is VLBI?
VLBI (Very Long Baseline Interferometry) is a geometric technique to measure the
time difference between the arrival, at two antennas, of a radio wave front emitted from
a distant quasar. The time difference measurements are so precise that they are
translated to determine the relative positions of the antennas to within a few millimetres
even if the antennas are separated by several thousands of kilometres.
The History of VLBI
Idea of Very Long Baseline Interferometry was introduced by Russian radio
astronomers in 1965 and realised technically in the USA in 1966. Since 1983 regular
geodetic and astrometrical VLBI sessions have been promoted internationally under
supervision of the Goddard Space Flight Centre (GSFC, Maryland, USA). Australia
participated in first intercontinental VLBI experiments in the late 1960's jointly with
American, Canadian and Russian scientists. Currently, two Australian radio telescopes –
Tidbinbilla and Hobart – participate in geodetic VLBI experiments.
Networks
About 50 VLBI antennas around the world are combined in networks for different
special programs. But there are 3 permanent networks that operate for geodetic research
in local regions – the European network, the Asia-Pacific network and the North
American network. Australia participates in activities of the Asia-Pacific network
having a 26-metre VLBI dish in Hobart, Tasmania.

10-6
 Hobart VLBI Antenna

The Role of Geodetic VLBI

- Geodetic VLBI observations provided the first direct confirmation of tectonic plate
motion at the end of the 1980's. Now VLBI observations measure the motions with
accuracy less than 1 mm/year jointly with other space geodetic techniques such as
GPS and SLR , etc.
- VLBI observations from stable distant quasars are an important component for the
establishment of a reference system of coordinates. The radio system has replaced an
optical reference system bearing on star positions. The optical system has been used
for the last 200 years and had an average accuracy of star positions to about 10
milliarcseconds (equivalent to 0.01"). The current average accuracy of quasar
positions observed by radio systems is about 0.1-0.2 milliarcseconds (50-100 times
better); therefore International Astronomical Union has recommended using radio
systems for broad VLBI application.
For information on the Australian Long Baseline Array (LBA), check out
http://www.atnf.csiro.au/vlbi/.

10-7
Satellite Altimetry:
Satellite altimetry is used to measure the height of the ocean surface and to monitor
variations in sea level, which is directly connected to studying the effects of global
warming. Ocean currents and gravity anomalies over the ocean can also be determined.
See CW’s lecture for details.

It is necessary to calibrate altimeters for absolute bias and bias drift before they can be
used for geodetic applications. This can be done, for example, using GPS buoys.

10-8
The following is a map of the sea level observed by TOPEX/Poseidon:

Major Satellite Missions

The major satellite missions used for geodetic applications are summarised here.

Gravity Field Modelling:

• CHAMP (Challenging Minisatellite Payload) – launched July 15, 2000. See
http://www.gfz-potsdam.de/pb1/op/champ/.

• GRACE (Gravity Recovery And Climate Experiment) – launched March 17,

2002. See http://www.csr.utexas.edu/grace/.

10-9
 • GOCE (Gravity Field and Steady-State Ocean Circulation Explorer) – to be
launched 2007. See http://www.esa.int/export/esaLP/goce.html.

Glaciology:
• ICESat (Ice, Cloud and Land Elevation Satellite) – launched Jan 13, 2003; to
measure changes in polar ice sheets (ice sheet mass balance, sea-level studies).
See http://icesat.gsfc.nasa.gov/.

Ocean Surface Topography Modelling:

TOPEX/POSEIDON and JASON-1: See http://topex-www.jpl.nasa.gov/.

Remote Sensing:
European Remote Sensing Satellites ERS-1/ERS-2: See http://earth.esa.int/ers/.

Orbit Determination:
DORIS is a Doppler satellite tracking system developed for precise orbit determination
and precise ground location. It is onboard the TOPEX/POSEIDON, Jason-1 and
ENVISAT altimetric satellites and the remote sensing satellites SPOT-2, SPOT-4 and
SPOT-5. It also flew with SPOT-3. For more info see http://ids.cls.fr/welcome.html and
http://www.jason.oceanobs.com/html/doris/welcome_uk.html.

10-10
Satellite Laser Ranging (SLR)

The following is taken from the NASA publication available at

http://ilrs.gsfc.nasa.gov/docs/slrover.pdf (also in the KGG350 directory on the
Geography server).

A global network of SLR stations measure the instantaneous round-trip flight time of
ultrashort pulses of light to satellites equipped with special reflectors. SLR was initiated
by NASA in 1964 and now has a ranging precision of a few mm. This is the most
accurate technique to determine the geocentric position of earth satellites (orbits). It is
used for the precise calibration of radar altimeters, which is essential for the
measurement of global MSL changes of a few mm/yr. SLR has the ability to measure
temporal variations in the Earth’s gravity field and monitor station motions with
reference to the geocentre in an absolute system. It is also used for verification of the
Theory of General Relativity.

SLR contributes to the modelling and evaluation of long-term climate change by:
 Providing a reference system for post-glacial rebound, sea level and ice volume
change.
 Determining temporal mass redistribution of the solid Earth, ocean and
atmosphere system.
 Monitoring response of atmosphere to seasonal variations in solar heating.

SLR provides precise orbit determination for:

 Spaceborne radar altimeter missions mapping the ocean surface (used to model
global ocean circulation).
 Mapping volumetric changes in continental ice masses.
 Land topography.

Check out the International Laser Ranging Service website: http://ilrs.gsfc.nasa.gov.

Lunar Laser Ranging (LLR)

LLR measures the round-trip travel times of light pulses between stations on Earth and
four retro-reflectors on the surface of the Moon. It contributes to lunar sciences, the
theory of gravitation, and is a key IERS technique for connecting reference frames. LLR
is used to determine the obliquity of the ecliptic, the orientation of the dynamical frame
of the Solar System in the extragalactic reference frame, and the long-period nutation
and precession. Its accuracy is a few mm in range.

10-11
Interferometric Synthetic Aperture Radar (InSAR)
Some of the following is based on information available on the web at
http://www.gmat.unsw.edu.au/snap/work/diff_insar.htm.
A synthetic aperture radar (SAR) is a coherent radar system that generates high
resolution remote sensing imagery. Signal processing uses magnitude and phase of the
received signals over successive pulses from elements of a 'synthetic aperture' to create
an image.

Interferometric SAR (InSAR) systems exploit the phase difference from two SAR
images acquired over the identical scene, thereby useful geodetic information such as
the digital elevation model (DEM) can be derived. Two sets of SAR image values are
registered by using correlation techniques in the InSAR processing software to an
accuracy of 1/8 of a pixel. The phase differences between two SAR images acquired
over the same region are then calculated to generate the interferogram. These phase
differences in the interferogram wrap around in cycles of 2pi, and need to be
'unwrapped' in order to obtain absolute phase values. Finally, the absolute phase can be
converted into height information.

If used in differential mode, InSAR can be facilitated for the measurement of small
deformations of the terrain (down to the millimetre level) that have occurred between
the two different image acquisitions. The topographic contribution to the phase
difference is carefully removed using a DEM and atmospheric disturbances can be
either accounted for using independent observations such as GPS, or neglected if the
tropospheric delay is considered to be homogeneous at the time of the radar image
acquisition. InSAR can be used for a wide range of applications such as the monitoring
of:
 Subsidence due to underground mining or oil/gas fields.
 Urban subsidence due to tunnel construction and/or groundwater extraction.
 Ice flow in polar regions.
 Seismic deformation and volcanic activity.
 Landslide or mass movement.
 Beach or coast erosion.

The figure below shows an InSAR interferogram of the “Loose Tooth” at the front of
the Amery Ice Shelf, Antarctica acquired from an ERS-1/2 tandem pair in 1996. The
developing east-west rift can easily be recognised.

10-12
 ↑
Transverse rift
has just started
to develop

A more recent satellite image shows the current situation at the “Loose Tooth” – it will
break off soon!

10-13