Mathematical Geodesy

Maa-6.3230

Martin Vermeer

16th February 2015

Figure: Cluster of galaxies Abell 2218, distance 2 billion light years, acts as a gravi-
tational lense, The geometry of space-time within the cluster is non-Euclidean.
Image credit: NASA/ESA
Course Description

Workload 7 cr
Teaching Period III
Learning Outcomes After completing the course, the student

◦ Is able to do simple computations on the sphere and understands the geometry of the
reference ellipsoid
◦ Is able to solve, using tools like MatlabTM , the geodetic forward and inverse problems
etc. on the reference ellipsoid
◦ Masters the basics of global and local reference systems and is able to execute trans-
formations
◦ Masters the basics of Gaussian and Riemannian surface theories and can derive the
metric tensor, Christoffel symbols and curvature tensor for simple surfaces
◦ Masters the basic math of map projections, esp. conformal ones, and the behaviour
of map scale, and is able to compute the isometric latitude.

Content Spherical trigonometry, geodetic co-ordinate computations in ellipsoidal and rectangu-

lar spatial co-ordinate systems, astronomical co-ordinates, co-ordinate system transforma-
tions, satellite orbits and computations, Gaussian and Riemannian surface theories, map
projection computations.
Foreknowledge
Equivalences Replaces course Maa-6.230.
Target Group
Completion Completion in full consists of the exam, the exercise works and the calculation
exercises.
Workload by Component

◦ Lectures 16 × 2 h = 32 h
◦ Independent study 55 h
◦ Exercise works 2 × 12 h = 24 h (independent work)
◦ Calculation exercises at home 16, of which must be done 12 ×6 h = 72 h (independent
work)
◦ Exam 3 h
◦ Total 186 h

Grading The grade of the exam becomes the grade of the course, 1-5
Study Materials Lecture notes. Background material Hirvonen: Matemaattinen geodesia; Torge:
Geodesy
Teaching Language Suomi
Course Staff and Contact Info Martin Vermeer, room Technopolis Innopoli 3D (Vaisalantie 8)
4.143, martin dot vermeer at aalto dot fi

i
Reception times: email martin dot vermeer at aalto dot fi, or phone 050 3574139
CEFR-taso
Lisätietoja

ii
Contents

1. Spherical trigonometry 1
1.1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2. Spherical excess . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.3. The surface area of a triangle on a sphere . . . . . . . . . . . . . . . . . . . . . . 2
1.4. A rectangular spherical triangle . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.5. A general spherical triangle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.6. Deriving the formulas with the aid of vectors in space . . . . . . . . . . . . . . . 5
1.7. Polarization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.8. Solving the spherical triangle by the method of additaments . . . . . . . . . . . . 7
1.9. Solving the spherical triangle by Legendre’s method . . . . . . . . . . . . . . . 8
1.10. The forward geodetic problem on the sphere . . . . . . . . . . . . . . . . . . . . . 8
1.11. The geodetic inverse problem on the sphere . . . . . . . . . . . . . . . . . . . . . 9
1.12. The half-angle cosine rule . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2. The geometry of the reference ellipsoid 11

2.1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.2. The geodesic as solution to a system of differential equations . . . . . . . . . . . 11
2.3. An invariant . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.4. The geodetic main problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.5. The geodetic inverse problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

3. Co-ordinates on the reference ellipsoid 15

3.1. Representations of the sphere and the ellipsoid . . . . . . . . . . . . . . . . . . . 15
3.2. Various latitude types . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.3. Measures for the flattening . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
3.4. Relationships between different types of latitude . . . . . . . . . . . . . . . . . . 16
3.5. Co-ordinates in the meridional ellipse . . . . . . . . . . . . . . . . . . . . . . . . . 17
3.6. Three-dimensional rectangular co-ordinates on the reference ellipsoid . . . . . . . 17
3.7. Computing geographic co-ordinates from rectangular ones . . . . . . . . . . . . . 18
3.8. Meridian arc length . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

4. Reference systems 21
4.1. The GRS80 system and geometric parameters . . . . . . . . . . . . . . . . . . . . 21
4.2. Gravimetric parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
4.3. Reference frames . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
4.4. The orientation of the Earth . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
4.5. Transformations between systems . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

5. Using rotation matrices 29

5.1. General . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
5.2. Chaining matrices in three dimensions . . . . . . . . . . . . . . . . . . . . . . . . 30
5.3. Orthogonal matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
5.4. Topocentric systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
5.5. From geocentric to topocentric and back . . . . . . . . . . . . . . . . . . . . . . . 33
5.6. The geodetic main and inverse problems with rotation matrices . . . . . . . . . . 34
5.6.1. Geodetic main problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

iii
Contents

5.6.2. Geodetic inverse problem . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

5.6.3. Comparison with ellipsoidal surface solution . . . . . . . . . . . . . . . . . 35

6. Co-ordinate systems and transformations 37

6.1. Geocentric terrestrial systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
6.2. Conventional Terrestrial System . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
6.3. Polar motion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
6.4. The Instantaneous Terrestial System . . . . . . . . . . . . . . . . . . . . . . . . . 38
6.5. The quasi-inertial geocentric system . . . . . . . . . . . . . . . . . . . . . . . . . 38

7. Reference systems and realizations 41

7.1. Old and new reference systems; ED50 vs. WGS84/GRS80 . . . . . . . . . . . . . 41
7.2. WGS84 and ITRS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
7.3. Co-ordinate system realizations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
7.4. Realization of WGS84 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
7.5. Realizations of ITRS/ETRS systems . . . . . . . . . . . . . . . . . . . . . . . . . 42
7.6. The three-dimensional Helmert transformation . . . . . . . . . . . . . . . . . . . 43
7.7. Transformations between ITRF realizations . . . . . . . . . . . . . . . . . . . . . 44

8. Celestial co-ordinate systems 47

8.1. Sidereal time . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
8.2. Trigonometry on the celestial sphere . . . . . . . . . . . . . . . . . . . . . . . . . 49
8.3. Using rotation matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
8.4. On satellite orbits . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
8.5. Crossing a given latitude in the inertial system . . . . . . . . . . . . . . . . . . . 52
8.6. The satellite’s topocentric co-ordinates . . . . . . . . . . . . . . . . . . . . . . . . 53
8.7. Crossing a given latitude in the terrestrial system . . . . . . . . . . . . . . . . . . 54
8.8. Determining the orbit from observations . . . . . . . . . . . . . . . . . . . . . . . 54

9. The surface theory of Gauss 59

9.1. A curve in space . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
9.2. The first fundamental form (metric) . . . . . . . . . . . . . . . . . . . . . . . . . 60
9.3. The second fundamental form . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
9.4. Principal curvatures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
9.5. A curve embedded in a surface . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
9.6. The geodesic . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

10.The surface theory of Riemann 73

10.1. What is a tensor? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
10.2. The metric tensor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
10.3. The inverse metric tensor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
10.3.1. Raising or lowering sub- or superscripts of a tensor . . . . . . . . . . . . . 79
10.3.2. The eigenvalues and -vectors of a tensor . . . . . . . . . . . . . . . . . . . 80
10.3.3. The graphic representation of a tensor . . . . . . . . . . . . . . . . . . . . 80
10.4. The Christoffel symbols . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
10.5. The geodesic revisited . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
10.6. The curvature tensor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
10.7. Gauss curvature and spherical excess . . . . . . . . . . . . . . . . . . . . . . . . 86
10.8. The curvature in quasi-Euclidean geometry . . . . . . . . . . . . . . . . . . . . . 88

11.Map projections 91
11.1. Map projections and scale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
11.1.1. On the Earth surface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
11.1.2. In the map plane . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
11.1.3. The scale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92

iv
 Contents

11.1.4. The Tissot-indicatrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94

11.2. The Lambert projection (LCC) . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
11.3. On the isometric latitude . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
11.4. The Mercator projection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
11.5. The stereographic projection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
11.6. The Gauss-Krüger projection . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
11.7. Curvature of the Earth surface and scale . . . . . . . . . . . . . . . . . . . . . . . 104

12.Map projections in Finland 107

12.1. Traditional map projections . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
12.2. Modern map projections . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
12.3. The triangulated affine transformation used in Finland . . . . . . . . . . . . . . . 108
12.3.1. Plane co-ordinates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108

A. Isometric latitude on the ellipsoid 111

B. Useful equations connecting the main radii of curvature 113

C. The Christoffel symbols from the metric 115

D. The Riemann tensor from the Christoffel symbols 117

v
 Chapter 1
Spherical trigonometry

1.1. Introduction

The formulas of spherical geometry are very useful in geodesy. The surface of the Earth, which
to first approximation is a plane, is in second approximation (i.e., in a small, but not so small,
area) a spherical surface. Even in case of the whole Earth, the deviation from spherical shape
is only 0.3%.
The starry sky again may be treated as a precise spherical surface, the radius of which is
undefined; in practical computations we often set R = 1.

1.2. Spherical excess

See figure 1.1. Let us assume that the radius of the sphere is 1. The “front half” of the sphere is
a semi-sphere, the surface area of which is 2π. A triangle is formed between three great circles.
The same great circle form, on the back surface of the sphere, a “antipode triangle” of the same
size and shape.
When the surface area of the whole semi-sphere is 2π, the area of the “orange slice” bounded by
α
two great circles will be · 2π, where α is the angle between the great circles. We obtain
π
A1 + A2 = 2α
A1 + A3 = 2β
A1 + A4 = 2γ

and
A1 + A2 + A3 + A4 = 2π.
By summing up the first three equations we obtain

2A1 + A1 + A2 + A3 + A4 = 2 (α + β + γ)

β α
A1
γ
A2

A3
A4

Figure 1.1.: Spherical triangles on a semi-sphere. The back side surface of the sphere has been
depicted in a lighter shade with its “antipode triangles”

1
Chapter 1. Spherical trigonometry

i.e.,
A1 = α + β + γ − π = ε,
where ε is called the spherical excess.
If the radius of the sphere is not 1 but R, we obtain
A1
A1 = εR2 ⇒ ε = .
R2
Here ε is expressed in radians. If ε is not in radians, we may write
ρunit A1
ε [unit] = ,
R2
where ρunit is the conversion factor of the unit considered, e.g., for degrees, 57.29577951308232087721
or for gons, 63.66197723675813430801.
As we see is the spherical excess inversely proportional to R2 , i.e., directly proportional to the
total curvature R−2 . It is also directly proportional to the surface area of the triangle.
This is a special case of a more general rule:
The directional closing error of a vector which is transported in a parallel way around
the closed edge of a surface is the same as the integral over the surface of the total
curvature.
As a formula: ˆ
ε= Kdσ,
A
where σ is the variable for surface integration, and K is the total curvature of the surface
according to K.F. Gauss, which thus can vary from place to place. E.g., on the surface of an
ellipsoid
1
K= ,
MN
where M is the meridional curvature (in the North-South direction) and N the so-called trans-
verse curvature in the East-West direction. Both depend on the latitude ϕ. In a smallish area,
the internal geometry of the√ellipsoidal surface does not differ noticeably from a spherical surface,
the radius of which is R = M N .
If the triangle of the surface of the sphere is small compared to the radius of the Earth, also
the spherical excess will be small. In the limit we have ε → 0 ja α + β + γ = π exactly. We say
that the a plane surface (or a very small part of a spherical surface) forms a Euclidean space,
whereas a spherical surface is non-euclidean.

1.3. The surface area of a triangle on a sphere

If the triangle isn’t very large – i.e., just as large as geodetic triangulation triangles generally
on, at most some 50 km –, we may calculate its surface area using the formula for the plane
triangle:
1 1
A = a · ha = ab sin γ,
2 2
where ha is the height of the triangle relative to the a side, i.e., the straight distance of corner
point A from side a.
sin β
Because according to the sine rule b = a , it also follows that
sin α
a2 sin β sin γ
A= .
2 sin α

2
 1.4. A rectangular spherical triangle

A
A D
α 1
sin b
cos b
sin a cos b
1 b a
O cos b
D O cos a cos b
E
b
b c
O c A
. D γ. A
a β 1 sin c
E . a
C sin c sin c sin β
β c E β
O E D
B cos c sin c cos β
Figure 1.2.: A rectangular spherical triangle. The partial triangles are lifted out for visibility

At least for computing the spherical excess, these approximate formulas are good enough:

A
ε= ,
R2

where also the approximate value for R, e.g., R ≈ a = 6378137 m, is completely sufficient.

1.4. A rectangular spherical triangle

This case is depicted in figure 1.2. Many simple formulas follow directly from the figure and the
separately drawn plane triangles:

EO = cos c = cos a cos b

DE = sin c cos β = sin a cos b (1.1)
AD = sin c sin β = sin b

By interchanging the roles of a and b (and thus of α and β) we obtain furthermore

sin c cos α = sin b cos a

sin c sin α = sin a (1.2)

of which the first yields

sin b
cos α = cos a = cos a sin β,
sin c
according to the last equation in group (1.1).

By dividing the first of group (1.2) by its second, we obtain

cot α = cot a sin b

and the second of group (1.1) by its third, correspondingly

cot β = cot b sin a.

3
Chapter 1. Spherical trigonometry

a
c

γ
α b1 b2
.
Figure 1.3.: General spherical triangle

1.5. A general spherical triangle

The formulae for a spherical triangle are obtained by dividing the triangle into two right triangles,
see figure 1.3. here, the third side is b = b1 + b2 .
If we apply to the sub-triangles the formulas derived above, we obtain:

cos a = cos h cos b2 ,

sin a cos γ = cos h sin b2 ,
sin a sin γ = sin h,

cos c = cos h cos b1 ,

sin c cos α = cos h sin b1 ,
sin c sin α = sin h.

By substituting

sin b1 = sin (b − b2 ) = sin b cos b2 − cos b sin b2 ,

cos b1 = cos (b − b2 ) = cos b cos b2 + sin b sin b2

we obtain

cos c = cos h (cos b cos b2 + sin b sin b2 ) =

= cos b (cos h cos b2 ) + sin b (cos h sin b2 ) =
= cos b cos a + sin b sin a cos γ, (1.3)

the so-called cosine rule of spherical trigonometry, and

sin c cos α = cos h (sin b cos b2 − cos b sin b2 ) =

= sin b (cos h cos b2 ) − cos b (cos h sin b2 ) =
= sin b cos a − cos b sin a cos γ.

From the two “sin h” formulas we obtain

sin c sin α = sin a sin γ,

or, more generally,

sin a sin b sin c
= = , (1.4)
sin α sin β sin γ
the so-called sine rule of spherical trigonometry.

4
 1.6. Deriving the formulas with the aid of vectors in space

C
γ
b
xC
a

A α
c
xA β
ϕA b B
a
c xB
ϕB
O y
x λB
Figure 1.4.: The spherical triangle ABC, for deriving the cosine and sine rules using vectors in
space

For comparison the corresponding formulas for a plane triangle:

c2 = a2 + b2 − 2ab cos γ

and
a b c
= = .
sin α sin β sin γ
At least in the case of the sine rule it is clear, that in the limit for a small triangle sin a −→ a
etc., in other words, the spherical sine rule morphs into the one for a plane triangle. For the
cosine rule this is not immediately clear.

1.6. Deriving the formulas with the aid of vectors in space

 π 
If on the sphere we look at a triangle consisting of two points A ϕA = − b, λA = 0 ja
2
 π   π 
B = ϕB = − a, λB = γ and a pole C = ϕC = , λC = arbitrary , see figure 1.4, we may
2 2
write two vectors:

       
xA cos ϕA xB cos ϕB cos λB
xA =  yA  =  0  , xB =  yB  =  cos ϕB sin λB  .
zB sin ϕA zB sin ϕB

The vectors’ dot product is

cos c = xA · xB = cos ϕA cos ϕB cos λB + sin ϕA sin ϕB =

= sin b sin a cos γ + cos b cos a,

the cosine rule for a spherical triangle.

The cross product of the vectors is
 
− sin ϕA cos ϕB sin λB
xA × x B =  sin ϕA cos ϕB cos λB − cos ϕA sin ϕB  =
cos ϕA cos ϕB sin λB
 
− cos a sin b sin γ
=  cos a sin b cos γ − sin a cos b  .
sin b sin a sin γ

5
Chapter 1. Spherical trigonometry

π−c
π
π−β 2 π−α
π
γ 2
π
2 b
a
π
2
α
c β
π−a
π π−b
π
2
2

π−γ

Figure 1.5.: Polarization of a spherical triangle

When the third vector is 


0
xC =  0  ,
1
we obtain the volume of the parallelepiped spanned by the three vectors (i.e., twice the volume
of the tetrahedron ABCO) as follows:

Vol {xA , xB , xC } = (xA × xB ) · xC = sin b sin a sin γ.

The volume contained by the three vectors does not however depend on the order of the vectors,
so also
Vol {xA , xB , xC } = sin b sin c sin α = sin a sin c sin β.
Division yields

sin a sin γ = sin c sin α,

sin b sin γ = sin c sin β,

i.e.,
sin a sin b sin c
= = ,
sin α sin β sin γ
the sine rule for a spherical triangle.

1.7. Polarization

For every corner of the triangle we may define an “equator” or great circle, one “pole” of which
is that corner point. If we do this, we obtain three “equators”, which themselves also form a
triangle. This procedure is called polarization.
Because the angular distance between the two corner points on the surface of the sphere is the
length of the side, the angle between the planes of two such great circles must be the same as
this length. And the “polarization triangle”’s angle is 180◦ minus this.

6
 1.8. Solving the spherical triangle by the method of additaments

The polarization method is symmetric: the original triangle is also the polarization of the
polarization triangle. The intersection point of two edges of the polarization triangle is at a
distance of 90◦ from both “poles”, i.e., the corners of the original triangle, and the edge between
them thus is the “equator” of the intersection point.
For symmetry reasons also the length of an edge of the polarization triangle equals 180◦ minus
the corresponding angle of the original triangle.
For an arbitrary angle α we have:

sin (180◦ − α) = sin α

cos (180◦ − α) = − cos α

Because of this we obtain of the spherical trigonometry cosine rule (1.3) the following polarized
version:
− cos γ = (− cos β) (− cos α) + sin β sin α (− cos c)
or
cos γ = − cos β cos α + sin β sin α cos c,
a formula with which one may calculate an angle if the two other angles and the side between
them are given.

1.8. Solving the spherical triangle by the method of additaments

In the additaments method we reduce a spherical triangle to a plane triangle by changing the
lengths of the sides. Generally all angles and one side of a triangle are given, and the problem
is to compute the other sides.
As the sides are small in comparison with the radius R of the Earth, we may write (series
expansion):
s2
 
1 3 s
sin ψ = ψ − ψ + . . . ≈ 1− .
6 R 6R2
Now the spherical sine rule is (s = a, b, c):

a (1 − ∂a) b (1 − ∂b) c (1 − ∂c)

= = ,
sin α sin β sin γ

s2 a2 b2 c2
jossa ∂s = : ∂a = , ∂b = ja ∂c = .
6R2 6R2 6R2 6R2
The method works now so, that
1. From the known side we subtract its additament ∂s;
2. The other sides are computed using the sine rule for a plane triangle and the known angle
values ;
3. To the computed sides are now added their additaments.
The additaments are computed using the best available approximate values; if they are initially
poor, we iterate.
The additament correction
s0 = s (1 − ∂s)
may be changed from a combination of a multiplication and a subtraction into a simple subtrac-
tion by taking logarithms:

ln s0 = ln s + ln (1 − ∂s) = ln s + (0 − ∂s) = ln s − ∂s,

7
Chapter 1. Spherical trigonometry

∆λ

ψ 12
A12
ϕ2
1
ψ12

ϕ1
.λ
. ekvaattori 2

λ1
Figure 1.6.: The triangle 1 − 2 − N on the globe. N is the North pole

or, in base-ten logarithms

10
log s0 =10 log s − M · ∂s,

where M =10 log e = 0.43429448. In the age of logarithm tables this made the practical
computations significantly easier.

1.9. Solving the spherical triangle by Legendre’s method

In the Legendre method the reduction from a spherical to a plane triangle is done by changing
the angles. If again all angles and one edge is given, we apply the following formula:

a b c
= = ,
sin (α − ε/3) sin (β − ε/3) sin (γ − ε/3)

i.e., from every angle we subtract one third of the spherical excess ε.
It is however important to understand that the further calculations must be made using the
original angles α, β, γ! The removal of the spherical excess is only done for the computation of
the unknown sides of the triangle.
Nowadays these approximate methods (additaments and Legendre) are no longer used. It is
easy to compute directly by computer using the spherical sine rule.

1.10. The forward geodetic problem on the sphere

The spherical trigonometry cosine and sine rules, suitably applied:

sin ϕ2 = sin ϕ1 cos ψ12 + cos ϕ1 sin ψ12 cos A12 ,

8
 1.11. The geodetic inverse problem on the sphere

and
sin (λ2 − λ1 ) sin A12
= ⇒
sin ψ12 cos ϕ2
 
sin ψ12 sin A12
λ2 = λ1 + arcsin .
cos ϕ2

1.11. The geodetic inverse problem on the sphere

The cosine and sine rules of spherical trigonometry, suitably applied:

cos ψ12 = sin ϕ1 sin ϕ2 + cos ϕ1 cos ϕ2 cos (λ2 − λ1 ) , (1.5)

sin (λ2 − λ1 )
sin A12 = cos ϕ2 . (1.6)
sin ψ12

1.12. The half-angle cosine rule

The above spherical cosine rule (1.3):

cos a = cos b cos c + sin b sin c cos α,

is numerically ill-behaved when the triangle is very small compared to the sphere, in other
words, if a, b, c are small. E.g., the triangle Helsinki-Tampere-Turku is very small compared to
the globe, some 200 km/6378 km ∼ 0.03. Then cos b cos c ∼ 0.999, but sin b sin c ∼ 0.0009! We
are adding two terms of which one is approx. 1 and the other about a thousand times smaller.
That is the way to lose numerical precision.
For solving this, we first remark, that
α
cos α = 1 − 2 sin2 ;
2
a
cos a = 1 − 2 sin2 ;
2
and

cos b cos c + sin b sin c cos α = (cos b cos c + sin b sin c) + sin b sin c (cos α − 1) =
α
= (cos b cos c + sin b sin c) − 2 sin b sin c sin2 ;
2
b − c
cos b cos c + sin b sin c = cos (b − c) = 1 − 2 sin2 ;
2
after a few rearrangements we obtain the half-angle cosine rule for a spherical
triangle, which also for small triangles is wel behaved1 :

a b−c α
sin2 = sin2 + sin b sin c sin2 .
2 2 2

1
. . . which, surprise, contains only sines!

9
 Chapter 2
The geometry of the reference ellipsoid

2.1. Introduction

The reference ellipsoid is already a pretty precise description of the true figure of the Earth. The
deviations of mean sea level from the GRS80 reference ellipsoid are of magnitude ±100 m.
Regrettably, the geometry of the reference ellipsoid is not as simple as that of the sphere.
Nevertheless, rotational symmetry brings in at least one beautiful invariant.
The geodetic forward and reverse problems have traditionally been solved by series expansions
of many terms, the coefficients of which also contain many terms. Here we rather offer numerical
methods, which are conceptually simpler and easier to implement in an error-free way.

2.2. The geodesic as solution to a system of differential equations

We may write in a small rectangular triangle (dy, dx the metric east and north shift, see Fig.
2.1):

dx = M (ϕ) dϕ = cos Ads,

dy = p (ϕ) dλ = sin Ads ja
dA = sin ϕdλ,

where p = N cos ϕ is the distance from the axis of rotation, and M and N are the meridional
and transversal curvatures, respectively.
The system of equations, normalized:

dϕ cos A
= ,
ds M (ϕ)
dλ sin A
= , (2.1)
ds p (ϕ)
dA dλ sin ϕ sin A
= sin ϕ = .
ds ds p (ϕ)

Group 2.1 is valid not only on the ellipsoid of revolution; it applies to all figures of revolution.
On a rotationally symmetric body we have p (ϕ) = N (ϕ) cos ϕ, i.e.,

dλ sin A
= ,
ds N (ϕ) cos ϕ
dA tan ϕ sin A
= .
ds N (ϕ)

If the initial condition is given as ϕ1 , λ1 , A12 , we may obtain the geodesic ϕ (s) , λ (s) , A (s)
as a solution parametrized by arc length s. Numerically computing the solution using the

11
Chapter 2. The geometry of the reference ellipsoid

dA

N
p − dp π dx

co
−ϕ

tϕ
2
p

. p − dp
dx
p= dλ
Nc
os ϕ A ds
. A dx
dy
ds
.
M

dϕ dy

Figure 2.1.: The geometry of integrating the geodesic

MatLab software is also fairly easy thanks to the ODE routines (“Ordinary Differential Equation”)
provided.
Transformation to rectangular form is easy:

X (s) = N (ϕ (s)) cos ϕ (s) cos λ (s) ,

Y (s) = N (ϕ (s)) cos ϕ (s) sin λ (s) ,
b2
Z (s) = N (ϕ) sin ϕ (s) .
a2

2.3. An invariant

Let us look closer at the quantity p (ϕ) = N (ϕ) cos ϕ, the distance of a point from the rotation
axis. let us compute the s derivative:

dp dp dx
= = − sin ϕ cos A.
ds dx ds

Now, with the equation 2.1

dA sin ϕ sin A
=
ds p

we obtain by division

dp cos A
=− p.
dA sin A

Now

12
 2.4. The geodetic main problem

d (p sin A) dp
= sin A + p cos A =
dA dA
cos A
= − p · sin A + p cos A =
sin A
= p (− cos A + cos A) = 0.

Result:

the expression p sin A is an invariant.1

This applies on all rotationally symmetric bodies, i.e., also on the ellipsoid of revolution – where
this is called the Clairaut equation –, and of course on the plane. This invariant can be used
to eliminate the differential equation in A from the system 2.1. This yields
dϕ cos A (ϕ)
= , (2.2)
ds M (ϕ)
dλ sin A (ϕ)
= , (2.3)
ds p (ϕ)
where A (ϕ) is obtained from the invariant formula
p (ϕ1 )
sin A (ϕ) = sin A12 . (2.4)
p (ϕ)
The shape of the object is defined by giving the function p (ϕ).

2.4. The geodetic main problem

Solving the forward geodetic problem now amounts simply to substituting the also given arc
length s12 into this solution.
Computing the solution is easiest in practice using numerical integration; the methods are found
in textbooks on numerical analysis, and the routines needed in many numerical libraries.
The “classical” alternative, series expansions found in many older textbooks, are more efficient
in theory but complicated.
For the ellipsoid we may specialize the equations 2.2, 2.3 ja 2.4 with the aid of the following
expressions:

−3/2
M (ϕ) = a 1 − e2 1 − e2 sin2 ϕ


,
−1/2
p (ϕ) = N (ϕ) cos ϕ = a 1 − e2 sin2 ϕ


· cos ϕ.

2.5. The geodetic inverse problem

A direct numerical method for solving the inverse problem is iteration.

 
(1) (1)
Let us have as given ϕ1 , λ1 and ϕ2 , λ2 . First we compute the approximate values2 A12 , s12 ,
 
(1) (1)
and solve the geodetic forward problem in order to compute ϕ2 , λ2 . Thus we obtain the
closing errors
(1) (1)
∆ϕ2 ≡ ϕ2 − ϕ2 ,
(1) (1)
∆λ2 ≡ λ2 − λ2 .
2
. . . e.g., using the closed formulas of spherical trigonometry.

13
Chapter 2. The geometry of the reference ellipsoid
 
(2) (2)
Next, these closing errors could be used for computing improved values A12 , s12 , and so on.

The nearly linear dependence between (A12 , s12 ) and (ϕ2 , λ2 ) can be approximated using spher-
ical geometry. The formulas needed (1.5, 1.6):

cos s12 = sin ϕ1 sin ϕ2 + cos ϕ1 cos ϕ2 cos (λ2 − λ1 ) ,

sin s12 sin A12 = cos ϕ2 sin (λ2 − λ1 ) .

From this, by differentiation:

− sin s12 ∆s12 = [sin ϕ1 cos ϕ2 − cos ϕ1 sin ϕ2 cos (λ2 − λ1 )] ∆ϕ2 −
− cos ϕ1 cos ϕ2 sin (λ2 − λ1 ) ∆λ2 ,
sin s12 cos A12 ∆A12 + cos s12 sin A12 ∆s12 = − sin ϕ2 sin (λ2 − λ1 ) ∆ϕ2 +
+ cos ϕ2 cos (λ2 − λ1 ) ∆λ2 .

So, if we write
 
sin ϕ1 cos ϕ2 − cos ϕ1 sin ϕ2 cos (λ2 − λ1 ) − cos ϕ1 cos ϕ2 sin (λ2 − λ1 )
A= ,
− sin ϕ2 sin (λ2 − λ1 ) cos ϕ2 cos (λ2 − λ1 )
 
0 − sin s12
B= ,
sin s12 cos A12 cos s12 sin A12
we obtain as iteration formulas:
" # " # " # " # " #
(i+1) (i) (i) (i) (i)
s12 s12 ∆s12 s12 −1 ∆ϕ2
(i+1) = (i) + (i) = (i) +B A (i) .
A12 A12 ∆A12 A12 ∆λ2
 
(i+1) (i+1)
Using the new values s12 , A12 we repeat the computation of the geodetic forward problem
 
(i+1) (i+1)
to obtain new values ϕ2 , λ2 until convergence. The matrices A, B can be recomputed
if needed using better approximate values.

14
 Chapter 3
Co-ordinates on the reference ellipsoid

3.1. Representations of the sphere and the ellipsoid

An implicit representation of the circle is

x2 + y 2 − a2 = 0,

where a is the radius (Pythagoras). The parametric representation is

x = a cos β,
y = a sin β.

From this we obtain an ellipse by “squeezing” the y axis by the factor b/a, i.e.,

x = a cos β,
y = b sin β,

from which again

 x 2  y 2
+ = sin2 β + cos2 β = 1
a b
or
x2 y 2
+ 2 −1=0
a2 b
is the implicit representation.

3.2. Various latitude types

The latitude on the ellipsoid of revolution can be defined in at least three different ways. Let us
consider a cross-section of the ellipsoid, itself an ellipse; the so-called meridian ellipse.
The figure shows the following three concepts of latitude:
1. Geographical latitude ϕ: the direction of the ellipsoidal normal relative to the plane of the
equator;
2. Geocentric latitude φ (tai ψ): the angle of the line connecting the point with the origin,
relative to the plane of the equator;
3. Reduced latitude β: point P is shifted straight in the y direction to the circle around
the merian ellipse, to become point Q. The geocentric latitude of point Q is the reduced
latitude of point P .
Reduced latitude is used only in theoretical contexts. In maps, geographical (i.e., geodetic, or
sometimes, ellipsoidal) latitude is used. Geocentric latitude appears in practice only in satellite
and space geodesy.
The longitude λ is the same whether we are considering geographical, geocentric or reduced
co-ordinates.

15
Chapter 3. Co-ordinates on the reference ellipsoid

.P
b

β x
φ ϕ .
a O S R U V

Figure 3.1.: The meridian ellipse and various types of latitude

3.3. Measures for the flattening

The flattening of the reference ellipsoid is described by a variety of measures:

a−b
1. The flattening f = ;
a
a2 − b2
2. The first eccentricity (square) e2 = ;
a2
2 a2 − b2
3. The second eccentricity (square) e0 = .
b2

3.4. Relationships between different types of latitude

From figure 3.1 can be seen that

PR b
= .
QR a
In the triangles ORQ and ORP we have

QR PR
tan β = ja tan φ = ;
OR OR
so
tan φ PR b
= = ,
tan β QR a
i.e.,
a
tan β = tan φ.
b

In the triangles RV Q, RV P :
π  QR π  PR
tan −β = ja tan −ϕ = ,
2 VR 2 VR

16
 3.5. Co-ordinates in the meridional ellipse

i.e.,
cot ϕ PR b
= = ,
cot β QR a
eli
a
tan ϕ = tan β.
b
By combining still
a2
tan ϕ = tan φ.
b2

3.5. Co-ordinates in the meridional ellipse

Let us compute the co-ordinates x, y of the point P on the surface of the ellipsoid as follows.
We mark the distance P T with the symbol N , the so-called transversal radius of curvature, i.e.,
the radius of curvature of the ellipsoid in the West-East direction.
Now we have
x = N cos ϕ. (3.1)
Also
b2 2


P R = OR tan φ = OR tan ϕ = N cos ϕ · 1 − e tan ϕ
a2
using OR = x = N cos ϕ. The end result is because y = P R:

y = N 1 − e2 sin ϕ.


(3.2)

The equations 3.1, 3.2 represent a description of the meridian ellipse as a function of geodetic
latitude ϕ. Remember that N is a function of ϕ too, so not a constant! In fact

2
x2 y 2 N2 2 N 2 1 − e2
+ 2 = cos ϕ + sin2 ϕ =
a2 b a2 b2
N2 2 N 2 b2
= cos ϕ + sin2 ϕ =
a2  a2 a2
N2 b2

2 2
= cos ϕ + 2 sin ϕ = 1;
a2 a

the latter condition yields

a
N=p ,
1 − e2 sin2 ϕ
a2 − b2
using the definition e2 ≡ .
a2

3.6. Three-dimensional rectangular co-ordinates on the reference

ellipsoid

The above formulas are easily generalized: if x and y are co-ordinates within the meridian
section, the rectangular co-ordinates are

X = x cos λ = N cos ϕ cos λ,

Y = x sin λ = N cos ϕ sin λ,
Z = y = N 1 − e2 sin ϕ.

17
Chapter 3. Co-ordinates on the reference ellipsoid

T
. P
. ϕ
Q
S
.
O φ

Figure 3.2.: Suorakulmaisista koordinaateista maantieteellisiin

If we look at points not on the surface of the reference ellipsoid but above or below it in space,
we may write

X = (N + h) cos ϕ cos λ,
Y = (N + h) cos ϕ sin λ, (3.3)
Z = N 1 − e2 + h sin ϕ.

Here, h is the straight distance of the point from the surface of the ellipsoid (“ellipsoidal height”).
This quantity is interesting because satellite positioning devices can be said to directly measure
precisely this quantity (more precisely, they measure X, Y, Z and compute from these h).

3.7. Computing geographic co-ordinates from rectangular ones

This, the reverse problem from that of equations (3.3) , isn’t quite simple to solve.
Closed solutions exist, but are complicated. Of course an iterative solution based directly on
equations (3.3) is certainly possible and often used.
Computing the longitude λ is extremely simple:

Y
tan λ = .
X
A possible stumbling block is identifying the correct quadrant for λ.
ϕ and h are more complicated. See figure 3.2, where point P ’s rectangular co-ordinates X, Y, Z
are known and the geographical ϕ, λ, h are to be computed.
Let us first compute the geocentric latitude using the formula:

Z
tan φ = √
X2 +Y2
and the geometric distance (radius) by the equation:
p
OP = X 2 + Y 2 + Z 2 .

If point P would be located on the reference ellipsoid, we could determine its geographical latitude
ϕ by the following formula:
Z
tan ϕP = √
(1 − e2 ) X 2 + Y 2

18
 3.8. Meridian arc length

(proof directly from eqs. (3.3).) Now that P is located above the ellipsoid, we obtain in this
way the latitude ϕQ of point Q , where Q is the intersection of the ellipsoid and the radius of
P , for which (geometrically obviously) the above ratio is the same as for P .
Now we compute

XQ = N cos ϕQ cos λ,
YQ = N cos ϕQ sin λ,
ZQ = N 1 − e2 sin ϕQ ,

from which q
OQ = 2 + Y 2 + Z2
XQ Q Q

and thus
P Q = OP − OQ.
Additionally, in the little triangle T QP we may compute (ϕQ abbreviated to ϕ):

∠T QP = ϕ − φ

and thus

TP = P Q sin (φ − ϕ) ,
T Q = P Q cos (ϕ − φ) .

Now
h = PS ≈ TQ
and
TP
ϕP ≈ ϕQ − , (3.4)
PO
1
or perhaps a hair’s width more precise

TP
ϕP ≈ ϕQ − cos (ϕ − φ) . (3.5)
PO
This procedure is in practice fairly precise. If h = 8000 m and ϕ = 45◦ , then T P ≈ 26 m,
between the ϕP solutions (3.4) and (3.5) there is a difference of 0.1 mm2 , which is also the order
of magnitude of the error that is possibly present in the different solutions. The approximation
P S ≈ T Q contains 0.05 mm of error.3 .

3.8. Meridian arc length

The length of a meridian arc, a quantity needed, e.g., with map projections (UTM, Gauss-
Krüger) is computed by integration.
Above we already defined the quantity N , the transverse radius of curvature. The other radius
of curvature of the Earth surface is the meridional radius of curvature M . If it is given as a
function of latitude ϕ, we compute an element of path ds as follows:

ds = M dϕ.
1
Or then, not. Instead of P O:n one should take M (ϕ) + h, where M the meridional radius of curvature.
2
Linearly T P (1 − cos (φ − ϕ)).
2
3 1 TP
.
2 OP

19
Chapter 3. Co-ordinates on the reference ellipsoid

Now we may calculate the length of a meridian arc as follows:

ˆ ϕ0
s (ϕ0 ) = M dϕ. (3.6)
0

On the reference ellipsoid

a 1 − e2


M=
3/2 ,
1 − e2 sin2 ϕ
i.e., ˆ ϕ0
−3/2
2
1 − e2 sin2 ϕ


s (ϕ0 ) = a 1 − e dϕ.
0

Here, the last factor can be expanded into a series – because e2 sin2 ϕ  1 – and integrated
termwise. See the literature. Of course also a numeric approach is possible, nowaday it may
even be the superior alternative. MatLab offers for this purpose its QUAD (quadrature) routines.

20
 Chapter 4
Reference systems

4.1. The GRS80 system and geometric parameters

Nowadays, the overwhelmingly most used global geodetic reference system is the Geodetic Ref-
erence System 1980, GRS80. The parameters defining it (e.g., (Heikkinen, 1981)) are given in
table 4.1.

Some of the parameters are geometric (a), some are dynamic (J2 ,ω). Other geometric and
dynamic parameters may be derived as follows ((Moritz, 1992), following (Heiskanen and Moritz,
1967, eqs. 2-90, 2-92)):
e2 2 me0
 
J2 = 1− ⇒
3 15 q0
2me0 e2
e2 = 3J2 + .
15q0
Here(Heiskanen and Moritz, 1967, eq. 2-70))
ω 2 a2 b
m=
GM
and (based on the definitions)
be0 = ae
with the aid of which
4 ω 2 a3 e3
e2 = 3J2 + . (4.1)
15 GM 2q0
Furthermore we know ((Heiskanen and Moritz, 1967, eq. 2-58))
 
3 3
2q0 e = 1 + 0 2 arctan e0 − 0
0


e e
ja
e
e0 (e) = √ .
1 − e2
Now we may compute e2 iteratively using eq. (4.1) , which computes q0 e0 (e) anew in every

step. The result is (table 4.2):

Quantity Value Explanation

a 6378137 m semi-major axis
GM 3986005 · 108 m3 s−2 Earth mass ×G
J2 108263 · 10−8 Dynamic form factor
ω 7292115 · 10−11 rad s−1 Angular rotation rate
Table 4.1.: GRS80 defining parameters

21
Chapter 4. Reference systems

Quantity Value Explanation

e2 0.00669438002290 Ensimmäinen eksentrisyys neliöön
2
e0 0.00673949677548 Toinen eksentrisyys neliöön
b 6356752.314140 m Lyhyt akselipuolikas
1/f 298.257222101 Käänteinen litistyssuhde
Table 4.2.: GRS80 derived parameters

Quantity Value (WGS84) Remark

a 6378137 m same
1/f 298.257223563 different!
b 6356752.314245 m diff. 0.1 mm
Table 4.3.: WGS84 ellipsoidal parameters

a2 − b2 b2 1 a−b b
Here we used 1 − e2 = 1 − 2
= 2
ja 1 − = 1 − = , i.e.,
a a f a a

1 2
 
2 1 p
1−e = 1− ⇒ = 1 − 1 − e2 .
f f

Often we use (a, f ) together to define the GRS80 reference ellipsoid.

The official reference system of the GPS system is the World Geodetic System 1984 (WGS84),
whose reference ellipsoid is almost identical with GRS80. However, not exactly, table 4.3:

Most often the difference, a bit over 0.1 mm at most, can be neglected. It is apparently due to
imprecise numerics.

4.2. Gravimetric parameters

Computing geodetic parameters is done as follows ((Heiskanen and Moritz, 1967, kaavat 2-73,
2-74)):

m e0 q00
 
GM
γe = 1−m− ,
ab 6 q0
m e0 q00
 
GM
γp = 1 + ,
a2 3 q0

where ((Heiskanen and Moritz, 1967, kaava 2-67)):

  
0 0

1 1 0
q0 e = 3 1 + 0 2 1 − 0 arctan e − 1.
e e
γp − γe
The solution is again obtained iteratively, yielding (f ∗ ≡ ):
γe

22
 4.3. Reference frames

Quantity Value Explanation

γe 9.7803267715 ms−2 Normal gravity at equator
γp 9.8321863685 ms−2 Normal gravity at poles
f∗ 0.00530244011229 “Gravity flattening”

As a check, we may still compute Clairaut’s equation in its precise form ((Heiskanen and
Moritz, 1967, eq. 2-75)):
ω2b 0
 
∗ 0 q0
f +f = 1+e .
γe 2q0

A simple closed, beautiful formula for normal gravity on the reference ellipsoid is Somigliana-
Pizzetti’s formula:
aγe cos2 ϕ + bγp sin2 ϕ
γ=p .
a2 cos2 ϕ + b2 sin2 ϕ

4.3. Reference frames

Every twenty-four hours, the Earth rotates around its axis relative to the heavens at what is
very nearly a constant rate, about what it very nearly a fixed axis. The direction of this rotation
axis will serve as the z axis of both the celestial and the terrestrial frame. In order to completely
define the orientation of our reference frame, we then need to conventionally fix two longitudes:
1. On the celestial sphere: we take for this the vernal equinox, where the Sun crosses the
equator S-N
2. On the Earth: the International Meridian Conference in Washington DC, 1884, chose
Greenwich as the prime meridian.
A bonus of this choice, which was realized after the conference, was that at the same time
was defined a single, unified global time system, comprising 15◦ broad hourly time zones,
so – especially in the United States, which was expanding Westward over many time zones
– the trains would run on time.
See figure 4.1.
Red denotes an ECEF (Earth-Centred, Earth-Fixed) reference frame, which co-rotates with the
solid Earth, so the x axis always lies in the plane of the Greenwich meridian. This is
also called a CT (Conventional Terrestrial) System. Locations on the Earth’s surface are
(almost) constant in this kind of frame, and can be published, e.g., on maps. However,
moving vehicles, ocean water and atmospheric air masses will sense “pseudo-forces” (like
the Coriolis force) due to the non-uniform motion of this reference frame
blue denotes a (quasi-)inertial system, which does not undergo any (rapid) rotations relative
to the fixed stars. Also called a celestial reference frame, as the co-ordinates of the fixed
stars are nearly constant in it and may be published. Also the equations of motion of,
e.g., satellites or gyroscopes apply strictly, without pseudo-forces induced by non-uniform
reference system motion.

A Conventional Terrestrial System (CTS) is defined as follows:

◦ the origin of the frame coincides with the centre of mass of the Earth

23
Chapter 4. Reference systems

Z
Rotational motion

Greenwich
Y

Y0
X0
Spring
Equinox
X

sidereal time θ
Figure 4.1.: Geocentric reference frames

◦ the Z axis is directed along the rotation axis of the Earth, more precisely the Conventional
International Origin (CIO) , i.e., the average direction of the axis over the time span 1900-
1905
◦ the XZ-plane is parallel to the zero meridian as defined by “Greenwich”, more precisely
by: earlier the BIH (Bureau International de l’Heure, International Time Bureau), today
the IERS (International Earth Rotation and Reference Systems Service), based on their
precise monitoring of the Earth rotation.
In figure 4.2 we see the Earth orbit or ecliptic, the Earth axis tilt relative to the ecliptic plane,
and how this tilt causes the most impressive climating variation observable to human beings:
the four seasons.

4.4. The orientation of the Earth

The orientation of the Earth’s rotation axis undergoes slow changes. Relative to the stars, i.e., in
inertial space, this motion consists of precession and nutation. It is caused by the gravitational
torque exerted by the Sun and the Moon, which are either in (Sun) or close to (Moon) the
ecliptic plane. See figure 4.3.

If we study the motion of the Earth’s axis, and Earth rotation in general, relative to a reference
frame connected to the solid Earth itself, we find different quantitities:
◦ Polar motion: this consists of an annual (forced) component and a 14-months component
called the Chandler wobble.
◦ Length of Day.

24
 4.5. Transformations between systems

Celestial
North Pole

Autumnal equinox

True Earth orbit Spring

Summer
Winter
Winter Summer
solstice solstice

Autumn

Apparent
Solar
path
Ecliptic (zodiac)

Celestial equator

Vernal equinox (”Aries”)

Hic sunt dracones

Figure 4.2.: Geometry of the Earth’s orbit and rotation axis. The seasons indicated are boreal

Together these are called Earth Orientation Parameters (EOP). They are nowadays monitored
routinely, and available after the fact from the International Earth Rotation Service over the
Internet. The variation of these parameters is geophysically well understood, e.g., for the Chan-
dler wobble it is mainly the pressure of the Earth’s oceans and atmosphere that is responsible
((Gross, 2000)).

4.5. Transformations between systems

See the following diagram, which depicts only the rotations between the various systems:

Φ, Λ xp , yp θ0
Local ⇐⇒ Conventional ⇐⇒ Instantaneous ⇐⇒ Real
Astronomical Terrestrial Terrestrial Astronomical

Here, Φ, Λ are local astronomical co-ordinates (direction of the plumbline), while xp , yp are the
co-ordinates of the pole in the CIO system. θ0 is Greenwich Apparent Sidereal Time (GAST).
In the sequel we shall show how these transformations are realized in practice.

25
Chapter 4. Reference systems

α Cyg (Deneb)

Lyra
α Umi
α Lyr (Vega)
Ursa α Umi
Minor (Polaris) Precession (26 000 yr)
Dragon
Nutation (18 yr)
Ecliptic
Pole Ecliptic plane
24 hr
Moon

Sun

Figure 4.3.: Precession, nutation and the torques from Sun and Moon

Polar motion 1970-2000

2000
1995
1990
1985
1980
1975 0.8
0.6
1970 0.4
-0.5 0.2 y
0 0
0.5
1 -0.2
x
Figure 4.4.: Polar motion

26
 4.5. Transformations between systems

Reference pole (CIO) Real (instantaneous) pole

90◦ − φB
o

Plumb-
line B
Polar motion
90◦ − φA (Exaggerated)
A

90◦ − φC
Earth centre of mass

Figure 4.5.: How to monitor polar motion using latitude observatories

27
 Chapter 5
Using rotation matrices

5.1. General

Always when we change the orientation of the axes of a co-ordinate system, we have, written in
rectangular co-ordinates, a multiplication with a rotation matrix.
Let us investigate the matter in the(x, y) plane, figure 5.1.
The new co-ordinate
x0P = OU = OR cos α,
where
OR = OS + SR = xP + P S tan α = xP + yP tan α.
By substituting
x0P = (xP + yP tan α) cos α = xP cos α + yP sin α.

In the same fashion

yP0 = OT = OV cos α,
where
OV = OQ − V Q = yP − P Q tan α = yP − xP tan α,
where
yP0 = (yP − xP tan α) cos α = −xP sin α + yP cos α.
Summarizingly in a matrix equation:
 0    
x cos α sin α x
0 = .
y − sin α cos α y

The place of the minus sign is the easiest to obtain by sketching both pairs of axes on paper,
mark the angle α, and infer graphically whether the for a point on the positive x axis (i.e.:
y = 0) the new y 0 co-ordinate is positive or negative: in the above case
 
0
  x
y = − sin α cos α = − sin α · x ,
0
i.e., y 0 < 0, i.e., the minus sign is indeed in the lower left corner of the matrix.

0 y
y P
Q 0
x
V
T U
α x
O S R
Figure 5.1.: Rotation in the plane

29
Chapter 5. Using rotation matrices

SR

z z = z0
S z0
R
z 00 y 00
0
y
y β y0
y
x α
x
x0 x00 = x0
Figure 5.2.: The associativity of rotations

5.2. Chaining matrices in three dimensions

In a three-dimensional co-ordinate system we may write a two-dimensional rotation matrix as

follows:  0    
x cos α sin α 0 x
 y 0  =  − sin α cos α 0   y  ,
z0 0 0 1 z
i.e., the z co-ordinate is copied as such z 0 = z, while x and y transform into each other according
to the above formula.
If there are several transformations in sequence, we obtain the final transformation by “chaining”
the transformation matrices. I.e., if

r00 = Sr0 , r0 = Rr,

in which    
cos α sin α 0 1 0 0
R =  − sin α cos α 0  , S =  0 cos β sin β  ,
0 0 1 0 − sin β cos β
then (associativity):
r00 = S (Rr) = (SR) r,
i.e., the matrices are multiplied with each other.
Remember that
RS 6= SR,
in other words, matrices and transformations are not commutative1 !
See figure 5.2.

5.3. Orthogonal matrices

Rotation matrices are orthogonal, i.e.

RRT = RT R = I; (5.1)

their inverse matrix is the same as the transpose.

1
Two dimensional rotations are in fact commutative; they can be described also by complex numbers.

30
 5.3. Orthogonal matrices

E.g.,
 −1     
cos α sin α cos α − sin α cos (−α) sin (−α)
= = ,
− sin α cos α sin α cos α − sin (−α) cos (−α)
completely understandable, because this is a rotation around the same axis, by the same amount,
but in the opposite direction.
The above formula written in the following way:
n 
X 1 if j = k
Rij Rik = .
0 if j 6= k
i=1

The columns of a rotation matrix are orthonormal, their norm (length) is 1 and they are mutually
orthogonal. This can be seen for the case of our example matrix:
p
cos2 α + sin2 α = 1,
cos α · sin α + (− sin α) · cos α = 0.

Often we encounter other orthogonal matrices:

1. The mirror matrix for an axis, e.g.:
 
1 0 0
M2 ≡  0 −1 0  ,
0 0 1

which inverts the direction, or algebraic sign, of the y co-ordinate.

2. The axes interchange matrix (permutation):
 
0 1 0
P12 ≡  1 0 0  .
0 0 1

3. Inversion of all axes:  

−1 0 0
X =  0 −1 0  .
0 0 −1

Both M and P differ from rotation matrices in this way, that their determinant is −1, when
for rotation matrices it is +1. The determinant of the X matrix is (−1)n , with n the number
of dimensions (i.e., 3 in the above case). A determinant of −1 means that the transformation
changes a right handed co-ordinate frame into a left handed one, and conversely.
If we multiply, e.g., M2 and P12 , we obtain
 
0 1 0
M2 P12 =  −1 0 0  .
0 0 1

The determinant of this is +1. However, it is again a rotation matrix:

cos 90◦ sin 90◦ 0

 

R3 (90◦ ) =  − sin 90◦ cos 90◦ 0 !

0 0 1

All orthogonal transformations having positive determinants are rotations.

(Without proof still, that all orthogonal transformations can be written as either a rotation
around a certain axis, or a mirroring through a certain plane.)

31
Chapter 5. Using rotation matrices

Z T
x
z
z T
y
K
s

y ζ
Y
A O
x
K
X

Figure 5.3.: Local astronomical co-ordinates

5.4. Topocentric systems

Also “local astronomical”. Note that, whereas the geocentric system is “unique”, i.e., there is
only one of a certain type, there are as many local systems as there are points on Earth, i.e., an
infinity of them.
The system’s axes:
1. The z axis points to the local zenith, straight up.
2. The x axis points to the local North.
3. The y axis is perpendicular to both others and points East.
In Figure 5.3 the situation of the local topocentric system in the global context is depicted.
The spherical co-ordinates of the system are:
◦ The azimuth A
◦ The zenith angle ζ, alternatively the elevation angle () η = 100 gon − ζ.
◦ Distance s
The transformation between a point P ’s topocentric spherical co-ordinates and rectangular co-
ordinates is:
   
x s cos A sin ζ
 y  =  s sin A sin ζ  .
z T s cos ζ T
The inverse transformation:
p
x2 + y 2
ζ = arctan ,
z
y
A = 2 arctan p .
x + x2 + y 2

The latter formula is known as the half-angle formula and avoids the problem of finding the
correct quadrant for A. The result is in the interval (−180◦ , 180◦ ] and negative values may be
incremented by 360◦ to make them positive.

32
 5.5. From geocentric to topocentric and back

Astronomical Φ, Λ
co-ordinates x z
Plumb-
z line y
x
π−Φ
x y
Geoid
profile −x
Z

Y z0 Λ
z0
X Λ y0 y0
x0 Z x0

Y
X
Figure 5.4.: From the geocentric to the topocentric system. The matrix R1 mentioned in the
text was left out here

5.5. From geocentric to topocentric and back

Let (X, Y, Z) be a geocentric co-ordinate system and (x, y, z) a topocentric instrument co-
ordinate system (i.e., the x axis points to the zero direction of the instrument instead of North;
in the case of a theodolite, the zero direction on the horizontal circle.)
In this case we can symbolically write:

x = R1 R2 R3 (X − X0 ) ,

where the rotation matrices R3 , R2 , R1 act in succession to transform X into x. See figure 5.4.
X0 denotes the co-ordinates of the local origin in the geocentric system.
The inverse transformation chain of this is

X = X0 + R3T R2T R1T x,

as can be easily dervied by multiplying the first equation from the left by the matrix R1T = R1−1 ,
then by the matrix R2T , and then by the matrix R3T , and finally by moving X0 to the other side.

R3 rotates the co-ordinate frame around the z axis from the Greenwich meridian to the
local meridian of the observation site, rotation angle Λ:
 
cos Λ + sin Λ 0
R3 =  − sin Λ cos Λ 0  . (5.2)
0 0 1

Seen from the direction of the z axis we see (figure), that

x0 = x cos Λ + y sin Λ,
y 0 = −x sin Λ + y cos Λ.

(The correct algebraic signs should always be established by the aid of a sketch! Also
the directional conventions of different countries may differ.)

33
Chapter 5. Using rotation matrices

y
y0

x0

λ x

R2 turns the co-ordinate frame around the y axis, in such a way that the z axis points
to the North celestial pole instead of to the zenoth. The rotation angle needed for
this is Φ − 90◦ .  
sin Φ 0 − cos Φ
R2 =  0 1 0 . (5.3)
+ cos Φ 0 sin Φ

R1 rotates the co-ordinate frame around the new z axis or vertical axis by the amount
A0 , after which the x axis points to the azimuth of the zero point of the instrument’s
horizontal circle:  
cos A0 + sin A0 0
R1 =  − sin A0 cos A0 0  .
0 0 1

5.6. The geodetic main and inverse problems with rotation matrices

It is possible to solve the geodetic main and inverse problems three-dimensionally, without using
the surface geometry of the ellipsoid of revolution.
The idea is based on that the three-dimensional co-ordinate of a point or points are given, e.g.,
in the form (ϕ, λ, h) relative to some reference ellipsoid; and that is given or to be computed
the azimuth, elevation angle and distance of a second point as seen from the first point. In the
forward problem one should compute the secont point’s ellipsoidal co-ordinates (ϕ, λ, h).

5.6.1. Geodetic main problem

Given the ellipsoidal co-ordinates (ϕA , λA , hA ) of point A and in point A, the azimuth AAB ,
of another point B, the distance sAB , and either the elevation ηAB or the zenith angle zAB ≡
90◦ − ηAB .
Now we have to compute the co-ordinates (ϕB , λB , hB ) of point B.
As follows:
1. Transform the local A-topocentric co-ordinates of B, (AAB , sAB , zAB ) to rectangular:
   
xAB cos AAB sin zAB
xAB ≡  yAB  = sAB  sin AAB sin zAB  .
zAB cos zAB

2. Transform, using rotation matrices, these rectangular co-ordinate differences into geocen-
tric2 :
XAB = R3T R2T xAB ,
in which R3 and R2 are already given, equations (5.2) ja (5.3).
2
More precisely, into co-ordinate differences in the geocentric orientation – as in this case the origin is not the
Earth’s centre of mass!

34
 5.6. The geodetic main and inverse problems with rotation matrices

3. Add to the result the geocentric co-ordinates of point A:

XB = XA + XAB ,

jossa  
(N (ϕA ) + hA ) cos ϕA cos λA
XA =  (N (ϕA ) + hA )
cos ϕA
sin λA  .
N (ϕA ) 1 − e2 + hA sin ϕA

4. Transform the geocentric co-ordinates of B obtained back to ellipsoidal form (3.7) in the
way depicted:
XB → (ϕB , λB , hB ) .

5.6.2. Geodetic inverse problem

Given the ellipsoidal co-ordinates of two points (ϕA , λA , hA ) and (ϕB , λB , hB ). To be computed
the topocentric spherical co-ordinates of point B, AAB , zenith angle zAB ja etäisyys sAB .
1. Transform the ellipsoidal co-ordinates of A and B into geocentric co-ordinates: XA , XB .
2. Compute the relative vector between the points

XAB = XB − XA .

3. In point A, transform this vector into the topocentric rectangular system

xAB = R2 R3 XAB ;

4. transform to spherical co-ordinates by translating the formula

 
cos AAB sin zAB
xAB = sAB  sin AAB sin zAB  ,
cos zAB

with the familiar arc tangent and Pythagoras formulas (and using the half-angle formula
to avoid quadrant problems):
q
sAB = x2AB + yAB
2 + z2 ,
AB
yAB y
tan AAB = = 2 arctan qAB ,
xAB xAB + x2AB + yAB
2
q
x2AB + yAB
2
tan zAB = .
zAB

5.6.3. Comparison with ellipsoidal surface solution

The solution thus obtained is, concerning the azimuths, very close to the one obtained by using
the geodesic between the projections of points A and B on the surface of the ellipsoid. However,
not identical. The azimuths are so-called normal plane azimuths, which differ by a fraction of an
arc second from the geodesic’s azimuths even on a distance of 100 km. A very small difference,
but not zero!
The distance is of course the straight distance in space, not the length of the geodesic. This
distance is substantial.

35
 Chapter 6
Co-ordinate systems and transformations

6.1. Geocentric terrestrial systems

Generally geocentric systems, like WGS84, are defined in the following way:
1. The origin of the co-ordinate system coincides with the centre of mass of the Earth.
2. The Z axis of the co-ordinate system points in the direction of the Earth’s rotation axis,
i.e., the direction of the North pole.
3. The X axis of the co-ordinate system is parallel to the Greenwich meridian.
Such “rotating along” systems are called terrestrial. Also ECEF (Earth Centred, Earth Fixed).

6.2. Conventional Terrestrial System

To this we must however make the following further restrictions:

1. As the direction of the Z axis we use the so-called CIO, i.e., the Conventional International
Origin, the average place of the pole over the years 1900-1905.
The instantaneous or true pole circles around the CIO in a quasi-periodic fashion: the polar
motion. The main periods are a year and approx. 435 days (the so-called “Chandler
wobble”), the amplitude being a few tenths of a second of arc — corresponding on the
Earth’s surface to a few metres.
2. Nowadays the zero meridian plane no is longer based on observation by the Greenwichin
observatory, but on worldwide VLBI observations. These are co-ordinated by the Inter-
national Earth Rotation Service (IERS). So, it is no longer precisely the meridian of the
Greenwich observatory.
In this way we obtain a system co-rotating with the solid Earth, i.e., an ECEF (Earth Cen-
tered, Earth Fixed) reference system, e.g., WGS84 or ITRSxx (ITRS = International Terrestrial
Reference System, xx year number; created by the IERS). Another name used is Conventional
Terrestrial System (CTS).

6.3. Polar motion

The direction of the Earth’s rotation axis is slightly varying over time. This polar motion has
two components called xP and yP , the offset of the instantaneous pole from the CIO pole in the
direction of Greenwich, and perpendicular to it in the West direction, respectively. The trans-
formation between the instantaneous and conventional terrestrial references is done as follows:

XIT = RY (xp ) RX (yp ) XCT .

Here, note that the matrix RY denotes a rotation by an amount xP about the Y axis, i.e., the
Y axis stays fixed, while the X and Z co-ordinates change. Similarly, the matrix RX denotes a
rotation yp about the X axis, which changes only the Y and Z co-ordinates. The matrices are

37
Chapter 6. Co-ordinate systems and transformations

North Z Greenwich
pole meridian Ellipsoidal
h P normal

.
ϕ
. Y

.
. ϕ
λ
.
.

Figure 6.1.: The ellipsoid and geocentric co-ordinates

   
cos xP 0 − sin xP 1 0 0
RY (xp ) =  0 1 0  and RX (yp )  0 cos yP sin yP  .
sin xP 0 cos xP 0 − sin yP cos yP
Because the angles xp and yp are very small, order of magnitude second of arc, we may approx-
imate sin xp ≈ xp and cos xp ≈ 1 (same for yp ), as well as xP yP ≈ 0, obtaining
    
1 0 −xp 1 0 0 1 0 −xp
RY (xp ) RX (yp ) ≈  0 1 0  0 1 yp  =  0 1 yp  .
xp 0 1 0 −yp 1 xp −yp 1

6.4. The Instantaneous Terrestial System

If we take, instead of the conventional pole, the instantaneous pole, i.e., the direction of the
Earth’s rotation axis, we obtain, instead of the conventional, the so-called instantaneous terres-
trial system (ITS). This is the system to use with, e.g., astronomical or satellite observations,
because it describes the true orientation of the Earth relative to the stars.
The transformation between the conventional and the instantaneous system happens with the
aid of polar motion parameters: if they are xP , yP — xp points to Greenwich and yp to the East
from the Greenwich meridian — we obtain (see above):
      
X X 1 0 −xp X
 Y  = RY (xp ) RX (yp )  Y  ≈  0 1 yp   Y  ,
Z IT Z CT xp −yp 1 Z CT

because the angles xp , yp are so extremely small.

6.5. The quasi-inertial geocentric system

The quasi-inertial, also celestial, or real astronomical (RA) reference frame, drawn in blue in
figure 4.1, is a geocentric system, like the conventional terrestrial system. It is however celestial
in nature and the positions of stars are approximately constant in it. It is defined as follows:

38
 6.5. The quasi-inertial geocentric system

1. The origin of the co-ordinate system again coincides with the Earth’s centre of mass.
2. The Z axis of the co-ordinate system is again oriented in the direction of the Earth’s
rotation axis, i.e., the North Pole.
3. But: The X axis is pointing to the vernal equinox.
Such a reference system doesn’t rotate along with the solid Earth. It is (to good approximation)
inertial. We also refer to it as an equatorial co-ordinate system. In this system the direction
co-ordinates are the astronomical right ascention and declination α, δ.
If the “celestial spherical co-ordinates” α, δ are known, we may compute the unit direction vector
as follows:    
X cos δ cos α
 Y  =  cos δ sin α  .
Z RA sin δ

Here we must consider, however, that, while δ is given in degrees, minutes and seconds, α is
given in time units. They must first be converted to degrees etc. One hour corresponds to 15
degrees, one minute to 15 minutes of arc, and one second to 15 seconds of arc.
Going from rectangular to spherical again requires the following formulae (note the use of the
half-angle formula for α, which is precise over the whole range and avoids the problem of iden-
tifying the right quadrant):

δ = arcsin Z
Y
α = 2 arctan √ .
X+ X2 + Y 2

Again, negative values for α can be made positive by adding 24 h .

The co-ordinates, or places, of stars are apparent, a technical term meaning “as they appear at
a certain point in time 1 ”. The places α, δ read from a celestial chart are not apparent. They
refer to a certain point in time, e.g., 1950.0 or 2000.0. Obtaining apparent places requires a long
reduction chain, taking into account precession, nutation, the variations in Earth rotation, and
also the annual parallax and the possible proper motion of the star.
The apparent places of stars are found from the reference work “Apparent Places of Fundamental
Stars” precalculated and tabulated according to date.
We can transform between RA and IT (instantaneous terrestrial) co-ordinates as follows:
      
X X cos (θ0 ) − sin (θ0 ) 0 X
 Y  = RZ (−θ0 )  Y  =  sin (θ0 ) cos (θ0 ) 0   Y  .
Z RA Z IT 0 0 1 Z IT

here, θ0 on Greenwich Apparent Sidereal Time, or GAST.

1
. . . from the centre of the Earth. The fixed stars are so far way, however, that the location of the observer
doesn’t matter.

39
 Chapter 7
Reference systems and realizations

7.1. Old and new reference systems; ED50 vs. WGS84/GRS80

In Finland like in many European countries, the traditional reference system is non-geocentric
and based on an old reference ellipsoid, the International Ellipsoid computed by John Fillmore
Hayford, and adopted by the International Union of Geodesy and Geophysics (IUGG) in 1924.
European Datum 1950 (ED50) was created in 1950 by unifying the geodetic networks of the
countries of Western Europe, and was computed on the Hayford ellipsoid.
The newer systems, both World Geodetic System (WGS84) and Geodetic Reference System 1980
(GRS80) are designed to be geocentric.
So, the differences can be summarized as:
1. Reference ellipsoid used: International Ellipsoid (Hayford) 1924 vs. GRS80
2. Realized by terrestrial measurements vs. based on satellite (and space geodetic) data
3. Non-geocentric (100 m level) vs. geocentric (cm level).
The figure of a reference ellipsoid is unambiguously fixed by two quantities: the semi-major exis
or equatorial radius a, and the flattening f .
◦ International ellipsoid 1924: a = 6378388 m, f = 1 : 297
◦ GRS80: a = 6378137 m, f = 1 : 298.257222101
The reference ellipsoid of GPS’s WGS84 system is in principle the same as GRS80, but due to
poor numerics it ended up with
◦ f = 1 : 298.257223563. The net result is that the semi-minor axis (polar radius) of WGS84
is longer by 0.1 mm1 compared to GRS80.
To complicate matters, as the basis of the ITRS family of co-ordinate systems, and their real-
izations ITRF, was chosen the GRS80 reference ellipsoid. The parameters of this system were
already given in Tables 4.1 and 4.2.

7.2. WGS84 and ITRS

Both WGS84 and the International Terrestrial Reference System (ITRS) are realized by com-
puting co-ordinates for polyhedra of points (stations) on the Earth’s surface. The properties of
these systems are:
◦ Geocentric, i.e., the co-ordinate origin and centre of reference ellipsoid is the Earth’s centre
of mass (and the Earth’s mass includes oceans and atmosphere, but not the Moon!). This
is realized by using observations to satellites, whose equations of motion are implicitly
geocentric.
1
Computation:
∆f ∆ (1/f ) 0.000001462
∆b = − (a − b) = 1 (a − b) = · (21384 m) = 0.000105 m.
f ( /f ) 298.25722

41
Chapter 7. Reference systems and realizations

◦ The scale derives from the SI system. This is realized by using range measurements by
propagation of electromagnetic waves. The velocity of these waves in vacuum is conven-
tionally fixed to 299 792 458 m s−1 . Thus, range measurement becomes time measurement
by atomic clock, which is very precise.
◦ Orientation: originally the Conventional International Origin (CIO) of the Earth axis,
i.e., the mean orientation over the years 1900-1905, and the direction of the Greenwich
Meridian, i.e., the plumbline of the Greenwich transit circle. Currently, as the orientation
is maintained by the International Earth Rotation and Reference Systems Service (IERS)
using VLBI and GPS, this is no longer the formal definition; but continuity is maintained.
◦ The current definition uses the BIH (Bureau International de l’Heure) 1984 definition of
the conventional pole, and their 1984 definition of the zero meridian plane. Together, X,
Y and Z form a right-handed system.

7.3. Co-ordinate system realizations

Internationally, somewhat varying terminology is in use concerning the realization of co-ordinate

systems or datums.
◦ ISO: Co-ordinate reference system / co-ordinate system
◦ IERS: Reference system / reference frame
◦ Finnish: koordinaattijärjestelmä / koordinaatisto ((Anon., 2008))
The latter of the names is used to describe a system that was implemented in the terrain, using
actual measurements, producing co-ordinate values for the stations concerned; i.e., a realization.
Then also, a datum was defined, with one or more datum points being kept fixed at their
conventional values.
The former refers to a more abstract definition of a co-ordinate system, involving the choice of
reference ellipsoid, origin (Earth center of mass, e.g.) and axes orientation.

7.4. Realization of WGS84

Because “WGS84” is often referred to as the system in which satellite positioning derived co-
ordinates are obtained, we shall elaborate a little on how this system has been actually realized
over time. Our source is (Kumar and Reilly, 2006). The first version of WGS84 was released in
1987 by the US Defense Mapping Agency, currently the NGA (National Geospatial-Intelligence
Agency). After that, it was updated in 1994 (G730), 1996 (G873) and 2001 (G1150).
As you will see, there are a number of problems even with the latest realization of WGS84.
For this reason it is better to consider WGS84 as an approximation at best, of the reference
frames of the ITRF/ETRF variety. The precision of this approximation is clearly sub-metre,
so using WGS84 for metre precision level applications should be OK. See the following note:
ftp://itrf.ensg.ign.fr/pub/itrf/WGS84.TXT.
If you want more confusion, read (Stevenson, 2008).

7.5. Realizations of ITRS/ETRS systems

All these systems are the responsibility of the international geodetic community, specifically the
IERS (International Earth Rotation and Reference Systems Service). “I” stands for International,
“E” for European. The “S” stands for “system”, meaning the principles for creating a reference

42
 7.6. The three-dimensional Helmert transformation

frame before actual realization. With every ITRS corresponds a number of ITRF’s (“Frames”),
which are realizations, i.e., co-ordinate solutions for networks of ground stations computed from
sets of actual measurements. Same for ETRS/ETRF, which are the corresponding things for
the European area, where the effect of the slow motion of the rigid part of the Eurasian tectonic
plate has been corrected out in order to obtain approximately constant co-ordinates.
Data used for realizing ITRF/ETRF frames: mostly GPS, but also Very long Baseline Interfer-
ometry (VLBI) providing a strong orientation; satellite and lunar laser ranging (SLR, LLR) con-
tributing to the right scale, and the French DORIS satellite system. Nowadays also GLONASS
is used.
Currently the following realizations exist for ITRS: ITRF88, 89, 90, 91, 92, 93, 94; 96, 97;
ITRF2000, ITRF2005 and ITRF2008.
The definition of an ITRFyy is as follows (McCarthy, 1996):
◦ The mean rotation of the Earth’s crust in the reference frame will vanish globally (cf.
for ETRF: on the Eurasian plate). Obviously then, co-ordinates of points on the Earth’s
surface will slowly change due to the motion of the plate that the point is on. Unfortunately
at the current level of geodetic precision, it is not possible to define a global co-ordinate
frame in which points are fixed.
◦ The Z-axis corresponds to the IERS Reference pole (IRP) which corresponds to the BIH
Conventional terrestrial Pole of 1984, with an uncertainty of 0.005”
◦ The X-axis, or IERS Reference Meridian, similarly corresponds to the BIH zero meridian
of 1984, with the same uncertainty.
Finally, note that the Precise Ephemeris which are computed by IGS (the International GPS
Geodynamics Service) and distributed over the Internet, are referred to the current (newest)
ITRF, and are computed using these co-ordinates for the tracking stations used.

7.6. The three-dimensional Helmert transformation

The form of the transformation, in the case of small rotation angles, is

    (1)  
X (2)
 
1 ez −ey X tx
 (2) 
 Y  = (1 + m)  −ez 1 ex  ·  Y (1)  +  ty  , (7.1)
 
Z (2) ey −ex 1 Z (1) tz
 T
where tx ty tz is the translation vector of the origin, m is the scale factor correction
from unity, and (ex , ey , ez ) are the (small) rotation angles about the respective axes. Together
we thus have seven parameters. The superscripts (1) and (2) refer to the old and new systems,
respectively.
Eq. (7.1) can be re-written and linearized as follows, using mex = mey = mez = 0, and replacing
h iT T
the vector X (1) Y (1) Z (1) by approximate values X 0 Y 0 Z 0

. This is allowed as
m and the e angles are all assumed small.

    (1)  
X (2) − X (1)
 
m ez −ey X tx
 (2)
 Y − Y (1)  ≈  −ez m ex   Y (1)  +  ty  ≈
  
Z (2) − Z (1) ey −ex m Z (1) tz
  0   
m ez −ey X tx
≈  −ez m ex   Y 0  +  ty  .
ey −ex m Z0 tz

43
Chapter 7. Reference systems and realizations

An elaborate rearranging yields

 
m
 
 ex 
(2) (1)
  0 
Xi − Xi −Zi0 +Yi0
 
Xi 0 1 0 0  ey 
 (2) (1)  0 0
−Xi0
 
 Yi − Yi  =  Yi +Zi 0 0 1 0 
 ez .
 (7.2)
(2)
Zi − Zi
(1) Zi0 −Yi0 +Xi0 0 0 0 1 

tx 
 
 ty 
tz

Here we have added for generality a point index i, i = 1, . . . , n. The number of points is then
n, the number of “observations” (available co-ordinate differences) is 3n. The full set of these
“observation equations” then becomes
 (2) (1) 
X1 − X1 
X10 0 −Z10 +Y10 1 0 0

(2) (1)
Y1 − Y1
  0 0
   Y1 +Z1 0 −X10 0 1 0 
 (2) (1)   

 Z1 − Z1  
  Z10 −Y10 +X10 0 0 0 1 
 

 ..
 
  .. .. .. .. .. .. ..

 m
 .   . . . . . . .  


 
 

 ex 
(2) (1)
Xi0 −Zi0 +Yi0 1 0 0


 Xi − Xi
 
  0 
 ey 
Yi0 +Zi0 −Xi0 0 1 0

 (2) (1) =
 
0

ez . (7.3)
 Yi − Yi  
Zi0 −Yi0 +Xi0
   
 (2)
Zi − Zi
(1)   0 0 0 1  
tx 
    
.. .. .. .. .. .. .. ty 
   
 ..  
. . . . . . .

.
   
    tz
Xn0 −Zn0 +Yn0 1 0 0
   
0
Xn(2) − X (1)
   
Yn0 +Zi0 −Xn0 0 1 0
   
0
Yn(2) − Yn(1)
   
 
0 0 0
Zn −Yn +Xn 0 0 0 1
Zn(2) − Zn(1)

This is a set of observation equations of form ` + v = Ab x (but without the residuals vector v
which is needed to make the equality true in the presence of observational uncertainty). There
are seven unknowns b x on the right. They can be solved in the least-squares sense if we have
co-ordinates (X, Y, Z) in both the old (1) and the new (2) system for at least three points, i.e.,
nine “observations” in the observation vector `. In fact, two points and one co-ordinate from a
third point would suffice. However, it is always good to have redundancy.

7.7. Transformations between ITRF realizations

For transformation parameters between the various ITRF realizations, see the IERS web page:
http://itrf.ensg.ign.fr/trans_para.php. As an example, the transformation parameters
from ITRF2008 to ITRF2005, at epoch 2005.0, http://itrf.ensg.ign.fr/ITRF_solutions/
2008/tp_08-05.php:

T1 T2 T3 D R1 R2 R3
mm mm mm ppb 0.001” 0.001” 0.001”
-0.5 -0.9 -4.7 0.94 0.000 0.000 0.000
± 0.2 0.2 0.2 0.03 0.008 0.008 0.008
Rate 0.3 0.0 0.0 0.00 0.000 0.000 0.000
± 0.2 0.2 0.2 0.03 0.008 0.008 0.008

44
 7.7. Transformations between ITRF realizations

These parameters2 are to be used as follows:

      
X  D −R3 R2  X
 Y  (t) = 1 +  R3 D −R1   Y  (t) +
Z IT RF 2005 −R2 R1 D Z IT RF 2008
 
  
D −R3 R2 X
d
+ (t − t0 )  R3 D −R1   Y  (t) +
dt
−R2 R1 D Z IT RF 2008
   
T1 T1
d
+  T2  + (t − t0 )  T2  =
dt
T3 T3
     
 0.94 0 0  X
= 1+ 0 0.94 0  · 10−9 ·  Y  (t) +
−0 0 0.94 Z IT RF 2008
 
 
−0.5 + 0.3 (t − 2005.0)
+  −0.9 
−4.7

d
with the numbers given, and forgetting about the uncertainties. Here refers to the rates, of
dt
which only that of T1 is non-vanishing in this example.

2
Note the change in parameter names compared to the previous. For (tx , ty , tz ) we now have (T 1, T 2, T 3); m is
now called D; and (ex , ey , ez ) became (R1, R2, R3).

45
 Chapter 8
Celestial co-ordinate systems

8.1. Sidereal time

The transformation from the inertial system to the terrestrial one goes through sidereal time.
◦ Greenwich Apparent Sidereal Time, GAST, symbol θ0
◦ The apparent sidereal time at the observation location, LAST, symbol θ.

θ = θ0 + Λhms ,

in which Λ is the astronomical longitude of the place of observation, converted to suitable

time units.
GAST is
◦ the transformation angle between the inertial and the terrestrial (“co-rotating”) systems,
i.e.
◦ the angle describing the Earth’s orientation in the inertial system, i.e.
◦ the difference in longitude between the Greenwich meridian and the vernal equinox.
Also the Greenwich apparent sidereal time is tabulated – afterwards, when the precise Earth
rotation is known. GAST can be computed to one second precision based on the calendar and
civil time.
If a few minutes of precision suffices, we may even tabulate GAST as a function of day of the
year only. The table for the annual part of sidereal time per month is:

Month θm Month θm Month θm Month θm

January 6 40 April 12 40 July 18 40 October 0 40
February 8 40 May 14 40 August 20 40 November 2 40
March 10 40 June 16 40 September 22 40 December 4 40

In constructing this table, the following knowledge was used: on March 21 at 12 UTC in Green-
wich, the hour angle of the Sun, i.e., of the vernal equinox, i.e., sidereal time, is 0h . This
Greenwich sidereal time θ0 consists of two parts: an annual part θa , and a clock time τ (UTC
or Greenwich Mean Time). So we obtain the annual part of sidereal time by subtracting:
θa = θ0 − τ = −12h , i.e., 12h after adding a full turn, 24h .
If on March 21, or more precisely at midnight after March 20, sidereal time is 12h 00m , then
sidereal time for March 0 is 12h 00m − 4 × 20m = 10h 40m ; Remember that one day corresponds
to about four minutes.
We may fill out the table using the rule that 1 month ≈ 2t . (In principle we could slightly
improve the table’s accuracy by taking into account the varying lengths of the months. However,
the cycle of leap years causes error of similar magnitude.)

47
Chapter 8. Celestial co-ordinate systems

Right
o
oW 180 W 165 o ascension
165 E
o W hGr 15
0 oE α
0
15 o
30 N
o W

13
5

5
o
13 h

E
o
45 N

0W

12
12 o

0E
o
o
60 N
W

105
105 o

oE
o
75 N

90oE
90oW

75 oE
75 W
o

vernal Greenwich

E
60

60 o
o

equinox
W

True

E
o
45
o

45
W

30 o o E
Mean W λ 30
15 oW o
0o 15 E Earth
LMST rotation
LAST
GMST
GAST
Figure 8.1.: Sidereal time

After this, local sidereal time is obtained as follows:

θ = θm + θd + τ + Λhms
= θa + τ + Λhms =
= θ0 + Λhms

where
θm the value taken from the above table
θd four times the day number within a month
θa = θm + θd annual part of sidereal time
τ time (UTC)
Λ the longitude of the Earth station converted to hours, minute and seconds (15◦ = 1h , 1◦ =
4m , 10 = 4s ).

See the pretty figure 8.1. We have the following quantities:

◦ GAST = Greenwich Apparent Sidereal Time (= θ0 )

48
 8.2. Trigonometry on the celestial sphere

apparent rotation
of the celestial sphere
Zenith
North Meridian
pole
h Equator

δ
Apparent
diurnal
motion
East

North Observer h

South
West
Figure 8.2.: Hour angle and other co-ordinates on the celestial sphere

◦ LMST = Local Mean Sidereal Time (= θ )

◦ the equinox varies irregularly with time due to precession and nutation. That’s why we
distinguish Mean and Apparent. The difference is called the equation of equinoxes, ee.
◦ h is the hour angle
◦ hGr is the Greenwich hour angle
◦ α is the right ascension (of a celestial object)
◦ Λ is the longitude (of a terrestrial object).

h = θ − α,
hGr = θ0 − α,
θ = θ0 + Λ,
θ = θ0 + Λ;
θ − θ = θ0 − θ0 = ee.

8.2. Trigonometry on the celestial sphere

On the celestial sphere we have at least two different kinds of co-ordinates: local and equatorial.
Local spherical co-ordinates are related to the local astronomical rectangular system as follows
(x to the North, z to the zenith):

     
x cos A sin ζ cos A cos η
 y  =  sin A sin ζ  =  sin A cos η  ,
z cos ζ sin η
where A is the azimut (from the North clockwise), ζ is the zenith angle and η is the height or
elevation angle.
Equatorial co-ordinates are α, δ, right ascension and declination; their advantage is, that the
co-ordinates (“places”) of stars are nearly constant. The disadvantage is that there is no simple
relationship to local astronomical co-ordinates.

49
Chapter 8. Celestial co-ordinate systems

Celestial pole

West h East

90◦ − δ
90◦ − Φ

Zenith A
ζ
Star
Figure 8.3.: Fundamental triangle of astronomy

An intermediate form of co-ordinates is: hour angle and declination, h, δ. The hour angle is
defined as shown in figure 8.3, the angular distance around the Earth’s rotation axis (or at the
celestial pole) from the meridian of the location.
Equation:
h = θ0 + Λ − α,

where θ0 is Greenwich apparent sidereal time (GAST), Λ is the local (astronomical) longitude,
and α is the right ascension of the star.
If the star is in the meridian, we have h = 0 and α = θ0 + Λ. On this is based the use of a transit
instrument: if, of the three quantities θ0 , α or Λ, two are knownthe third may be computed.
According to the application we speak of astronomical position determination, time keeping of
determination of the places of stars. “One man’s noise is another man’s signal ”.
On the celestial sphere there is a fundamental triangle of astronomy: it consists of the star, the
celestial North pole, and the zenith. Of the angles of the triangle we mention t (North pole) and
A (the zenith), of its sides, 90◦ − Φ (pole-zenith), 90◦ − δ (star-pole) and ζ (star-zenith).
The sine rule:

− sin A sin h
= .
cos δ sin ζ

The cosine rule:

cos ζ = sin δ sin Φ + cos δ cos Φ cos h,

sin δ = sin Φ cos ζ + cos Φ sin ζ cos A.

We compute first, using the cosine rule, either δ or ζ, and then, using the sine rule, either h or
A. Thus we obtain either (A, ζ) ↔ (h, δ) in both directions.

8.3. Using rotation matrices

The various transformations between celestial co-ordinate systems can be derived also in rect-
angular co-ordinates by using rotation matrices.

50
 8.3. Using rotation matrices

 T
Let the topocentric co-ordinate vector (length 1) be r = x y z . In spherical co-ordinates
this is    
x cos A sin ζ
r =  y  =  sin A sin ζ  .
z cos ζ
The same vector we may also write in local astronomical co-ordinates:
 0   
x cos α cos δ
r0 =  y 0  =  sin α cos δ  .
z0 sin δ
The transformation between them is:
1. the direction of the x axis is changed from North to South
2. the new xz axis pair is rotated by an amount 90◦ −Φ to the South, Φ being the astronomical
latitude.
3. the new xy axis pair is rotated to the West (clockwise) from the local meridian to the
vernal equinox by an amount θ, the local (apparent) sidereal time.
The matrices are:
 
−1 0 0
M1 =  0 1 0 ,
0 0 1
cos (90◦ − Φ) 0 sin (90◦ − Φ)
   
sin Φ 0 − cos Φ
R2 =  0 1 0 = 0 1 0 ,
◦ ◦
− sin (90 − Φ) 0 cos (90 − Φ) cos Φ 0 sin Φ
 
cos θ − sin θ 0
R3 =  sin θ cos θ 0  .
0 0 1

Let us compute
  
cos θ − sin θ 0 − sin Φ 0 cos Φ
R3 R2 M1 =  sin θ cos θ 0  0 1 0 =
0 0 1 cos Φ 0 sin Φ
 
− cos θ sin Φ − sin θ cos θ cos Φ
=  − sin θ sin Φ cos θ sin θ cos Φ  .
cos Φ 0 sin Φ

After this:
    
cos α cos δ − cos θ sin Φ − sin θ cos θ cos Φ cos A sin ζ
 sin α cos δ  =  − sin θ sin Φ cos θ sin θ cos Φ   sin A sin ζ  ,
sin δ cos Φ 0 sin Φ cos ζ
where we identify immediately
sin δ = sin Φ cos ζ + cos Φ sin ζ cos A,
the cosine rule in the triangle star-celestial pole-zenith.
The inverse transformation is, based on orthogonality (the transpose!)
    
cos A sin ζ − cos θ sin Φ − sin θ sin Φ cos Φ cos α cos δ
 sin A sin ζ  =  − sin θ cos θ 0   sin α cos δ  .
cos ζ cos θ cos Φ sin θ cos Φ sin Φ sin δ

With these, we can do the transformation of spherical co-ordinates with the aid of three-
dimensional “direction cosines”.

51
Chapter 8. Celestial co-ordinate systems

8.4. On satellite orbits

We describe here the computation of circular satellite orbits around a spherical Earth; this is
sufficient at least for enabling visual satellite observations and finding the satellites.
There are thousands of satellites orbiting the Earth, which for the most part are very small. A
few hundred however are so large, generally last stages of launcher rockets, that they can be seen
after dark in the light of the Sun even with the naked eye. With binoculars these can be observed
easily. The heights of their orbits vary from 400 km to over 1000 km; the inclination angle of
the orbital plane may vary a lot, but certain inclination values, like 56◦ , 65◦ , 72◦ , 74◦ , 81◦ , 90◦
and 98◦ are especially popular.
In a class of their own are the Iridium satellites, which have stayed in orbit after an ill-fated
mobile telephone project. Every Iridium satellite has a long metal antenna which reflects sunlight
in a suitable orientation extremely brightly. Predictions for the Iridium satellites are found from
the World Wide Web.
When we know the satellite’s equator crossing, i.e., the time of the “ascending node” of the orbit,
t0 , and the right ascension (”inertial “longitude”) α0 , we can compute the corrections to time
and longitude for various latitudes like we describe in the sequel.

8.5. Crossing a given latitude in the inertial system

If the target latitude ϕ is given, we can compute the distance ν from the ascending node
(“downrange-angle”) in angular units as follows:
sin φ
sin ν = ,
sin i
where i is the inclination of the satellite orbit. From this follows again the elapsed time using
Kepler’s third law; the period or time of completing one orbit is:
r
4π 2 3
P = a .
GM
From this
ν
∆τ = P,
2π
the flight time from the equator to latitude φ.
The azimut angle between satellite orbit and local meridian is obtained from the Clairaut
formula (luku 2.3):
cos φ sin A = cos (0) cos i,
π
because at the equator (φ = 0) A = − i, i.e.
2
cos i
sin A = .
cos φ

Now we obtain the difference in right ascension with the equator crossing using the sine rule for
a spherical triangle:
sin ∆α sin φ tan φ
= ⇒ sin ∆α = .
sin A sin i tan i
After this, we obtain the satellite’s right ascension and time when crossing latitude ϕ:
τ = τ0 + ∆τ,
α = Ω + ∆α.
Here, Ω is the right ascension of the ascending node of the satellite orbit.
Both longitude and right ascension α (ja ∆α) are reckoned positive to the East.

52
 8.6. The satellite’s topocentric co-ordinates

θ
r
R
A
Φ
ν ϕ

i
Kevät-
tasaus Ω ∆α

Figure 8.4.: The places of satellite and Earth station in the inertial system

8.6. The satellite’s topocentric co-ordinates

Generally, the satellite orbit’s height h is given; according to its definition, the height of the
satellite from the geocentre is
r = krk = ae + h,
where ae is the equatorial radius of the Earth. In spherical approximation ae = R. Then the
rectangular, inertial co-ordinates of the satellite are
 
r cos φ cos αG
r =  r cos φ sin αG  ,
r sin φ

where φ is the satellite’s geocentric latitude (or, if one wants to put it this way, the “geocentric
declination” δG ) and αG the geocentric right ascension.
Let the geocentric latitude of the ground station be Φ and its longitude Λ; then, at moment t,
the geocentric right ascension of the ground station is

θ = Λ + τ + θa ,

where τ is the time (UTC) and θa the annual part of sidereal time. θ is the same as the local
sidereal time, which thus represents the orientation of the local meridian in (inertial) space.
Now the Earth station’s rectangular, inertial co-ordinates are
 
R cos Φ cos θ
R =  R cos Φ sin θ  .
R sin Φ

Subtraction yields:  
d cos δT cos αT
d = r − R ≡  d cos δT sin αT  ,
d sin δT
from which we may solve the topocentric co-ordinates:

d3 d2
tan δT = p 2 2
, tan αT =
d1 + d2 d1

53
Chapter 8. Celestial co-ordinate systems

These are thus the satellite’s right ascension and declination as seen against the starry sky. With
the aid of a star chart and binoculars we may now wait for and observe the satellite.
We also obtain the satellite’s distance:
q
d= d21 + d22 + d23 .
Using the distance we may compute the satellite’s visual brightness of magnitude. The distances
are generally of order 500-1000 km and the magnitudes 2-5.

8.7. Crossing a given latitude in the terrestrial system

If we do the computation in the terrestrial system, we must have been given the longitude of the
equator crossing λ0 (λ0 = Ω − θ0 (τ0 ), where θ0 is GAST at the time τ0 of the equator crossing.)
The time difference ∆τ from the equator crossing to latitude Φ is obtained in the same way;
however, now we compute the longitude difference as follows:
∆λ = ∆α − ∆τ · ω,
where ω is the rotation rate of the Earth, some 0.25◦ per minute. Thus we obtain the satellite’s
longitude:
λ = λ0 + ∆λ.

The geocentric co-ordinates of both the ground station and the satellite can also be described
in the terrestrial system. In this system the co-ordinates of the ground station are
 
cos Φ cos Λ
R = R  cos Φ sin Λ 
sin Φ
and the satellite co-ordinates  
cos ϕ cos λ
r = r  cos ϕ sin λ  .
sin ϕ
From this again we obtain the terrestial topocentric vector:
 
d cos δT cos (Λ − hT )
d = r − R ≡  d cos δT sin (Λ − hT )  ,
d sin δT
from which we may solve the topocentric declination δT and hour angle hT .
If we then want the right ascension for use with a celestial chart, all we need to do is subtract
it from the local sidereal time:
αT = θ − hT = (θ0 + Λ) + hT = θ0 + (Λ − hT .) .
In this, θ0 is Greenwich sidereal time, and (Λ − hT ) the ”hour angle of Greenwich”.

8.8. Determining the orbit from observations

If the satellite’s place on the sky has been observed at two different points in time τ1 and τ2 –
i.e., we have as given (α (τ1 ) , δ (τ1 )) and (α (τ2 ) , δ (τ2 )) – we may compute firstly the topocentric
direction vectors (unit vectors):
   
cos δ (τ1 ) cos α (τ1 ) cos δ (τ2 ) cos α (τ2 )
d (τ1 ) d (τ2 )
e1 = =  cos δ (τ1 ) sin α (τ1 )  , e2 = =  cos δ (τ2 ) sin α (τ2 )  .
d1 (τ1 ) d2 (τ2 )
sin δ (τ1 ) sin δ (τ2 )

54
 8.8. Determining the orbit from observations

Zenith
eU Moment τ1

Moment τ2

Satellite orbit α, δ
e

d
r
d

Horizon η
R Celestial
sphere
Observer
r

O
Figure 8.5.: Satellite orbit determination from two visual observations

See Figure 8.5. When also the ground station vector R has been computed, we may compute
d1 , d2 by the cosine rule if a suitable value1 for the satellite height h — or equivalently, for the
satellite orbit’s radius r = R + h – has been given:

r2 = R2 + d2 + 2Rd cos (∠e, R) ⇒

⇒ d2 + 2Rd sin η + R2 − r2 = 0 ⇒
p
−2R sin η ± 4R2 sin2 η − 4 (R2 − r2 )
⇒ d(1,2) = =
q 2
= −R sin η ± r2 − R2 1 − sin2 η ,

in which sin η = cos (∠e, R) = he · eU i is the projection of the satellite’s direction vector onto
the local vertical, and  
cos Φ cos θ
R 
eU = = cos Φ sin θ  ,
R
sin Φ
the “zenith vector” of the observer, a unit vector pointing straight upward. The value sin η, the
sine of the elevation angle, may be calculated directly as the dot product he · eU i when both
1
An experienced observer can guess a satellite’s height from the observed motion in the sky surprisingly precisely.

55
Chapter 8. Celestial co-ordinate systems

vectors are given in rectangular form.

Only the positive solution makes physical sense:
q
d = r2 − R2 1 − sin2 η − R sin η.

Thus we obtain d1 , d2 and thus d1 = d1 e1 , d2 = d2 e2 . For the satellite velocity we obtain

kd2 − d1 k
v= .
τ2 − τ1
We know however what the velocity for a circular orbit should be at height h: according to
Kepler’s third law
r r
4π 2 2π (R + h) GM
P = (R + h)3 ⇒ vk = = .
GM P R+h
Let us prepare the following little table:

Height (km) 500 750 1000 1500

Speed (m/s) 7612.609 7477.921 7350.139 7113.071

We see that the satellite’s linear speed of flight decreases only slowly with height. Therefore we
may use the “observed velocity” v for correcting the height h according to the following formula:
vk
h0 = h ,
v
where vk is the velocity according to Kepler (from the table for the value h), v the calculated
velocity, and h0 the improved value for the satellite height. This process converges already in
one step to almost the correct height.
Note that the height thus obtained is, in the case of an elliptical orbit, only (approximately) the
“ overflight height”! The real height will vary along the orbit. For the same reason, the height
obtained will not be good enough to calculate the period P (or equivalently: the semi-major
axis a) from. The real period can be inferred only of the satellite has been observed for at least
a couple of successive days.
The computed vector for the change in position between the two moments of observation, d2 −
d1 = r2 − r1 also tells about the inclination of the orbit, by calculating its vectorial product
with, e.g., the satellite’s geocentric location vector r1 , as follows:

(d2 − d1 ) × r1 = kd2 − d1 k kr1 k cos i,

from which i may be solved. And if i and some “footpoint” (ϕ, λ) of the satellite is known, also
Ω and (with P ) τ0 may be computed. Then we can already generate predictions!
When the height of the satellite (and its approximate period) as well as the inclination are
known, we can also compute the rapid precessional rate of the ascending node2 :

dΩ
Ω (τ ) = Ω (τ0 ) − (τ − τ0 ) =
dτ r
3 GM  ae 2
= Ω (τ0 ) − (τ − τ0 ) · J2 cos i,
2 a3 a
2
The formula applies for circular orbits.

56
 8.8. Determining the orbit from observations

in which J2 is the dynamic flattening of the Earth, value J2 = 1082.62 · 10−6 . One of the first
achievements of the satellite era was the precise determination of J2 3 . As a numerical formula:

dΩ cos i 
= −6.52927 · 1024 3.5 m3.5 degrees/day

dτ a
(note the unit!) Thus we obtain the following table (unit degree/days):

Ht./Incl. 500 750 1000 1500

0◦ -7.651 -6.752 -5.985 -4.758
56◦ -4.278 -3.776 -3.347 -2.661
65◦ -3.233 -2.854 -2.529 -2.011
74◦ -2.109 -1.861 -1.650 -1.311
81◦ -1.197 -1.056 -0.936 -0.744
90◦ 0 0 0 0
0.9856 97◦ .401 98◦ .394 99◦ .478 101◦ .955

In the table the last row is special. The value 0.9856 degree/day is the Sun’s apparent angular
velocity relative to the background of the stars. If the precessional velocity of the satellite’s
orbital plane is set to this value, the satellite will always fly over the same area at the same
local solar time. In this way we achieve a so-called heliostationary orbit, also known as “no-
shadow” orbit, the advantages of which are continuous light on the solar cells, and in remote
sensing, always the same elevation angle of the Sun when imaging the Earth surface. Given is
the inclination angle which, for each value of the mean height for each column of the table, gives
precisely such a sun-stationary orbit. The precession rate of the satellite’s orbital plane relative
to the sun is now
dΩ
q≡ − 0.9856 degree/day
dτ
E.g., for a satellite whose height is 500 km and inclination 56◦ :

q = −4.278 − 0.9856 = −5.2636 degree/day,

and the period of one cycle is 360/q = 68 days ≈ 2.2 months. This is the time span after which
the satellite appears again, e.g., in the evening sky of the same latitude.

3
In fact, this motion of the ascending node is so rapid, that without considering it it is not possible to generate
usable orbit predictions.

57
 Chapter 9
The surface theory of Gauss

Carl Friedrich Gauss1 was among the first to develop the theory of curved surfaces. The theory
he developed assumed still, that the two-dimensional surface is inside three-dimensional space
– many derivations are then simple and nevertheless the theory applies as such to the curved
surface of the Earth: note that Gauss was a geodesist who measured and calculated the geodetic
networks of Hannover and Brunswick using the method of least squares.
Let a curved surface S be given in three-dimensional space R3 . The surface is parametrized by
the parameters (u, v). Example: the surface of the Earth, parametrization (ϕ, λ).
In three-dimensional space, we may describe point positions by creating an orthonormal triad
of base vectors
{e1 , e2 , e3 } .
On this base, a point, or location vector, is

x = x1 e1 + x2 e2 + x3 e3 = xe1 + ye2 + ze3 .

We will now often choose to write x by its representation on this basis:

 1   
x x
x =  x2  =  y  .
x3 z

These three parameters thus form a parametrization of R3 .

9.1. A curve in space

A curve running through space C can be parametrized by a parameter s. Then, the points on
the curve are
x (s) .
If it holds for the parameter, that

ds2 = dx2 + dy 2 + dz 2 ,

we say that C is parametrized according to distance (or arc length).

Examples:

1. The numbers on a measurement tape form a parametrization by distance.

2. If you drive along a road, the trip meter readings form a parametrization of the road
according to distance.
3. Along Mannerheimintie, the kkj co-ordinate x constitutes a parametrization, however not
according to distance (because the direction of the road is variable).
1
Carl Friedrich Gauss, (1777 – 1855), was also among the first to speculate on the possibility of non-Euclidean
geometry. Others were János Bolyai and Nikolai I. Lobachevsky.

59
Chapter 9. The surface theory of Gauss

Let us assume in the sequel, that the parametrization used is unambiguous and differentiable
(and thus continuous).
The tangent of the curve C is obtained by differentiation:

dx (s)
t (s) = = xs ;
ds
the length of the tangent is
s   2  2 r 2
dx 2 dy dz dx + dy 2 + dz 2 √
ktk = + + = = 1 = 1.
ds ds ds ds2

This applies only if s is a parametrization by distance.

An arbitrary parametrization t can always be converted into a parametrization by distance in
the following way:
ˆ t ˆ tp 2
ds dx + dy 2 + dz 2
s(t) = dτ = dτ = (9.1)
0 dτ 0 dτ
ˆ t  2  2  2
s
dx dy dz
= + + dτ,
0 dτ dτ dτ

i.e., in differential form

s 2  2  2
dx (t) dy (t) dz (t)
ds = dt + + ,
dt dt dt

from which s (t) can always be computed by integration 9.1.

Another differentiation yields the curvature vector :

dt (s) d2
k (s) = = 2 x (s) .
ds ds

9.2. The first fundamental form (metric)

The Gauss first fundamental form:

ds2 = Edu2 + 2F dudv + Gdv 2 . (9.2)

ds is a path element within the surface.

Later we shall see that this fundamental form is the same as the metric of the surface under
consideration, and an alternative way of writing is

ds2 = g11 du2 + g12 dudv + g21 dvdu + g22 dv 2 .

If the point x is on the surface S, we may take the derivative of its co-ordinates:
   
∂x/∂u ∂x/∂v
∂x  ∂x 
xu = = ∂y/∂u  , xv = = ∂y/∂v  .
∂u ∂v
∂z/∂u ∂z/∂v

These vectors are called the tangent vectors of surface S and parametrization (u, v).

60
 9.3. The second fundamental form

v + dv

xv

v
xu
G /2
1

1/2
x E

u + du
u

O
Figure 9.1.: The Gauss first fundamental form

From these we we obtain, as “dot products” of two vectors in space, the elements of the funda-
mental form:

E = hxu · xu i ,
F = hxu · xv i ,
G = hxv · xv i .

In figure 9.1 we have depicted the Gauss first fundamental form and the tangent vectors xu , xv .
It is easy to show (chain rule in three dimensions (x, y, z)), that

ds2 = dx2 + dy 2 + dz 2 =
"   2  2 #
∂x 2 ∂y ∂z
= + + du2 +
∂u ∂u ∂u
     
∂x ∂x ∂y ∂y ∂z ∂z
+2 + + dudv +
∂u ∂v ∂u ∂v ∂u ∂v
"   2  2 #
∂x 2 ∂y ∂z
+ + + dv 2 =
∂v ∂v ∂v
= hxu · xu i du2 + 2 hxu · xv i dudv + hxv · xv i dv 2 ,

from which (9.2) follows directly.

We see, e.g., that in the direction of the v curves (dv = 0):

ds2 = Edu2 ,

in other words, E represents the metric distance between two successive curves (u, u + 1).
Similarly G represents the distance between two successive v curves. The closer the curves
are to each other, the larger xu or xv and also the larger E or G. F again represents the angle
between the u and v curves: it vanishes if the angle is straight.

9.3. The second fundamental form

The normal on a surface is the vector which is orthogonal to every curve running within the
surface, also the u and v curves. We write
 
nx
n =  ny  .
nz

61
Chapter 9. The surface theory of Gauss

Apparently
hn · xu i = hn · xv i = 0. (9.3)
Let us also require
knk = 1,
or
n2x + n2y + n2z = 1.

We differentiate x a second time:

∂2x
xuu = ,
∂u2
∂x
xuv = ,
∂u∂v
∂x
xvv = .
∂v 2
Now, Gauss’s second fundamental form is

edu2 + 2f dudv + gdv 2 ,

where

e = hn · xuu i ,
f = hn · xuv i , (9.4)
g = hn · xvv i .

Based on condition (9.3) we also have

∂
0= hn · xu i = hnu · xu i + hn · xuu i ⇒ e = −nu · xu ,
∂u
∂
0= hn · xu i = hnv · xu i + hn · xuv i ⇒ f = −nu · xv ,
∂v
∂
0= hn · xv i = hnu · xv i + hn · xuv i ⇒ f = −nv · xu , (9.5)
∂u
∂
0= hn · xv i = hnv · xv i + hn · xvv i ⇒ g = −nv · xv .
∂v

See figure 9.2. When moving from location to location, the normal vector’s direction changes:
when we move from point x (u, v) to point x0 (u + du, v), the normal changes n → n0 = n + dn’.
In the same way, when moving from point x to x00 (u, v + dv), the normal changes n → n00 =
n + dn00 .
As a formula:
∂n ∂n
dn = du + dv = nu du + nv dv.
∂u ∂v

The norm of the normal vector, i.e., its length, is always 1, like we saw above; that’s why the
vector can change only in two directions, either in the direction of tangent vector xu , or in the
direction of tangent vector xv .
Let us separate them by projection:

hdn · xu i = hnu · xu i du + hnv · xu i dv,

hdn · xv i = hnu · xv i du + hnv · xv i dv,

62
 9.4. Principal curvatures

−gdv
−f du −edu
n + nv dv −f dv
n + nu du

n
xv
xu

v + dv
x u + du

O v u
Figure 9.2.: Explaining geometrically the Gauss second fundamental form

in which we may directly identify the elements of the second fundamental form e, f, g with the
aid of (9.5):
− hdn · xu i = edu + f dv,
− hdn · xv i = f du + gdv.

Contrary to the first fundamental form, the second fundamental form has no corresponding
object in the Riemann surface theory. It exists only for surfaces embedded in a surrounding
(Euclidean) space.
Often we can also find a “tensorial” notation:
β11 = e,
β12 = β21 = f,
β22 = g.

9.4. Principal curvatures

The Gauss second fundamental form describes in a way the curvature of a surface in space, by
depicting how the direction of the normal vector changes, when we travel either in the u or in the
v co-ordinate curve direction. Unfortunately this is not enough for an absolute characterization
of the curvature, because the parametrization (u, v) is
1. an arbitrary choice, and
2. not metrically scaled.
The latter means that if the direction of the normal vector n changes by an amount dn when
we travel a distance du along the v co-ordinate curve, we still don’t know how many metres the
distance du corresponds to. If it is a long distance, then the same change in the normal vector
dn means only a small curvature of the surface; if it is a short distance, then the same change
dn corresponds to a large surface curvature.
The fact that the parameter curves u = constant and v = constant generally not are perpendic-
ular to each other, makes this problem even trickier.
Write the first and second fundamental forms in matrix form:
   
E F e f
H= ,B= .
F G f g

63
Chapter 9. The surface theory of Gauss

Form the matrix:

   
−1 G −F 1 e f
C≡H B= =.
−F E EG − F 2 f g
 
1 Ge − F f Gf − F g
=
EG − F 2 Ef − F e Eg − F f

This is called the shape operator, and the above, the Weingarten2 equations, see http://en.
wikipedia.org/wiki/Differential_geometry_of_surfaces#Shape_operator
We can say that multiplication by the inverse of the H matrix, i.e., the first fundamental form,
which describes the length of a distance element, performs a metric scaling of the B matrix 3 .

Principal curvatures:
The matrix C has two eigenvalues: the values κ1,2 for which

(C − κI) x = 0 (9.6)
 T
for suitable value pairs x = du dv . The solutions κ1,2 are called the principal curvatures
of surface S. They are invariant with respect to the chosen parametrization (u, v). The cor-
 T  T
responding value pairs x1 = du1 dv1 and x2 = du2 dv2 define the local principal
directions of curvature on the surface.
Other invariants:
det B eg − f 2
1. The product κ1 κ2 = det C = = is the Gaussian curvature.
det H EG − F 2
1 1 1 eG + gE − 2f F
2. The half-sum (κ1 + κ2 ) = (C11 + C22 ) = is the mean curvature.
2 2 2 EG − F 2

Principal directions of curvature:

We may also study the eigenvectors of C, which are called principal directions of curvature. For
this, write

0 = H (C − κ1 I) x1 = (B − κ1 H) x1 ,
0 = H (C − κ2 I) x2 = (B − κ2 H) x2 .

Multiply the first from the left with xT2 , and the second with xT1 and transposing:

0 =xT2 Bx1 − κ1 xT2 Hx1 ,

T
0 = xT1 Bx2 − κ2 xT1 Hx2 = xT2 Bx1 − κ2 xT2 Hx1 ,

as, both B and H being symmetric matrices, we have

T
xT1 Bx2 = xT2 Bx1 and

T
xT1 Hx2 = xT2 Hx1 .
2
Julius Weingarten, German mathematician 1836-1910
3
A more technical description: B is a covariant tensor βij , and H the covariant metric tensor gij . H −1 again
corresponds to the contravariant tensor g ij . C is now the “mixed tensor” βki = g ij βjk , whose tensorial
eigenvalue problem is  
βji − κδji xj = 0,
the same equation as (9.6).
See chapter 10.

64
 9.4. Principal curvatures

Subtraction of the two yields

(κ2 − κ1 ) xT2 Hx1 = 0.
This shows that if the principal radii of curvature are different 4 , then the expression xT2 Hx1
vanishes. This expression5 can be interpreted as an inner product: in fact, in Cartesian plane
co-ordinates in the tangent plane, the matrix H = I, and we have xT2 x1 = 0, or x1 ⊥x2 in the
Euclidean sense.
The principal directions of curvature are mutually perpendicular.
This is just a special case of self-adjoint (symmetric) operators having mutually orthogonal
eigenvectors, e.g., the eigenfunctions of Sturm-Liouville theory (http://en.wikipedia.org/
wiki/Sturm%E2%80%93Liouville_theory).

Example:
The co-ordinates of a point on the surface of the ellipsoidal Earth are
 
N (ϕ) cos ϕ cos λ
x =  N (ϕ) cos ϕ sin
λ
.
2
N (ϕ) 1 − e sin ϕ

From this  d 
cos λ (N (ϕ) cos ϕ)
 dϕ 
∂x  d 
xϕ = = sin λ (N (ϕ) cos ϕ) ;
 
∂ϕ  dϕ 

d 
1 − e2 (N (ϕ) sin ϕ)
dϕ
Using the formulas derived in the Appendix B we obtain:
 
− sin ϕ cos λ
xϕ = M (ϕ)  − sin ϕ sin λ  .
+ cos ϕ
 
− cos ϕ sin λ
∂x
xλ = = N (ϕ)  + cos ϕ cos λ  .
∂λ
0
The surface normal is obtained as the vectorial product, normalized:
hxϕ × xλ i
n= ,
kxϕ × xλ k

where
− cos2 ϕ cos λ
 

hxϕ × xλ i = N M  − cos2 ϕ sin λ =

2 2
− sin ϕ cos ϕ cos λ − sin ϕ cos ϕ sin λ
cos2 ϕ cos λ
   
cos λ
= −N M  cos2 ϕ sin λ  = −N M cos2 ϕ  sin λ  ,
sin ϕ cos ϕ tan ϕ

the norm of which is

p
kxϕ × xλ k = N M cos2 ϕ 1 + tan2 ϕ = N M cos ϕ.
4
And if they are not, any linear combination of x1 and x2 will again be an eigenvector, and we can always choose
two that are mutually perpendicular.
5
In index notation: gij xi2 xj1 .

65
Chapter 9. The surface theory of Gauss

Thus  
cos ϕ cos λ
n = −  cos ϕ sin λ  ,
sin ϕ
not a surprising result. . .
Let us compute the first fundamental form:

E = hxϕ · xϕ i = M 2 sin2 ϕ sin2 λ + cos2 λ + cos2 ϕ = M 2 ,

F = hxϕ · xλ i = 0,
G = hxλ · xλ i = N 2 cos2 ϕ = p2 .

We calculate for use in calculating the second fundamental form

   
+ sin ϕ cos λ + cos ϕ sin λ
nϕ =  + sin ϕ sin λ  ; nλ =  − cos ϕ cos λ 
− cos ϕ 0

and thus (equations 9.5)

e = − hnϕ · xϕ i = +M,
f = − hnϕ · xλ i = − hnλ · xϕ i = 0,
g = − hnλ · xλ i = +N cos2 ϕ.

In other words:

M2
   
E F 0
H = = ,
F G 0 N 2 cos2 ϕ
   
e f M 0
B = = , and
f g 0 N cos2 ϕ
1
 
0
C = H −1 B =  M 1  .
 
0
N

The principal curvatures κ1,2 are C’s characteristic values or eigenvalues, solutions of the eigen-
value problem
(C − κI) x = 0;

we obtain the values by solving

1
 
 M −κ 0
det (C − κI) = 0 ⇒ det  =0 ⇒

1
0 −κ
   N
1 1 1 1
⇒ −κ − κ = 0 ⇒ κ1 = , κ2 = .
M N M N

As could be expected. . .

9.5. A curve embedded in a surface

If a curve C runs inside a curved surface, we may study other interesting things.

66
 9.5. A curve embedded in a surface

The tangent vector

 1   
6 i t du/ds
If we call t = = “the t vector’s components in the (u, v) co-ordinate system”,
t2 dv/ds
dx dx du dx dv
this is t = = + = t1 xu + t2 xv . In other words, the tangent to the curve is also
ds du ds dv ds
one of the tangents to the surface, and lies inside the tangent plane.

The curvature vector

dt d 1
t xu + t 2 xv =


k= =
ds ds
dt1 dt2
= xu + xv +
ds ds

2
2
+ xuu t1 + 2xuv t1 t2 + xvv t2 .

In other words, the “curvature vector’s components in the (u, v) co-ordinate system” ontain other
things besides just the derivatives of the component values dt1 /ds ja dt2 /ds.

We write

xuu = Γ111 xu + Γ211 xv + hxuu · ni n,

xuv = Γ112 xu + Γ212 xv + hxuv · ni n, (9.7)
xvv = Γ122 xu + Γ222 xv + hxvv · ni n;

i.e., we develop the three-dimensional vectors on the base (xu , xv , n) of the space.

Here appear — or naturally arise — the Γ symbols or “Christoffel symbols” which we discuss later
(Chapter 10). They illustrate the reality that differentiating a vector (also known as “parallel
transport”, see Chapter 10) in curvilinear co-ordinates on a curved surface is not a trivial thing.

The third term on the right hand side however represents, according to definition (9.4) , the
elements of the second fundamental form e, f, g.

We obtain

dt1
 
1 1 2 1 1 2 1 2 2



k = + Γ11 t + 2Γ12 t t + Γ22 t xu +
ds
 2 
dt
2
2
+ + Γ211 t1 + 2Γ212 t1 t2 + Γ222 t2 xv +
ds

2
2 
+ e t1 + 2f t1 t2 + g t2 n.

Here, the first two terms represent the internal curvature kint of the curve C, the curvature
inside surface S; the latter term represents the exterior curvature kext , the curvature of the
curve “along with” the itself curved surface.

The internal curvature again has two components in “(u, v) co-ordinates”, which can be read
from the above equation. We write

kint = k 1 xu + k 2 xv ,

6
Don’t get nervous about the use of a superscript (upper index). It’s just a writing convention, the good sense
of which we will argue later on

67
Chapter 9. The surface theory of Gauss

where
"
2
2 #
k1 dt1 /ds + Γ111 t1 + 2Γ112 t1 t2 + Γ122 t2
 
i
k = = =
k2

2
2
dt2 /ds + Γ211 t1 + 2Γ212 t1 t2 + Γ222 t2
 
2 X
X 2
1
 dt /ds + Γ1ij ti tj 
 
i=1 j=1
=  ,
 
 X2 X 2 
 dt2 /ds + Γ2ij ti tj 
i=1 j=1

where we have “economized” the formulas by using summation signs. The external curvature
again is
kext = hk · ni n,
where
2
2
hk · ni = e t1 + 2f t1 t2 + g t2 .

Example:
a latitude circle on the Earth surface.
 
N (ϕ) cos ϕ cos λ
x =  N (ϕ) cos ϕ sin

λ
;
2
N (ϕ) 1 − e sin ϕ
   
−N cos ϕ sin λ − sin λ
dx dx dλ  1
t= = = +N cos ϕ cos λ  =  + cos λ  ;
ds dλ ds N cos ϕ
0 0
 
− cos λ
dt dt dλ  1
k= = = − sin λ  .
ds dλ ds N cos ϕ
0

External curvature:  
cos ϕ cos λ
n =  cos ϕ sin λ  ,
sin ϕ
so
1 n
cos2 λ + sin2 λ n = − .


kext = hk · ni n = −
N N
Because the vector n is a unit vector, we may infer, that the exterior curvature of the curve is
precisely the inverse of the transverse radius of curvature, and directed inward. This is the same
as the curvature of the surface taken in the direction of the curve.
The internal curvature is
 
cos λ
− + cos ϕ cos λ  
cos ϕ  tan ϕ − cos λ sin ϕ
 
1 
kint = k − kext =  sin λ =  − sin λ sin ϕ  .
N − + cos ϕ sin λ  N
 cos ϕ  + cos ϕ
+ sin ϕ

We obtain for the length of this vector

tan ϕ
kkint k = .
N

68
 9.6. The geodesic

p/sin ϕ = N cot ϕ

kint

p
k ϕ

kext

N = p/cos ϕ

Figure 9.3.: The curvatures and radii of curvature of a latitude circle

1 1
In Figure 9.3 we see the curvature vector k, length k = = ; its internal part kint ,
N cos ϕ p (ϕ)
1 tan ϕ 1
length k sin ϕ = sin ϕ = ; and its external part kext , length k cos ϕ = cos ϕ =
p (ϕ) N (ϕ) p (ϕ)
1
. In the figure one also sees how the distance from the Earth’s rotation axis is in every
N (ϕ)
direction (k, kint ja kext ) the inverse of the curvature. This is also intuitively clear from rotational
symmetry.

9.6. The geodesic

Internally, in surface co-ordinates

We obtain the formula for the geodesic by requiring kint = k i = 0, i.e., the curve has no interior
curvature (the exterior curvature cannot be eliminated as the curve is on a curved surface):
2 2
dti X X i j k
+ Γjk t t = 0. (9.8)
ds
j=1 k=1

This approach is developed further later on in chapter 10.5.

Externally, using vectors in space

An alternative, three-dimensional (“exterior”) form is obtained by observing, that the geodesic
is only externally curved, i.e., by writing
    1 
dt 
1 2

1 2


2 2
 1 2  e f t
= kext = e t + 2f t t + g t n= t t n.
ds f g t2

Because

ht · xu i = t1 hxu · xu i + t2 hxu · xv i = Et1 + F t2 ,

ht · xv i = t1 hxv · xu i + t2 hxv · xv i = F t1 + Gt2 ,

69
Chapter 9. The surface theory of Gauss

it follows that −1 

t1
   
E F ht · xu i
= .
t2 F G ht · xv i
If we define    
t1 ht · xu i
≡
t2 ht · xv i
we may write
 −1  
 E F 
ht · xu i ht · xv i ·
dt  F G

=   n.
 
−1 
ds 
 
e f E F ht · xu i 
·
f g F G ht · xv i

Here  −1   −1

E F e f E F
= H −1 BH −1
F G f g F G
following the earlier used notation7 .
Example: on the ellipsoidal surface of the Earth, we have

M −2 M −2
   
0 M 0
H −1 BH −1 = =
0 N −2 cos−2 ϕ N cos2 ϕ 0 N −2 cos−2 ϕ
M −3 1/M 3
   
0 0
= =
0 N −3 cos−2 ϕ 0 cos ϕ/p3

and    
ht · xu i M ht · eN i
= ,
ht · xv i p ht · eE i
from which we obtain
  
dt  1/M 0  ht · eN i
= ht · eN i ht · eE i n=
ds 0 cos ϕ/p ht · eE i
 
1 1
= ht · eN i2 + ht · eE i2 n.
M N
Here, the unit vectors
   
− sin ϕ cos λ − cos ϕ sin λ
xu xv
eN = =  − sin ϕ sin λ  , eE = =  + cos ϕ cos λ  .
kxu k kxv k
cos ϕ 0
1
The approach is geometrically intuitive. The expression in wave brackets ht · eN i2 +
M
1 1 1
ht · eE i2 = cos2 A+ sin2 A is precisely the curvature of the surface in the direction
N M N
of the curve.
dx
In addition we have the equations = t, altogether 2 × 3 = 6 ordinary differential
ds
equations.
The method of space vectors has both advantages and disadvantages compared to the surface
co-ordinate method (e.g., equations (2.1)).
Advantage: in surface co-ordinates (ϕ, λ) there are inevitably two poles, singularities where the
curvature of the latitude circles goes to infinity and numerical methods can be ill-behaved.
This will not happen in rectangular space co-ordinates.
7
In index notation: g ij βjk g kl , see chapter 10. A logical notation for this would be β il .

70
 9.6. The geodesic

Disadvantages:

1. More equations, more computational effort.

Z
2. At every point, we must compute M and N and for this ϕ = arctan , where
p (1 − e2 ) p
p = X 2 + Y 2.
3. The “roll-in” and “roll-out” of the computation of the geodesic presupposes the transfor-
mation of (ϕ, λ) to (X, Y, Z) co-ordinates in the starting point, and back in the end point.

71
 Chapter 10
The surface theory of Riemann

An important theoretical frame in which curved surfaces and curves are often described, is
Riemann’s1 surface theory. In this theory, we study the curved surface intrinsically, i.e., without
taking into account that the curved surface of the Earth is “embedded” in a complete three-
dimensional (Euclidean) space.
This makes the surface theory of Riemann useful also in situations, where this embedding higher-
dimensional space doesn’t necessarily even exist. For example, in General Relativity we describe
spacetime (x, y, z, t) as a curved manifold according to Riemann’s theory, the curvature param-
eters of which relate to the densities of mass and current through Einstein2 ’s field equations.
The gravitational field is the expression of this curvature.

10.1. What is a tensor?

In physics, more than in geodesy, we often encounter the concept of “tensor”. What is a tensor?

Vectors

First of all, a vector. Initially we shall look at the situation in two-dimensional Euclidean space,
i.e., the plane.
A pair of co-ordinate differences between two adjacent points is a vector:
 
i ∆x
v=v = (10.1)
∆y

is a vector, written with a superscript.

In the other co-ordinate system we write this vector as

∆x0
 
i0
v=v = .
∆y 0

Note that here we have both a symbolic notation (v) and a component notation (v i , i = 1, 2).
The components depend on the chosen co-ordinate system (∆x0 6= ∆x, ∆y 0 6= ∆y, although

2
2
always ∆x2 + ∆y 2 = ∆x0 + ∆y 0 = constant — an invariant), but the symbolic notation
does not depend on this. A vector is always the same thing, even if its co-ordinates transform
with the co-ordinate system. A vector is always the same “arrow in space”.
There exists the following transformation formula for changing the components of a vector:
0
X 0
vi = αii v i , (10.2)
i
1
Georg Friedrich Bernhard Riemann (1826 – 1866), German mathematician.
2
Albert Einstein, 1879 – 1955, German-born theoretical physicist, discoverer of relativity theory, and icon of
science.

73
Chapter 10. The surface theory of Riemann

or in “matrix language”
v0 = Av,
where
v1
   
i ∆x
v=v = =
v2 ∆y
is the column matrix formed by the components, and
 1
α1 α21
  
0 cos θ sin θ
A = αii = =
α12 α22 − sin θ cos θ
is the component matrix of the transformation operator.
Every quantity that transforms according to equation (10.2) is called a vector. Familiar vector
quantities are velocity, acceleration, force, . . . what they have in common is, that they can be
graphically presented as arrows.

Tensors

A tensor is just a square object — a matrix — which has the same transformation property
from on co-ordinate system to another, but for every index :
0 0
X 0 j0
Ti j = αii αj T ij .
i,j
3
In “matrix language” this is
T 0 = AT AT ,
where A is the same as defined above.
We have already encountered many tensors in geodetic theory:
1. The inertial tensor of the Earth.
2. The gravity gradient or Marussi tensor
  ∂2W
∂2W ∂2W
 
∂g
 ∂x   ∂x2
∂x∂y ∂x∂y 
 ∂g   
  ∂2W
∂2W ∂2W 
M = = .

 ∂y   ∂y 2
∂y∂x ∂x∂z 
 ∂g   
 ∂2W ∂2W 
∂2W
∂z ∂z∂x
∂z∂y ∂z 2
" #
σx2 σxy
3. The variance “matrix” is really a tensor: Var (x) = , where σx and σy are the
σxy σy2
mean errors of the co-ordinates, i.e., the components of x, and σxy is the covariance.
 
E F
4. In the earlier discussed surface theory of Gauss, the first, H = , and the second,
F G
 
e f
B= fundamental form, as well as C = H −1 B;
f g
5. Also the elements of the map
" projection
# fundamental form (to be discussed later) E, F , G
e e e
Ee Fe
constitute a tensor He = , a metric tensor of sorts. The object H −1 H
e =
Fe Ge
 
Ee Fe
 M2 Mp 
  may again be called the scale tensor.
 Fe Ge 
Mp p2
3
Verify that the matrix equation produces the same index summations as the index equation!

74
 10.1. What is a tensor?

The geometric representation of a tensor and invariants

In the same way that the geometric depiction of a vector is an arrow, is the geometric depiction
of a tensor an ellipse (in two dimensions) or an ellipsoid (in three dimensions).The lengths of
the principal axes of the ellipse/ellipsoid depict the eigenvalues4 of the tensor; the directions of
the principal axes again are the directions of the tensor’s eigenvectors.
In Euclidean space and in rectangular co-ordinates, tensors T ij are typically symmetric; there-
fore, to different eigenvalues λi , λj , i 6= j belong eigenvectors xi , xj that are mutually perpendic-
ular, as proven in mathematics textbooks.
In an n-dimensional space, a tensor T has n independent invariants. The eigenvalues λi , i =
1, . . . , n are of course invariants. So are their sum and product.
X X
λi = T ii ,
i i

on the two-dimensional plane

λ1 + λ2 = T 11 + T 22 ,
the sum of the diagonal elements or trace 5 ; ja
Y
λi = det (T ) ,
i

on the two-dimensional plane again

λ1 λ2 = det (T ) = T 11 T 22 − T 12 T 21 ,

the determinant of the tensor.

The trace of the variance matrix σx2 +σy2 we already know as the point variance σP2 ; it was chosen
precisely because it is an invariant, independent of the direction of the xy axes. Also the trace
of the gravity gradient tensor is known:
X ∂2W ∂2W ∂2W
Mii = + + ≡ ∆W,
∂x2 ∂y 2 ∂z 2
i

the Laplace operator.

Tensors in general co-ordinates

And what about non-Euclidean spaces, with non-rectangular co-ordinate systems? Well, in that
case the difference between sub- and superscripts becomes meaningful, and the reason we write
superscripts may finally be appreciated.
A contravariant vector transforms as follows:
0 0
X
vi = αii v i (10.3)
i

and a covariant vector as follows: X

vi0 = αii0 vi . (10.4)
i

(These look very similar, but are not the same!)

4
More precisely, the lengths of the semi-axes are the square roots of the eigenvalues λi .
5
germ. Spur.

75
Chapter 10. The surface theory of Riemann

Where the “prototype” of a contravariant vector is the co-ordinate differences if two (nearby)
points, as above:  
i ∆x
v = ,
∆y
is the prototype of a covariant vector the gradient operator :
 
∂V ∂V /∂x
vj = = . (10.5)
∂xj ∂V /∂y

V (x, y) is some scalar field in space.

If we take for a model vector  
i dx
v = ,
dy
we obtain
dx0 ∂x0/∂x ∂x0/∂y
    
i0 dx
v = = ;
dy 0 ∂y 0/∂x ∂y 0/∂y dy
also with the chain rule

∂V /∂x0
    
∂x/∂x0 ∂x/∂y 0 ∂V /∂x
= .
∂V /∂y 0 ∂y/∂x0 ∂y/∂y 0 ∂V /∂y
0
∂xi i0 ∂xi
The coefficient matrices = α i and = αii0 are apparently each other’s inverse matrices.
∂xi ∂xi0
So, if the matrix form of the covariant transformation equation’s (10.4) transformation param-
eters αii0 is A (i row and i0 column index), then the matrix of the contravariant transformation
0
parameters (equation 10.3) αii must be A−1 (i0 row and i column index). From this circumstance,
the names “covariant” and “contravariant” derive.
A tensor may have both super- and subscripts. Under a co-ordinate transformation, each index
changes according to its “character”. There may even be more than two indices.

“Trivial” tensors

1. If we compute the gradient of a contravariant vector xi we obtain

∂xi
 
1 0
= ≡ δji ,
∂xj 0 1

the so-called Kronecker6 delta, in practice a unit matrix. It is a tensor:

0 0
X 0
δji 0 = αjj 0 αii δji = αii δji αjj 0
i,j

or in matrix language
0 0
i 0↓ ↓ i ↓
 i −1 j↓
I = A [I] [A] = A−1 A = I,
0
j → i → j → j0 →
0 0
because if the matrix form of αjj 0 is A, then that of αii – or correspondingly, that of αjj is
A−1 , as shown above.
6
Leopold Kronecker (1823 – 1891), German mathematician

76
 10.2. The metric tensor

2. The “corkscrew tensor” of Tullio Levi-Civita7 in three dimensions:


 0 if, of ijk, two are the same
ijk = 1 if ijk is an even permutation of the numbers (123)
−1 if ijk is an odd permutation of the numbers (123)


That is, 123 = 231 = 312 = 1, 132 = 321 = 213 = −1, all others = 0.
Ks. http://folk.uio.no/patricg/teaching/a112/levi-civita/.

10.2. The metric tensor

The metric tensor, or the metric, describes the Pythagoras theorem in a curved space. It is
identical to the earlier discussed Gaussian First Fundamental Form. In the ordinary plane (R2 )
we may choose rectangular co-ordinates (x, y), after which we may write the distance s between
to points 1, 2 as:
s2 = ∆x2 + ∆y 2 ,
where ∆x = x2 − x1 , ∆y = y2 − y1 are the co-ordinate differences between the points. This
equation applies throughout the plane. Its differential version looks the same:

ds2 = dx2 + dy 2 .

This is now written into the following form:

ds2 = gij dxi dxj , (10.6)

where    
1 0 i j dx
gij = ja dx = dx = .
0 1 dy
This is called the metric of the rectangular co-ordinate system in the Euclidean plane. In
equation (10.6) it is assumed that we sum over the indices i and j; we call this the Einstein
summation convention. Always when in an equation we have a subscript and a superscript with
the same name, we sum over it. I.e., in this case
2 X
X 2
2
ds = gij dxi dxj .
i=1 j=1

An alternative way of writing using matrix notation is a quadratic form:

ds2 = hdx · Hdxi = dxT Hdx =

H dx
dxT z }|   }|1 { .
{ z
= z }| { g11 g12 dx
1 2
dx dx
g21 g22 dx2

If the surface is not curved, one may always find a co-ordinate system which is everywhere
rectangular and both co-ordinates “scaled” correctly such that to a co-ordinate difference of 1
m corresponds also a location difference of 1 m. Then the gij matrix or metric tensor has the
form of the unit matrix, like above.
If the surface is curved, we can find a unit matrix only in some places. E.g., on the surface of
the Earth, only within a strip in the vicinity of the equator. It won’t be possible in a unified
manner over the whole surface of the Earth.

77
Chapter 10. The surface theory of Riemann

v + ∆v C B

v
q∆
v α p∆u
A
u u
∆
+
u
Figure 10.1.: Skewed metric in the plane

Nevertheless we an choose also on a non-curved surface a co-ordinate system that isn’t rectan-
gular but skewed, and where the co-ordinates are scaled arbitrarily. See figure 10.1.
In this case we find with the aid of the cosine rule:

s2 = p2 ∆u2 + q 2 ∆v 2 + 2p∆uq∆v cos α,

or, differentially
ds2 = p2 du2 + q 2 dv 2 + 2pq cos αdudv,
or, in index notation
ds2 = gij dxi dxj
where
p2
   
pq cos α i du
gij = , dx = .
pq cos α q2 dv

The matrix representation:

gij , ↓ i, → j dxj
dxi
2
z }| { z
  }| { .
ds = z }| {  p2 pq cos α du
du dv
pq cos α q2 dv

On a curved surface gij will depend on place, gij xi , where xi is a “vector” of parameters


T
describing the surface, xi = u v

. Also on a non-curved surface can gij depend on place,
e.g., when choosing curvilinear, e.g., polar, co-ordinates.
Examples:
1. The surface of a spherical Earth:

ds2 = R2 dϕ2 + R2 cos2 ϕdλ2 ,

i.e.,
R2
   
0 i dϕ
gij (ϕ, λ) = , dx = ,
0 R cos2 ϕ
2 dλ
from which in “matrix language”
gij , ↓ i, → j dxj
dxi
ds2 = z
z }| {
  }| { .
z
{  R2 0 dϕ
}|
dϕ dλ
0 R2 cos2 ϕ dλ
7
Tullio Levi-Civita (1873 – 1941), Italian mathematician.

78
 10.3. The inverse metric tensor

2. Polar co-ordinates (ρ, θ) in the plane:

ds2 = dρ2 + ρ2 dθ2 ,
i.e.,    
1 0 i dρ
gij = , dx = ,
0 ρ2 dθ
from which
gij , ↓ i, → j dxj
dxi
2
z }| {  z }| {
ds = z }| { 1 0 dρ .
dρ dθ 2
0 ρ dθ

3. In three dimensions in the air space using “aviation co-ordinates”: nautical miles
North dN , nautical miles East dE, feet above sea level dH; in metres
ds2 = (1852)2 dN 2 + (1852)2 dE 2 + (0.3048)2 dH 2 ,
eli
18522
   
0 0 dN
gij =  0 18522 0  dxi =  dE  ,
0 0 0.30482 dH
ja
gij , ↓ i, → j dxj
dxi
z }| {
 z }| {
18522

ds2 = z }| { 0 0 dN .
dN dE dH  0 18522 0   dE 
0 0 0.30482 dH

10.3. The inverse metric tensor

The inverse tensor of the metric tensor gij is written g ij . As a matrix, it is the inverse matrix
of gij :
g ij = (gij )−1 ,
or, in index notation
g ij gjk = δki .
This is nothing but the definition of the inverse matrix:
H −1 H = I,
where I is the unit matrix, the matrix representation of the Kronecker tensor δki . Written
out:
g ij , ↓ i, → j gjk , ↓ j, → k δki , ↓ i, → k
z }| { z }| {
 =  z }| {
g 11 g 12
 
g11 g12 1 0 .
21 22 g21 g22
g g 0 1

10.3.1. Raising or lowering sub- or superscripts of a tensor

The sub- or superscript of an arbitrary tensor can be “raised” or “lowered” by multiplying with
the metric tensor gij or inverse tensor g ij :
Tji = g ik Tkj = gjk T ki .

All forms Tij , T ij , Tji designate the same tensor, written in different ways. As a special case,
δji = gjk g ki = g ik gkj , i.e., the Kronecker delta tensor is a “mixed form” of the metric tensor
and could be written gji . The delta style of writing has established itself, however.

79
Chapter 10. The surface theory of Riemann

0
v 01 v
v 02

yi
2 2
u +y
v1 v
v2

u2

u1 u1 + y 1
Figure 10.2.: Christoffel symbols and parallel transport

10.3.2. The eigenvalues and -vectors of a tensor

The eigenvalue problem for a square tensor has the following form:

T ij − λg ij xj = 0,

(Tij − λgij ) xj = 0,
Tji − λδji xj = Tji − λδji xi = 0.

All three forms are equivalent, which is easily proved. For eigenvectors we have xi = gij xj . If
the tensor T ij (or Tji , or Tij ) is symmetric (meaning T ij = T ji etc.), then the eigenvalues λ are
real and the eigenvectors mutually orthogonal: if xi , y i are different eigenvectors, then

gij xi y j = 0.

There are as many eigenvalues as there are dimensions in the space, i.e., in the plane R2 two.

10.3.3. The graphic representation of a tensor

The quadratic form

Tij xi xj = 1
defines an ellipsoid (in R2 an ellipse) that may be considered the graphic of Tij . E.g., the inertial
ellipsoid, the variance ellipsoid.

10.4. The Christoffel symbols

The metric tensor isn’t yet the same thing as curvature. It isn’t even the same thing as the
curvature of the co-ordinate curves; to study this, we need the Christoffel 8 symbol s.
The Christoffel symbols describe what happens to a vector when it is transported parallelly along
the surface; more precisely, what happens to its components.
See figure 10.2. When a vector v, the components of which are v i , is transported parallelly from
i
one point to another over a distance y i , its components will change by amounts ∆v i = v 0 − v i ,
8
Elwin Bruno Christoffel (1829–1900), German mathematician-physicist

80
 10.4. The Christoffel symbols

although for the “vector itself” v0 = v. When both the transported vector and the distance of
transport are small, we may assume that the change depends linearly on both the vector and
the direction of transport, in the following way:

∆v i = Γijk v j y k .

The Christoffel symbols Γijk do not form a tensor like the metric tensor does. They describe the
“curvilinearity” of a curvilinear co-ordinate system, i.e., a property of the co-ordinate system. A
co-ordinate transformation that removes – at least locally if not everywhere – the curvature of
the co-ordinate curves, also makes the elements of Γijk vanish in that point (for tensors, such a
local “transforming away” will never succeed!).
E.g., on the Earth surface in the (ϕ, λ) co-ordinate system on the equator the co-ordinate
curves are locally non-curved and later we shall show, theat there indeed all Γijk vanish. An-
other example is from the general theory of relativity, where the components of acceleration are
Christoffel symbols: acceleration can be transformed away by moving to a “falling along”
reference system. In Einstein’s falling elevator the people inside the elevator accelerate with
respect the Earth surface but are weightless (i.e., the acceleration vanishes) in a reference system
connected to the elevator.
The Christoffel symbols can be computed from the metric tensor; the equation is (complicated
proof; see appendix C):  
1 ∂gk` ∂g`j ∂gjk
Γijk = g i` + − . (10.7)
2 ∂xj ∂xk ∂x`
(g ij ≡ (gij )−1 .) From this is can be seen, that the Christoffel symbols are, like, the first
derivatives of the metric. In a straight-lined co-ordinate system on a non-curved surface gij is a
constant and thus all Γijk vanish. Also, because gij is symmetric (gij = gji ), we obtain

Γijk = Γikj .

The Christoffel symbols are useful when writing the equation of the geodesic 9 in this formalism
(cf. equation (9.8)):
d i
t + Γijk tj tk = 0.
ds
Here, ti is the tangent vector of the geodesic ti = dxi/ds.
Examples:

1. On the surface of the spherical Earth (ϕ, λ). We use a notation where ϕ and λ symbolize
the index values 1 and 2 (for these symbols Einstein’s summation convention thus fails
to work!):
   
∂gk` 0 0 ∂gk` 0 0
= , = ,
∂ϕ 0 −R2 sin 2ϕ ∂λ 0 0
so only
∂gλλ ∂g22
= = −R2 sin 2ϕ
∂ϕ ∂ϕ
9
In fact we write
Dti dti
≡ + Γijk tj tk ,
ds ds
the so-called absolute or covariant derivative which is a tensor; and

Dti
=0
ds
is then the equation for the geodesic, which thus applies in curvilinear co-ordinates.

81
Chapter 10. The surface theory of Riemann

is non-vanishing. Then, remembering that

 2 −1  −2 
i` −1 R 0 R 0
g = (gi` ) = = ,
0 R2 cos2 ϕ 0 R−2 cos−2 ϕ

the only non-zero elements are:

 
ϕ 1 −1 ∂gλλ 1 1
Γλλ = (gϕϕ ) − == − R−2 · R2 sin 2ϕ = + sin 2ϕ = sin ϕ cos ϕ,
2 ∂ϕ 2 2
 
1 ∂gλλ 1 sin 2ϕ
Γλϕλ = (gλλ )−1 R−2 cos−2 ϕ −R2 sin 2ϕ = −



= = − tan ϕ,
2 ∂ϕ 2 2 cos2 ϕ
 
1 ∂gλλ
Γλλϕ = (gλλ )−1 = − tan ϕ.
2 ∂ϕ
Note that at the equator ϕ = 0 all Γijk = 0.
2. Polar co-ordinates in the plane (ρ, θ):
   
∂gk` 0 0 ∂gk` 0 0
= , = ,
∂ρ 0 2ρ ∂θ 0 0

yielding
∂g22 ∂gθθ
= = 2ρ,
∂ρ ∂ρ
so
 
1 −1 ∂gθθ 1
Γρθθ = (gρρ ) − = · −2ρ = −ρ,
2 ∂ρ 2
 
1 −1 ∂gθθ 1 1 1
Γθρθ = Γθθρ = (gθθ ) + = · · 2ρ = + .
2 ∂ρ 2 ρ2 ρ

10.5. The geodesic revisited

Here are alternative formulas for integrating the geodesic on the sphere (generalization to the
ellipsoid of revolution is complicated but possible):
dξ
+ η 2 sin ϕ cos ϕ = 0,
ds
dη
− 2ηξ tan ϕ = 0.
ds

Here we have already used the formulas derived above for the elements of Γijk .
Simultaneous integration of the differential equations would give
cos A
 
 
ξ R
=  sin A ,
 
η
R cos ϕ
as a function of s, from which A and ϕ may be computed.
To this set we add the definition equations of the tangent
dϕ
= ξ = t1 ,
ds
dλ
= η = t2 ,
ds

82
 10.6. The curvature tensor

C
D
u` + ∆u` 0i
∆u`
v v i

∆uk

A B
`
u

uk uk + ∆uk
Figure 10.3.: The curvature tensor and parallel transport around a closed grid path

In this way also λ comes along.

Perhaps this approach seems overly complicated; its major theoretical advantage is, that it will
work for all curvilinear co-ordinate systems on the surface considered, also, e.g., for stereographic
map projection co-ordinates used in the polar areas, as long as we first manage to write down
the metric tensor of the co-ordinate curves.
 T
The length of the tangent vector ξ η is computed as follows:
 2  2
2 i j 2 cos A 2 2 sin A
ds = gij t t = R + R cos ϕ = 1,
R R cos ϕ

if we remember that
R2
 
0
gij =
0 R cos2 ϕ
2

on the surface of the sphere, The lengths of the tangent vectors have always to be 1; this
requirement is fulfilled if s is a parametrization of the curve “by distance”.

10.6. The curvature tensor

The curvature is described again in a slightly more complicated way by means of parallel transport
around a closed path — a small rectangle10 .

i
As seen from figure 10.3 we obtain the change ∆v i ≡ v 0 − v i of the vector v i , which depends
on
1. the orientation of the closed rectangular path, the side indices k and `;
2. the size of the rectangle, measured in curvilinear co-ordinates: let the length of the uk side
be∆uk and the length of the u` side, ∆u` ;
3. the transported vector v i itself.
In the following way:
∆v i = Rjk`
i
v j ∆uk ∆u` . (10.8)
i
Here, the creature Rjk` is called the Riemann curvature tensor. In two-dimensional space (i.e.,
4
on a surface) it has 2 = 16 elements.
10
More precisely, a parallellogram

83
Chapter 10. The surface theory of Riemann

We can compute Riemann from the Christoffels11 :

i ∂ i ∂
Rjk` = k
Γj` − ` Γijk + Γikm Γm i m
j` − Γ`m Γjk . (10.9)
∂x ∂x
From this we spot immediately the following antisymmetry:
i i
Rjk` = −Rj`k .

There exist many more such symmetries; the number of independent components is actually
small.
Esimerkkejä:
1. Surface of a spherical Earth (ϕ, λ):
Based on the antisymmetry property we may say that only elements for which k 6= l can
differ from zero. Additionally the Christoffel symbols depend only on ϕ. Thus we
obtain the auxiliary terms:

 
∂ ϕ ∂ 1
Γ = sin 2ϕ = + cos 2ϕ,
∂ϕ λλ ∂ϕ 2
∂ λ ∂ λ ∂ 1
Γ = Γ = (− tan ϕ) = − 2 ,
∂ϕ ϕλ ∂ϕ λϕ ∂ϕ cos ϕ
the others zero.
The terms Γikm Γm
jl are obtained in the following way:

Γϕ λ λ ϕ ϕ λ 2
λλ Γϕλ = Γλϕ Γλλ = Γλλ Γλϕ = − sin ϕ,
Γλϕλ Γλλϕ = Γλϕλ Γλϕλ = + tan2 ϕ.

By combining

ϕ ϕ ∂ ϕ ∂ ϕ ϕ
Rϕϕλ = −Rϕλϕ = Γϕλ − Γϕϕ + Γϕ m m
ϕm Γϕλ − Γλm Γϕϕ = 0;
∂ϕ ∂λ
λ λ ∂ λ ∂ λ 1
Rϕϕλ = −Rϕλϕ = Γϕλ − Γϕϕ + Γλϕm Γm λ m
ϕλ − Γλm Γϕϕ = − + tan2 ϕ = −1;
∂ϕ ∂λ cos2 ϕ
ϕ ϕ ∂ ϕ ∂ ϕ ϕ
Rλϕλ = −Rλλϕ = Γλλ − Γ + Γϕ m m 2 2
ϕm Γλλ − Γλm Γλϕ = cos 2ϕ + sin ϕ = cos ϕ;
∂ϕ ∂λ λϕ
λ λ ∂ λ ∂ λ
Rλϕλ = −Rλλϕ = Γλλ − Γ + Γλϕm Γm λ m
λλ − Γλm Γλϕ = 0.
∂ϕ ∂λ λϕ

2. Polar co-ordinates (ρ, θ) plane:

∂ ρ
Γ = −1,
∂ρ θθ
∂ θ ∂ θ 1
Γ = Γθρ = − 2 ,
∂ρ ρθ ∂ρ ρ
the others again vanish. Then:

Γρθθ Γθρθ = Γθθρ Γρθθ = Γρθθ Γθθρ = −1,

1
Γθρθ Γθθρ = Γθρθ Γθρθ = + 2 .
ρ
11
The derivation is found in appendix D.

84
 10.6. The curvature tensor

After this:
ρ θ
Rρρθ = −Rρθρ = 0;
θ θ ∂ θ ∂ θ 1 1
Rρρθ = −Rρθρ = Γ − Γ + Γθρm Γm θ m
ρθ − Γθm Γρρ = − 2 + 2 = 0;
∂ρ ρθ ∂θ ρρ ρ ρ
ρ ρ ∂ ρ ∂ ρ ρ
Rθρθ = −Rθθρ = Γ − Γθρ + Γρρm Γm m
θθ − Γθm Γθρ = −1 + 1 = 0;
∂ρ θθ ∂θ
θ θ
Rθρθ = −Rθθρ = 0.
So the whole Riemann tensor vanishes, as it should, because the surface is not curved.
From the Riemann tensor we obtain the smaller Ricci12 -tensor in the following way:
2
X
i i
Rjk = Rjik = Rjik . (10.10)
i=1

This is a symmetric tensor, Rij = Rji .

Examples
1. Let us continue with the computation of Rij for the case of a spherical surface:
ϕ λ
Rϕϕ = Rϕϕϕ + Rϕλϕ = +1;
ϕ λ
Rϕλ = Rϕϕλ + Rϕλλ = 0 + 0 = 0;
λ ϕ
Rλλ = Rλλλ + Rλϕλ = cos2 ϕ.

2. For polar co-ordinates in the plane, all Rij = 0.

We may continue this process to obtain the curvature scalar 13 :
2
(gij )−1 Rji .
X
R = g ij Rji = (10.11)
j=1

For the example case:

1. Spherical surface:
1
 
0    −2 
 R2  1 0 R 0
Rki ij
= g Rjk = 1 = ,
0 + cos2 ϕ 0 R−2

0 2 2
R cos ϕ
i.e.14
X 1 cos2 ϕ 2
R= Rii = 2
+ 2 2
= 2.
R R cos ϕ R
i

2. In polar co-ordinates Rij = 0 , i.e., R = 0.

More generally we have (without proof):
2
R = 2K = 2κ1 κ2 = , (10.12)
R1 R2
twice the Gauss total curvature K, the inverse of the product of the two principal radii of
curvature.
12
Gregorio Ricci-Curbastro (1853 – 1925), Italian mathematician, inventor of tensor calculus
13
This is in fact the “trace” of Rij , more precisely, of Rji ≡ g ik Rkj , written Rii , see above. This also serves as an
example of how in general curvilinear co-ordinates an index may be “raised” from covariant to contravariant,
or “lowered”, using the metric tensor gij or its inverse g ij . In rectangular co-ordinates these are unit matrices
and the difference between super- and subscripts is without consequence.
14
Here we use the symbol R both for the radius of the Earth and for the curvature scalar. We hope it doesn’t
cause confusion.

85
Chapter 10. The surface theory of Riemann

du2
K

du1
θ = Kdu1 du2
Figure 10.4.: The curvature of a two-dimensional surface may be characterized by just one pa-
rameter: θ.

10.7. Gauss curvature and spherical excess

Contrary to the eigenvalues of the Gauss second fundamental form βki = g ij βjk (chapter 9.4)
are the eigenvalues of the tensor Rji unrelated to the principal curvatures of the surface15 κ1 and
i
κ2 . In the two-dimensional case both the Rij tensor and the Rjkl have only one essentially
independent element, which is related to the Gauss total curvature K = κ1 κ2 .
This is not hard to show:
1. in the equation (10.8) ∆u∆v describe the sides of a small diamond shape, around which the
vector v i is transported in a parallel fashion. In two dimensions there are only two choices
for the sides in the uk , u` directions: ∆uk ∆u` and ∆u` ∆uk . One represents clockwise
i i
transport, the other, counterclockwise. The corresponding elements Rjk` = −Rj`k thus
are essentially the same.
2. In the same equation v j and ∆v i are mutually perpendicular. ∆v i represents a small
rotation of vector v j . On a surface, in two dimensions, the rotation is described by one
angle. Because the angle between two parallelly transported vectors doesn’t change, we
may infer that this rotation angle is the same for all vectors v i . Again we find only one
independent parameter. See figure 10.4.
In other words, when give a ∆u1 ∆u2 diamond shape, the expression Rj12 i
∆u1 ∆u2 =
 
i cos θ sin θ
−Rj21 ∆u2 ∆u1 is a 2 × 2 rotation matrix , containing only one free pa-
− sin θ cos θ
rameter.
However, for higher dimensionalities, the number of independent elements in the tensors of
Riemann and Ricci is larger. This case is interesting because of General Relativity (four
dimensions!) though not for geodesy.
We can remark here, that the spherical excess already discussed in section 1.2 is a special case
of the change in direction of a vector that is parallelly transported around: the spherical excess
is the small change of direction of a vector transported around a closed triangle!
When we transport a vector around a larger surface area, it is the same as if we had transported
it successively around all of the little patches making up the surface area. This is depicted in
figure 10.5. In every little patch du1 du2 , area dS, the change in direction of the vector amounts
to KdS, where K is the Gauss total curvature at the patch location.
By generalization one obtains from this the following integral equation (Gauß-Bonnet16 for
the triangle, Theorema elegantissimum):
¨
ε= K (ϕ, λ) dS,
∆
15
The radii of curvature exist only in the three-dimensional space surrounding the Earth surface, in which it is
embedded. R on the other hand is an intrinsic property of the surface.
16
Pierre Ossian Bonnet (1819 – 1892), French mathematician

86
 10.7. Gauss curvature and spherical excess

3
D
4

ABCD

2
1
B

A
Figure 10.5.: Transport of a vector around a larger surface area. θABCD = θ1 + θ2 + θ3 + θ4 .

where K is a function of place.

On the reference ellipsoid K = (M N )−1 and
¨ ¨
1 1
ε= dS = M N cos ϕdϕdλ =
MN MN
∆ ∆
¨
= dσ = σ∆ ,
∆

which is the surface area of the corresponding triangle, but on the unit sphere.17 .
From this follows that
on the reference ellipsoid the spherical excess depends only on the directions of the
normals in the corner points (ϕi , λi ) , i = 1, 2, 3, not on the shape of the ellipsoidal
surface. Shorter: spherical excesses are computed on the sphere, even while all other
computations are done on the ellipsoid.
It suffices that the geographical co-ordinates of the corner points, which describe the ellipsoidal
normal, are known
A fast geometric way of computing the spherical excess of a triangle is the following:
Let the corner points  
cos ϕi cos λi
xi =  cos ϕi sin λi  , i = 1, 2, 3;
sin ϕi
we use the polarization method (1.7) in three dimensions, in order to find the “poles” of the
triangle’s sides:

hx2 × x3 i hx3 × x1 i hx1 × x2 i

y1 = , y2 = , y3 = .
kx2 × x3 k kx3 × x1 k kx1 × x2 k
17
This is not quite exactly true. . . the triangle symbols under the integral sign dϕdλ and under the integral sign
dσ are not exactly corresponding. On the reference elipsoid the triangle consists of geodesics, on the unit
sphere, of great circle segments; however, a geodesic on the ellipsoid does not map to a great circle!
This may be suspected already from the observation, that a long geodesic around the ellipsoid does not
generally close upon itself, while a great circle around a sphere always does.
Furthermore, the angles of an ellipsoidal and a spherical triangle are not individually equal; their sums (and
with that, their spherical excesses) however are.

87
Chapter 10. The surface theory of Riemann

Now the inter-pole distances correspond to the angles of the spherical triangle (more precisely,
π minus those angles):

ky2 − y3 k
α1 = π − 2 arctan ,
ky2 + y3 k
ky1 − y3 k
α2 = π − 2 arctan ,
ky1 + y3 k
ky − y2 k
α3 = π − 2 arctan 1 ,
ky1 + y2 k

numerically strong equations18 . Here αi is the angle at corner point i. The spherical excess is
now
X3
ε= αi − π.
i=1

10.8. The curvature in quasi-Euclidean geometry

We mention without proof that in a Riemann space we may transform curvilinear co-ordinated
always in such a way, that in a certain point P
1. the metric tensor is the unit matrix,

1 if i = j
gij = g ij =
0 otherwise,

2. the metric tensor is (locally) stationary, i.e.


∂gij 
= 0.
∂xk xk =xk
P

In this case we speak of a quasi-Euclidean neighbourhood around P .

In this case the Christoffel symbols (equation 10.7) are
 
i 1 i` ∂gk` ∂g`j ∂gjk
Γjk = g + − =
2 ∂xj ∂xk ∂x`
 
1 ∂gki ∂gij ∂gjk
= j
+ k
− ≈0
2 ∂x ∂x ∂xi

based on the above assumption of stationarity.

Of the Riemann curvature tensor (equation 10.9) the last two terms vanish:

i ∂ i ∂
Rjk` = k
Γj` − ` Γijk =
∂x  ∂x 
1 ∂ ∂g`i ∂gij ∂gj`
= + − −
2 ∂xk ∂xj ∂x` ∂xi
 
1 ∂ ∂gki ∂gij ∂gjk
− + − =
2 ∂x` ∂xj ∂xk ∂xi
   
1 ∂ ∂g`i ∂gki 1 ∂ ∂gj` ∂gjk
= − − − .
2 ∂xj ∂xk ∂x` 2 ∂xi ∂xk ∂x`
18
In some programming languages, we have for arctan (x/y) the symmetrical alternative form atan2 (x, y), where
the zerodivide problem does not occur.

88
 10.8. The curvature in quasi-Euclidean geometry

If we now derive the Ricci tensor (equation 10.10):

i 1 X ∂ 2 gj`
Rj` = Rji` =− +
2 (∂xi )2
i
∂ 2 g`i ∂ 2 gii ∂ 2 gji
 
1X
+ − ` j + ` i =
2 ∂xi ∂xj ∂x ∂x ∂x ∂x
i
1 X ∂ 2 g`i ∂ 2 gij ∂ 2 gii
 
1
= − ∆gj` + + i `− ` j .
2 2 ∂xi ∂xj ∂x ∂x ∂x ∂x
i

This already looks a lot more symmetric. The symbol ∆ here means the Laplace operator
X ∂2
∆= .
i (∂xi )2

Finally we still derive the curvature scalar

X
R = g jk Rkj = Rjj =
j
1X X X ∂ 2 gij 1 X X ∂ 2 gii
= − ∆gjj + − =
2
j i
∂xi ∂xj
j
2 (∂xj )2 i j
1X X X ∂ 2 gij 1X
= − ∆gjj + − ∆gii =
2 ∂xi ∂xj 2
j i j i
X XX ∂2g ij
= − ∆gii + .
∂xi ∂xj
i i j

In the special case that the co-ordinate curves are (everywhere, not just in point P ) orthogonal,
we obtain gij = 0 if i 6= j with its derivatives of place, i.e.
X X ∂ 2 gii
R=− ∆gii + =
(∂x i )2
i i
∂ 2 g11 ∂ 2 g22
 2
∂ 2 g11 ∂ 2 g22 ∂ 2 g22

∂ g11
= + − + + + =
(∂x1 )2 (∂x2 )2 (∂x1 )2 (∂x2 )2 (∂x1 )2 (∂x2 )2
 2
∂ 2 g22

∂ g11
=− + . (10.13)
(∂x2 )2 (∂x1 )2

This equation will be of use in the sequel.

89
 Chapter 11
Map projections

Map projections are needed because the depiction of the curved Earth surface on a plane is
not possible without error at least for larger areas. In this chapter we discuss the deformations
introduced by a map projection in terms of its scale error, using the tools developed in the
previous chapters.

11.1. Map projections and scale

11.1.1. On the Earth surface

On the surface of the Earth a distance element dS may consist of an element of latitude dϕ
and an element of longitude dλ. These correspond to linear distances M (ϕ) dϕ and p (ϕ) dλ,
respectively.
According to Pythagoras, the length of the diagonal of the postage stamp (dϕ, dλ) is

dS 2 = M 2 dϕ2 + p2 dλ2 . (11.1)

Now dS 2 defines a metric,

X
dS 2 = gij dxi dxj = xT Hx,
i,j

where
M2 0
   
E F
gij = = =H (11.2)
0 p2 F G

and
dx1
   
dϕ
x= = .
dx2 dλ

Here we have suitably defined two ways of writing, the index notation gij , dxi and the matrix-
vector notation H, x. The elements of the matrix are also the elements of the Gauss First
Fundamental Form on the Earth surface E, F, G. We see that E = M 2 , F = 0 ja G = p2 .
For the azimuth A again we obtain the following equation:

pdλ
tan A = , (11.3)
M dϕ

from which, with the above,

dS sin A = pdλ, (11.4)

dS cos A = M dϕ. (11.5)

91
Chapter 11. Map projections

11.1.2. In the map plane

When we project this little rectangle into the map plane, we obtain the sides dx and dy, and
according to Pythagoras the diagonal is

ds2 = dx2 + dy 2 .

Now we may calculate an element of distance in the map plane as follows:

∂x ∂x
dx = dϕ + dλ,
∂ϕ ∂λ
∂y ∂y
dy = dϕ + dλ,
∂ϕ ∂λ
i.e.,
 2  2
2 ∂x ∂x ∂y ∂y
ds = dϕ + dλ + dϕ + dλ =
∂ϕ ∂λ ∂ϕ ∂λ
e 2 + 2Fedϕdλ + Gdλ
= Edϕ e 2, (11.6)

where

∂x 2
  2
∂y
E
e = + ,
∂ϕ ∂ϕ
∂x ∂x ∂y ∂y
Fe = + ,
∂ϕ ∂λ ∂ϕ ∂λ
 2  2
∂x ∂y
G =
e + .
∂λ ∂λ

E,
e Fe ja Ge are the Gauss First Fundamental Form in the map plane. If we interpret the element
of distance in the map plane ds2 as a metric of the Earth surface, then we obtain also here a
metric tensor, " #
Ee Fe
gf
ij = ≡ H.
e (11.7)
Fe Ge

The corresponding metric is X

ds2 = gf i j T e
ij dx dx = x Hx,
i,j
T
where dx1 = dϕ and dx2 = dλ, , i.e., dxi = dxj = x = dϕ dλ . Also here we see as


alternatives the index notation and the matrix-vector notation, which describe the same thing.

11.1.3. The scale

The scale is now the ratio

ds
, m=
dS
which apparently depends on the direction of the distance element, i.e., the azimuth A.
Let us write the eigenvalue problem:
i j
P
2 ds2 i,j g
fij dx dx
m = = ⇒
dS 2 i j
P
i,j gij dx dx
X X
i j 2
⇒ ij dx dx − m
gf gij dxi dxj = 0. (11.8)
i,j i,j

92
 11.1. Map projections and scale

T  T
The same equation in “matrix language”, if x = dx1 dx2

= dϕ dλ :
 
xT H e − m2 H x = 0.

Equations (11.4, 11.5) The read:

dS sin A = pdλ,
dS cos A = M dϕ.
Let us now look at all vectors that are of form
cos A
 
 
dϕ  M 
x= =  sin A . (11.9)
dλ
p
The lengths of these vectors on the Earth surface are:
dS 2 = M 2 dϕ2 + p2 dλ2 =
 2  2
2 cos A 2 sin A
= M +p =
M p
= cos2 A + sin2 A = 1.
So: on the surface of the Earth these vectors form a circle of unit radius.
Substitute (11.9) into equation (11.8), using (11.6):
cos A 2 e sin A 2
   
cos A sin A
E
e +G + 2Fe −
M p M p
"  2  2 #
2 2 cos A 2 sin A
−m M +p =0
M p

or, after cleaning up,

Ee Ge Fe
m2 = cos 2
A + sin2
A + 2 sin A cos A = (11.10)
M2 p2 Mp
! !
1 E
e Ge 1 E
e Ge Fe
= 2
+ 2 + 2
− 2 cos 2A + sin 2A.
2 M p 2 M p Mp

From this we obtain the stationary values:

!
d 2 Ge Ee Fe
0= m = 2
− 2 sin 2A + cos 2A,
dA p M Mp

in other words  
E G
−
e e
M2 p2 e 2 − GM
Ep e 2
tan 2A = = .
Fe FeM p
Mp
This yields two maximum and two minumum values, which all four are at a distance of 90◦ from
each other1 .
These eigenvalues are obtained by writing the equation (11.10) as follows:
 
E
e
2 Fe
  M2 − m
 
 M p  cos A
cos A sin A  
 sin A = 0;

Fe Ge
2
−m
Mp p2
1
If F = 0, we obtain the condition sin 2A = 0, which is fulfilled when A = k · 90◦ , k = 0, 1, 2, 3.

93
Chapter 11. Map projections

This presupposes that the determinant of the matrix in the middle vanishes:
 
0 = det H −1 H e − m2 I =
! !
Ee
2 Ge
2 Fe2
= − m − m − =
M2 p2 M 2 p2
!
Ee Ge 1  e e e2 
= m4 + − 2 − 2 m2 + 2 2 E G−F .
M p M p

From this r
  2  
E G E G
+ ± + − 4 M 12 p2 E e − Fe2
eG
e e e e
M2 p2 M2 p2
m21,2 = .
2
These two solutions are called the principal scale factors m1 , m2 .
If the H
e matrix is a diagonal matrix:
" #
E
e 0
H
e = ,
0 G
e

we obtain
! !
  Ee Ge
det H −1 H
e − m2 I = − m2 − m2 =0 ⇒
M2 p2
Ee G e
m21,2 = , .
M 2 p2
s s
E
e Ge
m1 = is the meridional scale factor, m 2 = is the scale factor in the direction of the
M2 p2
parallel. In intermediate directions A (azimuth) the magnification factor is then
q
m = m21 cos2 A + m22 sin2 A.

11.1.4. The Tissot-indicatrix

The matrix (“scale tensor”)

 
Ee Fe
2
e = M
H −1 H Mp 
 = g ik gekj = geji

 Fe G 
e
Mp p2

is often visualized as an ellipse on the Earth surface (or, correspondingly, in the map plane).
The eigenvalues of the matrix are m21 and m22 ; the axes of the visualizing ellipse are m1 in the
meridional direction and m2 in the direction of the parallel. This ellipse is called the indicatrix
of Tissot2 . See subsection 10.3.3.
Of the many map projections used, we need to mention especially those that are conformal, i.e.,
the scale in a point is the same in all directions:

m1 = m2 = m,
2
Nicolas Auguste Tissot (1824-1897), French cartographer

94
 11.2. The Lambert projection (LCC)

M dϕ = 1 m1
λ
y

pdλ = 1
m2

Figure 11.1.: The Tissot indicatrix

ρ ρ
0

dt
ds dt
ds
ϕ0
θ

Figure 11.2.: Lambert projection

The Tissot indicatrix is a circle. With a conformal projection, small circles, squares and local
angles and length ratios are mapped true. Surface areas are distorted, however.
Conformal projections have the useful property, that the projection formulas in one co-ordinate
direction can be derived when the formulas for the other co-ordinate direction are given.
3
A sometimes useful property is, that surface areas are mapped true , even
 though the shapes
−1 e
of small circles or squares are distorted. This requirement is det H H = m21 m22 = constant.

11.2. The Lambert projection (LCC)

The Lambert4 projection (LCC, Lambert Conformal Conical) is a conformal conical projec-
tion. The scale for all latitude circles is constant; it is however different for different latitudes,
and attains its maximum value for a certain latitude ϕ0 in the middle of the area mapped.
Generally this value is m > 1; in that case we have two standard parallels for which the scale
m = 1.
3
E.g. when depicting the surface density of population or some other phenomenon.
4
Johann Heinrich Lambert (1728 – 1777) Swiss methematician, physicist and astronomer.

95
Chapter 11. Map projections

75
N o
60
N o
45
N o
30
N o
15
N o
0 60
W o
o

E
90 o
45 o
W
o E
30 o 75
W
oE
15 o
W 60
o
0o o o 45 E
15 E 30 E

Figure 11.3.: Example of a Lambert projektion (software used: m_map (http://www2.ocgy.

ubc.ca/~rich/map.html))

The images of the meridians in the map are straight lines, which intersect in one point, which
is also the common centre of all latitude circle images.
Polar co-ordinates in the map plane are

x = ρ sin θ,
y = ρ0 − ρ cos θ,

where ρ = ρ0 describes the latitude circle ϕ = ϕ0 .

Also
θ = nλ. (11.11)

Let ds be the distance element on the Earth surface and dt the corresponding element in the
map plane; then along the latitude circle

ds = p (ϕ) dλ
dt = nρdλ

and by dividing
dt nρ
= .
ds p
Based on conformality this will also be valid along the meridian.
Now we have ds = M (ϕ) dϕ, i.e.
dt dt ds nρ
= = M.
dϕ ds dϕ p

If we reckon ρ positive toward the South, i.e., dρ = −dt, we obtain:

dρ M (ϕ)
= −nρ (ϕ) . (11.12)
dϕ p (ϕ)

Using equations (11.11, 11.12) we may compute (θ, ρ) if given (ϕ, λ). However, let us first move
ρ to the left hand side:
ˆ ϕ
M (ϕ0 ) 0
 
d M
ln ρ = −n ⇒ ρ = ρ0 exp −n 0
dϕ .
dϕ p ϕ0 p (ϕ )

If we define ˆ ϕ
M (ϕ0 ) 0
ψ (ϕ) ≡ dϕ ,
0 p (ϕ0 )

96
 11.3. On the isometric latitude

the so-called isometric latitude, we obtain

exp {−nψ (ϕ)}

ρ = ρ0 exp {−n (ψ (ϕ) − ψ (ϕ0 ))} = ρ0 . (11.13)
exp {−nψ (ϕ0 )}

Furthermore we must choose a value n. We do this so, that the scale is stationary at the reference
latitude ϕ0 :
dm
= 0 jos ϕ = ϕ0 .
dϕ

The scale is
dρ dt nρ
m≡ =− =− ,
ds ds p

i.e.,
dm n dρ n dp
= −ρ 2 =
dϕ p dϕ p dϕ
n2 ρ nρ nρ
= − 2 M + 2 M sin ϕ = − 2 M (n − sin ϕ) ,
p p p

(using 11.12 and Appendix B) which ought to vanish for ϕ = ϕ0 . It does by choosing

n = sin ϕ0 .

As an initial condition we must still choose ρ0 ; it yields for the scale at the reference latitude
ρ0
m (ϕ0 ) = n .
p (ϕ0 )

Alternatively one may choose a value ϕ1 where m (ϕ1 ) = 1 (a standard latitude): then
ρ1
n = 1,
p (ϕ1 )

from which ρ1 ≡ ρ (ϕ1 ) follows. Now ρ (ϕ) is obtained by integrating (11.12) from a starting
value, either ρ0 or ρ1 , with the integration interval being either [ϕ0 , ϕ] or [ϕ1 , ϕ].
The inverse operation is easy for θ → λ; solving in reverse the differential equation (11.12)
(computing ϕ when ρ is given) can be done as follows:

1. invert analytically equation (11.13):

ψ (ϕ0 )
ψ (ϕ) = − (ln ρ − ln ρ0 ) ;
n

2. perform the inverse computation ψ → ϕ (see below).

11.3. On the isometric latitude

The isometric latitude, ˆ ϕ

M (ϕ)
ψ= dϕ,
0 p (ϕ)
can be computed numerically (quadrature; the QUAD routines of Matlab). However, for the
spherical case a closed solution exists
π ϕ
ψ = ln tan + , (11.14)
4 2

97
Chapter 11. Map projections

easily proven by differentiating using the chain rule.

Even for the ellipsoid of revolution a closed solution exists:
 π ϕ   1 + e sin ϕ  2e
" #
ψ = ln tan + . (11.15)
4 2 1 − e sin ϕ

See appendix A.
Reverse operation: if ψ is given and ϕ to be computed, we may take the equation
dψ M (ϕ)
=
dϕ p (ϕ)
and turn it upside down:
dϕ p (ϕ)
= .
dψ M (ϕ)
This equation has the general form
dy
= f (y, t)
dt
and may be numerically solved using Matlab’s ODE routines.
An alternative way for the sphere is to analytically invert equation (11.14):
 π
ϕ = 2 arctan exp ψ − .
4
This is also useful as the first iteration step for the ellipsoidal case:
 π
ϕ(0) = 2 arctan exp ψ − ,
4
after which   !e  
 1 − e sin ϕ(i) 2 π 
ϕ(i+1) = 2 arctan  exp ψ − .
 1 + e sin ϕ(i) 4

This converges rapidly.

11.4. The Mercator projection

The classical Mercator projection is obtained as a limiting case of Lambert by choosing the
limit n → 0, ρ0 → ∞, but nevertheless nρ0 = 1, and also ϕ0 = 0. Then
ρ = ρ0 exp {−n (ψ (ϕ) − ψ (ϕ0 ))} ≈
≈ ρ0 − nρ0 (ψ (ϕ) − ψ (ϕ0 )) =
= ρ0 − ψ (ϕ) .
Let us choose y ≡ − (ρ − ρ0 ) and x ≡ λ and we have the projection formulas of Mercator
x = λ,
ˆ ϕ
M (ϕ0 ) 0
y = ψ (ϕ) = dϕ .
0 p (ϕ0 )

Here we see the isometric latitude in its most simple glory: in the case of spherical geometry
π ϕ
y = ln tan + .
4 2
Mercator is not a “lamp projection”: There is no projection centre from which the “light”
emanates.

98
 11.5. The stereographic projection

y
dt
ds dt
ds
x
λ ϕ

Figure 11.4.: The Mercator projection

o
75 N

60oN

o
45 N

o
30 N

15oN

o
0

o o o o o o o o o o o
60 W 45 W 30 W 15 W 0 15 E 30 E 45 E 60 E 75 E 90 E

Figure 11.5.: An example of the Mercator projection

11.5. The stereographic projection

This so-called azimuthal projection is also conformal and also a limiting case of the Lambert
projection.
Let us choose in the equation (11.13) n = 1 and choose ϕ0 = 0 (the equator) and ρ0 correspond-
ingly, i.e., ρ0 ≡ ρ (ϕ0 ). In that case

ρ = ρ0 exp {−ψ (ϕ)} .

π
Unfortunately in the limit ϕ → the function ψ (ϕ) diverges; we may define
2
π
ρ = ρ0 exp {−ψ (ϕ)} if ϕ < ,
2
π
ρ=0 if ϕ = ,
2
after which ρ (ϕ) is continuous.

99
Chapter 11. Map projections

ρ ρ0
Projection plane

ρ0
π/2 −ϕ
Equator ϕ

( π2 −ϕ)/2

Projection centre
Figure 11.6.: The stereographic projection

In the spherical case we obtain

ρ = ρ0/tan( π + ϕ )
=
4 π 2 ϕ 
= ρ0 cot + =
 4π 2ϕ 
= ρ0 tan − =
 4  2 
1 π
= ρ0 tan −ϕ .
2 2

This case is depicted in figure 11.6. ρ0 is the distance from the central point (“South pole”) to the
projection plane. In this case the projection is truly a “lamp projection”. . . but unfortunately
only in the case of spherical geometry.
Let us derive the scale of the stereographical projection:

dρ dρ dψ dϕ M 1 ρ exp {−ψ (ϕ)}

= = −ρ0 exp {−ψ} · · = − = −ρ0 .
dx dψ dϕ dx p M p p (ϕ)

The value is negative because x grows to the North and ρ to the South. The co-ordinate x is
the metric “northing”.
By using equation 11.15 we obtain
 − e
dρ ρ0 π ϕ 1 + e sin ϕ 2
=− cot + ,
dx N cos (ϕ) 4 2 1 − e sin ϕ

or close to the pole

 − e
dρ ρ0 1+e 2
≈− .
dx 2 1−e
In case of the Earth where (GRS80) e = 0.08181919104281097693, we obtain for the last term
in e 0.99331307907268199009. If we want unity for the scale at the pole, we must set

ρ0 = 2.01346387371353381594.

100
 11.6. The Gauss-Krüger projection

o
oW 180 W 165 o
165 E
oW 15 o
0 0
15 E
o
30 N
o W

13
5

5
o
13

E
o
45 N

12
12 o

0
0

oE
o
60 N

105
105 o

oE
o
75 N

90 E
90 W

o
o

75 oE
75 W
o

E
60

60 o
oW

E
45

45 o
W o

30 o oE
W 30
15 oW o
0
o 15 E

Figure 11.7.: An example of the stereographic projection

Reference Latitude Number Computation Number

latitude interval of support points point (ϕ, ∆λ) of terms
1◦ 1◦ − 20◦ 20 (19◦ .333, 10◦ .0) 12
61◦ 61◦ − 80◦ 20 (79◦ .333, 20◦ .0) 16
70◦ 61◦ − 80◦ 20 (79◦ .333, 20◦ .0) 15
70◦ 61◦ − 80◦ 39 (79◦ .333, 20◦ .0) 15
65◦ 61◦ − 70◦ 19 (60◦ .333, 10◦ .0) 14
Table 11.1.: Results of Gauss-Krüger test computations. “Number of terms”: the number of
terms in the polynomial expansion needed to achieve a change under ±1 mm in the
computed Gauss-Krüger co-ordinates

11.6. The Gauss-Krüger projection

A good practical example of a map projection is the Gauss-Krüger projection, in use in

Finland. It is a transversal cylindrical projection which is conformal.
Like has been traditionally the habit, we will present projection formulas in the form of series
expansions. But, differently from tradition, we will determine the coefficients of the expansion
numerically.
Let us proceed in the following way. First we choose a suitable starting projection, e.g., an
ordinary Mercator, the equations of which are simple. So, we map the surface of the Earth
ellipsoid onto the map plane of an ordinary Mercator:

v = λ − λ0 ,
ˆ ϕ
M (ϕ0 ) 0
u = dϕ .
0 p (ϕ0 )
Next we construct, in the Mercator plane, an analytical mapping

u + iv → x + iy,

one property of which is, that x is of true length along the central meridian λ − λ0 = y = v = 0.

dx = M (ϕ) dϕ,

101
Chapter 11. Map projections

i.e., ˆ ϕ
M ϕ0 dϕ0 .


x= (11.16)
0
Furthermore we have on the central meridian:

y = 0.

Now we have a boundary value problem: the sought for complex map co-ordinate z ≡ x + iy is
given as a function of the starting projection’s map co-ordinate w ≡ u+iv on the edge y = v = 0
i.e., the real axis. The problem is to determine the function in the whole complex plane. See
figure 11.8.
Intermezzo. In complex analysis we talk of analytical functions. An analytical function is such
a mapping
z = f (w)
which is differentiable. Not just once, but infinitely many times. In this case, the Cauchy-
Riemann conditions apply:
∂x ∂y ∂y ∂x
= , =− ,
∂u ∂v ∂u ∂v
if z = x + iy, w = u + iv.
 T  T
This means that the small vector du dv is mapped to a vector dx dy according
to the following equation:
       
dx a −b du cos θ − sin θ du
= =K ,
dy b a dv sin θ cos θ dv
∂x ∂y ∂y ∂x
where a = = and b = = − , precisely as the Cauchy-Riemann conditions
∂u ∂v ∂u ∂v
require.
From this it is seen that the mapping of local vectors (du, dv) → (dx, dy) is a scaling and
rotation; i.e., we have a conformal mapping.
Almost all familiar functions are analytical in the complex plane: the exponent, the log-
arithm, trigonometric and hyperbolic functions, and especially power series expansions.
The sum and product of two analytical functions is again analytic.

From the Cauchy-Riemann conditions we still derive the Laplace equations:

∂2x ∂2x ∂2y ∂2y
+ = 0, + = 0.
∂u2 ∂v 2 ∂u2 ∂v 2

Let us try, as a general solution, a series expansion:

∞
X
z = a0 + a1 w + a2 w2 + a3 w3 + . . . = ak w k . (11.17)
k=0

Here we define for the sake of generality

z ≡ (x − x0 ) + iy,
w ≡ (u − u0 ) + iv.

The values
ˆ ϕ0
x0 ≡ M (ϕ) dϕ,
ˆ0
ϕ0
M (ϕ)
u0 ≡ dϕ
0 p (ϕ)

102
 11.6. The Gauss-Krüger projection

Figure 11.8.: The Gauss-Krüger projection as a boundary value problem. The boundary
values at the central meridian are marked with crosses, the directions of integration
with arrows

have been chosen to be compatible with a suitable reference latitude ϕ0 . These series expansions
are thus meant to be used only in a relatively small area, not the whole projection zone. Thus
the number of terms needed also remains smaller.
The meridian conditions (11.16) now form a set of observation equations:

xi = a0 + a1 ui + a2 u2i + a3 u3i + . . .

from which the coefficients aj can be solved, if a sufficient number of support points (ui , xi ) is
given (the crosses in figure 11.8).
Applying the found coefficients aj to equation (11.17) yields the solution for the whole complex
plane. This solution may be “squeezed” into the following, more computationally suitable, form
(so-called Clenshaw summation):

z = a0 + w (a1 + w (a2 + w (a3 + w (· · · + wan )))) .

The computing sequence is multiplication-addition-multiplication-addition. . . the powers of1 w

need not be computed separately. Remember that also the intermediate results are complex
numbers!
Fortunately complex numbers belong to the popular programming languages either as a fixed
part (Matlab, Fortran) or a standard library (C++).
From the same equation we obtain also easily the scale and the meridian convergence. At said,
a complex analytical map is a scaling plus a rotation. When the equation of the map is (11.17),
its differential version is
dz
∆z = ∆w, (11.18)
dw
where
∞
dz X
= a1 + 2a2 w + 3a3 w2 + 4a4 w4 + · · · = kak wk−1 .
dw
k=1

Which also can easily we written into the Clenshaw form.

103
Chapter 11. Map projections

Now we may write the complex equation (11.18) in real-valued matrix form
     
dz dz
 
 Re dw −Im  
dx dw  du
=     
 dv ,
dy  dz dz
Im Re
dw dw
 T  T
where dx = dx dy , dw = du dv ,
 
dz ∂x ∂y
Re = = = K cos θ,
dw ∂u ∂w
the scale; and  
dz ∂y ∂x
Im = =− = K sin θ,
dw ∂u ∂v
from which θ, the meridian convergence, may be resolved.
The choice of the classical Mercator as a starting projection is not the only alternative.
Probably the number of terms of the series expansion (11.17) could be reduced by using a
Lambert projection for the reference latitude ϕ0 . Whether the saving achievable is worth the
trouble, is a so-called Good Question.

11.7. Curvature of the Earth surface and scale

Recall the equation 10.13 derived above:
∂ 2 g11 ∂ 2 g22
 
R=− + .
(∂x2 )2 (∂x1 )2
Because according to equation 10.12 we have R = 2K, we obtain
1 ∂ 2 g11 ∂ 2 g22
 
K=− + .
2 (∂x2 )2 (∂x1 )2

Let us now take a rectangular co-ordinate system in the map plane (x, y) and transfer its co-
ordinate lines back to the curved surface of the Earth, forming a curvilinear co-ordinate system
(ξ, η). At the origin or central meridian in the map plane the metric of this co-ordinate system
gij is, at least for ordinary map projections, quasi-Euclidean, in other words, this metric is the
unit matrix and it is stationary at the origin. In this case the theory in the previous chapter
(section 10.8) applies.
In the case of a conformal projection the form of this metric is
 
−2 1 0
gij = m ,
0 1
where m is the scale of the map projection, which thus depends on location xi .
Derive
∂ 2 gξξ ∂ 2 gηη
 
1 1
= − ∆ m−2 .


K=− +
2 ∂η 2 ∂ξ 2 2
Now
∂2 ∂2
 
−2
m−2 (ξ, η) =


∆ m = +
∂ξ 2 ∂η 2
   
∂ −3 ∂m ∂ −3 ∂m
= −2 m −2 m =
∂ξ ∂ξ ∂η ∂η
∂m−3 ∂m ∂2m ∂m−3 ∂m ∂2m
   
= −2 · + m−3 2 − 2 · + m−3 2
∂ξ ∂ξ ∂ξ ∂η ∂η ∂η
≈ −2∆m,

104
 11.7. Curvature of the Earth surface and scale

if we consider the stationarity of m (and thus of m−3 ) and m ≈ 1.

So
K = ∆m.

So, there is a simple relationship between the Gauss curvature radius and the second derivative
of the scale (more precisely, the Laplace operator ∆) As a consequence, there also is a rela-
tionship between the scale’s second derivatives in the North-South and East-West directions: if,
e.g.,
∂2m ∂2m
= 0 ⇒ =K
∂ξ 2 ∂η 2
etc. In the below table we give some examples – remember that the area considered is always a
neighbourhood of the central point or central meridian or standard parallel, where m ≈ 1 and
stationary!

Projection ∂ 2 m/∂x2 ∂ 2 m/∂y 2

Mercator K 0
Lambert conical K 0
K K
Oblique stereographic
2 2
Gauss-Krüger, UTM 0 K

In the table we used again, instead of (ξ, η),(x, y), i.e., we substituted

∂2m ∂2m ∂2m ∂2m

−→ and −→ ,
∂ξ 2 ∂x2 ∂η 2 ∂y 2
which is allowed based on quasi-Euclidicity.
The scale of Mercator and Lambert projections is a constant in the direction of the y co-
ordinate, in other words, from map West to map East. In transversal cylindrical projections
the scale again is constant in the direction of the x co-ordinate, i.e., the direction of the central
meridian. The oblique stereographic projection again is symmetric and the scale behaves in the
same way in all compass directions from the central point.
Thus, we may use the second derivatives of place of the scale for classifying map projections.
E.g., Lambert is most suited for countries extending in the West-East direction (Estonia),
whereas Gauss-Krüger again is best for countries extending in the North-South direction
(New Zealand). The conformal azimuthal projection called oblique stereographic is suitable for
“square” countries (The Netherlands).

105
Chapter 11. Map projections

N N N

(a) Gauss-Krüger (b) Mercator, Lambert (c) Stereographic

Figure 11.9.: The classification of map projections as “vertical”, “horizontal” and “square” pro-
jections

106
 Chapter 12
Map projections in Finland

12.1. Traditional map projections

In Finland, the traditional map projection has been Gauß-Krüger with zone width 3◦ . System
name: kkj (“National Map Grid Co-ordinate System”) created in 1970 (Parm, 1988). Following
characteristics:
◦ Based on International (Hayford) reference ellipsoid of 1924; datum was taken from the
European datum of 1950 by keeping fixed the triangulation point Simpsiö (nr. 90), at
. latitude and longitude values from the ED50 European adjustment, and
. geoidal undulation from the Bomford astro-geodetic geoid (Bomford, 1963).
◦ Map plane co-ordinates were obtained using Gauß-Krüger for central meridians of 19◦ , 21◦ , 24◦ , 27◦ , 30◦ ;
for small-scale all-Finland maps, 27◦ is used (the ykj system).
◦ These co-ordinates (x, y) were further transformed in the plane using a four-parameter
similarity (“Helmert”) transformation in order to achieve agreement with the pre-existing
provisional vvj (“Old State System”, also “Helsinki System”) co-ordinates, cf. (Ollikainen,
1993).
Equation:
      
x 1.00000075 −0.00000439 x −61.571 m
= +
y kkj
0.00000439 1.00000075 y ED50
95.693 m

12.2. Modern map projections

The modern Finnish system is quite different:

◦ Based on the GRS80 reference ellipsoid of 1980; datum is called EUREF-FIN, created by
keeping fixed four stations fixed to their ITRF96 values at epoch 1997.0: the permanent
GPS stations Metsähovi, Vaasa, Joensuu and Sodankylä (Ollikainen et al., 2000). Then,
a transformation (Boucher and Altamimi, 1995) was applied to obtain co-ordinates in
ETRF89. Thus, the datum is correctly described as ETRF89, but the epoch remains
1997.0, as no correction for individual station motion (mostly, glacial isostatic adjustment)
was made in the transformation.
Equation:
 
0 −Ṙ3 Ṙ2
XE (tC ) = XIyy (tC ) + Tyy +  Ṙ3 0 −Ṙ1  XIyy (tC ) · (tC − 1989.0)
 
−Ṙ2 Ṙ1 0
yy

with tC observations central epoch, yy = (19)96. The values T96 and Ṙi,96 are tabulated in
(Boucher and Altamimi, ) Tables 3 and 4.

107
Chapter 12. Map projections in Finland

◦ For small-scale and topographic maps, the UTM projection is used with a central meridian
of 27◦ (zone 35) for the whole country, producing the ETRS-TM35FIN plane co-ordinate
system. This also defines the map sheet division. However, on maps in parts of Finland
where another central meridian would be more appropriate (like zone 34, central meridian
21◦ ), the corresponding co-ordinate grid is also printed on the map, in purple (Anon.,
2003).

◦ For large scale maps, such as used for planning and cadastral work, the Gauß-Krüger
projection continues to be used (but based on the above reference ellipsoid and datum),
with a central meridian interval of only one degree: ETRS-GKn, where n designates the
central meridian longitude. This avoids the problem of significant scale distortions.

12.3. The triangulated affine transformation used in Finland

12.3.1. Plane co-ordinates

The National Land Survey offers a facility to convert kkj co-ordinates to the new ETRS89-
TM35FIN system of projection co-ordinates. The method is described in the publication
(Anon., 2003), where it is proposed to use for the plane co-ordinate transformation between
the projection co-ordinates of ETRS-89 and the ykj co-ordinate system, a triangle-wise affine
transformation.
Inside each triangle we may write the affine transformation can be written like

x(2) = ∆x + a1 x(1) + a2 y (1)

y (2) = ∆y + b1 x(1) + b2 y (1)
   
where x(1) , y (1) are the point co-ordinates in ETRS-GK27, and x(2) , y (2) are the co-ordinates
of the same point in ykj. This transformation formula has six parameters:
 ∆x,
 ∆y, a1, a2 , b1 ja
(1) (1)
b2 . If, in the three corners of the triangle, are given both x , y and x(2) , y (2) , we can
solve for these uniquely.
The transformation formula obtained is inside the triangles linear and continuous across the
edges, but not differentiable: the scale is discontinuous across triangle edges. Because the
mapping is not conformal either, the scale will also be dependent upon the direction considered.
A useful property of triangulation is, that it can be locally “patched”:
 if better data is available
(i) (i)
in the local area – a denser point set, whose co-ordinate pairs x , y , i = 1, 2 are known –
then we can take away only the triangles of that area and replace them by a larger number of
smaller triangle, inside which the transformation will become more precise. This is precisely the
procedure that local players, like municipalities, can use to advantage.
Write these equations in vector form:
" #   " (1) #
x(2)
 
∆x a1 a2 x
(2) = + .
y ∆y b1 b2 y (1)

 T
Most often the co-ordinates in the (1)and (2) datums are close to each other, i.e., ∆x ∆y
are small. In that case we may write the shifts

δx ≡ x(2) − x(1) = ∆x + (a1 − 1) x(1) + a2 y (1) ,

δy ≡ y (2) − y (1) = ∆y + b1 x(1) + (b2 − 1) y (1) .

108
 12.3. The triangulated affine transformation used in Finland

819
230
31

44 2
85
301 203
22
441

233
415 503
440
5 510
4
118
453
The Lappeenranta base network
Delaunay triangulation (blue)
Voronoi diagram (red)

Figure 12.1.: Lappeenranta densification of the national triangular grid

If we now define
     
∆x a11 a12 a1 − 1 a2
∆x ≡ , A= ≡ ,
∆y a21 a22 b1 b2 − 1

we obtain the short form

δx = ∆x + Ax(1) .

Also in this generally, if the co-ordinates are close together, the elements of A will be numerically
small. Let there be a triangle ABC. Then we have given the shift vectors of the corners
(1)
δxA = ∆x + AxA ,
(1)
δxB = ∆x + AxB ,
(1)
δxC = ∆x + AxC .

Write this out in components, with ∆x, A on the right hand side:
(1) (1)
δxA = ∆x + a11 xA + a12 yA
(1) (1)
δyA = ∆y + a21 xA + a22 yA
(1) (1)
δxB = ∆x + a11 xB + a12 yB
(1) (1)
δyB = ∆y + a12 xB + a22 yB
(1) (1)
δxC = ∆x + a11 xC + a12 yC
(1) (1)
δyC = ∆y + a21 xC + a22 yC

109
Chapter 12. Map projections in Finland

111111111111111
000000000000000
000000000000000
111111111111111 00000000
11111111
APC ) 11111111
C

p B = 111
000
000000000000000
111111111111111 00000000
ω(∆ ABC ) 11111111
00000000
000= ω(∆
111
000000000000000
111111111111111
000000000000000
111111111111111
000000000000000
111111111111111 00000000
11111111
000000000000000
111111111111111 00000000
11111111
p A = 000
111 =
000000000000000
111111111111111
000000000000000
111111111111111 00000000
11111111
00000000
11111111
000000000000000
111111111111111 P 11111111
00000000
0000000000000000000000000000000
1111111111111111111111111111111
0000000000000000000000000000000
1111111111111111111111111111111
000000000000000
111111111111111
0000000000000000000000000000000
1111111111111111111111111111111
000000000000000
111111111111111 00000000
11111111
= ω(∆ PBC )
0000000000000000000000000000000
1111111111111111111111111111111
0000000000000000000000000000000
1111111111111111111111111111111
0000000000000000000000000000000
1111111111111111111111111111111
0000000000000000000000000000000
1111111111111111111111111111111
0000000000000000000000000000000
1111111111111111111111111111111
A
0000000000000000000000000000000
1111111111111111111111111111111 11111111
00000000
ω(∆ ABC )
1111111111111111111111111111111
0000000000000000000000000000000
0000000000000000000000000000000
1111111111111111111111111111111
0000000000000000000000000000000
1111111111111111111111111111111 00000000
11111111
00000000
11111111
0000000000000000000000000000000
1111111111111111111111111111111
0000000000000000000000000000000
1111111111111111111111111111111
0000000000000000000000000000000
1111111111111111111111111111111

00000000
11111111
0000000000000000000000000000000
1111111111111111111111111111111
0000000000000000000000000000000
1111111111111111111111111111111
0000000000000000000000000000000
1111111111111111111111111111111

00000000
11111111
0000000000000000000000000000000
1111111111111111111111111111111
0000000000000000000000000000000
1111111111111111111111111111111
0000000000000000000000000000000
1111111111111111111111111111111

00000000
11111111
0000000000000000000000000000000
1111111111111111111111111111111
0000000000000000000000000000000
1111111111111111111111111111111
0000000000000000000000000000000
1111111111111111111111111111111
0000000000000000000000000000000
1111111111111111111111111111111
0000
1111
0000000000000000000000000000000
1111111111111111111111111111111
0000
1111
0000000000000000000000000000000
1111111111111111111111111111111
0000
1111 00000000
11111111
0000 = ω(∆ ABP )
00000000
11111111
C
0000000000000000000000000000000
1111111111111111111111111111111
1111
p =
0000000000000000000000000000000
1111111111111111111111111111111
0000000000000000000000000000000
1111111111111111111111111111111

00000000
11111111
0000000000000000000000000000000
1111111111111111111111111111111
ω(∆ ABC )
0000000000000000000000000000000
1111111111111111111111111111111
0000000000000000000000000000000
1111111111111111111111111111111
0000000000000000000000000000000 B
1111111111111111111111111111111
Figure 12.2.: Computing barycentric co-ordinates as the ratio of the surface areas of two triangles

or in matrix form
 
(1) (1)
  1 0 xA 0 yA 0  
δxA  (1) (1) 
∆x
 δyA   0 1 0 xA 0 yA   ∆y 
   (1) (1)  

 δxB  
= 1 0 xB 0 yB 0   a11
 
,
δyB   (1) (1)   a21
0 1 0 xB 0 yB 
 
    
 δxC   (1) (1)   a12 
1 0 xC 0 yC 0 

δyC a22

(1) (1)
0 1 0 xC 0 yC

from which they can all be solved.

Let us write the coordinates (x, y) as follows:

x = pA xA + pB xB + pC xC ,
y = pA yA + pB yB + pC yC,

with the further condition pA + pB + pC = 1. Then also

δx = pA δxA + pB δxB + pC δxC , (12.1)

A B C
δy = p δyA + p δyB + p δyC. (12.2)

The set of three numbers pA , pB , pC is called the barycentric co-ordinates of point P See figure

12.2.
ω (∆BCP )
They can be found as follows (geometrically pA = etc., where ω is the surface area
ω (∆ABC)
of the triangle) using determinants:
     
 xB xC x   xC xA x   xA xB x 
    
 yB yC y   yC yA y   yA yB y 
    
 1 1 1   1 1 1   1 1 1 
A B C
p = , p =  , p =  .
 xA xB xC   xA xB xC   xA xB xC 
   
 yA yB yC   yA yB yC   yA yB yC 
   
 1 1 1   1 1 1   1 1 1 
These equations are very suitable for coding.

110
 Appendix A
Isometric latitude on the ellipsoid

We follow the presentation from the book (Grossman, 1964).

The starting formula is ˆ ϕ
M (ϕ0 ) 0
ψ (ϕ) = dϕ .
0 p (ϕ0 )
As a differential equation

M 1 − e2
dψ = dϕ = dϕ.
1 − e2 sin2 ϕ cos ϕ


p

The integrand is decomposed into partial fractions:

1 − e2 1 − e2 sin2 ϕ − e2 cos2 ϕ
= =
1 − e2 sin2 ϕ cos ϕ 1 − e2 sin2 ϕ cos ϕ

1 e2 cos ϕ
= − =
cos ϕ 1 − e2 sin2 ϕ
1 e2 cos ϕ (1 − e sin ϕ) + e2 cos ϕ (1 + e sin ϕ)
= − =
cos ϕ 2 (1 + e sin ϕ) (1 − e sin ϕ)
 
1 e e cos ϕ e cos ϕ
= + − − .
cos ϕ 2 1 + e sin ϕ 1 − e sin ϕ

The integral of the first term is

ˆ ϕ
1 0
π ϕ
ψ= dϕ = ln tan + .
0 cos ϕ0 4 2
Proof by using the chain rule:

d ln tan π4 + ϕ2 d tan π4 + ϕ2 d π4 + ϕ2




dψ
= =
d tan π4 + ϕ2 d π4 + ϕ2



dϕ dϕ
1 1 1
= π ϕ
2 π ϕ
=
tan 4 + 2 cos 4 + 2 2
1 1 1
= π ϕ
π ϕ
= π

= .
2 sin 4 + 2 cos 4 + 2 sin 2 + ϕ cos ϕ

This is the full solution in the case that e = 0 (solution for the sphere).
In the case of the ellipsoid the second integral
ˆ   ˆ 0
e cos ϕ f (ϕ)
− dϕ = dϕ = ln f (ϕ) = ln (1 + e sin ϕ) ,
1 + e sin ϕ f (ϕ)
where we designate f ≡ 1 + e sin ϕ. In the same way
ˆ  
e cos ϕ
− dϕ = − ln (1 − e sin ϕ) ,
1 − e sin ϕ

111
Appendix A. Isometric latitude on the ellipsoid

and the end result is

π ϕ e
ψ = ln tan + + (ln (1 + e sin ϕ) − ln (1 − e sin ϕ)) =
4 2 2
 π ϕ   1 + e sin ϕ  2e
!
= ln tan + .
4 2 1 − e sin ϕ

112
 Appendix B
Useful equations connecting the main radii of
curvature

When we have given the principal radii of curvature of the ellipsoid of revolution:

−1/2
N (ϕ) = a 1 − e2 sin2 ϕ ,

−3/2
M (ϕ) = a 1 − e2 1 − e2 sin2 ϕ


,

we can calculate by brute-force derivation:

d
(N (ϕ) cos ϕ) = −M (ϕ) sin ϕ,
dϕ
d M (ϕ)
(N (ϕ) sin ϕ) = + cos ϕ.
dϕ 1 − e2
Furthermore
d
2 d
−3
M 2 (ϕ) = a2 1 − e2 1 − e2 sin2 ϕ =
dϕ dϕ
e2 M 2 N 2
= 3 sin 2ϕ,
a2
d d
N 2 (ϕ) cos2 ϕ = 2N


(N (ϕ) cos ϕ) =
dϕ dϕ
= −M N sin 2ϕ.

113
 Appendix C
The Christoffel symbols from the metric
Let us start from the definition of the metric in Chapter 9.2:
     
∂x ∂x ∂x ∂x
 ∂u1 · ∂u1 · 2 
 
∂x ∂x 1
gij = · =     ∂u ∂u   ,
∂ui ∂uj  ∂x ∂x ∂x ∂x 
· ·
∂u2 ∂u1 ∂u2 ∂u2
where ui = u1 , u2 is the parametrization (“co-ordinate frame”) of a curved surface in space (or

more generally a sub-space). Differentiate:

 
∂ ∂ ∂x ∂x
gjk = · =
∂ui ∂ui ∂uj ∂uk
 2   2 
∂ x ∂x ∂ x ∂x
= · + · . (C.1)
∂ui ∂uj ∂uk ∂ui ∂uk ∂uj
Correspondingly, by interchanging indices:
 2   2 
∂ ∂ x ∂x ∂ x ∂x
gki = · + · (C.2)
∂uj ∂uj ∂uk ∂ui ∂uj ∂ui ∂uk
 2   2 
∂ ∂ x ∂x ∂ x ∂x
gij = · + · (C.3)
∂uk ∂uk ∂ui ∂uj ∂uk ∂uj ∂ui

Compute equation (C.1) plus equation (C.2) minus equation (C.3):

 2 
∂ ∂ ∂ ∂ x ∂x
gjk + j gki − k gij = 2 · . (C.4)
∂ui ∂u ∂u ∂ui ∂uj ∂uk
 
∂x ∂x
Let’s write the second derivatives of x on the local base , , n as we did in equation
∂u1 ∂u2
(9.7), even if in a slightly different notation:
∂2x ∂x
i j
≡ Γ`ij ` + βij n, (C.5)
∂u ∂u ∂u
which implicitly defines the Γ symbols. Substitution into equation (C.4) yields
     
` ∂x ∂x ∂x 1 ∂ ∂ ∂
Γij · + βij n · k = gjk + j gki − k gij .
∂u` ∂uk ∂u 2 ∂ui ∂u ∂u
Here we recognise    
∂x ∂x ∂x
· = g`k ja n · k = 0,
∂u` ∂uk ∂u
or
 
1 ∂ ∂ ∂
Γ`ij g`k = gjk + j gki − k gij ⇒
2 ∂ui ∂u ∂u
 
1 `k ∂ ∂ ∂
⇒ Γ`ij = g gjk + j gki − k gij ,
2 ∂ui ∂u ∂u
equation (10.7). Here, g ij gjk = δki , i.e., g ij is the inverse matrix of gij .

115
 Appendix D
The Riemann tensor from the Christoffel
symbols

The Riemann tensor equation is derived with the age of parallel transport of a vector around
a closed, small co-ordinate rectangle. Of the spatial vector v is transported parallelly inside a
ui -parametrized surface S , the following holds

∂v
= 0.
∂ui
Let us write v on the basis of the tangent vectors:

∂x
v = vi .
∂ui
Then
∂v ∂v i ∂x i ∂ x
2
0= = + v =
∂uj ∂uj ∂ui ∂ui ∂uj
∂v i ∂x ∂x
= + Γijk v k i + βij n,
∂uj ∂ui ∂u
∂v i ∂x ∂x
using (C.5). Here, the v derivative consists of two parts: the “interior” part, j i
+Γijk v k i ,
∂u ∂u ∂u
embedded in the surface, and the “exterior” part, βij n, aperpendicular to the surface. When the
surface S has been given, we can only zero the internal part, i.e.

∂v i
+ Γijk v k = 0 (D.1)
∂uj

describes the parallel transport of the vector v i within the surface.

Let us now consider a small rectangle ABCD, side lengths ∆uk and ∆u` (see Fig. 10.3), along
co-ordinate curves. The sides AB and CD are on opposite sides, the running co-ordinate being
uk . Similarly BC and AD are opposite, the running co-ordinate being u` .
Derive the change in v i over the distance AB :

∂v i
∆AB v i = ∆uk = −Γikm v m ∆uk .
∂uk
In the same way
∆CD v i = +Γikm v m ∆uk .
For the side BC we obtain
∂v i
∆BC v i = ∆u` = −Γi`m v m ∆u`
∂u`
and
∆DA v i = +Γi`m v m ∆u` .

117
Appendix D. The Riemann tensor from the Christoffel symbols

sum together these four terms:

∆ABCD v i =
 i m

Γkm v CD − Γikm v m AB ∆uk − Γi`m v m DA − Γi`m v m BC ∆u` =




   
∂ i m ` k ∂
Γ v ∆u ∆u` =
i m k



= Γ v ∆u ∆u −
∂u` km ∂uk `m
 i
∂Γkm ∂Γi`m m m
 
m i ∂v i ∂v
= − v + Γ km − Γ `m ∆u` ∆uk .
∂u` ∂uk ∂u` ∂uk

Equation (D.1) gives

∂v m m h ∂v
m
= −Γ`h v , = −Γm h
kh v ;
∂u` ∂uk
substituting:

∂Γikm ∂Γilm
  
∆ABCD v i = v m + Γilm Γm i m

h
− kh − Γ Γ
km lh v ∆u` ∆uk =
∂ul ∂uk
!
∂Γikj ∂Γi`j
= − + Γi`m Γm i m
kj − Γkm Γ`j v j ∆u` ∆uk ,
∂u` ∂uk

where we have changed the names of the indices m → j (in the first two terms) ja h → j (in the
last two terms).
Here we see in ready form the Riemann curvature tensor :

i
∂Γikj ∂Γi`j
Rj`k = − + Γi`m Γm i m
kj − Γkm Γ`j ,
∂ul ∂uk
aapart from the names of the indices and interchanges of type Γijk = Γikj , just what already was
given in eq. (10.9).

118
Bibliography

Anon. (2003). JHS154. ETRS89 -järjestelmään liittyvät karttaprojektiot, tasokoordinaatistot

ja karttalehtijako (Map projections, plane co-ordinates and map sheet division in relation
to the ETRS89 system). Web site, Finnish National Land Survey. URL: http://www.
jhs-suositukset.fi/suomi/jhs154, accessed August 30, 2005.
Anon. (2008). JHS153. ETRS89-jr̈jestelmn̈ mukaiset koordinaatit Suomessa. Web site, Finnish
National Land Survey. URL: http://www.jhs-suositukset.fi/suomi/jhs153, accessed
Oct. 19, 2010.
Bomford, G. (1963). Report of Study Group No. 10. The Geoid in Europe and Connected
Countries. In Traveaux de l’Association Internationale de Géodésie, number 22.
Boucher, C. and Altamimi, Z. Specifications for reference frame fixing in the analysis of a
EUREF GPS campaign. Version 4, 24-10-2008, http://etrs89.ensg.ign.fr/memo-V7.
pdf.
Boucher, C. and Altamimi, Z. (1995). Specifications for Reference Frame Fixing in the Analysis
of a EUREF GPS Campaign. In VI Meeting of the TWG in Bern, March 9 and 10, Veröff.
der Bayerischen Kommission für die Internationale Erdmessung, number 56, pages 265–267,
München.
Gross, R. S. (2000). The excitation of the Chandler wobble. Geophys. Res. Lett., 27:2329–2332.
Grossman, W. (1964). Geodätische Rechnungen und Abbildungen. Verlag Konrad Wittwer,
Stuttgart.
Heikkinen, M. (1981). Solving the shape of the Earth by using digital density models. Report
81:2, Finnish Geodetic Institute, Helsinki.
Heiskanen, W. A. and Moritz, H. (1967). Physical Geodesy. W.H. Freeman and Company, San
Francisco.
Kumar, M. and Reilly, J. P. (2006). Is definition of WGS 84 correct. URL:
http://mycoordinates.org/is-definition-of-wgs-84-correct/.
McCarthy, D. D. (1996). IERS Conventions (1996). (IERS Technical Note ; 21).
Technical report, Central Bureau of IERS - Observatoire de Paris. URL:
http://www.iers.org/nn 11216/IERS/EN/Publications/TechnicalNotes/tn21.html.
Moritz, H. (1992). Geodetic Reference System 1980. In Tscherning, C. C., editor, Geodesist’s
Handbook 1992, pages 187–192. International Association of Geodesy, Copenhagen.
Ollikainen, M. (1993). GPS-koordinaattien muuntaminen Kartastokoordinaateiksi (Transform-
ing GPS co-ordinates to the National Grid Co-ordinate System). Tiedote 8, Finnish Geode-
tic Institute, Helsinki.
Ollikainen, M., Koivula, H., and Poutanen, M. (2000). The densification of the EUREF network
in Finland. Publication 129, Finnish Geodetic Institute, Kirkkonummi.
Parm, T. (1988). Kansallisen koordinaattijärjestelmän luominen Suomessa (The creation of a
national co-ordinate system in Finland). Maanmittaus, (1).
Stevenson, P. (2008). The Truth About WGS84. URL: http://www.ashgps.com/out3/General

119