Mathematical Geodesy Maa-6.3230
Martin Vermeer 21st March 2012
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 3 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, solve, using tools to ols 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 transformations
◦ Masters the basics of Gaussian and Riemannian surface theories and can derive the metric tensor, Christoffel symbols and curvature tensor for simple surfaces
◦ Masters Masters the basic math of map projections, projections, esp. conformal conformal ones, and the behaviour behaviour of map scale, and is able to compute the isometric latitude.
Content Spherical trigonometry, geodetic co-ordinate computations in ellipsoidal and rectangular spatial co-ordinate systems, astronomical co-ordinates, co-ordinate system transformations, 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 Materials Lecture notes. Background material Hirvonen: Matemaattinen geodesia; Torge: Geodesy. Teaching Language Suomi Course Staff and Contact Info Martin Vermeer, Vermeer, room Gentti Gentti 4.143, name@aalto.fi name@aalto.fi Reception Reception times Wednesdays at 12:15 CEFR-taso Lis¨ Lis¨ atietoja atietoja
i
Contents 1. Spherical Spherical trigonom trigonometry etry 1.1. Introdu Introducti ction on . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2. Spheri Spherical cal exce excess ss . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3. The surfac surfacee area of a triangle triangle on a sphere sphere . . . . . . . . . . . 1.4. A rectangular rectangular spherical spherical triangle triangle . . . . . . . . . . . . . . . . . 1.5. A general general spheric spherical al triangle triangle . . . . . . . . . . . . . . . . . . . 1.6. Deriving Deriving the formulas formulas with the the aid of vectors vectors in space . . . . 1.7. Polarizatio Polarization n . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.8. Solving Solving the spherical spherical triangle by the method method of additaments additaments . 1.9. Solving Solving the spherica sphericall triangle triangle by Legendre’s method . . . . 1.10. The forward geodetic problem on the sphere . . . . . . . . . . 1.11. The geodetic inverse problem on the sphere . . . . . . . . . . 1.12. The half-angle cosine rule . . . . . . . . . . . . . . . . . . . . 2. The 2.1. 2.2. 2.3. 2.4. 2.5.
geometry geometry of the reference reference ellipsoid ellipsoid Introdu Introducti ction on . . . . . . . . . . . . . . . . . . . . . . . . . . . . The geodesic geodesic as solution solution to a system of differentia differentiall equations equations An invarian invariantt . . . . . . . . . . . . . . . . . . . . . . . . . . . . The geodet geodetic ic main main problem problem . . . . . . . . . . . . . . . . . . . The geodetic geodetic inv inverse erse problem problem . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . .
. . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . .
3. Co-ordina Co-ordinates tes on the reference ellipsoid ellipsoid 3.1. Representa Representations tions of the sphere sphere and the ellipsoid ellipsoid . . . . . . . . . . . . 3.2. Various latitude latitude types types . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3. Mea Measur sures es for the flatte flattenin ningg . . . . . . . . . . . . . . . . . . . . . . . 3.4. Relationship Relationshipss between between different different types of latitude . . . . . . . . . . . 3.5. Co-ordinates Co-ordinates in the the meridional meridional ellipse ellipse . . . . . . . . . . . . . . . . . . 3.6. Three-dimen Three-dimensional sional rectangular rectangular co-ordinates co-ordinates on the reference reference ellipsoid ellipsoid 3.7. Computing Computing geographic geographic co-ordinates co-ordinates from rectangula rectangularr ones . . . . . . 3.8. Meridi Meridian an arc leng length th . . . . . . . . . . . . . . . . . . . . . . . . . . . 4. Reference Reference systems systems 4.1. The GRS80 GRS80 system system and geometric geometric parameters parameters . 4.2. Gravimetri Gravimetricc parameters parameters . . . . . . . . . . . . . 4.3. Refere Reference nce fram frames es . . . . . . . . . . . . . . . . . 4.4. The orien orientati tation on of the Earth Earth . . . . . . . . . . . 4.5. Transformations ransformations betwee between n systems . . . . . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . . . . . . . . . . . . . . . . . . . . .
. . . . .
5. Using rotation rotation matrices matrices 5.1. Genera Generall . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2. Chaining Chaining matrices matrices in three three dimensions dimensions . . . . . . . . . . . . . . . . . 5.3. Orthogonal Orthogonal matrices matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4. Topocentric opocentric system systemss . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.5. From geocentric geocentric to topocentric topocentric and back back . . . . . . . . . . . . . . . . 5.6. The geodetic main and inver inverse se problems problems with rotation matrices matrices . . . 5.6.1. 5.6.1. Geodeti Geodeticc main main probl problem em . . . . . . . . . . . . . . . . . . . . .
. . . . . . . .
. . . . .
. . . . . . . . . . . .
. . . . . . . . . .
. . . . . . . .
. . . . .
. . . . . . . . . . . . . .
. . . . . . . .
. . . . .
. . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . .
1 1 1 2 3 4 5 6 7 8 8 9 9
11 . . . . 11 . . . . 11 . . . . 12 . . . . 13 . . . . 13 . . . . . . . .
. . . . .
. . . . . . . . . . . . . . . .
. . . . .
. . . . .
. . . . . . . .
15 15 15 16 16 17 17 18 19
. . . . .
21 21 22 23 24 25
29 . . . . . 29 . . . . . 30 . . . . . 30 . . . . . 32 . . . . . 33 . . . . . 34 . . . . . 34
iii
Contents 5.6.2. Geodetic inverse inverse problem problem . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.6.3. Comparison Comparison with with ellipsoid ellipsoidal al surface surface solution solution . . . . . . . . . . . . . . . . .
6. Co-ordina Co-ordinate te systems and transform transformations ations 6.1. Geocentric Geocentric terrest terrestrial rial system systemss . . . . . 6.2. Conventi Conventional onal Terrest Terrestrial rial System . . . . 6.3. Polar Polar motion motion . . . . . . . . . . . . . . 6.4. The Instantaneou Instantaneouss Terres Terrestial tial System System . 6.5. The quasi-ine quasi-inertial rtial geocentri geocentricc system system .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
7. Reference Reference systems and and realizations realizations 7.1. Old and new referenc referencee systems; systems; ED50 vs. vs. WGS84/GRS80 WGS84/GRS80 7.2. WGS84 WGS84 and ITRS ITRS . . . . . . . . . . . . . . . . . . . . . . 7.3. Co-ordinate Co-ordinate system system realizations realizations . . . . . . . . . . . . . . . 7.4. Realiz Realizati ation on of WGS84 WGS84 . . . . . . . . . . . . . . . . . . . . 7.5. Realizations Realizations of ITRS/ETRS ITRS/ETRS systems systems . . . . . . . . . . . . 7.6. The three-dimens three-dimensional ional Helmert Helmert transformation transformation . . . . . . 7.7. Transformations ransformations betwee between n ITRF realizations realizations . . . . . . . . 8. Celestial Celestial co-ordin co-ordinate ate systems 8.1. Sidere Sidereal al time time . . . . . . . . . . . . . . . . . . . . . 8.2. Trigonometry rigonometry on the celestial celestial sphere sphere . . . . . . . . 8.3. Using rotation rotation matrice matricess . . . . . . . . . . . . . . . 8.4. On satell satellite ite orbi orbits ts . . . . . . . . . . . . . . . . . . 8.5. Crossing Crossing a given latitude latitude in the inertial inertial system system . . 8.6. The satellite satellite’s ’s topocentric topocentric co-ordinat co-ordinates es . . . . . . . 8.7. Crossing Crossing a given latitude latitude in the terrestri terrestrial al system system . 8.8. Determining Determining the the orbit from observ observations ations . . . . . . 9. The 9.1. 9.2. 9.3. 9.4. 9.5. 9.6.
surfa surface ce theory theory of Gauss Gauss A curve curve in spac spacee . . . . . . . . . . . The first first fundament fundamental al form (metric) (metric) The second second fundamen fundamental tal form . . . . Principal Principal curvatures curvatures . . . . . . . . . A curve curve embedde embedded d in a surface surface . . . The geodesic geodesic . . . . . . . . . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . . . . . . . . . .
. . . . . . . . . . . . . .
. . . . . . . . . . . . . .
. . . . . . . . . . . . . .
10.The surface theory of Riemann 10.1. What is a tensor? . . . . . . . . . . . . . . . . . . . . . . . 10.2. The metric tensor . . . . . . . . . . . . . . . . . . . . . . . 10.3. The inverse metric tensor . . . . . . . . . . . . . . . . . . 10.3.1. Raising or lowering sub- or superscripts of a tensor 10.3.2. The eigenvalues and -vectors of a tensor . . . . . . 10.3.3. The graphic representation of a tensor . . . . . . . 10.4. The Christoffel symbols . . . . . . . . . . . . . . . . . 10.5. The geodesic revisited . . . . . . . . . . . . . . . . . . . . 10.6. The curvature tensor . . . . . . . . . . . . . . . . . . . . . 10.7. Gauss curvature and spherical excess . . . . . . . . . . . 10.8. The curvature in quasi-Euclidean geometry . . . . . . . . 11.Map projectio projections ns 11.1. Map projections and scale . . 11.1.1. On the Earth surface . 11.1.2. In the map plane . . . 11.1.3. The scale . . . . . . .
iv
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . .
. . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . .
. . . .
. . . . . . . . . . . . . . . . . . . .
. . . . .
. . . . . . .
41 . . . . . . . 41 . . . . . . . 41 . . . . . . . 42 . . . . . . . 42 . . . . . . . 42 . . . . . . . 43 . . . . . . . 44
. . . . . . . .
. . . . . . . .
47 . . . . 47 . . . . 49 . . . . 50 . . . . 52 . . . . 52 . . . . 53 . . . . 54 . . . . 54
. . . . . .
57 57 58 60 61 65 67
. . . .
89 89 89 90 90
. . . .
. . . .
. . . . . .
. . . . .
. . . . . . . . . . .
. . . .
. . . . . .
. . . . .
71 . . . . . . . . 71 . . . . . . . . 75 . . . . . . . . 77 . . . . . . . . 78 . . . . . . . . 78 . . . . . . . . 78 . . . . . . . . 78 . . . . . . . . 80 . . . . . . . . 81 . . . . . . . . 84 . . . . . . . . 86 . . . .
. . . . . .
. . . . .
. . . . . .
. . . .
. . . . . .
. . . . .
. . . . . .
. . . .
. . . . . .
. . . . . . . .
. . . . .
37 37 37 37 38 38
. . . . .
. . . . . . . .
. . . . .
35 35
. . . .
. . . . . .
. . . .
Contents 11.1.4. The Tissot-indicatrix . . . . . . 11.2. The Lambert projection (LCC) . . . . 11.3. On the isometric latitude . . . . . . . . 11.4. The Mercator projection . . . . . . . 11.5. The stereographic stereographic projecti pro jection on . . . . . . ¨ger projection . . . . . 11.6. The Gauss-Kr¨ Gauss-Kruger u 11.7. Curvature of the Earth surface and scale
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
12.Map projections in Finland 12.1. Traditional map projections . . . . . . . . . . . . . . . . 12.2. Modern map projections . . . . . . . . . . . . . . . . . . 12.3. The triangulated affine transformation used in Finland . 12.3.1. Plane co-ordinates co-ordinates . . . . . . . . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . .
. . . . . . .
. . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . .
. 92 . 93 . 95 . 96 . 97 . 99 . 102
. . . .
105 105 105 106 106
. . . .
A. Isometric latitude on the ellipsoid
109
B. Useful equations connecting the main radii of curvature
111
C. The Christoffel Christoffel symbols from from the metric
113
D. The Riemann tensor from the Christoffel symbols
115
v
Chapter 1 Spherical trigonometry 1.1. Introduction 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 spherical surface. surface. Even Even in case of the whole Earth, the deviation deviation from spherical spherical shape is only 0.3%. The The starr starry y sky sky aga again in ma may y be trea treate ted d as a prec precis isee sphe spheri rica call surfa surface ce,, the the radiu radiuss of whic which h is undefined; in practical computations we often set R = 1.
1.2. Spherical Spherical excess excess See figure 1.1 1.1.. Let us assume assume that that the radius radius of the sphere sphere is 1. The “fron “frontt half” of the sphere sphere is a semi-sphere, the surface area of which is 2π 2 π . A triangle is formed between three great circles. The same great circle circle form, on the back surface surface of the sphere, sphere, a “antipode “antipode triangle” of the same size and shape. When the surface area of the whole semi-sphere is 2 π, the area of the “orange “orange slice” slice” bounded 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.: 1.1.: Spherical Spherical triangles on a semi-sphere semi-sphere.. The back side surface surface of the sphere has been depicted depicted in a lighter lighter shade with its “antipode “antipode triangles”
1
Chapter Chapter 1. Spherical Spherical trigonometry trigonometry i.e., A1 = α + β + β + γ
− π = ε,
where ε is called the spherical excess . If the radius of the sphere is not 1 but R, we obtain A1 = εR2
⇒ ε = RA12 .
Here ε is expressed in radians . If ε is not in radians, we may write ε [unit] =
ρunit A1 , 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 variabl ariablee for surfac surfacee integ integrat ration ion,, and K is the total curvature of the surface accordi according ng to K.F. K.F. Gauss, Gauss, which which thus thus can vary from place place to place. place. E.g., E.g., on the surfac surfacee of an ellipsoid 1 K = , M N where M is the meridional curvature (in the North-South direction) and N the so-called so-called transverse 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. 1.3. The surfac surface e area of a triang triangle le on a sphere 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.
·
Because according to the sine rule b = a
sin β , it also follows that sin α
a2 sin β sin β sin γ A= . 2sin α
2
1.4. A rectangular rectangular spherical triangle triangle
A A
α
1
1 b O c
cos b
sin b
D
a O cos cos a
sin a cos b b
E
c
b
.
a β
E .
b cos b
O
D
β
a
A
D γ .
1
C
c cos c
O
A sin c
sin c
E
E
sin c sin β
β
D sin c cos β B 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: ε= where also the approximate value for R, e.g., R
A , R2
≈ a = 6378137 m, is completely completely sufficient. sufficient.
1.4. A rectangular rectangular spherical spherical triangle 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 cos α = cos a
sin b = cos a sin β, sin c
according to the last equation in group (1.1 ( 1.1). ). By dividing the first of group (1.2 ( 1.2)) by its second, we obtain cot α = cot a sin b and the second of group (1.1 ( 1.1)) by its third, correspondingly cot β = cot b sin a.
3
Chapter Chapter 1. Spherical Spherical trigonometry trigonometry
a
c h
γ
α
b2
b1
. Figure 1.3.: General spherical triangle
1.5. A general spherica sphericall triangle triangle The formulae for a spherical triangle are obtained by dividing the triangle into two right triangles, see figure 1.3 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 sin(b
− b2) cos b1 = cos cos (b − b2 )
= sin b cos b2
− cos b sin 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” formula formulass we obtain obtain sin c sin α = sin a sin γ, or, more generally,
sin a sin b sin c = = , sin α sin β sin γ
the so-called sine rule of spherical trigonometry .
4
(1.4)
1.6. Deriving Deriving the formulas with with the aid of vectors vectors in space C
γ b
a
α
A
ϕA
xC
c
xA
β
b a
c
B
xB
O
ϕB y
λB
x
Figure 1.4.: The spherical triangle ABC , ABC , for deriving the cosine and sine rules using vectors in space
For comparison the corresponding formulas for a plane triangle: c2 = a 2 + b 2 and
− 2ab cos γ
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., etc., in other words, words, the spheri spherical cal sine rule morphs morphs into the one for a plane plane triangle. triangle. For the cosine rule this is not immediately clear.
−→
1.6. 1.6. Derivi Deriving ng the formul formulas as with with the aid of vectors vectors in space space
π If on the spher sphere e 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 =
xA yA zB
=
cos ϕA 0 sin ϕA
, xB =
xB yB zB
=
−
cos ϕB cos λB cos ϕB sin λB 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 xA
× xB
=
=
− sin ϕA cos ϕB sin λ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 Chapter 1. Spherical Spherical trigonometry trigonometry
π
π
π 2
− β
−c
γ π 2
b
a
β
c
−a
π 2
π 2 π
−α
π 2
α π
π
π 2
π
−b
− γ
Figure 1.5.: Polarization of a spherical triangle
When the third vector is xC =
0 0 1
,
we obtain the volume of the parallelepiped spanned by the three vectors (i.e., twice the volume of the tetrahedron ABCO) ABCO) as follows: Vol xA , xB , xC = (xA
{
}
× xB ) · xC = sin b sin a sin γ.
The volume contained contained by the three vectors vectors does not however however depend on the order of the vectors, 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. Pola Polarizat rization ion For every corner of the triangle we may define an “equator” or great circle, one “pole” of which is that corner corner point. point. If we do this, this, we obtain obtain three three “equ “equator ators” s”,, which which themselv themselves es 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. length. And the “polarization “polarization triangle” triangle”’s ’s angle is 180◦ minus this.
6
1.8. Solving Solving the spherical spherical triangle by the method of additament additaments s The polarization method is symmetric : the the origi origina nall trian triangle gle is also also the the polari polariza zati tion on of the the polariza polarization tion triangle triangle.. The interse intersecti ction on point point of two two edges edges of the polarizatio polarization n triangl trianglee is at a ◦ distance of 90 from both “poles”, “poles”, i.e., the corners corners of the original triangle, triangle, and the edge between between them thus thus is the “equator” “equator” of the intersection intersection point. 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◦ cos (180◦
− α) − α)
= sin α =
− cos α
Because of this we obtain of the spherical trigonometry cosine rule ( 1.3 1.3)) the following polarized version : cos γ = ( cos β ) ( cos α) + sin β sin β sin α ( cos c)
−
−
−
−
or cos γ =
− cos β cos β cos α + sin β 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 Solving the spherical triangle triangle by the method of additament additamentss In the additaments method we reduce a spherical triangle to a plane triangle by changing the lengths lengths of the sides. sides. Generally Generally all angles and one side of a triangle 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): 1 3 s s2 sin ψ = ψ ψ +... 1 . 6 R 6R2
−
≈
−
Now the spherical sine rule is (s ( s = a,b,c): a,b,c): a (1 ∂a) ∂a ) b (1 ∂b) ∂b ) c (1 ∂c) ∂c ) = = , sin α sin β sin γ
−
−
−
s2 a2 b2 c2 jossa ∂s = : ∂a = , ∂b = ja ∂c = . 6 R2 6R2 6 R2 6R2 The method works now so, that 1. From the known side we subtract its additament additament ∂s; ∂s ; 2. The other sides are computed computed using the sine rule for a plane triangle triangle and the known angle values ; 3. To the computed sides are now added their additaments. The additaments additaments are computed using the best availabl availablee approximate approximate values; values; if they are initially initially poor, we iterate. The additament correction s = s (1
− ∂s) ∂s )
may be changed from a combination of a multiplication and a subtraction into a simple subtraction by taking logarithms: ln s = ln s + ln ln (1
− ∂s) ∂s ) = ln s + (0 − ∂s) ∂s ) = ln s − ∂s, 7
Chapter Chapter 1. Spherical Spherical trigonometry trigonometry
N ∆λ
2
ψ 1 2
A12
ϕ2
1
ψ12 ϕ1
.
.
o r i t to e k vaa t t
λ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 s =10 log s
− M · ∂s,
where M =10 log e = 0.434 434294 29448. 48. In the the age of log logari arith thm m tabl tables es this made the prac practi tica call computations significantly easier.
1.9. Solving Solving the spherical triangle triangle by Legendre’s 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 sin(α sin(α
− ε/3) ε/3)
=
b c = , sin(β sin(β ε/3) ε/3) sin(γ sin(γ ε/3) ε/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. Nowada Nowa days ys these these appro approxim ximate ate methods (additam (additamen ents ts and Legendre Legendre)) are no longer longer used. used. It is easy to compute directly by computer using the spherical sine rule.
1.10. The forwa forward rd geodetic problem problem on the sphere 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 invers inversee problem on the sphere and sin(λ sin(λ2 λ1 ) sin ψ12
−
=
sin A12 cos ϕ2
⇒
λ2 = λ1 + arcsin
sin ψ12 sin A12 cos ϕ2
.
1.11. 1.11. The geodetic geodetic invers inverse e probl problem em on the sphere sphere The cosine and sine rules of spherical trigonometry, suitably applied:
cos ψ12 = sin ϕ1 sin ϕ2 + cos ϕ1 cos ϕ2 cos(λ cos(λ2 sin(λ sin(λ2 λ1 ) sin A12 = cos ϕ2 . sin ψ12
−
− λ1) ,
(1.5) (1.6)
1.12. The half-angle half-angle cosine cosine rule The above spherical cosine rule (1.3 ( 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 triangle HelsinkiHelsinki-T Tampere-Turku ampere-Turku is very small compared to the globe, globe, some 200 km km/ /637 63788 km 0.03. Then Then cos b cos c 0.999, but sin b sin c 0.0009 0009!! We are adding two two terms of which which one is approx. approx. 1 and the other about a thousand times smaller. smaller. That is the way to lose numerical precision.
∼
∼
∼
For solving this, we first remark, that
− 2sin2 α2 ; a 1 − 2sin2 ; 2
cos α = 1 cos a = and
cos b cos c + sin b sin c cos α = (cos b cos c + sin b sin c) + sin b sin c (cos α
− 1) =
− 2sin b sin c sin2 α2 ; b−c cos b cos c + sin b sin c = cos(b cos(b − c) = 1 − 2sin2 ; = (cos b cos c + sin b sin c)
2
after after a few rearra rearrange ngemen ments ts we obtai obtain n the half-a half-angl ngle e cosine cosine rule rule for a spheri spherica call 1 triangle, which also for small triangles is wel behaved :
sin2
1
a b c α = sin2 + sin b sin c sin2 . 2 2 2
−
. . . which, which, surprise, surprise, contains contains only only sines! sines!
9
Chapter 2 The geometry of the reference ellipsoid 2.1. Introduction Introduction The reference reference ellipsoid ellipsoid is already a pretty pretty precise description description of the true figure of the Earth. The deviations of mean sea level from the GRS80 reference ellipsoid are of magnitude 100m.
±
Regret Regrettabl tably y, the geometry geometry of the referenc referencee ellips ellipsoid oid is not as simple simple as that that of the sphere. sphere. Neverthel Nevertheless, ess, rotational symmetry brings in at least one beautiful beautiful invarian invariant. t. 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. 2.2. The geodesic geodesic as soluti solution on to a sys system tem of differe differenti ntial al equati equations ons We may write in a small rectangular triangle (dy,dx ( dy,dx the metric east and north shift, see Fig. 2.1): 2.1 ): dx = M ( M (ϕ) dϕ = cos Ads, dy = p (ϕ) dλ = sin Ads ja Ads ja dA = sin ϕdλ, where p = N cos 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ϕ ds dλ ds dA ds
cos A , M ( M (ϕ) sin A = , p (ϕ) dλ sin ϕ sin A = sin ϕ = . ds p (ϕ) =
(2.1)
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 ( N (ϕ)cos ϕ, i.e., dλ ds dA ds
= =
sin A , N ( N (ϕ)cos ϕ tan ϕ sin A . N ( 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. Numeri Numerical cally ly computin computingg the solutio solution n using using the
11
Chapter Chapter 2. The geometry geometry of the reference reference ellipsoid ellipsoid
ϕ dA
N
p
c o t ϕ
−
d p
p
.
π 2
dx
−ϕ
p − d p dx p =
N
dλ
c o s ϕ
A
.
ds dy
M
A
ds
.
dx
dy
dϕ
Figure 2.1.: The geometry of integrating the geodesic
MatLab software is also fairly easy thanks to the ODE routines (“Ordinary Differential Equation” E quation”)) provided. Transformation to rectangular form is easy: X ( X (s) = N ( N (ϕ (s))cos ϕ (s)cos λ (s) , Y (s) = N ( N (ϕ (s))cos ϕ (s)sin λ (s) , b2 Z ( Z (s) = N ( N (ϕ)sin ϕ (s) . a2
2.3. 2.3. An invaria invariant nt Let us look closer at the quantity p (ϕ) = N ( N (ϕ)cos ϕ, the distance of a point from the rotation axis. let us compute the s derivative: dp dp dx = = ds dx ds
− sin ϕ cos A.
Now, with the equation 2.1 dA sin ϕ sin A = ds p we obtain by division
dp = dA
A − cos p. sin A
Now
12
2.4. The geodetic geodetic main problem problem d ( p sin A) dA
dp sin A + p cos A = 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ϕ ds dλ ds
= =
cos A (ϕ) , M ( M (ϕ) sin A (ϕ) , p (ϕ)
(2.2) (2.3)
where A (ϕ) is obtained from the invariant formula p (ϕ1 ) sin A (ϕ) = sin A12 . p (ϕ)
(2.4)
The shape of the object is defined by giving the function p (ϕ).
2.4. 2.4. The geodetic geodetic main probl problem em Solving the forward geodetic problem now amounts simply to substituting the also given arc length s12 into 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:
− −
M ( M (ϕ) = a 1
e2
1
e2 sin2 ϕ
−
p (ϕ) = N ( N (ϕ)cos ϕ = a 1
3/2
−
,
e2 sin2 ϕ
1/2
−
· cos ϕ.
2.5. The geodetic inverse inverse proble problem m A direct numerical method for solving the inverse problem is iteration . Let us have as given ϕ1 , λ1 and ϕ2 , λ2 . First we compute the approximate approximate values values and solve the geodetic forward problem in order to compute closing errors (1)
∆ϕ2
(1)
∆λ2 2
≡ ≡
(1)
(1)
ϕ2 , λ2
2
(1) (1) A12 , s12
,
. Thus Thus we obtain obtain the
(1)
ϕ2
− ϕ2, (1) λ2 − λ2 .
. . . e.g., using the closed formulas formulas of spherical spherical trigonometry trigonometry..
13
Chapter Chapter 2. The geometry geometry of the reference reference ellipsoid ellipsoid Next, these closing errors could be used for computing improved values
(2) (2) A12 , s12
, and so on.
The nearly linear linear dependence dependence between between ( A12, s12 ) and (ϕ (ϕ2 , λ2 ) can be approximated using spherical geometry . The formulas needed (1.5, ( 1.5, 1.6): 1.6):
cos s12 = sin ϕ1 sin ϕ2 + cos ϕ1 cos ϕ2 cos(λ cos(λ2 sin s12 sin A12 = cos ϕ2 sin(λ sin(λ2
− λ1) .
− λ1) ,
From this, by differentiation:
− sin s12∆s12
= [sin ϕ1 cos ϕ2
sin s12 cos A12 ∆A12 + cos s12 sin A12 ∆s12 =
− cos ϕ1 sin ϕ2 cos(λ cos(λ2 − λ1 )] ∆ϕ2 − − cos ϕ1 cos ϕ2 sin(λ sin(λ2 − λ1 ) ∆λ2 , − sin ϕ2 sin(λ sin(λ2 − λ1 ) ∆ϕ2 + +cos ϕ2 cos(λ cos(λ2 − λ1 ) ∆λ2 .
So, if we write A=
sin ϕ1 cos ϕ2 cos ϕ1 sin ϕ2 cos(λ cos(λ2 sin ϕ2 sin(λ sin(λ2 λ1 )
−
−
B=
−
− λ1) − cos ϕ1 cos ϕ2 sin(λ sin(λ2 − λ1 ) cos ϕ2 cos(λ cos(λ2 − λ1 ) − sin s12 ,
0 sin s12 cos A12 cos s12 sin A12
,
we obtain as iteration formulas :
(i+1) s12 (i+1) A12
(i) s12 (i) A12
=
(i+1)
Using the new values s12
(i+1)
, A12
(i+1)
(i+1)
+
(i) ∆s12 (i) ∆A12
=
(i) s12 (i) A12
+B
1
−
(i)
A
∆ϕ2 (i) ∆λ2
.
we repeat the computation of the geodetic forward problem
to obtain new values ϕ2 , λ2 until until convergenc convergence. e. The matrices A, B can be recomputed if needed using better approximate values.
14
Chapter 3 Co-ordinates on the reference ellipsoid 3.1. Represe Representati ntations ons 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 obtain an ellipse by “squeezing” “squeezing” the y axis by the factor b/a, b/a, i.e., x = a cos β, y = b sin β, from which again or
x a
2
+
y b
2
= sin2 β + β + cos2 β = 1
x2 y2 + 2 a2 b
−1 =0
is the implicit representation.
3.2. Various Various latitude latitude types 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 Geographical latitude latitude ϕ: the direction of the ellipsoidal normal relative to the plane of the equator; 2. Geocentric Geocentric latitude latitude φ (tai ψ ): the angle angle of the line connect connecting ing the point with the origin, origin, relative to the plane of the equator; 3. Reduced Reduced latitude β : poin point P is shifted straight in the y direction to the circle around the merian ellipse, to become point Q. The geocentric geocentric latitude latitude of point p oint Q is the reduced latitude of point P . P . Reduced Reduced latitude is used only in theoretical theoretical contexts. contexts. In maps, geographical 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 Chapter 3. Co-ordinates Co-ordinates on the reference ellipsoid ellipsoid y
Q
.P b
φ
β
ϕ S
O
a
x
. R
U
V
T
Figure 3.1.: The meridian ellipse and various types of latitude
3.3. Measures Measures for for the flattening flattening The flattening of the reference ellipsoid is described by a variety of measures: 1. The flattenin flatteningg f =
a
− b; a
a2
2. The first eccentricit eccentricity y (square) e2 =
− b2 ;
a2 2
3. The second eccentricit eccentricity y (square) (square) e =
a2
− b2 .
b2
3.4. Relationsh Relationships ips between between different different types of latitude From figure 3.1 can be seen that PR b = . QR a In the triangles ORQ and ORP we have tan β = so
QR PR ja tan φ = ; OR OR
tan φ PR b = = , tan β QR a
i.e., tan β =
a tan φ. b
In the triangles RV Q, RV P : P :
−
π tan 2
16
−
QR π β = ja tan VR 2
ϕ =
PR , VR
3.5. Co-ordinates Co-ordinates in the meridional ellipse ellipse i.e., cot ϕ PR b = = , cot β QR a eli tan ϕ =
a tan β. b
By combining still a2 tan ϕ = 2 tan φ. b
3.5. Co-ordinat Co-ordinates es in the meridional meridional ellipse 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 , 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 N cos ϕ.
(3.1)
Also
b2 P R = OR tan φ = OR 2 tan ϕ = N cos N cos ϕ a using OR = x = N cos N cos ϕ. The end result is because y = P R:
· 1 − e2
tan ϕ
−
y = N 1
e2 sin ϕ.
(3.2)
The equations 3.1 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 y2 + 2 a2 b
= = =
N 2 1 e2 N 2 2 cos ϕ + sin2 ϕ = 2 2 a b 2 2 2 N N b 2 cos ϕ + sin2 ϕ = 2 2 2 a a a 2 N b2 2 cos ϕ + sin2 ϕ = 1; 2 2 a a
the latter condition yields N = using the definition e
2
a2
− b2 ≡ a2 .
− 1
a e2 sin2 ϕ
,
3.6. Three-dime Three-dimensiona nsionall rectangula rectangularr co-ordina co-ordinates tes on the reference reference ellipsoid The above formula formulass are easily easily generaliz generalized: ed: if x and y are co-ordinates within the meridian section, section, the rectangular rectangular co-ordinates co-ordinates are X = x cos λ = N cos N cos ϕ cos λ, Y = x sin λ = N cos N cos ϕ sin λ,
−
Z = y = N 1
e2 sin ϕ.
17
Chapter Chapter 3. Co-ordinates Co-ordinates on the reference ellipsoid ellipsoid T
. Q
.
P
ϕ
. 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 λ, Z =
− 2
N 1
e
(3.3)
+ h sin ϕ.
Here, h is the straight distance of the point from the surface of the ellipsoid (“ellipsoidal height”). This quantity quantity is interesting interesting becaus b ecausee satellite satellite positioning devices can be said to directly directly measure precisely this quantity (more precisely, they measure X , Y , Z and and compute from these h).
3.7. Computing Computing geographic geographic co-ordina co-ordinates tes from rectangula rectangularr ones This, the reverse problem from that of equations ( 3.3 3.3)) , isn’t quite simple to solve. Closed solutions solutions exist, exist, but are complicated. complicated. Of course an iterative solution solution based directly directly on equations equations (3.3 3.3)) is certainly possible and often used. Computing Computing the longitude longitude λ is extremely simple: tan λ =
Y . X
A possible stumbling block is identifying the correct quadrant for λ. ϕ and h are more complicated. See figure 3.2 3.2,, where point P ’s 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:
√ X 2Z + Y 2
tan φ =
and the geometric distance (radius) by the equation: 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 Meridian arc length (proof (proof directly directly from eqs. eqs. (3.3 3.3). ).)) No Now w that 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 , P , for which (geometrically obviously) the above ratio is the same as for P . P . Now we compute X Q = N cos N cos ϕQ cos λ, Y Q = N cos N cos ϕQ sin λ,
−
e2 sin ϕQ ,
Z Q = N 1 from which
2 + Y 2 + Z 2 X Q Q Q
OQ = and thus
P Q = OP
− OQ.
Additionally, in the little triangle T QP we may compute (ϕ ( ϕQ abbreviated abbreviated to ϕ): ∠T QP =
ϕ
−φ
and thus T P = P Q sin(φ sin(φ TQ =
− ϕ) , P Q cos(ϕ cos(ϕ − φ) .
Now h = P S and
≈ ϕQ − PT P O ,
(3.4)
≈ ϕQ − PT P O cos(ϕ cos(ϕ − φ) .
(3.5)
ϕP or perhaps a hair’s width more precise ϕP
≈ TQ
1
This This procedu procedure re is in practice practice fairly fairly precise. precise. If h = 8000m and ϕ = 45◦ , then T P 26m, 2 between the ϕP solutions (3.4 (3.4)) and (3.5 (3.5)) there is a difference of 0.1 mm , 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. error.3 .
≈
≈
3.8. 3.8. Meridi Meridian an arc length length The length of a meridian arc, a quantity needed, e.g., with map projections (UTM, GaussKruger) u ¨ ger) is computed by integration. Above we already defined the quantity N , N , the transverse radius of curvature . The other radius of curvature of the Earth surface is the meridional radius of curvature M . M . If it is is given given as a function of latitude ϕ, we compute an element of path ds as follows: ds = Mdϕ. 1
Or then, not. Instead of P O:n one should take M ( M (ϕ) + h, where M the meridional meridional radius of curvature. curvature. Linearly T P (1 − cos(φ cos(φ − ϕ)). 2 1 T P 3 . 2 OP 2
19
Chapter Chapter 3. Co-ordinates Co-ordinates on the reference ellipsoid ellipsoid Now we may calculate the length of a meridian arc as follows: ϕ0
s (ϕ0 ) =
ˆ
Mdϕ.
(3.6)
0
On the reference ellipsoid
− − − ˆ − M =
i.e.,
s (ϕ0 ) = a 1
e2
a 1
3/2 e2 sin2 ϕ
1
ϕ0
e2
1
,
e2 sin2 ϕ
3/2
−
dϕ.
0
Here, the last factor can be expanded into a series – because e2 sin2 ϕ 1 – and integrated termw termwise ise.. See the litera literature ture.. Of course course also a numer numeric ic approach approach is possibl possible, e, nowa nowada day y 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 system and geometric geometric parame parameters ters Nowadays, the overwhelmingly most used global geodetic reference system is the Geodetic Reference System 1980 , GRS80. GRS80. The paramete parameters rs defining defining it (e.g., (e.g., (Heikk ( Heikkinen, inen, 1981 1981)) )) are given in table 4.1.
Some of the parameters are geometric (a ( a), some are dynamic (J ( J 2 ,ω ). Othe Otherr geomet geometri ricc and dynamic parameters may be derived as follows ((Moritz, (( Moritz, 1992), 1992), following (Heiskanen (Heiskanen and Moritz, 1967,, eqs. 2-90, 2-92)): 1967 J 2 e2
e2 3
− ⇒
2 me = 1 15 q 0 2me e2 = 3J 2 + . 15q 15 q 0
Here(Heisk Here(Heiskanen anen and Moritz, Moritz, 1967, eq. 2-70)) ω 2 a2 b m= GM and (based on the definitions) be = ae with the aid of which
4 ω 2 a3 e 3 e = 3J 2 + . 15 GM 2q 0 Furthermore we know ((Heisk ((Heiskanen anen and Moritz, Moritz, 1967, eq. 2-58)) 2
2q 0 e = ja
1+
e (e) =
3 e 2
arctan e
(4.1)
− e3
√ 1 e− e2 .
Now we may compute e2 iteratively usin using g eq. (4.1 4.1)) , which computes q 0 e (e) anew in every step. The result is (table 4.2 4.2): ): Quantity
Value
Explanation
a GM J 2 ω
6378137 m 3986005 108 m3 s−2 108263 10−8 7292115 10−11 rads−1
semi-ma jor axis Earth mass G Dynamic form factor Angular rotation rate
·
·
·
×
Table 4.1.: GRS80 defining parameters
21
Chapter Chapter 4. Reference Reference systems systems Quantity
Value
Explanation
e2 2 e b 1/f
0.006694380 0.00669 43800229 022900 0.00 0.00673 673949 949677 677548 548 6356 635675 752. 2.31 3141 4140 40 m 298.257222101
Ensimm Ensimmainen ¨ainen eksentrisyys neli o¨on ¨on Toine oinen n ekse eksen ntrisy trisyys ys neli nelio¨on ¨on Lyh Ly hyt akse akseli lipu puol olik ikas as Ka¨anteinen ¨anteinen litistyssuhde
Table 4.2.: GRS80 derived parameters Quan Quanttity
Value (WGS84)
Remark
a 1/f b
6378137 m 298. 298.25 2572 7223 2356 5633 635675 635 6752. 2.314 314245 245 m
same diffe differe ren nt! diff. diff. 0.1 mm
Table 4.3.: WGS84 ellipsoidal parameters
− b2 b2 1 a−b b Here we used 1 − e = 1 − = 2 ja 1 − = 1 − = , i.e., 2 a a f a a 2
a2
1
2
−e
=
− ⇒ 1 f
1
2
1 =1 f
−
− 1
e2 .
Often we use (a, (a, f ) 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 4.3::
Most Most often often the differen difference, ce, a bit over over 0.1 mm at most, most, can be neglec neglected ted.. It is apparen apparently tly due to imprecise numerics.
4.2. Gravimetr Gravimetric ic parameters parameters Computing geodetic parameters is done as follows ((Heiskanen (( Heiskanen and Moritz, 1967, kaavat 2-73, 2-74)):
− − − − γ e =
γ p =
GM m e q 0 1 m ab 6 q 0 GM m e q 0 1 + , a2 3 q 0
,
where ((Heisk ((Heiskanen anen and Moritz, Moritz, 1967 1967,, kaava 2-67)): q 0 e = 3 1 +
1 e 2
1
1 arctan e e
The solution is again obtained iteratively, yielding (f ( f ∗
22
≡ γ p γ −e γ e ):
1.
4.3. Reference Reference frames frames Quantity
Value
Explanation
γ e γ p f ∗
9.7803 78032677 267715 15 ms−2 9.8321 83218636 863685 85 ms−2 0.00 0.00530 530244 244011 011229 229
Normal gravity at equator Normal gravity at poles “Gra “G ravi vity ty flatt flatten enin ing” g”
As a check, we may still compute Clairaut ’s equation in its precise form ((Heiskanen (( Heiskanen and Moritz, Moritz, 1967, eq. 2-75)): ω2b q f + f ∗ = 1 + e 0 . γ e 2q 0
A simple closed, beautiful formula for normal gravity on the reference ellipsoid is SomiglianaPizzetti’s formula: aγ e cos2 ϕ + bγ p sin2 ϕ γ = . a2 cos2 ϕ + b2 sin2 ϕ
4.3. Reference Reference frame framess Every twenty-four hours, the Earth rotates around its axis relative to the heavens at what is very nearly nearly a constant constant rate, about what it very nearly a fixed axis. The direction direction of this rotation axis will serve as the z axis of both the celestial celestial and the terrestrial terrestrial frame. In order to completely completely define the orientation of our reference frame, we then need to conventionally fix two longitudes : 1. On the celesti celestial al sphere: sphere: we take take for this this the vernal equinox , where the Sun crosses the equator equator S-N 2. On the Earth: the Internat Internation ional al Meridi Meridian an Confer Conferenc encee in Washing ashington ton DC, 1884 1884,, chose 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 axis alway alwayss lies lies in the plane plane of the Greenw Greenwic ich h meridi meridian. an. This This is also called a CT (Conventi (Conventional onal Terrestr Terrestrial) ial) System. Locations Locations on the Earth’s surface surface are (almos (almost) t) consta constant nt in this this kind kind of frame, frame, and can be publis published hed,, e.g., e.g., on maps. maps. Ho Howe weve ver, r, 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 celestial reference frame, frame, as the co-ordinates co-ordinates of the fixed stars stars are nearly nearly constan constantt in it and may be publis published hed.. Also Also the equatio equations ns of motion 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 Chapter 4. Reference Reference systems systems Z
Rotational motion
Greenwich Y
Y
X
Spring Equinox
X
sidereal time θ Figure 4.1.: Geocentric Geocentric reference frames
◦ the Z axis is directed along the rotation axis of the Earth, more precisely the Conventional
Internation International al Origin (CIO) , i.e., the average direction direction of the axis over the time span 19001905
◦ the X Z -plane -plane is parallel to the zero meridian as defined by “Greenwich”, more precisely
by: earlier earlier the BIH (Bureau (Bureau International International de l’Heure, l’Heure, International 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. 4.4. The orient orientati ation on of the Earth 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 systems
Celestial North Pole Autumnal Autumnal equinox True Earth orbit
Spring
Summer
Winter
Winter solstice
Summer 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 nowaday nowadayss monitored monitored routinely, and available after the fact from the International Earth Rotation Service over the Internet. Internet. The variation variation of these parameters parameters is geophysica geophysically lly well understood, e.g., for the Chandler wobble it is mainly the pressure of the Earth’s oceans and atmosphere that is responsible ((Gross, ((Gross, 2000)). 2000)).
4.5. Transformations ransformations between systems systems See the following diagram, which depicts only the rotations between the various systems: Φ, Λ Local Astronomical
⇐⇒
x p , y p Conventional Terrestrial
⇐⇒
θ0 Instantaneous Terrestrial
⇐⇒
Real Astronomical
Here, Φ, Φ, Λ are local astronomical co-ordinates (direction of the plumbline), while x p , y p 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 Chapter 4. Reference Reference systems systems
α Cyg (Deneb)
Lyra α Umi
α Lyr (Vega)
Ursa Minor Dragon
α Umi (Polaris)
Precession (26 000 yr) 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 1970
-0.5
0
0.5
1 -0.2
x
Figure 4.4.: Polar motion
26
0
0.2
0.4
0.6 y
0.8
4.5. Transformations between systems systems
Reference pole (CIO) Real (instantaneous) pole o
◦
90
Plumbline 90
◦
−φ
−φ
B
B
Polar motion (Exaggerated)
A
A
◦
90
−φ
C
Earth centre of mass C
Figure 4.5.: How to monitor polar motion using latitude observatories
27
Chapter 5 Using rotation matrices 5.1. 5.1. Genera Generall Always when we change the orientation of the axes of a co-ordinate system, we have, written in rectangular rectangular co-ordinates co-ordinates,, a multiplication with a rotation matrix . Let us investigate the matter in the(x, the( x, y ) plane, figure 5.1. The new co-ordinate xP = OU = OR cos α,
where OR = OS + SR = xP + P S tan S tan α = xP + yP tan α. By substituting xP = (xP + yP tan α)cos α = xP cos α + yP sin α.
In the same fashion yP = OT = OV cos α,
where OV = OQ where yP = (yP
− V Q = yP − P Q tan α = yP − xP tan α, − xP tan α)cos α = −xP sin α + yP cos α.
Summarizingly in a matrix equation:
x y
=
cos α sin α sin α cos α
−
x 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 co-ordinate is positive or negative: in the above case
y =
−
sin α cos α
x 0
=
− sin α · x ,
i.e., y < 0, i.e., the minus sign is indeed in the lower left corner of the matrix.
y
y
P
Q
x
V
U
T α O
x S
R
Figure 5.1.: Rotation in the plane
29
Chapter Chapter 5. Using rotation rotation matrices
SR z
z = z
S
z
R
y
y
y
z
β
y
x
y
α
x
x =x
x
Figure 5.2.: The associativity of rotations
5.2. Chaining Chaining matrices matrices in three dimensions dimensions In a three-dimensional co-ordinate system we may write a two-dimensional rotation matrix as follows: x cos α sin α 0 x y = sin α cos α 0 y , z 0 0 1 z
−
i.e., the z co-ordinate is copied as such z = 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 r = S r , r = Rr,
in which R= then (associativity):
−
cos α sin α 0 sin α cos α 0 0 0 1
, S =
1 0 0
0 0 cos β sin β sin β cos β
−
,
S (Rr) = (SR) SR ) r, r = S ( i.e., the matrices are multiplied with each other. Remember that RS = SR,
in other words, matrices and transformations are not commutative 1 ! See figure 5.2.
5.3. Orthogonal Orthogonal matrice matricess Rotation matrices are orthogonal , i.e. RRT = RT R = I ; 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.1)
5.3. Orthogonal Orthogonal matrices matrices E.g.,
cos α sin α sin α cos α
−
1
−
cos α sin α
=
− sin α cos α
cos( α) sin ( α) 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
Rij Rik =
i=1
1 if j = k . 0 if j = k
The columns of a rotation matrix are orthonormal, orthonormal, their norm (length) (length) is 1 and they are mutually orthogonal. This can be seen for the case of our example matrix:
cos2 α + sin2 α = 1,
cos α sin α + ( sin α) cos α = 0.
·
−
·
Often we encounter other orthogonal matrices: 1. The mirror mirror matrix for an axis, e.g.: e.g.: M 2
≡
1 0 0
0 0 1 0 0 1
−
,
which inverts the direction, or algebraic sign, of the y co-ordinate. 2. The axes interchange interchange matrix (permutation): (permutation):
≡ −
P 12 12
3. Inversion Inversion of all axes : X =
1 0 0
0 1 0 1 0 0 0 0 1
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 rotation matrices matrices it is +1. The determi determinan nantt of the X matrix is ( 1)n , with n the number of dimens dimension ionss (i.e., (i.e., 3 in the above above case). case). A determ determina inant nt 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., M 2 and P 12 12 , we obtain M 2 P 12 12 =
−
0 1 0 1 0 0 0 0 1
.
The determinant of this is +1. However, it is again a rotation matrix: R3 (90◦ ) =
−
cos90◦ sin90◦ 0 sin sin 90◦ cos90◦ 0 0 0 1
!
All orthogonal orthogonal transformatio transformations ns having having positive positive determinan determinants ts are rotations. 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 Chapter 5. Using rotation rotation matrices
Z z
z
T K
s y
T
x
y
ζ
Y O
A
x
K
X
Figure 5.3.: Local astronomical astronomical co-ordinates co-ordinates
5.4. Topocentric systems systems Also “local astronomical” astronomical”.. Note that, whereas whereas the geocentric system system is “uniqu “unique” e”,, 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 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 gon − ζ . ◦ Distance s
The transformation between a point P ’s P ’s topocentric spherical co-ordinates and rectangular coordinates is:
x y z
The inverse transformation:
=
T
s cos A sin ζ s sin A sin ζ s cos ζ
ζ = arc arctan tan
A = 2 arct rctan
. T
x2 + y 2 , z y
x+
x2 + y2
.
The latter formula is known as the half-angle formula and avoids the problem of finding the correct quadrant for A. The result result is in the inter interv val ( 180◦ , 180◦ ] and negative values may be incremented by 360◦ to make them positive.
−
32
5.5. From geocentric geocentric to topocentric topocentric and back
Astronomical Φ, Λ co-ordinates Plumbz line x
y
π Geoid profile
y
x
Z
z
x
−Φ −x Λ z
z
Y X Λ
y
y
x
Z
x
Y X Figure 5.4.: From From the geocentric geocentric to the topocentric system. system. The matrix R1 mentioned in the text was left out here
5.5. 5.5. From geocentr geocentric ic to topocentr topocentric ic and back Let (X (X , Y , Z ) be a geoce geocent ntri ricc co-o co-ord rdin inat atee syst system em and (x,y,z) x,y,z) a topocen topocentri tricc instrument instrument cocoordinate system (i.e., (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 RT 1 x, −1 as can be easily dervied by multiplying the first equation from the left by the matrix RT 1 = R1 , then by the matrix R2T , and then by the matrix R3T , and finally by moving X 0 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 Λ: R3 =
−
cos Λ + si sin Λ 0 sin Λ cos Λ 0 0 0 1
.
(5.2)
Seen from the direction of the z axis we see (figure), that x = x cos cos Λ + y sinΛ, sinΛ, y =
−x sin sin Λ + y cosΛ. 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 Chapter 5. Using rotation rotation matrices
y
y
x z
λ R2
x
turns the co-ordinate frame around the y axis, in such a way that the z axis points to the North North celest celestial ial pole instead instead of to the zenoth zenoth.. The rotation rotation angle needed needed for ◦ this is Φ 90 . sin Φ 0 cosΦ R2 = 0 1 0 . (5.3) + co 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 problem problemss with rotation matrices 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 (ϕ,λ,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). ϕ,λ,h).
5.6.1. 5.6.1. Geodetic main main problem 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 ransform the local A-topocentric -topocentric co-ordinates co-ordinates of B , (AAB , sAB , zAB ) to rectangular: xAB
≡ xAB yAB zAB
= sAB
cos AAB sin zAB sin AAB sin zAB cos zAB
.
2. Transform, ransform, using rotation rotation matrices, matrices, these rectangular rectangular co-ordinate co-ordinate differences differences into geocen2 tric : T XAB = RT 3 R2 xAB , in which R3 and R2 are already given, equations (5.2 ( 5.2)) ja (5.3 (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 inverse problems with rotation matrices matrices 3. Add to the result the geocentric co-ordinate co-ordinatess of point A: XB = XA + XAB ,
jossa XA =
− (N ( N (ϕA ) + hA )cos ϕA cos λA (N ( N (ϕA ) + hA )cos ϕA sin λA N ( N (ϕA ) 1 e2 + hA sin ϕA
.
4. Transform ransform the geocentric geocentric co-ordinates co-ordinates of B obtained back to ellipsoidal form (3.7 (3.7)) in the way depicted: (ϕB , λB , hB ) . XB
→
5.6.2. 5.6.2. Geodetic inverse inverse problem 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¨ etaisyys ¨aisyys sAB . 1. Transform ransform the ellipsoidal ellipsoidal co-ordinates co-ordinates of A and B into geocentric co-ordinates: XA , XB . 2. Compute Compute the relative 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 transform to spherical spherical co-ordinates co-ordinates by translating translating the formula xAB = sAB
cos AAB sin zAB sin AAB sin zAB cos zAB
,
with the familiar arc tangent and Pythagoras formulas (and using the half-angle formula to avoid quadrant problems): sAB = tan AAB =
tan zAB =
2 + y2 + z 2 , xAB AB AB
yAB = 2 arct arctan an xAB x 2 + y2 xAB AB
zAB
yAB
AB +
2 + y2 xAB AB
,
.
5.6.3. 5.6.3. Compari Comparison son with ellipsoidal ellipsoidal surface solution The solution thus obtained is, concerning the azimuths, very close to the one obtained by using the geodesic between 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 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 distanc distancee is of course course the straight straight distanc distancee in space, space, not the length length of the geodesic geodesic.. This This distance is substantial.
35
Chapter 6 Co-ordinate systems and transformations 6.1. Geocentric Geocentric terrestrial terrestrial systems systems Generally geocentric systems, like WGS84, are defined in the following way: 1. The origin of the co-ordinate co-ordinate system coincides coincides with the centre 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 “rotating along” along” systems systems are called called terrestrial . Also ECEF (Earth Centred, Earth Fixed).
6.2. Conventio Conventional nal Terrestria errestriall System System To this we must however make the following further restrictions : 1. As the directi direction on 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 motion.. The main periods periods are a year year and approx. approx. 435 days days (the so-call so-called ed “Chandler wobble”), the amplitude being a few tenths of a second of arc — corresponding on the Earth’s surface to a few metres. 2. Nowaday Nowadayss the zero meridian meridian plane no is longer based on observ observation by the Greenwic Greenwichin hin observatory observatory,, but on worldwide worldwide VLBI observations observations.. These are co-ordinated co-ordinated by the International 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 Centered, Earth Fixed) reference reference system, system, e.g., WGS84 WGS84 or ITRSxx ITRSxx (ITRS = International Terrestrial Reference System, xx year number; created by the IERS). Another name used is Conventional Terrestrial System (CTS).
6.3. Pola Polarr motion motion The direction of the Earth’s Earth’s rotation axis is slightly slightly varying 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 direction of Greenwic Greenwich, h, and perpendicular perpendicular to it in the West West direction, respectivel respectively y. The transformation formation between between the instantaneous instantaneous and conventi conventional onal terrestrial terrestrial references is done as follows: follows: XIT = RY (x p ) RX (y p ) 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 y p about the X axis, which changes only the Y and Z co-ordinates. The matrices are
37
Chapter Chapter 6. Co-ordinate Co-ordinate systems systems and transformatio transformations ns
Greenwich meridian
Z
North pole
P
h
.
ϕ
.
Ellipsoidal normal
Y
.
.
λ
ϕ . .
X
Figure 6.1.: The ellipsoid and geocentric co-ordinates
RY (x p ) =
cos xP 0 0 1 sin xP 0
− sin xP 0 cos xP
and RX (y p )
1 0 0
0 0 cos yP sin yP sin yP cos yP
−
.
Because the angles x p and y p are very small, order of magnitude second of arc, we may approximate sin x p x p and cos x p 1 (same for y p ), as well as xP yP 0, obtaining obtaining
≈
≈
RY (x p ) RX (y p )
≈
1 0 0 1 x p 0
−x p 0 1
≈
1 0 0
0 0 1 y p y p 1
−
=
1 0 x p
0 1 y p
−
−x p y p 1
.
6.4. The Instantaneous Instantaneous Terre Terrestial stial 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, conventional, the so-called so-called instantaneous terrestrial 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 of polar motion parameters : if they are xP , yP — x p points to Greenwich and y p to the East from the Greenwich meridian — we obtain (see above):
X Y Z
IT
= RY (x p ) RX (y p )
≈ X Y Z
CT
because the angles x p , y p are so extremely small.
1 0 x p
0 1 y p
−
−x p y p 1
X Y Z
,
CT
6.5. The quasi-inerti quasi-inertial al geocentric system system The quasi-inertial, also celestial, or real astronomical (RA) reference frame, drawn in blue in figure 4.1 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 quasi-inertial geocentric geocentric system 1. The origin of the co-ordinate co-ordinate system again coincides coincides with the Earth’s centre 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: But: The X axis is pointing to the vernal equinox. equinox. Such a reference reference system doesn’t rotate along with the solid Earth. It is (to good approximat approximation) ion) inertial . We also also refer refer to it as an equatorial co-ordinate co-ordinate system. system. In this system the direction direction co-ordinates co-ordinates are the astronomical astronomical right ascention ascention and declination declination α, δ . If the “celestial “celestial spherical co-ordinates” co-ordinates” α, δ are known, we may compute the unit direction vector as follows: X cos δ cos δ cos α Y = cos δ sin δ sin α . Z RA sin δ
Here we must consider, however, that, while δ is given in degrees, minutes and seconds, α is given in time units . They They must must first be conve converte rted d to degree degreess etc. etc. One hour correspon corresponds ds 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 half-angle formula formula for α, which is precise over the whole range and avoids the problem of identifying tifying the right right quadrant): quadrant): δ = arcs arcsin in Z α = 2 arc arctan
Y √ . X + X + X 2 + 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 places α, δ read from a celestial chart are not apparent. apparent. They refer to a certain certain point in time, e.g., 1950.0 or 2000.0. Obtaining Obtaining apparent 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 f rom the reference work “Apparent Places of Fundamental Fundamental Stars” precalculate precalculated d and tabulated tabulated according to date. We can transform between RA and IT (instantaneous terrestrial) co-ordinates as follows:
X Y Z
= RZ ( θ0 )
−
RA
X Y Z
=
IT
cos(θ cos(θ0 ) sin(θ sin(θ0 ) 0
here, θ0 on Greenwich Apparent Sidereal Time, or GAST.
1
− sin(θ sin(θ0 ) cos (θ (θ0 ) 0
0 0 1
X Y Z
.
IT
. . . from the centre centre of the Earth. Earth. The fixed stars stars are so far way, way, howe howeve ver, r, that that the location location of the observ observer er doesn’t matter.
39
Chapter 7 Reference systems and realizations 7.1. Old and new reference reference systems; systems; ED50 vs. WGS84/GRS8 GS84/GRS80 0 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 Reference ellipsoid ellipsoid used: Internatio International nal Ellipsoid Ellipsoid (Hayford) 1924 vs. GRS80 2. Realized Realized by terrestrial terrestrial measurements measurements vs. based on satellite (and space geodetic) data 3. Non-geocent Non-geocentric ric (100 m level) level) vs. geocentric geocentric (cm level). 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 . f .
◦ International ellipsoid 1924: a = 6378388 6378388 m, f = 1 : 297 ◦ GRS80: a = 6378137 6378137 m, f = 1 : 298. 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. 298.25722356 257223563. 3. The net result is that the semi-minor semi-minor axis (polar radius) of WGS84 WGS84 1 is longer by 0.1 mm compared to GRS80.
To complicate matters, as the basis of the ITRS family of co-ordinate systems, and their realizations ITRF, was chosen the GRS80 reference ellipsoid . The paramete parameters rs of this this system system were were already given in Tables 4.1 and 4.2.
7.2. 7.2. WGS84 GS84 and ITRS Both WGS84 and the International Terrestrial Reference System (ITRS) are realized by computing 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: ∆b = −
∆f ∆ (1/f ) 0.000001462 · (21384m) (a − b) = 1 (a − b) = (21384m) = 0.000105m 000105m.. f ( /f ) 298. 298.25722
41
Chapter Chapter 7. Reference Reference systems systems and realizations realizations
◦ The scale deriv derives es from the SI system system..
This This is realize realized d by using range mea measur sureme ement ntss by propagation propagation of electromagne electromagnetic tic waves. waves. The velocity velocity of these waves waves in vacuum is conven conven-−1 tion tional ally ly fixed fixed to 299 299 792 458 m s . Thus, range measurement becomes time measurement by atomic clock, which is very precise.
◦ Orientation :
origin originally ally the Conve Convent ntion ional al Inter Internat nationa ionall Origin Origin (CIO) (CIO) of the Earth axis, axis, i.e., the mean orientation over the years 1900-1905, and the direction of the Greenwich Meridian, Meridian, i.e., the plumbline plumbline of the Greenwich Greenwich transit circle. Currently Currently,, as the orientation orientation is maintained maintained by the International International Earth Rotation 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-ordinat Co-ordinate e system system realizatio realizations ns Internationally, somewhat varying terminology is in use concerning the realization of co-ordinate systems systems or datums .
◦ ISO: Co-ordinate reference system / co-ordinate system ◦ IERS: Reference system / reference frame ◦ Finnish: koordinaattijarjestelm ((Anon., 2008)) 2008)) ¨arjestelma¨ / koordinaatisto ((Anon., 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 Then also also,, a datum was defined, with one or more datum points points being kept fixed at their conventional conventional values. 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 Realization of WGS84 WGS84 Because “WGS84” is often referred to as the system in which satellite positioning derived coordinates are obtained, we shall elaborate a little on how this system has been actually realized over time. Our source is (Kumar ( Kumar and Reilly Reilly, 2006 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 variet variety y. The precision precision of this approximation approximation is clearly sub-metre, 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, ( Stevenson, 2008). 2008).
7.5. Realization Realizationss of ITRS/ETRS ITRS/ETRS systems 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 stands for “system” “system”,, meaning meaning the principles for creating creating a reference reference
42
7.6. The three-dimension three-dimensional al Helmert Helmert transformation transformation frame before b efore actual realization. realization. With every ITRS corresponds a number number of ITRF’s ITRF’s (“Frames” “Frames”), ), which are realizations , i.e., co-ordinate solutions for networks of ground stations computed from sets of actual measurements measurements.. Same for ETRS/ETRF, which which are the corresponding 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 Interferometry (VLBI) providing a strong orientation; satellite and lunar laser ranging (SLR, LLR) contributing to the right scale, and the French DORIS satellite system. Nowadays also GLONASS is used. Curren Currently tly the follow following ing realizat realization ionss exist exist for ITRS: ITRF88, ITRF88, 89, 90, 91, 92, 93, 94; 96, 97; ITRF2000, ITRF2000, ITRF2005 ITRF2005 and ITRF2008. ITRF2008. The definition of an ITRFyy is as follows (McCarthy, ( McCarthy, 1996): 1996):
◦ The mean rotation of the Earth’s crust in the reference frame will vanish globally (cf.
for ETRF: on the Eurasian plate). plate). Obviously Obviously then, then, co-ordinates 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 -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, -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-dime three-dimension nsional al Helmert Helmert transfo transformat rmation ion The form of the transformation, in the case of small rotation angles, is
(2)
X Y (2) Z (2)
= (1 + m)
T
−
1 ez ey
ez 1 ex
−
−ey ex 1
· X (1) Y (1) Z (1)
tx ty tz
+
,
(7.1)
where tx ty tz is the translation vector of the origin, m is the scale factor correction from unity, and (e (ex , ey , ez ) are the (small) rotation angles about the respective axes. axes. Together ogether we thus have seven parameters. The superscripts (1) and (2) refer to the old and new systems, respectively. Eq. (7.1 (7.1)) can be re-written and linearized as follows, using mex = mey = mez = 0, and replacing replacing
≈ − ≈ − T
the vector X (1) Y (1) Z (1) by approximate values m and the e angles are all assumed small.
(2)
X − X (1) Y (2) − Y (1) Z (2) − Z (1)
m ez ey m ez ey
ez m ex
−ey
ez m ex
−ey
− −
ex m ex m
X 0 Y 0 Z 0
(1)
X Y (1) Z (1) X 0 Y 0 Z 0
+
+
T
. This is allowed as
≈
tx ty tz
tx ty tz
.
43
Chapter Chapter 7. Reference Reference systems systems and realizations realizations An elaborate rearranging yields
(2) X i (2) Y i (2) Z i
− − −
(1) X i (1) Y i (1) Z i
=
X i0 Y i0 Z i0
Z i0
0 +Z i0 0 0 Y i +X i0
−
−
+Y i0 X i0
−
0
1 0 0 0 1 0 0 0 1
m ex ey ez tx ty tz
.
(7.2)
Here we have added for generality a point index i, i = 1, . . . , n. n. The numbe numberr of points points is then n, the number of “observations” (available co-ordinate differences) is 3n 3 n. The The full set set of thes thesee “observ “observation equations” then become b ecomess
(2)
X 1 (2) Y 1 (2) Z 1
− X 1(1) − Y 1(1) − Z 1(1) .. .
(2) X i (2) Y i (2) Z i
− − −
(1) X i (1) Y i (1) Z i
.. .
X n(2) Y n(2) Z n(2)
− X (1) − Y n(1) − Z n(1)
=
X 10 0 Z 10 Y 10 +Z 10 0 0 0 Z 1 Y 1 +X 10
−
−
.. . X i0 Y i0 Z i0 .. .
.. .
−
.. .
.. .
Z i0
0 +Z i0 0 0 Y i +X i0
−
−
+Y 10 1 0 0 X 10 0 1 0 0 0 0 1
.. .
+Y i0 X i0
−
.. .
X n0 0 Z n0 Y n0 +Z i0 0 0 0 Z n Y n +X n0
−
−
.. .
.. .
.. .
0
1 0 0 0 1 0 0 0 1
.. .
.. .
.. .
.. .
+Y n0 1 0 0 X n0 0 1 0 0 0 0 1
−
m ex ey ez tx ty tz
.
(7.3)
This is a set of observation equations of form + v = Ax (but without the residuals vector v which is needed to make the equality true in the presence of observational uncertainty). There are seven unknowns x on the right right.. They They can be solve solved d in the leastleast-squ square aress sense sense if we have have co-ordinates 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 points and one co-ordi co-ordinat natee from a third point would suffice. However, it is always good to have redundancy.
7.7. Transformations ransformations between ITRF realizations 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:
± Rate
± 44
T 1 T 1 mm
T 2 T 2 mm
T 3 T 3 mm
D ppb
R1 R2 R3 0.001” 0.001” 0.001”
-0.5 0.2
-0.9 0.2
-4.7 0 .2
0.94 0.03
0.000 0.008
0.000 0.008
0.000 0.008
0.3 0.2
0.0 0.2
0 .0 0 .2
0.00 0.03
0.000 0.008
0.000 0.008
0.000 0.008
7.7. Transformations between ITRF realizations These parameters2 are to be used as follows:
X Y Z
(t) =
ITRF 2005 2005
− − − − − − − − · · − − − −−
+ (t
+
=
+
1+
D R3 R2
d t0 ) dt
R3 D R1
D R3 R2
R2 R1 D
R3 D R1
T 1 T 2 T 3
+ (t (t
1+
0.94 0 0 0 0.94 0 0 0 0.94
d t0 ) dt
T 1 T 2 T 3
X Y Z
R2 R1 D
(t) +
ITRF 2008 2008
X Y Z
(t) +
ITRF 2008 2008
=
10
9
−
X Y Z
(t) + ITRF 2008 2008
0.5 + 0. 0 .3 (t 2005 2005..0) 0.9 4.7
d with the numbers given, given, and forgetting about the uncertaint uncertainties. ies. Here refers to the rates, of dt which only that of T 1 is non-vanishi non-vanishing ng in this example.
2
Note the change in parameter names compared to the previous. For ( tx , ty , tz ) we now have (T ( T 1 1, T 2 T 2, T 3); T 3); m is now called D; and (e ( ex , ey , ez ) became (R (R1, R2, R3).
45
Chapter 8 Celestial co-ordinate systems 8.1. 8.1. Sidere Sidereal al 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. known. GAST can be computed computed to one second precision 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
January February March
θm 6 40 8 40 10 40
θm
Month
April May June
12 40 14 40 16 40
θm
Month
July August September
18 40 20 40 22 40
Month
October Novemb er Decemb er
θm 0 40 2 40 4 40
In constructing constructing this table, the following following knowledge knowledge was used: on March 21 at 12 UTC in GreenGreenh wic wich, the the hour hour angle angle of the the Sun, Sun, i.e. i.e.,, of the the verna vernall equi equino nox, x, i.e. i.e.,, side sidere real al time time,, is 0 . This Greenwich sidereal time θ0 consists of two parts : an annual part θa , and a clock time τ (UTC or Green Greenwic wich h Mea Mean n Time). Time). So we obtain obtain the annual annual part of sidere sidereal al time by subtract subtracting ing:: h h h θa = θ0 τ = 12 , i.e., 12 after adding a full turn, 24 .
−
−
If on March 21, or more precisely at midnight after March 20, sidereal time is 12 h 00m , then sidereal time for March 0 is 12 h 00m 4 20m = 10h 40m ; Remember that one day corresponds to about four minutes.
− ×
We ma may y fill out the table table using using the rule that 1 month month 2t . (In princip principle le we could slight slightly ly 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 Chapter 8. Celestial Celestial co-ordinate co-ordinate systems
o
180 W
o
o W 1 5 0
W 1 6 5
16 5 5 o E
hGr
1 5 o 0 E
o
30 N
o W W
5 1 3
Right ascension α 1 3 5 o o E
h o
45 N
W o W
0 1 2
1 2 0 o E
o
60 N
o W
0 5 1
1 0 5 o E
o
75 N
9 0 o W
E
o
0 9
o E 5 7
7 5 o W
vernal equinox True Mean
E o E 0 6
Greenwich
6 0 o o W
o E E 5 4
4 5 o o W 3 0 o W
λ 15 o W
LMST LAST GMST GAST
o
0
o
1 5 E
o E 3 0
Earth rotation
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 the day number within a month θa = θm + θd annual part of sidereal time τ time (UTC) Λ the longitude longitude of the Earth station converted converted to hours, minute and seconds (15◦ = 1h , 1◦ = 4m , 1 = 4s ).
See the pretty figure 8.1 8.1.. We have the following quantities:
◦ GAST = Greenwich Apparent Sidereal Time (= θ0) 48
8.2. Trigonometry rigonometry on the celestial celestial sphere sphere
apparent rotation of the celestial sphere Zenith North pole h
Meridian Equator
δ
Apparent diurnal motion
East North
t
Observer
South
West
Figure 8.2.: Hour angle and other co-ordinates on the celestial sphere
◦ LMST = Local Mean Sidereal Time (= θ ) ◦ the equinox equinox varies irregularly irregularly with time due to precession and nutation. nutation.
That’s why we distinguish Mean and Apparent. The difference is called the equation of equinoxes, ee. 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. Trigonomet rigonometry ry on the celestial celestial sphere On the celestial celestial sphere we have have at least two two different kinds of co-ordinates: co-ordinates: local and equatorial. equatorial. Local spherical co-ordinates are related to the local astronomical rectangular system as follows (x to the North, z to the zenith):
x y z
=
cos A sin ζ sin A sin ζ cos ζ
=
cos A cos η sin A 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 co-ordinates (“place (“places” s”)) of stars are nearly constant. constant. The disadvantag disadvantagee is that there is no simple simple relationship relationship to local astronomical astronomical co-ordinates. co-ordinates.
49
Chapter Chapter 8. Celestial Celestial co-ordinate co-ordinate systems
Celestial pole West
90
◦
$h$
− δ
East
◦
90
Zenith
ζ
−Φ
A
Star Figure 8.3.: Fundamental triangle of astronomy
An intermed intermediate iate form of co-ordinates is: hour angle and declination, h, δ . The The hour hour angle 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 instrumen instrument: t: if, of the three quantities 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 “One man’s noise is another man’s signal ” signal ”. On the celestial sphere there is a fundamental fundamental triangle triangle of astronomy astronomy : it consists of the star, the celestial celestial North pole, and the zenith. Of the angles of the triangle we mention mention t (North pole) and ◦ ◦ A (the zenith), of its sides, 90 Φ (pole-zenith), 90 δ (star-pole) (star-pole) and ζ (star-zenith).
−
−
The sine rule:
− sin A = sin h . cos δ
sin ζ
The cosine rule: cos ζ = sin δ sin δ sin Φ + cos cos δ cosΦcos δ cosΦcos h, sin δ = sinΦ cos z + cos cos Φ sin sin z cos A. We compute first, using the cosine rule, either δ or z , and then, using the sine rule, either h or A. Thus we obtain either (A, ( A, ζ ) (h, δ ) in both directions.
↔
8.3. Using Using rotation matrice matricess The various transformations between celestial co-ordinate systems can be derived also in rectangular co-ordinates co-ordinates by using rotation rotation matrices. matrices.
50
8.3. Using rotation rotation matrices
T
Let the topocentric topocentric co-ordinate co-ordinate vector vector (length 1) be r = x y z this is x cos A sin ζ y = sin A sin ζ . r= z cos ζ
. In spherical co-ordinates
The same vector we may also write in local astronomical co-ordinates:
r =
x y z
The transformation between them is:
=
cos α cos δ sin α cos δ sin δ
.
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: M1 =
R2 =
R3 = Let us compute
− −
1 0 0 0 1 0 0 0 1
cos θ sin θ 0
− −
− sin θ
cos α cos δ sin α cos δ sin δ
0 0 1
cos θ 0
−
cos θ sin θ 0
=
,
cos cos (90 (90◦ Φ) 0 sin (9 (90◦ 0 1 0 ◦ sin(90 Φ) 0 cos ( 90 90◦
R3 R2 M1 =
After this:
−
− − =
− Φ) − Φ)
=
sin Φ 0 0 1 cos Φ 0
− cosΦ 0 sin Φ
,
.
− sin θ cos θ 0
cos θ sinΦ sin θ sin Φ cos Φ
cos θ sinΦ sin θ sin Φ cos Φ
0 0 1
−
sin Φ 0 cos Φ 0 1 0 cos Φ 0 sin Φ
− sin θ cos θ 0
− sin θ cos θ 0
cos θ cosΦ sin θ cosΦ sin Φ
cos θ cosΦ sin θ cosΦ sin Φ
where we identify immediately
=
.
cos A sin ζ sin A sin ζ cos ζ
,
cos α cos δ sin α cos δ sin δ
.
sin δ = sin ϕ cos ζ + cos cos Φ sin sin ζ cos ζ cos A, the cosine rule in the triangle star-celestial pole-zenith. The inverse transformation is, based on orthogonality (the transpose!)
cos A sin ζ sin A sin ζ cos ζ
− − =
cos θ sinΦ sin θ cos θ cos Φ
− sin θ sinΦ cos θ sin θ cos Φ
cosΦ 0 sin Φ
With With these, these, we can do the transform transformatio ation n of spheri spherical cal co-ordinat co-ordinates es with with the aid of threethreedimensional “direction cosines”. cosines”.
51
Chapter Chapter 8. Celestial Celestial co-ordinate co-ordinate systems
8.4. On satellite satellite orbits 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. easily. The heights heights of their orbits vary vary from 400 km to over over 1000 km; the inclination 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 know the satellite’s satellite’s equator crossing, crossing, i.e., the time of the “ascending “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. 8.5. Crossi Crossing ng a given given latitude latitude in the inertia inertiall sys system tem If the target target latitud latitudee ϕ is given, given, we can com comput putee the distance distance ν from the ascend ascending ing node (“downrange-angle”) in angular units as follows: sin ϕ , sin i where i is the inclination of the satellite satellite orbit. From this follows again the elapsed time using Kepler’s third law; the period or time of completing one orbit is: sin ν =
P =
4π 2 3 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 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 spherical triangle: triangle: sin∆α sin∆α sin ϕ tan ϕ = sin∆α 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 satellite’s topocentric co-ordinates co-ordinates
θ
r
R
A
Φ
ϕ
ν i
Kev¨ Kevat¨attasaus
∆α
Ω
Figure 8.4.: The places of satellite and Earth station in the inertial system
8.6. The satellite’s satellite’s topocentric co-ordina co-ordinates tes Generally, the satellite orbit’s height h is given; according to its definition, the height of the satellite from the geocentre is r = r = ae + h,
where ae is the equatorial radius radius of the Earth. In spherical spherical approximation approximation ae = R. Then the rectangular, inertial co-ordinates of the satellite are r=
r cos ϕ cos αG r cos ϕ sin αG r sin ϕ
,
where φ is the satellite’s satellite’s geocentric geocentric latitude (or, if one wants wants to put it this way way, the “geocentric “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 sidereal time, which thus represent representss the orientation of the local meridian in (inertial) (inertial) space. space. Now the Earth station’s rectangular, inertial co-ordinates are R=
Subtraction yields: d=r
−R
R cosΦcos θ R cosΦsin θ R sinΦ
≡
.
d cos δ T T cos αT d cos δ T T sin αT d sin δ T T
from which we may solve the topocentric co-ordinates : tan δ T T =
,
d3
d2 tan αT = , d1 d21 + d22
53
Chapter Chapter 8. Celestial Celestial co-ordinate 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 : distance : 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 Crossing a given latitude in the terrestria terrestriall 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. 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
R=R
and the satellite satellite co-ordinates co-ordinates r=r
cosΦcosΛ cosΦsinΛ sinΦ
cos ϕ cos λ cos ϕ sin λ sin ϕ
.
From this again we obtain the terrestial topocentric vector: d=r
−R
≡
d cos δ T T cos(Λ d cos δ T sin (Λ T sin d sin δ T T
− tT ) − tT )
,
from which we may solve the topocentric declination δ T T and hour angle tT . 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 = θ
− tT = (θ0 + Λ) + tT = θ0 + (Λ − tT .) . In this, θ0 is Greenwich sidereal time, and (Λ − tT ) the ”hour angle of Greenwich”. 8.8. Determini Determining ng the orbit orbit from observations 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): e1 =
54
d1
d1
=
cos δ 1 cos α1 cos δ 1 sin α1 sin δ 1
, e2 =
d2
d2
=
cos δ 2 cos α2 cos δ 2 sin α2 sin δ 2
.
8.8. Determining Determining the orbit from observations observations When also the ground station vector R has been computed, we may compute d1 , d2 by the cosine rule if a suitable value 1 satelliitin korkeudelle h – tai vastaavasti, satelliitin radan s¨ s ateelle ¨ateelle r = R + h – has been given: r2
=
⇒ ⇒
R2 + d2 + 2Rd 2Rd cos(∠e, R) 2
2
d + 2Rd 2Rd sin η + R
d1,2 =
2
⇒
r =0
− ⇒ ± − − ± − − 4R2 sin2 η 2
−2R sin η
4 (R2
r2 )
=
−R sin η r2 R2 1 sin2 η , in which sin η = cos cos (∠e, R) = e · eU is the projection of the satellite’s direction vector onto =
the local vertical, and
eU =
R
R
cosΦcos θ cosΦsin θ sinΦ
=
,
the “zenith “zenith vector” of the observer, a unit vector pointing pointing straight upward. upward. The value value sin η , the sine of the elevation angle, may be calculated directly as the dot product e eU when both vectors are given in rectangular form.
·
Only the positive solution makes physical sense: d=
− − − r2
sin2 η
R2 1
R sin η.
Thus we obtain d1 , d2 and thus d1 = d1 e1 , d2 = d2 e2 . For the satellite velocity we obtain v=
d2 − d1 . τ 2 − τ 1
We know however what the velocity for a circular orbit should be at height h: accor accordi ding ng to Kepler’s third law
4π 2 P = (R + h)3 GM Let us prepare the following little table:
⇒
2π (R + h) vk = = P
GM . R+h
Height (km)
500
750
1000
1500
Spee Speed d (m/s (m/s))
7612 7612.6 .609 09
7477 7477.9 .921 21
7350 7350.1 .139 39
7113 7113.0 .071 71
We see that the satellite’s linear speed of flight decreases only slowly with height. Therefore we may use the “observed “observed velocity” velocity” v for correcting the height height h according to the following formula: h = h
vk , v
where vk is the velocity according to Kepler (from the table for the value h), v the calculated velocity, and h the improved improved value for the satellite satellite height. height. This process converges 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 “ ov “ overfligh erflightt height”! 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 equivalen equivalently: tly: the semi-major 1
An experienced observer can guess a satellite’s height from the observed motion in the sky surprisingly precisely.
55
Chapter Chapter 8. Celestial Celestial co-ordinate co-ordinate systems 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 = d2 − d1 r1 cos i,
from which which i may be solved. And if i and some some “footpoin “footpoint” t” (ϕ, λ) of the satellite is known, also Ω and (with P ) 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 node 2 :
Ω (τ ) τ ) = Ω (τ 0 ) = Ω (τ 0 )
− (τ − τ 0) ddτ Ω = − (τ −
·
3 τ 0 ) 2
GM ae a3 a
2
J 2 cos i,
in which J 2 is the dynamic flattening of the Earth , value J 2 = 1082. 1082.62 10−6 . One One of the first first 3 achievements of the satellite era was the precise determination of J 2 . 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◦ 56◦ 65◦ 74◦ 81◦ 90◦
-7.651 -4.278 -3.233 -2.109 -1.197 0
-6.752 -3.776 -2.854 -1.861 -1.056 0
-5.985 -3.347 -2.529 -1.650 -0.936 0
-4.758 -2.661 -2.011 -1.311 -0.744 0
0.9856
97◦ .401
98◦ .394
99◦ .478
101◦.955
In the table table the the last last row is speci special al.. The The value alue 0 .9856 degree/day is the Sun’s apparent angular velocity relati relative ve to the backgro background und of the stars. stars. If the precessi precessional onal velocit velocity y of the satellit satellite’s e’s orbital plane is set to this value, the satellite will always fly over the same area at the same local solar time. time. In this this way way we achiev achievee a so-cal so-called led heliostati heliostationa onary ry orbit, orbit, also also known known as “no“noshadow” orbit, the advantages of which are continuous light on the solar cells, and in remote sensing, sensing, always always the same elevation elevation angle of the Sun when imaging imaging the Earth surface. Given 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 360/q = 68 days days ≈ 2.2 months months.. This is the time span after which the satellite appears again, e.g., in the evening sky of the same latitude. 2 3
The formula applies for circular orbits. In fact, this motion of the ascending node is so rapid, that without considering it it is not possible to generate usable orbit predictions.
56
Chapter 9 The surface theory of Gauss Carl Friedrich Gauss1 was among the first to develop develop the theory of curved curved surfaces. 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 surface of the Earth: note that Gauss was a geodesist geodesist who measured 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, (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 + y e2 + z e3 .
We will now often choose to write x by its representation on this basis: x=
x1 x2 x3
=
x y z
.
These three parameters thus form a parametrization of R3 .
9.1. 9.1. A curve curve in spac space e A curve running through space C can be parametrized by a parameter s. Then, Then, the points 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 numbers on a measurement measurement tape form a parametrizat parametrization ion by distance. 2. If you drive along a road, the trip meter readings form a parametrization of the road according according to distance. distance. 3. Along Mannerheimi Mannerheimintie ntie,, 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 geometry.. Others were J´ anos Bolyai and Nikolai I. Lobachevsky. anos
57
Chapter Chapter 9. The surface 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: t (s) =
dx (s) = xs ; ds
the length of the tangent is
t =
dx ds
2
dy ds
+
2
2
dz ds
+
=
dx2 + dy2 + dz 2 = ds2
√
1 = 1.
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
t
ˆ ds ˆ dτ = dτ ˆ dx
s(t) =
0
0
t
=
0
2
+
dτ
dx2 + dy2 + dz 2 dτ = dτ
dy dτ
2
dz dτ
+
(9.1)
2
dτ,
i.e., in differential form
dx (t) dt
ds = dt
2
+
dy (t) dt
2
dz (t) dt
+
2
,
from which s (t) can always be computed by integration 9.1. Another differentiation yields the curvature vector : dt (s) d2 = 2 x (s) . k (s) = ds ds
9.2. The first fundamental fundamental form form (metric) (metric) The Gauss first fundamental form : ds2 = Edu 2 + 2Fdudv 2Fdudv + Gdv 2 .
(9.2)
ds is a path element within within the surface. 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 co-ordinates:: ∂ x = xu = ∂u
∂x/∂u ∂y/∂u ∂z/∂u
∂ x , xv = = ∂v
∂x/∂v ∂y/∂v ∂z/∂v
.
These vectors are called the tangent vectors of surface S and parametrization (u, (u, v ).
58
9.2. The first fundamental fundamental form (metric) (metric)
v + dv
xv
v
G/
1 2
E /
1 2
xu
u + du
x
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 fundamental form : E =
xu · xu , xu · xv , xv · xv .
F = G =
In figure 9.1 we have depicted depicted the Gauss first fundamental fundamental form and the tangent tangent vectors vectors xu , xv . It is easy to show (chain rule in three dimensions ( x,y,z)), x,y,z)), that ds2
= = +2 + =
dx2 + dy 2 + dz 2 =
∂x ∂u
2
+
∂x ∂x ∂u ∂v ∂x ∂v
∂y ∂u
+
2
+
2
+
∂y ∂y ∂u ∂v
∂y ∂v
∂z ∂u
2
+
+
∂z ∂v
2
du2 +
∂z ∂z ∂u ∂v 2
dudv +
dv 2 =
xu · xu du2 + 2 xu · xv dudv + xv · xv dv2,
from which (9.2 (9.2)) follows directly. We see, e.g., that in the direction of the v curves (dv (dv = 0): ds2 = Edu 2 , in other other words words,, E represent representss the metric distance distance between between two two successiv successivee curves curves ( u, u + 1). Similarly G represents the distance between two successive v curve curves. s. The closer closer the curves curves are to each other, the larger xu or xv and also the larger E or G. F again represents represents the angle between the u and v curves: it vanishes if the angle is straight.
59
Chapter Chapter 9. The surface surface theory of Gauss
9.3. The second second fundamental fundamental form 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 n=
Apparently
nx ny nz
.
n · xu = n · xv = 0.
(9.3)
Let us also require
n = 1, or nx2 + ny2 + n2z = 1. We differentiate x a second time: xuu
=
xuv
=
xvv
=
∂ 2 x , ∂u 2 ∂ x , ∂u∂v ∂ x . ∂v 2
Now, Gauss’s second fundamental form is edu2 + 2fdudv 2fdudv + gdv 2 , where e = f = g =
n · xuu , n · xuv , n · xvv .
(9.4)
Based on condition (9.3 (9.3)) we also have ∂ ∂u ∂ 0= ∂v ∂ 0= ∂u ∂ 0= ∂v
0=
n · xu = nu · xu + n · xuu n · xu = nv · xu + n · xuv n · xv = nu · xv + n · xuv n · xv = nv · xv + n · xvv
⇒ ⇒ ⇒ ⇒
e=
−nu · xu, f = −nu · xv , f = −nv · xu , g = −nv · xv .
(9.5)
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 x (u + du,v) du,v), the normal changes n n = n + dn’. In the same way, when moving from point x to x (u, v + dv), dv), the normal changes n n = n + dn .
→
As a formula: dn =
60
∂ n ∂ n du + dv = nu du + nv dv. ∂u ∂v
→
9.4. Principal Principal curvatures curvatures
−gdv
−edu
−f du
n + nv dv
−f dv
n + nu du n xv
xu
v + dv u + du
x
O
u
v
Figure 9.2.: Explaining Explaining geometricall geometrically y the Gauss second fundamental form
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 direction of tangent tangent vector vector xv . Let us separate them by projection :
dn · xu dn · xv
= =
nu · xu du + nv · xu dv, nu · xv du + nv · xv dv,
in which we may directly identify the elements of the second fundamental form e,f,g with the aid of (9.5 (9.5): ):
− dn · xu − dn · xv
= edu + fdv, = f du + gdv.
Contrary to the first fundamental form, the second fundamental form has no corresponding object in the Riemann Riemann surface surface theory. theory. It exists exists only only for surfaces surfaces embedded in a surrounding (Euclidean) space. Often Often we can also find a “tenso “tensorial rial”” notatio notation: n: β 11 11 = e, β 12 12 = β 21 21 = f, β 22 22 = g.
9.4. Principal Principal curvature curvaturess 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, ( u, v ) is 1. an arbitrary arbitrary choice, choice, and 2. not metrically metrically scaled. scaled.
61
Chapter Chapter 9. The surface surface theory of Gauss 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 corresponds to. 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 perpendicular to each other, makes this problem even trickier. Write the first and second fundamental forms in matrix form: H =
E F F G
e f f g
,B=
.
Form the matrix:
C
G F
≡ H 1B = − −
=
−F E
1
EG
−
F 2
1 EG
Ge Ef
− F 2
−Ff −Fe
− e f f g
=.
Gf F g E g 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 which describes describes the length of a distance element , performs a metric scaling of the B matrix 3 .
Principal curvatures: The matrix C has two eigenvalues 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, ( u, v ). The The corcorT 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 = is the Gaussian curvature . det H EG F 2 1 1 1 eG + gE 2f F 2. The half-su half-sum m (κ1 + κ2 ) = (C 11 is the mean curvature . 11 + C 22 22 ) = 2 2 2 EG F 2
− −
1. The product κ1 κ2 = det C =
− −
Principal directions of curvature: We may also study the eigenvectors of C , which are called principal directions of curvature . For this, write 0 = H ( H (C
− κ1I ) x1 = (B − κ1H ) x1, 0 = H ( H (C − κ2 I ) x2 = (B − κ2 H ) x2 . 2 3
Julius Weingarten, German mathematician mathematician 1836-1910 1 A more technical description: B is a covariant tensor β ij again ij , and H the covariant metric tensor gij . H corresponds corresponds to the contrav contravarian ariantt tensor gij . C is now now the “mixe “mixed d tensor tensor” ” β ki = gij β jk whose tensorial jk , whose eigenvalue problem is β ji − κδ ji xj = 0, −
the same equation as (9.6 ( 9.6). ). See chapter 10.
62
9.4. Principal Principal curvatures curvatures T Multiply the first from the left with x T transposing: 2 , and the second with x 1 and transposing:
0 =xT 2 B x1 0=
− κ1xT 2 H x1, T T T xT = xT 1 B x2 − κ2 x1 H x2 2 B x1 − κ2 x2 H x1 ,
as, both B and H being symmetric symmetric matrices, matrices, we have
Subtraction of the two yields
xT 1 B x2
T
= xT 2 B x1 and
xT 1 H x2
T
= xT 2 H x1 .
(κ2
− κ1) xT 2 H x1 = 0.
This shows that if the principal radii of curvature are different 4 , then the expression x T x 1 2 H x 5 vanishes. vanishes. This expression can be interpreted as an inner product : product : in fact, in Cartesia Cartesian n plane T co-ordinates in the tangent plane, the matrix H = I , and we have x 2 x1 = 0, or x1 x2 in the Euclidean sense.
⊥
The principal directions directions of curvature curvature are mutually mutually perpendi p erpendicular. cular. This is just a special case of self-adjoint (symmetric) operators having mutually orthogonal eigenvec eigenvectors, tors, e.g., the eigenfuncti eigenfunctions ons of Sturm-Liouvi Sturm-Liouville lle 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 x=
From this
xϕ =
− − − − −
∂ x = ∂ϕ
N ( N (ϕ)cos ϕ cos λ N ( N (ϕ)cos ϕ sin λ N ( N (ϕ) 1 e2 sin ϕ
.
d cos λ (N ( N (ϕ)cos ϕ) dϕ d sin λ (N ( N (ϕ)cos ϕ) dϕ d 1 e2 (N ( N (ϕ)sin ϕ) dϕ
Using the formulas derived in the Appendix B we obtain: sin ϕ cos λ sin ϕ sin λ +cos ϕ
M (ϕ) xϕ = M (
xλ =
∂ x = N ( N (ϕ) ∂λ
;
.
cos ϕ sin λ +cos ϕ cos λ 0
.
The surface normal is obtained as the vectorial product, normalized: n=
xϕ × xλ , xϕ × xλ
4
And if they are not, any linear combination of x 1 and x2 will again be an eigenvector, and we can always choose two that are mutually perpendicular. perpendicular. 5 In index notation: gij xi2 xj1 .
63
Chapter Chapter 9. The surface surface theory of Gauss where = N M
xϕ × xλ
=
−
−N M
the norm of which is
− cos2 ϕ cos λ − cos2 ϕ sin λ sin ϕ cos ϕ cos2 λ − sin ϕ cos ϕ sin2 λ
cos2 ϕ cos λ cos2 ϕ sin λ sin ϕ cos ϕ
xϕ × xλ = N M cos M cos Thus n=
not a surprisi surprising ng result. result. . .
2
−
=
−N M cos M cos
2
ϕ
=
cos λ sin λ tan ϕ
,
ϕ 1 + tan2 ϕ = N M cos M cos ϕ. cos ϕ cos λ cos ϕ sin λ sin ϕ
,
Let us compute the first fundamental form: E = F = G =
xϕ · xϕ = M 2 sin2 ϕ sin2 λ + cos2 λ xϕ · xλ = 0, xλ · xλ = N 2 cos2 ϕ = p2.
+ cos2 ϕ = M 2 ,
We calculate for use in calculating the second fundamental form nϕ =
and thus (equations (equations 9.5 9.5))
+sin ϕ cos λ +sin ϕ sin λ cos ϕ
−
e = f = g = In other words: H = B =
E F F G e f f g
det(C det(C
64
M 2 0 2 0 N cos2 ϕ
=
M 0 0 N cos N cos2 ϕ
=
,
, and
1 M 0
1 N
− κI ) x = 0;0;
− κI ) = 0 ⇒
⇒ − −
As could could be expecte expected. d. . .
0 . 1 N are C ’s ’s characteristic values or eigenvalues , solutions of the eigen-
we obtain the values by solving
κ
; nλ =
(C
1 M
+cos ϕ sin λ cos ϕ cos λ 0
− nϕ · xϕ = +M, − nϕ · xλ = − nλ · xϕ = 0, − nλ · xλ = +N cos N cos2 ϕ.
C = H −1 B = The principal curvatures κ1,2 value problem
−
κ
=0
⇒
det
1 M
−κ
0
1 0 κ N 1 1 κ1 = , κ2 = . M N
−
=0
⇒
9.5. A curve embedded embedded in a surface
9.5. 9.5. A curve embedde embedded d in a surface surface If a curve C runs inside a curved surface, we may study other interesting things.
The tangent vector
t1 du/ds If we call t = = “the t vector’s components in the (u, ( u, v) co-ordinate system”, 2 dv/ds t 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. 6 i
The curvature vector k=
dt ds
d 1 t xu + t2 xv = ds dt1 dt2 = xu + xv + ds ds 2 + xuu t1 + 2xuv t1 t2 + xvv t2 =
2
.
In other words, the“curvature vector’s components in the (u, ( u, v ) co-ordinate co-ordinate system”ontain other 1 things besides just the derivatives of the component values dt /ds ja dt2 /ds. /ds. We write xuu
1 2 = Γ11 xu + Γ 11 xv + xuu n n,
xuv
=
xvv
=
· 1 2 Γ12 xu + Γ 12 xv + xuv · n n, 1 2 Γ22 xu + Γ 22 xv + xvv · n n;
(9.7)
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 10). ). They illustrate illustrate the reality reality that that differentiating a vector (also known as “parallel transport”, see Chapter 10 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 9.4)) , the elements of the second fundamental form e,f,g. e,f,g . We obtain k =
+ +
dt1 1 + Γ11 t1 ds dt2 2 + Γ11 t1 ds
e t1
2
2
1 1 2 1 + 2Γ12 t t + Γ22 t2
2
2
2 1 2 2 + 2Γ12 t t + Γ22 t2
2
+ 2f 2 f t 1 t2 + g t2
2
xu + xv +
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 exterior curvature curvature kext , the curvature of the curve curve “along with” with” the itself itself curved curved surface. surface. The internal curvature again has two components in “(u, “( u, v ) co-ordinates”, which can be read from the above equation. We write kint = k1 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
65
Chapter Chapter 9. The surface surface theory of Gauss where i
k =
1
k k2
=
1
dt /ds dt2 /ds
1 + Γ11 2 + Γ11 2
dt1 /ds +
1 2
1 1 2 t + 2Γ12 t t 2 2 1 2 t1 + 2Γ12 t t 2
1 i j Γij tt
i=1 j=1 j =1 2 2
=
dt2 /ds +
2 i j Γij tt
i=1 j=1 j =1
1 + Γ22 2 + Γ22
t2 t2
2 2
=
,
where we have have “economi “economized” zed” the formulas by using summation signs. The external curvature curvature again is kext = k n n,
·
where
k · n = e
t1
2
+ 2f 2 f t 1 t2 + g t2
2
.
Example: a latitude circle on the Earth surface.
− − − − − − · − N ( N (ϕ)cos ϕ cos λ N ( N (ϕ)cos ϕ sin λ N ( N (ϕ) 1 e2 sin ϕ
x=
dx dx dλ = = t= ds dλ ds
N cos N cos ϕ sin λ +N cos N cos ϕ cos λ 0
dt dt dλ = = k= ds dλ ds External curvature:
n=
so
;
sin λ +cos λ 0
1 = N cos N cos ϕ
cos λ sin λ 0
cos ϕ cos λ cos ϕ sin λ sin ϕ
;
1 . N cos N cos ϕ
,
1 n cos2 λ + sin2 λ 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. kext = k n n =
The internal curvature is
kint = k
− kext = N 1
− −
cos λ + cos ϕ cos λ cos ϕ sin λ + cos ϕ sin λ cos ϕ +sin ϕ
We obtain for the length of this vector
kint = tanN ϕ . 66
=
tan ϕ N
− −
cos λ sin ϕ sin λ sin ϕ +cos ϕ
.
9.6. The geodesic geodesic
p/sin ϕ
= N cot N cot ϕ kint
p ϕ
k
kext
N = p/cos ϕ
Figure 9.3.: The curvatures and radii of curvature of a latitude circle 1 1 = ; its internal part kint , N cos N cos ϕ p (ϕ) 1 tan ϕ 1 length k sin ϕ = sin ϕ = ; and its external part kext, length k cos ϕ = cos ϕ = p (ϕ) N ( N (ϕ) p (ϕ) 1 . In the figure figure one also also sees sees how how the distanc distancee from from the Earth’s Earth’s rotation rotation axis is in every every N ( N (ϕ) direction (k, kint ja kext ) the inverse of the curvature. This is also intuitively clear from rotational symmetry. In Figure 9.3 we see the curvature vector k, length k =
9.6. 9.6. The geodesi geodesicc Internally, in surface co-ordinates We obtain the formula 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): dti + ds
2
2
i j k Γ jk t t = 0.
(9.8)
j=1 j =1 k=1
This approach is developed further later on in chapter 10.5 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
dt = kext = e t1 ds Because
2
t · xu t · xv
1 2
+ 2f 2f t t + g t
2 2
n=
t
1
t
2
e f f g
t1 t2
n.
= t1 xu xu + t2 xu xv = Et 1 + F t2 , =
· · t1 xv · xu + t2 xv · xv = F t1 + Gt2 ,
67
Chapter Chapter 9. The surface surface theory of Gauss it follows that
If we define
we may write dt = ds Here
· · ≡ · · · · · · t1 t2
1
−
E F F G
=
t1 t2
t xu
E F F G
1
−
.
t xu t xv
E F F G
t xv
t xu t xv
1
−
e f f g
e f f g
E F F G
E F F G
1
−
t xu t xv
n.
1
−
= H −1 BH −1
following the earlier used notation7 .
Example: on the ellipsoidal surface of the Earth, we have 1
−
1
−
H BH
= =
and
· · · · · · · · · · − − − M −2 0 −2 0 N cos−2 ϕ M −3 0 −3 0 N cos−2 ϕ t xu t xv
M −2 0 −2 0 N cos−2 ϕ
M
N cos N cos2 ϕ
1/M 3 0 0 cos ϕ/p3
=
=
=
M t eN p t eE
,
from which we obtain dt ds
= =
t eN
t eE
1 t eN M
2
+
1/M 0 0 cos ϕ/p
1 t eE N
2
t eN t eE
n=
n.
Here, the unit vectors xu eN = xu
=
sin ϕ cos λ sin ϕ sin λ cos ϕ
xv , eE = xv
cos ϕ sin λ +cos ϕ cos λ 0
=
.
1 The approach is geometricall geometrically y intuitiv intuitive. e. The expression expression in wave wave brackets brackets t eN 2 + M 1 1 1 2 2 2 cos A + sin A is precisely the curvature of the surface in the direction t eE = 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 ( 2.1)). )).
Advantage: in surface surface co-ordinate co-ordinatess (ϕ, λ) there are inevitably two poles , singularities where the curvature curvature of the latitude circles goes to infinity and numerical numerical methods can be ill-behav ill-behaved. ed. This will not happen in rectangular space co-ordinates. Disadvantages: 7
kl il In index notation: g ij β jk jk g , see chapter 10. A logical notation for this would be β .
68
9.6. The geodesic geodesic 1. More equations, equations, more computational computational effort. 2. At every point, we must compute M and N and for this ϕ = arctan p =
X 2 + Y 2 .
Z , where (1 e2 ) p
−
3. The “roll-i “roll-in” n” and “roll-out” “roll-out” of of the computation of the geodesic geodesic presupposes presupposes the transfortransformation of (ϕ, (ϕ, λ) to (X,Y,Z (X,Y,Z ) co-ordinates in the starting point, and back in the end point.
69
Chapter 10 The surface theory of Riemann An importa important nt theoret theoretica icall frame frame in which which curve curved d surfac surfaces es and curve curvess are often often descri described, bed, is 1 Riemann’s surface theory. theory. In this theory, theory, we study the curved curved surface surface intrinsically intrinsically,, i.e., without taking taking into into acc account ount that the curve curved d surfac surfacee of the Earth Earth is “emb “embedded” dded” in a complet ompletee thre threedimensional (Euclidean) space. This makes the surface theory of Riemann useful also in situations, where this embedding higherdimensional space doesn’t necessarily even exist. For example, in General Relativity we describe spacetime (x,y,z,t (x,y,z,t)) as a curved manifold according to Riemann’s theory, the curvature parameters 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. 10.1. What What is a tenso tensor? r? In physics, physics, more than in geodesy, geodesy, we often encounter encounter the concept of “tensor” “tensor”.. What is a tensor? tensor?
Vectors First of all, a vector . Initially Initially we shall look at the situation situation in two-dime two-dimension nsional al Euclidean space, i.e., the plane . A pair of co-ordinate differences between two adjacent points is a vector: i
v=v =
∆x ∆y
(10.1)
is a vector, written with a superscript. In the other co-ordinate system we write this vector as
v = vi =
∆x ∆y
.
Note that here we have both a symbolic notation (v) and a component notation (v (v i , i = 1, 2). The components depend on the chosen co-ordinate system (∆x (∆ x = ∆x, ∆y = ∆y, although 2 2 always ∆x ∆x2 + ∆y ∆ y 2 = ∆x + ∆y = constant — an invariant ), ), but the symbolic notation does not depend on this. A vector is always always the same thing, even if its co-ordinates co-ordinates transform with the co-ordinate co-ordinate system. A vector is always the same “arrow “arrow in space”. space”.
There exists the following transformation formula for changing the components of a vector:
vi =
αii v i ,
(10.2)
i
1 2
Georg Friedrich Bernhard Riemann (1826 – 1866), 1866), German mathematician mathematician. Albert Einstein, 1879 – 1955, German-born theoretical physicist, discoverer of relativity theory, and icon of science.
71
Chapter Chapter 10. The surface surface theory of Riemann or in “ma “matri trix x languag language” e” v = Av,
where
v1 v2
i
v=v =
∆x ∆y
=
is the column matrix formed by the components, and A=
αii
=
α11 α21 α21 α22
=
cos θ sin θ sin θ cos θ
−
is the component matrix of the transformation operator. Every quantity that transforms according to equation (10.2 ( 10.2)) is called a vector . Familiar amiliar vector quant quantiti ities es are veloci velocity ty,, accele accelerati ration, on, force, force, . . . what what they they have have in com common mon 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 : index :
T i j =
αii α j j T ij .
i,j
In “matrix “matrix languag language” e” this this is3 T = AT AT , where A is the same as defined above. We have already encountered many tensors in geodetic theory: 1. The inertial inertial tensor tensor of the Earth. 2. The gravity gravity gradient gradient or Marussi tensor
M =
∂ g ∂x ∂ g ∂y ∂ g ∂z
=
∂ 2 W ∂x 2 ∂ 2 W ∂y∂x ∂ 2 W ∂z∂x
∂ 2 W ∂x∂y ∂ 2 W ∂y 2 ∂ 2 W ∂z∂y
∂ 2 W ∂x∂y ∂ 2 W ∂x∂z ∂ 2 W ∂z 2
σx2
.
σxy , 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.
3. The varia variance nce “ma “matri trix” x” is really really a tensor: tensor: Var (x) =
4. In the earlier discusse discussed d surface theory theory of Gauss, the first, H = B=
e f f g
E F , and the second, F G
fundamental form, as well as C = H −1 B ;
5. Also the elements elements of the map projectio pro jection n fundamental fundamental form (to be discussed discussed later) E, E, F , G E E F F consti constitut tutee a tensor tensor H H = , a me metr tric ic tenso tensorr of sort sorts. s. The The object object H −1 H H = F F G
E E M 2 F F Mp
3
F F Mp G p2
may again be called the scale tensor .
Verify that the matrix equation produces the same index summations as the index equation!
72
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 eigenvalues 4 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 ; therefore, to different eigenvalues λi , λ j , i = j belong eigenvectors xi , x j that are mutually mutually perpendicular, as proven in mathematics textbooks.
In an n-dimensional space, a tensor T has n independent invariants. invariants. The eigenv eigenvalue alues s λi , i = 1, . . . , n are of course invariants. So are their sum and product .
λi =
i
T ii ,
i
on the two-dimensional plane λ1 + λ2 = T 11 + T 22 , the sum of the diagonal elements or trace 5 ; ja
λi = det(T det(T )) ,
i
on the two-dimansional plane again λ1 λ2 = det det (T ) T ) = T 11 T 22
− T 12T 21,
the determinant of the tensor. 2 The trace of the variance matrix σx2 + σy2 we already know as the point variance σP ; it was chosen precisely because it is an invariant, independent of the direction of the xy axes. Also the trace trace of the gravity gradient tensor is known:
M ii ii =
i
∂ 2 W ∂ 2 W ∂ 2 W + + ∂x 2 ∂y 2 ∂z 2
≡ ∆W,
the Laplace operator.
Tensors in general co-ordinates And what about non-Euclidean spaces, with non-rectangular co-ordinate systems ? systems ? Well, in that case the difference between sub- and superscripts becomes meaningful, and the reason we write superscripts superscripts may finally be appreciated. appreciated. A contravariant vector transforms as follows:
vi =
αii v i
(10.3)
αii vi .
(10.4)
i
and a covariant vector as follows: vi =
i
(These look very similar, but are not the same!) 4 5
More precisely, the lengths of the semi-axes are the square roots of the eigenvalues λi . germ. Spur .
73
Chapter Chapter 10. The surface surface theory of Riemann Where the “prototype” of a contravariant vector is the co-ordinate differences if two (nearby) points, as above: ∆x vi = , ∆y
is the prototype of a covariant vector the gradient operator : operator : v j =
∂V = ∂x j
∂V/∂x ∂V/∂y
.
(10.5)
V (x, y) is some scalar field in space. If we take for a model vector i
v = we obtain
dx dy
vi =
also with the chain rule
∂V/∂x ∂V/∂y
=
=
dx dy
,
∂x /∂x
∂x /∂y
∂y /∂x
∂y /∂y
∂x/∂x
∂x/∂y
∂y/∂x
∂y/∂y
dx dy
∂V/∂x ∂V/∂y
;
.
∂x i ∂x i i The coefficient matrices = α and = αii are apparently each other’s inverse matrices . i ∂x i ∂x i So, if the matrix form of the covariant transformation equation’s (10.4 (10.4)) transformation param i eters αi is A (i row and i column index), then the matrix of the contravariant transformation parameters (equation 10.3 10.3)) αii must be A−1 (i row and i column index). From this circumstance, the names “covarian “covariant” t” and “contrav “contravarian ariant” t” derive. derive.
A tensor may have both super- and subscripts. Under a co-ordinate transformation, each index changes changes according to its “character “character””. There may even be more than two indices.
“Trivial” “Trivial” tensors 1. If we compute the gradient of a contrav contravarian ariantt vector xi we obtain ∂x i = ∂x j
≡ 1 0 0 1
δ ji ,
the so-called Kronecker6 delta, in practice a unit matrix. It is a tensor: δ ji = α j j αii δ ji =
αii δ ji α j j
i,j
or in matrix language i I j
i A−1 i
↓ ↓ =
→
→
i [I ] j
↓ →
j [A] j
↓ →
= A−1 A = I ,
because if the matrix form of α j j is A, then that of αii – or correspondingly, that of α j j is A−1 , as shown above.
6
Leopold Kronecker (1823 – 1891), German mathematician
74
10.2. The metric metric tensor 2. The “corksc “corkscrew rew tensor tensor”” of Tull Tullio io Levi-Civita7 in three dimensions: dimensions: ijk =
−
0 1 1
if, of ijk,twoarethesame ijk ,twoarethesame if ijk if ijk is an even even permutat permutation ion of the numbe numbers rs (123) (123) if ijk if ijk is an odd perm permutati utation on of the the num numbers bers (123) (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. 10.2. The metric metric tensor tensor The metric tensor, or the metric , describes the Pythagoras theore theorem m in a curve curved d space. It is identical identical to the earlier discussed discussed Gaussian Gaussian First Fundamental undamental Form. Form. In the ordinary ordinary plane ( R2 ) we may choose rectangular co-ordinates (x, (x, y ), after which we may write the distance s between to points 1, 1, 2 as: s2 = ∆x2 + ∆y ∆y2 , where ∆x ∆x = x2 x1 , ∆y = y2 y1 are the co-ordinate co-ordinate differences differences between 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 dx j , where gij =
1 0 0 1
i
j
ja dx = dx =
(10.6)
dx dy
.
This is called the metric of the rectangu rectangular lar co-ordin co-ordinate ate system system in the Euclidea Euclidean n plane. plane. In equation (10.6 (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
2
ds =
2
gij dxi dx j .
i=1 j=1 j =1
An alternative way of writing using matrix notation is a quadratic form : ds2 = =
dx · H dx = dxT H dx = dxT
H
dx
g11 g12 g21 g22
dx1 dx2
dx
1
2
dx
.
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 correspo corresponds nds also a location location differen difference ce of 1 m. Then 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 places.. E.g., E.g., on the surface surface of the Earth, Earth, only only within within a strip strip in the vicini vicinity ty of the equato equator. r. It won’t won’t be possible possible in a unified unified manner over the whole surface of the Earth.
75
Chapter Chapter 10. The surface surface theory of Riemann v
+ ∆
v
s v ∆ q
v
α
p∆ p∆u
u
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 rectangular 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 + 2 p 2 p∆ ∆uq ∆v cos α, or, differentially ds2 = p2 du2 + q 2 dv 2 + 2 pq 2 pq cos cos αdudv, or, in index notation ds2 = gij dxi dx j where gij =
p2 pq cos pq cos α pq cos pq cos α q 2
, dxi =
du dv
.
The matrix representation:
ds2 =
gij ,
dxi
du dv
dx j
↓ i, → j
2
p pq cos pq cos α pq cos pq cos α q 2
.
du 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 Also on a non-cu non-curv rved ed surface surface can gij depend on place, e.g., when choosing choosing curvilinear curvilinear,, e.g., polar, p olar, co-ordinates. co-ordinates.
Examples: 1. The surface surface of a spherical spherical Earth: ds2 = R2 dϕ2 + R2 cos2 ϕdλ2 , i.e., gij (ϕ, λ) = 7
R2 0 0 R2 cos2 ϕ
Tullio Levi-Civita (1873 – 1941), Italian mathematician. mathematician.
76
, dxi =
dϕ dλ
,
10.3. The inverse inverse metric tensor tensor from which which in “matrix language” language” gij ,
dxi
dx j
↓ i, → j
ds2 =
R2 0 2 0 R cos2 ϕ
dϕ dλ
dϕ dλ
.
2. Polar Polar co-ordinates co-ordinates (ρ, (ρ, θ) in the plane: ds2 = dρ2 + ρ2 dθ2 , i.e., gij =
1 0 0 ρ2
dρ dθ
, dxi =
from which gij ,
i
dx
,
dx j
↓ i, → j
ds2 =
1 0 0 ρ2
dρ dθ
.
dρ dθ
3. In three dimensions in the air spac spacee using “aviatio “aviation n co-ordinates o-ordinates” ”: nautic nautical al miles miles North dN , dN , nautical miles East dE , feet above sea level dH ; in metres ds2 = (1852)2 dN 2 + (1852)2 dE 2 + (0. (0.3048)2 dH 2 , eli gij = ja 2
ds =
dxi
dN
↓ →
18522 0 0 2 0 1852 0 0 0 0.30482
dE dH
gij ,
i
dx =
i,
dN dE dH
,
dx j
j
18522 0 0 2 0 1852 0 0 0 0.30482
dN dE dH
.
10.3. The inverse inverse metric tensor tensor The inverse tensor of the metric tensor gij is written g ij . As a matrix matrix,, it is the invers inversee matrix of gij : gij = (gij )−1 , or, in index notation gij g jk = δ 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 ritten out: g ij , i, j g jk , j, k δ ki , i, k
↓ → ↓ → ↓ → g11 g12 g21 g22
g11 g12 g21 g22
=
1 0 0 1
.
77
Chapter Chapter 10. The surface surface theory of Riemann
10.3.1. 10.3.1. Raising Raising or lowering lowering sub- or superscripts of a tensor 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 : ki T ji = gik T kj kj = g jk T . ij i All forms T ij designate ate the same same tensor tensor,, writte written n in different different ways. ways. As a special special case, ij , T , T j design i ki ik δ j = g jk g = g gkj , i.e., the Kronecker delta tensor is a “mixed form” of the metric tensor and could be written g ji . The delta style of writing has established itself, however.
10.3.2. 10.3.2. The eigenvalues eigenvalues and -vectors of a tensor tensor The eigenvalue problem for a square tensor has the following form:
T ij
− λgij x j = 0, j (T ij ij − λgij ) x = 0, T ji − λδ ji x j = T ji − λδ ji
xi = 0.
All three forms are equivalen equivalent, t, which which is easily proved. proved. For eigenvec eigenvectors tors we have have xi = gij x j . If ij ji the tensor T ij (or T ji , or T ij ij ) is symmetric (meaning T = T etc.), then the eigenvalues λ are real and the eigenvectors mutually orthogonal: if xi , yi 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
2
R
two.
10.3.3. 10.3.3. The graphic representat representation ion of a tensor The quadratic form i j T ij ij x x = 1
defines an ellipsoid (in R2 an ellipse) that may be considered considered the graphic of T ij ij . E.g., the inertial ellipsoid, the variance ellipsoid.
10.4. The Christoff Christoffel el symbols The metric metric tensor tensor isn’t isn’t yet the same same thing thing as curv curvature. ature. It isn’t even even the same thing thing as the 8 curvature of the co-ordinate curves; to study this, we need the Christoffel symbol s. symbol s. The Christoffel symbols describe what happens to a vector when it is transported transported parallel parallel ly 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 ∆ vi = v vi , although for the “vector itself” v = v. When When both the transpor transported ted vector vector and the distanc distancee 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:
−
i j k ∆v i = Γ jk v y . 8
Elwin Bruno Christoffel (1829–1900), German mathematician-physicist
78
10.4. The Christoffel symbols
v1 v
v2
2
u +y
2
yi v1
v
v2
u2 u1
u1 + y1
Figure 10.2.: Christoffel symbols and parallel transport i The Christoffel symbols Γ jk do not form a tensor like the metric tensor does. They describe the “curviline “curvilinearity” arity” of a curviline curvilinear ar co-or co-ordinate dinate system, i.e., a prop property erty of the co-or co-ordinate dinate system. A co-ordinate transformation that removes – at least locally if not everywhere – the curvature of i the co-ordinate curves, also makes the elements of Γ jk vanish in that point (for tensors, such a local “transforming “transforming away” away” will never never succeed!). succeed!).
E.g., E.g., on the Earth Earth surfac surfacee in the (ϕ, λ) coco-ord ordina inate te system system on the equato equatorr the coco-ordi ordinat natee i curves are locally non-curved and later we shall show, theat there indeed all Γ jk vanish. anish. Another example is from the general theory of relativity, where the components of acceleration are symbols: acceleration acceleration can be b e transformed away away by moving moving to a “falling “falling along” Christoffel symbols: reference reference system. 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 connected to the elevator. elevator. The Christoffel symbols can be computed from the metric tensor; the equation is (complicated proof; see appendix C): ∂g j ∂g jk 1 ∂g k i Γ jk = gi + . (10.7) 2 ∂x j ∂x k ∂x
−
(g ij (gij )−1 .) From this this is can be seen seen,, that that the Chri Christ stoff offel el symbol symbolss are, are, like like,, the the first first derivativ derivatives es of the metric. In a straight-lin straight-lined ed co-ordinate system system on a non-curve non-curved d surface surface gij is a i constant and thus all Γ jk vanish. Also, because gij is symmetric (g (gij = g ji ), we obtain
≡
i i Γ jk = Γkj .
The Christoffel symbols are useful when writing the equation of the geodesic 9 in this formalism (cf. equation (9.8 (9.8)): )): d i i j k t + Γ jk t t = 0. ds Here, ti is the tangent vector of the geodesic ti = dxi/ds. 9
In fact we write
Dt i dti i ≡ + Γ jk tj tk , ds ds the so-called absolute or covariant derivative which is a tensor; and Dt i =0 ds is then the equation for the geodesic, which thus applies in curvilinear co-ordinates.
79
Chapter Chapter 10. The surface surface theory of Riemann
Examples: 1. On the surface surface of the spherical spherical Earth 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!): work!): ∂g k ∂g k 0 0 0 0 = , = , 2 0 0 0 R sin2ϕ sin2ϕ ∂ϕ ∂λ
−
so only ∂g λλ ∂g 22 = = ∂ϕ ∂ϕ
−R2 sin2ϕ sin2ϕ
is non-vanish non-vanishing. ing. Then, remembering remembering that i
1
−
g = (gi )
=
R2 0 2 0 R cos2 ϕ
1
−
R−2 0 −2 0 R cos−2 ϕ
=
,
the only non-zero elements are: Γϕ λλ λ Γϕλ =
−
1 = (gϕϕ )−1 2 1 (gλλ )−1 2
∂g λλ ∂ϕ
∂g λλ ∂ϕ
=
==
1 − 12 R 2 · R2 sin2ϕ sin2ϕ = + sin2ϕ sin2ϕ = sin ϕ cos ϕ, 2 −
−
1 −2 R cos−2 ϕ 2
λ Γλϕ =
1 (gλλ )−1 2
R2 sin2ϕ sin2ϕ =
∂g λλ ∂ϕ
=
sin2ϕ sin2ϕ − 2cos = − tan ϕ, 2ϕ
− tan ϕ.
i Note that at the equator ϕ = 0 all Γ jk = 0.
2. Polar co-ordinat co-ordinates es in the plane (ρ, (ρ, θ): ∂g k = ∂ρ
0 0 0 2ρ
∂g k , = ∂θ
0 0 0 0
,
yielding ∂g 22 ∂g = θθ = 2ρ, ∂ρ ∂ρ so Γρθθ
=
θ θ Γρθ = Γθρ =
−
1 ∂g θθ (gρρ )−1 2 ∂ρ 1 ∂g (gθθ )−1 + θθ 2 ∂ρ
1 2ρ = ρ, 2 1 1 1 = 2 ρ = + . 2 ρ2 ρ =
·−
−
· ·
10.5. The geodesic geodesic revisited revisited Here are alternative formulas for integrating the geodesic on the sphere (generalization to the ellipsoid of revolution is complicated complicated but possible): possible): dξ + η 2 sin ϕ cos ϕ = 0, ds dη 2ηξ tan ηξ tan ϕ = 0. ds
−
i Here we have already used the formulas derived above for the elements of Γ jk .
80
10.6. The curvature curvature tensor tensor 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 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. The length of the tangent vector 2
i j
ds = gij t t = R if we remember that
T
ξ η 2
is computed as follows:
cos A R
gij =
2
2
2
+ R cos ϕ
R2 0 2 0 R cos2 ϕ
sin A R cos ϕ
2
= 1,
on the the surfa surface ce of the the sphe sphere re,, The The leng length thss of the the tang tangen entt vect vectors ors have have alwa always ys to be 1; this this requiremen requirementt is fulfilled fulfilled if s is a parametrizati parametrization on of the curve “by distance” distance”..
10.6. The curvature curvature tensor tensor The curvature is curvature is described again in a slightly more complicated way by means of parallel transport around a closed path — a small rectangle10 . As seen from figure 10.3 we obtain the change ∆v ∆v i on
≡
− v
i
v i of the vector v i, which depends
1. the orientation orientation of the closed rectangular rectangular path, the side indices k and ; 2. the size of the rectangle, rectangle, measured measured in curvilinear curvilinear co-ordinates: co-ordinates: let the length of the uk side be∆u be∆uk and the length of the u side, ∆u ∆u ; 3. the transported transported vector vector v i itself. In the following way: i ∆v i = R jk v j ∆uk ∆u .
(10.8)
i Here, the creature R jk is called the Riemann curvature tensor . In two-dimensional space (i.e., on a surface) it has 24 = 16 elements. elements.
We can compute Riemann from the Christoffels11 : i R jk = 10 11
∂ i Γ ∂x k j
i m i m − ∂x∂ Γ jki + Γkm − Γm Γ j Γ jk .
(10.9)
More precisely, a parallellogram The derivation is found in appendix D.
81
Chapter Chapter 10. The surface surface theory of Riemann
C
D ∆ u u + ∆u
v
i
v
∆u
i
∆uk A
B
u uk + ∆u ∆uk
uk
Figure 10.3.: The curvature tensor and parallel transport around a closed grid path From this we spot immediately the following antisymmetry: i R jk =
i −R jk .
There exist many more such symmetries; the number of independent components is actually small.
Esimerkkej¨ Esimerkkej¨ a: a: 1. Surfac Surfacee of a spherica sphericall Earth (ϕ, (ϕ, λ): Based on the antisymmetry property we may say that only elements for which k = l can differ differ from zero. Additi Additiona onally lly the Christoffel symbols depend only on ϕ. Thus Thus we obtain the auxiliary terms:
∂ λ Γ ∂ϕ ϕλ
∂ ϕ Γ = ∂ϕ λλ ∂ λ = Γ = ∂ϕ λϕ
∂ 1 sin2ϕ sin2ϕ = +cos +cos 2ϕ, ∂ϕ 2 ∂ 1 ( tan ϕ) = , ∂ϕ cos2 ϕ
−
−
the others zero. i m The terms Γkm Γ jl are obtained in the following way: ϕ ϕ λ λ λ Γϕ λλ Γϕλ = Γλϕ Γλλ = Γλλ Γλϕ =
− sin2 ϕ,
λ λ λ λ Γϕλ Γλϕ = Γϕλ Γϕλ = + ta t an2 ϕ.
By combining
82
ϕ Rϕϕλ =
ϕ −Rϕλϕ
=
λ Rϕϕλ =
λ −Rϕλϕ
=
ϕ Rλϕλ =
ϕ −Rλλϕ
=
λ Rλϕλ =
λ −Rλλϕ
=
∂ ϕ Γ ∂ϕ ϕλ ∂ λ Γ ∂ϕ ϕλ ∂ ϕ Γ ∂ϕ λλ ∂ λ Γ ∂ϕ λλ
ϕ ϕ ϕ m m − ∂λ∂ Γϕϕ − Γλm + Γϕm Γϕλ Γϕϕ = 0;
1 λ λ m λ m − ∂λ∂ Γϕϕ + Γϕm Γϕλ − Γλm Γϕϕ = − 2 + tan2 ϕ = −1; cos ϕ ϕ ϕ ϕ m − ∂λ∂ Γλϕ + Γϕm Γm cos 2ϕ + sin2 ϕ = cos2 ϕ; λλ − Γλm Γλϕ = cos λ λ λ m − ∂λ∂ Γλϕ + Γϕm Γm λλ − Γλm Γλϕ = 0.
10.6. The curvature curvature tensor tensor 2. Polar co-ordin co-ordinates ates (ρ, (ρ, θ) plane:
∂ θ Γ ∂ρ ρθ
∂ ρ Γ = ∂ρ θθ ∂ θ = Γ = ∂ρ θρ
−1, − ρ12 ,
the others again vanish. Then: θ θ ρ θ Γρθθ Γρθ = Γθρ Γθθ = Γρθθ Γθρ =
−1,
θ θ θ θ Γρθ Γθρ = Γρθ Γρθ = +
1 . ρ2
After this: ρ Rρρθ =
θ −Rρθρ θ θ Rρρθ = −Rρθρ ρ Rθρθ =
ρ −Rθθρ θ θ Rθρθ = −Rθθρ
= 0; ∂ θ Γ ∂ρ ρθ ∂ ρ = Γ ∂ρ θθ = 0. =
1 1 θ m θ m − ∂θ∂ Γρρθ + Γρm − Γθm Γρθ Γρρ = − 2 + 2 = 0; ρ ρ ρ ρ m − ∂θ∂ Γθρρ + Γρm Γm θθ − Γθm Γθρ = −1 + 1 = 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
R jk =
i R jik
=
i R jik .
(10.10)
i=1
This is a symmetric tensor, Rij = R ji .
Examples 1. Let us continue continue with the computation 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-ordinate co-ordinatess in the plane, all Rij = 0. We may continue this process to obtain the curvature scalar 13 : 2
ij
R = g R ji =
(gij )−1 R ji .
(10.11)
j=1 j =1
For the example case: 12 13
Gregorio Ricci-Curbastro (1853 – 1925), 1925), Italian Italian mathematician mathematician,, invento inventorr of tensor tensor calculus calculus i ik i This is in fact the “trace” of Rij , more precisely, of Rj ≡ g Rkj , written Ri , see above. This also serves as an example example of how in general general curvilinear curvilinear co-ordinates co-ordinates an index may be “raised” “raised” from covariant covariant to contrav contravarian ariant, t, ij or “lowered” “lowered”,, using the metric tensor gij or its inverse g . In rectangular co-ordinates these are unit matrices and the difference difference between super- and subscripts subscripts is without without consequenc consequence. e.
83
Chapter Chapter 10. The surface surface theory of Riemann 1. Spherical Spherical surface: surface:
Rki = gij R jk =
i.e.14
1 R2 0
R=
0 1 R2 cos2 ϕ
i
Rii =
1 0 0 + co c os2 ϕ
=
R−2 0 0 R−2
,
1 cos2 ϕ 2 + = . R2 R2 cos2 ϕ R2
2. In polar co-ordinate co-ordinatess Rij = 0 , i.e., R = 0. More generally we have (without proof): R = 2K = 2κ1 κ2 =
2 , R1 R2
(10.12)
twice the Gauss total curvature K , the inverse of the product of the two principal radii of curvature.
10.7. Gauss curvature curvature and spherical spherical excess excess Contrary to the eigenvalues of the Gauss second fundamental form β ki = gij β jk (chapter 9.4 9.4)) i 15 are the eigenvalues of the tensor R j unrelated to the principal curvatures of the surface κ1 and i κ2 . In th the two-dimensional case both the Rij tensor and the R jkl 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 equation (10.8 (10.8)) ∆u ∆ u∆v describe the sides of a small diamond shape, around which the i vector v is transported in a parallel parallel fashion. In two dimensions dimensions there are only two choices choices k k k for the sides in the u , u directi directions ons:: ∆u ∆u and ∆u ∆u ∆u . One represe represent ntss clockwis clockwisee i i transport, transport, the other, countercl counterclockwis ockwise. e. The corresponding corresponding elements elements R jk = R jk thus are essentially the same.
−
2. In the same equation v j and ∆v ∆v i are mutual mutually ly perpendic perpendicula ular. r. ∆v i represents a small rotation of vector v j . On a surfac surface, e, in two dimensio dimensions, ns, the rotatio rotation n is described described by one angle . Because Because the angle between between two parallelly parallelly transported vectors vectors doesn’t change, change, we i may infer that this rotation angle is the same for all vectors v . Aga Again in we find find only one one independent parameter. See figure 10.4. i 1 2 In othe otherr word words, s, when when give give a ∆u1 ∆u2 diamon diamond d shape, shape, the expres expressio sion n R j12 j 12 ∆u ∆u = cos θ sin θ i 2 1 2 rotation matrix R j21 , containing only one free pa j 21 ∆u ∆u is a 2 sin θ cos θ rameter.
−
×
−
Howe Ho weve ver, r, for higher higher dimens dimension ionali alitie ties, s, the numbe numberr of indepen independen dentt elemen elements ts in the tensor tensorss of larger. This This case case is inter interest esting ing because because of General General Relativ Relativit ity y (four (four Riemann and Ricci is larger. 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 ! 14
Here we use the symbol R both for the radius radius of the Earth Earth and for the curv curvature ature scalar scalar.. We hope it doesn’t doesn’t cause confusion. 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.
84
10.7. Gauss curvature and spherical excess
du2 K du1 θ = K du1 du2 Figure 10.4.: The curvature of a two-dimensional surface may be characterized by just one parameter: θ.
C 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 . When we transport a vector around a larger surface area, it is the same as if we had transported it successive successively ly around all of the little patches patches making up the surface surface area. This is depicted in 1 2 figure 10.5. In every little patch du du , area dS , the change in direction of the vector amounts to K dS , 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, triangle, Theorema Theorema elegantissimum ): ): ε=
¨
K ( K (ϕ, λ) dS,
∆
where K is a function of place. On the reference ellipsoid K = (M N )−1 and ε=
¨
1 dS = M N
∆
¨ 1 M N cos N cos ϕdϕdλ = M N ¨ ∆
=
dσ = σ∆ ,
∆
which is the surface area of the corresponding triangle, but on the unit sphere .17 . From this follows that 16 17
Pierre Ossian Bonnet (1819 – 1892), French mathematician This is not quite exactly true. . . the triangle triangle symbols symbols under under the the integral integral sign sign dϕdλ and under the integral sign
85
Chapter Chapter 10. The surface surface theory of Riemann 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 xi =
cos ϕi cos λi cos ϕi sin λi sin ϕi
, i = 1, 2, 3;
we use the polarization method (1.7 1.7)) in three dimensions, in order to find the “poles” of the triangle’s sides: y1 =
x2 × x3 , y2 = x3 × x1 , y3 = x1 × x2 . x2 × x3 x3 × x1 x1 × x2
Now the inter-pole distances correspond to the angles of the spherical triangle (more precisely, π minus those angles):
y2 − y3 , − 2 arct arctan an y2 + y3 y1 − y3 , π − 2 arct arctan an y1 + y3 y − y2 , π − 2 arct arctan an 1 y + y2
α1 = π α2 = α3 =
1
numerically strong equations18 . He Here re αi is the angle at corner point i. The spheri spherical cal excess excess is now 3
ε=
αi
i=1
− π.
10.8. The curvature curvature in quasi-Eucl quasi-Euclidean idean geometry 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 tensor is the unit matrix, gij = gij =
1 if i = j 0 otherwise, otherwise,
2. the metric metric tensor is (locally) (locally) stationary , i.e. ∂g ij ∂x k
= 0. k xk =xP
In this case we speak of a quasi-Euclidean neighbourhood neighbourhood around P . P . dσ are not exactly exactly correspon correspondin ding. g. On the reference reference elipsoid elipsoid the triang triangle le consis consists ts of geodesi geodesics, cs, 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 ellipsoidal and a spherical spherical triangle are not individually individually equal; their sums (and with that, their spherical excesses) excesses) however however are. 18 In some programming languages, languages, we have for arctan ( x/y) x/y ) the symmetrical symmetrical alternative alternative form atan2 (x, y ), where the zerodivide problem does not occur.
86
10.8. The curvature curvature in quasi-Eucli quasi-Euclidean dean geometry In this case the Christoffel symbols (equation 10.7 10.7)) are
∂g jk 1 i ∂g k ∂g j g + 2 ∂x j ∂x k ∂x ∂g jk 1 ∂g ki ∂g ij + 2 ∂x j ∂x k ∂x i
i Γ jk =
=
−
≈
−
based on the above assumption of stationarity.
= 0
Of the Riemann curvature tensor (equation 10.9 10.9)) the last two terms vanish: i R jk
∂ i Γ ∂x k j 1 ∂ 2 ∂x k 1 ∂ 2 ∂x 1 ∂ 2 ∂x j
=
− ∂x∂ Γ jki =
=
− =
∂g i ∂g ij + ∂x j ∂x ∂g ki ∂g ij + ∂x j ∂x k ∂g i ∂g ki ∂x k ∂x
−
− − − −
∂g j ∂x i ∂g jk = ∂x i 1 ∂ ∂g j 2 ∂x i ∂x k
−
∂g jk ∂x
If we now derive the Ricci tensor (equation 10.10 10.10)): i = R j = R ji
+ =
1 2
−
− 12
i
∂ 2 g j (∂x i )2
i
∂ 2 gi ∂x i ∂x j
1 1 ∆g j + 2 2
∂ 2 gi ∂ 2 gij + i ∂x i ∂x j ∂x ∂x
i
.
+
∂ 2 gii ∂ 2 g ji + ∂x ∂x j ∂x ∂x i
−
=
−
∂ 2 gii ∂x ∂x j
.
This already looks a lot more symmetric. The symbol ∆ here means the Laplace operator ∆=
i
∂ 2 (∂x i )2
.
Finally we still derive the curvature scalar jk
R = g Rkj =
R jj =
j
= = =
− − − 1 2
∆g jj +
j
1 2
i
j
∆g jj +
j
i
j
∆gii +
i
i
j
∂ 2 gij ∂x i ∂x j ∂ 2 gij ∂x i ∂x j
∂ 2 gij . ∂x i ∂x j
− − 1 2 1 2
i
j
∂ 2 gii (∂x j )2
=
∆gii =
i
In the special case that the co-ordinate curves are (everywhere, not just in point P ) P ) orthogonal orthogonal,, we obtain gij = 0 if i if i = j with its derivatives of place, i.e.
R=− =
− −
∂ 2 gii
∆gii +
i
∂g 11
i
(∂x i )2
∂g 22
=
∂g 11
∂g 11
∂g 22
∂g 22
+ + + + = (∂x 1 )2 (∂x 2 )2 (∂x 1 )2 (∂x 2 )2 (∂x 1 )2 (∂x 2 )2 ∂g 11 ∂g 22 = + . (∂x 2 )2 (∂x 1 )2
(10.13)
This equation will be of use in the sequel.
87
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. 11.1. Map projec projectio tions ns and scale scale 11.1.1. 11.1.1. On the Earth 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λ. dλ. These These correspond correspond to linear linear distance distancess M ( M (ϕ) dϕ and p (ϕ) dλ, dλ, respectively. According According to Pythagoras , the length of the diagonal of the postage stamp ( dϕ,dλ) dϕ,dλ) is dS 2 = M 2 dϕ2 + p2 dλ2 . Now dS 2 defines a metric , dS 2 =
(11.1)
gij dxi dx j = xT H x,
i,j
where gij =
M 2 0 0 p2
and x=
=
dx1 dx2
=
E F F G dϕ dλ
= H
(11.2)
.
Here we have suitably defined two ways of writing, the index notation gij , dxi and the matrixvector vector notation 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: tan A =
pdλ , M dϕ
(11.3)
from which, with the above, dS sin dS sin A = pdλ,
(11.4)
dS cos dS cos A = Mdϕ.
(11.5)
89
Chapter Chapter 11. Map projections projections
11.1.2 11.1.2.. In the map plane When we project this little rectangle into the map plane, we obtain the sides dx and dy, dy , and according 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 dϕ + ∂ϕ ∂y dϕ + ∂ϕ
dx = dy =
∂x dλ, ∂λ ∂y dλ, ∂λ
i.e., ds
2
=
∂x ∂x dϕ + dλ ∂ϕ ∂λ
2
+
∂y ∂y dϕ + dλ ∂ϕ ∂λ
= Edϕ Edϕ 2 + 2Fdϕdλ F dϕdλ + Gdλ Gdλ2 , where
= (11.6)
∂x 2 ∂y 2 + , ∂ϕ ∂ϕ ∂x ∂x ∂y ∂y + , ∂ϕ ∂λ ∂ϕ ∂λ
E E = F F = G =
2
∂x ∂λ
2
+
∂y ∂λ
2
.
E, E, F F ja G 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, E E F F gij = H. (11.7) F F G
≡
The corresponding metric is
ds2 =
gij dxi dx j = xT H H x,
i,j
T
where dx1 = dϕ and dx2 = dλ, dλ, , i.e., dxi = dx j = x = dϕ dλ . Also Also here here we we see as alternatives the index notation and the matrix-vector notation, which describe the same thing.
11.1.3 11.1.3.. The scale scale The scale is now the ratio
ds , dS which apparently depends on the direction of the distance element, i.e., the azimuth A. m=
Let us write the eigenvalue problem : 2
m
⇒ − ⇒ =
ds2 = dS 2
i j i,j gij dx dx i j i,j gij dx dx
gij dxi dx j
i,j
90
m2
gij dxi dx j = 0.
i,j
(11.8)
11.1. Map projections and and scale
The same equation in “matrix language” language”,, if x = T
x
Equations (11.4 (11.4,, 11.5 11.5)) The read:
T
dx1 dx2
−
m2 H x = 0.
H H
=
dϕ dλ
T
:
dS sin dS sin A = pdλ, dS cos dS cos A = Mdϕ. Let us now look at all vectors that are of form x=
dϕ dλ
cos A M sin A p
=
.
(11.9)
The lengths of these vectors on the Earth surface are:
dS 2 = M 2 dϕ2 + p2 dλ2 = 2
= M
2
cos A M
+ p
sin A p
2
2
=
= cos2 A + sin2 A = 1.
So: on the surface of the Earth these vectors form a circle of unit radius . Substitute Substitute (11.9 11.9)) into equation (11.8 (11.8), ), using (11.6 (11.6): ):
−
cos A E E M 2
m
2
sin A p
+G
2
M
2
+ 2F F
2
cos A M
2
+ p
cos A sin A M p
sin A p
−
2
=0
or, after cleaning up, m2 = =
− − − −
E E G F F cos2 A + 2 sin2 A + 2 sin A cos A = 2 M p Mp 1 2
E E G + M 2 p2
1 2
+
E E M 2
G p2
cos2A cos2A +
(11.10)
F F sin2A. sin2A. Mp
From this we obtain the stationary stationary values : 0= in other words
d 2 G m = dA p2
tan2A tan2A =
E E M 2
sin2A sin2A +
G p2
Ep Ep 2
E M 2
F Mp
=
F F cos2A, cos2A, Mp GM GM 2
F 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 ( 11.10)) as follows:
1
cos A sin A
− − E E M 2
F F Mp
m
F F Mp
2
G p2
2
m
cos A sin A
= 0;
◦
If F = 0, we obtain the condition condition sin 2 A = 0, which is fulfilled when A = k · 90 , k = 0, 1, 2, 3.
91
Chapter Chapter 11. Map projections projections This presupposes that the determinant of the matrix in the middle vanishes :
− − − − − − − ± − −
0 = det H −1 H H =
E E M 2
m2
m21,2 =
E M 2
F F 2 = M 2 p2
G p2
m2
G p2
m2 +
E E M 2
= m4 + From this
m2 I =
E M 2
+ pG
2
1 E E G M 2 p2
2
+ pG
4 M 1 p
2
2
F F 2 .
E E G
2
F F 2
2
.
These two solutions are called the principal scale factors m1 , m2 .
If the H H matrix is a diagonal matrix:
we obtain
− − − H H =
det H −1 H H
m2 I
E E M 2
=
m21,2 =
E E 0 0 G
,
m2
G p2
m2
=0
⇒
E E G , . M 2 p2
E E G is the meridional scale factor , m2 = is the scale factor in the direction of the 2 M p2 parallel . In intermediate directions A (azimuth) the magnification factor is then m1 =
m=
m21 cos2 A + m22 sin2 A.
11.1.4. 11.1.4. The Tissot-indicatr Tissot-indicatrix ix The matrix (“scale tensor”)
H −1 H H =
E E M 2 F F Mp
F F Mp G p2
= gik gkj = g ji
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 meridional direction direction and m2 in the direction of the parallel. This ellipse is called the indicatrix 2 of Tissot . See subsection 10.3.3. 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
92
11.2. The Lambert projection (LCC) x ϕ m1
M dϕ = 1 λ
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 Conformal projections have have the useful property property, that the projection formulas in one co-ordinate direction can be derived when the formulas for the other co-ordinate direction are given. A sometimes useful property is, that surface areas are mapped true 3 , even though the shapes of small circles or squares are distorted. This requirement is det H −1 H H = m21 m22 = constant.
11.2. The Lambert project projection ion (LCC) The Lambert4 projection (LCC, Lambert Conformal Conical) is a conformal conical projection. The scale for all latitude circles is constant; constant; it is howeve howeverr different different for different different latitudes, 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 4
E.g. when depicting the surface density of population or some other phenomenon. Johann Heinrich Lambert (1728 – 1777) Swiss methematician, methematician, physicist physicist and astronomer. astronomer.
93
Chapter Chapter 11. Map projections projections 7 5 o o N 6 0 o o N 4 5 o o N 3 0 o o N 1 5 o o N 0 o o
6 0 o o W o E E 9 0
4 5 o o W o E 7 5
3 0 o o W
o E
1 5 o W 0 o
o E
4 5
o
o
30 30 E
15 E
6 0
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 ( M (ϕ) dϕ, dϕ, i.e. dt dt ds nρ = = M. dϕ ds dϕ p If we reckon ρ positive toward the South, i.e., dρ = dρ = dϕ
−dt, dt, we obtain:
M (ϕ) −nρ (ϕ) M ( . p (ϕ)
(11.12)
Using equations (11.11 (11.11,, 11.12 11.12)) we may compute (θ, ( θ, ρ) if given (ϕ, ( ϕ, λ). However, let us first move ρ to the left hand side: d ln ρ = dϕ
M n p
−
ϕ
⇒
ρ = ρ0 exp
If we define
ϕ
ψ (ϕ)
94
− ˆ
≡
ˆ 0
n
ϕ0
M ( M (ϕ ) dϕ , p (ϕ )
M ( M (ϕ ) dϕ . p (ϕ )
11.3. On the isometric isometric latitude the so-called isometric isometric latitude , we obtain ρ = ρ0 exp
exp {−nψ (ϕ)} {−n (ψ (ϕ) − ψ (ϕ0))} = ρ0 exp {−nψ (ϕ0)} .
(11.13)
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
dt nρ ≡ dρ =− =− , ds ds p
m i.e., dm dϕ
n dρ n dp ρ 2 = p dϕ p dϕ n2 ρ nρ = + M M sin M sin ϕ = p2 p2 =
−
−
− pnρ2 M ( M (n − sin ϕ) ,
(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 ) follow follows. s. No Now w ρ (ϕ) is obtained by integrating (11.12 ( 11.12)) from a starting value, either ρ0 or ρ1 , with the integration interval being either [ ϕ0 , ϕ] or [ϕ [ϕ1 , ϕ].
≡
λ; solving in reverse the differential equation (11.12 ( 11.12)) The inverse operation is easy for θ (computing ϕ when ρ is given) can be done as follows:
→
1. invert invert analytically analytically equation equation (11.13 ( 11.13): ): ψ (ϕ) =
− ψ (nϕ0) (ln ρ − ln ρ0) ;
2. perform the inverse inverse computation computation ψ
→ ϕ (see below).
11.3. 11.3. On the isomet isometric ric latitud latitude e The isometric isometric latitude, latitude,
ϕ
ψ=
ˆ 0
M ( M (ϕ) dϕ, p (ϕ)
can be b e computed computed numerically numerically (quadrature; (quadrature; the QUAD routin routines es of Matlab Matlab). ). Howev However, er, for the spherical case a closed solution exists ψ = ln tan tan
π ϕ + , 4 2
(11.14)
95
Chapter Chapter 11. Map projections projections easily proven by differentiating using the chain rule. Even for the ellipsoid of revolution a closed solution exists:
−
π ϕ ψ = ln tan + 4 2
e
1 + e sin ϕ 1 e sin ϕ
2
.
(11.15)
See appendix A.
Reverse operation: operation: if ψ is given and ϕ to be computed, we may take the equation dψ M ( M (ϕ) = dϕ p (ϕ) and turn it upside down: dϕ p (ϕ) = . dψ M ( M (ϕ) This equation has the general form dy = f ( 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 11.14): ):
−
ϕ = 2 arctan arctan exp ψ
− − − π . 4
This is also useful as the first iteration step for the ellipsoidal case: π ϕ(0) = 2 arctan arctan exp ψ , 4 after which ϕ(i+1) = 2 This converges rapidly.
arctan
e
e sin ϕ(i)
1 1 + e sin ϕ(i)
2
exp ψ
π 4
.
11.4. The Mercator Mercator projectio projection n 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 = λ, ϕ
y = ψ (ϕ) =
ˆ 0
M ( M (ϕ ) dϕ . p (ϕ )
Here we see the isometric latitude in its most simple glory: in the case of spherical geometry π ϕ y = ln tan tan + . 4 2
Mercator is not a “lamp “lamp projectio projection” n”:: There There is no projectio projection n centre centre from which which the “light “light””
emanates.
96
11.5. The stereographi stereographicc projection
y dt ds
dt
ds x
ϕ
λ
Figure 11.4.: The Mercator projection
o
75 N
o
60 N
o
45 N
o
30 N
o
15 N
o
0
o
60 W
o
45 W
o
30 W
o
15 W
0
o
o
15 E
o
30 E
o
45 E
o
60 E
o
75 E
o
90 E
Figure 11.5.: An example of the Mercator projection
11.5. The stereographi stereographicc project projection ion 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 11.13)) n = 1 and choose ϕ0 = 0 (the equator) and ρ0 correspondingly, i.e., ρ0 ρ (ϕ0 ). In that case
≡
ρ = ρ0 exp Unfortunately in the limit ϕ
{−ψ (ϕ)} .
π the function function ψ (ϕ) diverges; we may define 2 π ρ = ρ0 exp ψ (ϕ) if ϕ if ϕ < , 2 π ρ=0 if ϕ = , 2
→
{−
}
after which ρ (ϕ) is continuous.
97
Chapter Chapter 11. Map projections projections
ρ
ρ0
Projection plane
π/2
ρ0
−ϕ ϕ
Equator
(π 2
ϕ)/2
−
Projection centre Figure 11.6.: The stereographic projection
In the spherical case we obtain ρ =
ρ0/tan( π + ϕ )
= π ϕ = ρ0 cot + = 4 2 π ϕ = ρ0 tan = 4 2 1 π = ρ0 tan ϕ . 2 2 4
2
− −
This case is depicted in figure 11.6 11.6.. ρ0 is the distance from the central point (“South pole”) to the projectio projection n plane. plane. In this case the projection projection is truly truly a “lamp “lamp projection projection””. . . but unfortu unfortunat nately ely only in the case of spherical geometry. Let us derive the scale of the stereographical projection: dρ dρ dψ dϕ = = dx dψ dϕ dx
1 ρ exp {−ψ (ϕ)} −ρ0 exp {−ψ} · M · = − = −ρ0 . p M p p (ϕ)
The value is negative because x grows to the North and ρ to the South. South. The co-ordi co-ordinat natee x is the metric metric “northing” “northing”.. By using equation 11.15 we obtain dρ = dx
−
− ≈−
ρ0 π ϕ cot + N cos( N cos(ϕ ϕ) 4 2
or close to the pole dρ dx
ρ0 2
1+e 1 e
−
1 + e sin ϕ 1 e sin ϕ
−
−
e 2
,
e
2
.
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 01346387371353381594..
98
¨ ger projection 11.6.. The Gauss-Kruger 11.6 u
o
180 W
o
o W
W 1 6 5
1 5 0
16 5 5 o E
1 5 o 0 E
o
30 N
o W W 5 1 3
1 3 5 o o E
o
45 N
W o W 0 1 2
1 2 0 o E
o
60 N
o W
5 1 0
1 0 5 o E
o
75 N
9 0 o W
E
o
0 9
o E 5 7
7 5 o W
E o E 0 6
6 0 o W o E E 5 4
4 5 o o W 3 0 o W 15 o
W
0
o E
o E 3 0
1 5
o
Figure 11.7.: An example of the stereographic projection Reference lati latitu tud de 1◦ 61◦ 70◦ 70◦ 65◦
Latitude inte interv rval al 1◦ 61◦ 61◦ 61◦ 61◦
Number of supp suppor ortt poin points ts
Computation poin pointt (ϕ, ∆λ)
Number of terms
20 20 20 39 19
(19◦ .333 333,, 10◦ .0) (79◦ .333 333,, 20◦ .0) (79◦ .333 333,, 20◦ .0) (79◦ .333 333,, 20◦ .0) (60◦ .333 333,, 10◦ .0)
12 16 15 15 14
◦
− 20 − 80 − 80 − 80 − 70
◦ ◦ ◦ ◦
¨ger test computations. Gauss-Kruger u Table 11.1.: Results of Gauss-Kr¨ computations. “Number of terms” terms”: the number of terms in the polynomial expansion needed to achieve a change under 1mm in the ¨ger co-ordinates Gauss-Kruger u computed Gauss-Kr¨
±
11.6. 11.6. The GaussGauss-Kr Kr¨ uger u ¨ger projection ¨ger projectio A good practical example of a map projection is the Gauss-Kr¨ projection, n, in use in Gauss-Kruger u 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. expansions. But, differently differently from tradition, tradition, we will determine determine the coefficients coefficients of the expansion numerically . Let us proceed proceed in the followin followingg way way. First First we choose choose a suitab suitable le starting projection projection , e.g., an ordinary Mercator, the equati equations ons of which which are simple. simple. So, we map the surfac surfacee of the Earth Earth ellipsoid onto the map plane of an ordinary Mercator: v = λ
, ˆ − λM ( M (ϕ ) 0
ϕ
u =
0
p (ϕ )
dϕ .
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 ( M (ϕ) dϕ,
99
Chapter Chapter 11. Map projections projections i.e.,
ϕ
x=
ˆ
M ϕ dϕ .
0
Furthermore we have on the central meridian:
(11.16)
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., i.e., the real axis. axis. The problem problem is to determ determine ine the function function in the whole complex plane . See figure 11.8 11.8..
≡
≡
Intermezzo. In complex analysis we talk of analytical functions . An analytical function is such a mapping f (w) x = f ( which is differentiable. Not just once, but infinitely many times. In this case, the CauchyRiemann conditions apply: ∂x ∂y ∂y ∂x = , = , ∂u ∂v ∂u ∂v if x = x + iy, w = u + iv. iv .
−
−
This means that the small vector to the following equation: dx dy
=
a b
du dv
b a
du dv
∂x ∂y ∂y where a = = and b = = ∂u ∂v ∂u require.
T
is mapped to a vector
= K
cos θ sin θ
− sin θ cos θ
dx dy
du dv
T
according
,
− ∂x , precisely as the Cauchy-Riemann conditions ∂v
From this it is seen that the mapping of local vectors ( du, du, dv) dv) rotation ; i.e., we have a conformal mapping .
→ (dx, dx, dy) dy) is a scaling and
Almost all familiar functions functions are analytical in the complex plane: the exponent, exponent, the logarithm, 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 : ∂ 2 x ∂ 2 x ∂ 2 y ∂ 2 y + 2 = 0, + = 0. ∂u 2 ∂v ∂u 2 ∂v 2 Let us try, as a general solution, a series expansion :
∞
2
3
z = a0 + a1 w + a2 w + a3 w + . . . =
k=0
Here we define for the sake of generality z w
≡ ≡
(x
− x0) + iy, (u − u0 ) + iv.
The values ϕ0
x0 u0
100
≡ ≡
ˆ ˆ 0
0
ϕ0
M ( M (ϕ) dϕ, M ( M (ϕ) dϕ p (ϕ)
ak w k .
(11.17)
¨ ger projection 11.6.. The Gauss-Kruger 11.6 u
y
x
¨ger projection Figure 11.8.: The Gauss-Kr¨ projection as a boundar boundary y value value problem. problem. The boundary boundary Gauss-Kruger u values at the central meridian are marked with crosses, the directions of integration with arrows
have have been chosen chosen to be compatible with a suitable suitable reference reference latitude ϕ0 . These series series expansions expansions are thus meant to be used only in a relatively small area , not the whole projection zone. Thus Thus the number of terms needed also remains smaller. The meridian conditions (11.16 (11.16)) now form a set of observation equations : xi = a0 + a1 ui + a2 ui2 + a3 ui3 + . . . from which the coefficients a j can be solved, if a sufficient number of support points ( ui , xi ) is given (the crosses in figure 11.8 11.8). ). Applying the found coefficients a j to equation (11.17 (11.17)) yields the solution for the whole complex plane. This solution solution may be “squeezed” “squeezed” into into the following, following, more computationall computationally y suitable, form (so-called Clenshaw summation): z = a0 + w (a1 + w (a2 + w (a3 + w (
· · · + wan)))) .
The computing sequence sequence is multiplic multiplicationation-additi addition-m on-multip ultiplicati lication-add on-addition. ition. . . 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 When the equatio equation n of the map is (11.17 ( 11.17), ), its differentia differentiall version version is dz ∆z = ∆w, (11.18) dw
where dz dw
∞
2
4
= a1 + 2a 2a2 w + 3a 3a3 w + 4a 4a4 w +
··· =
Which also can easily we written into the Clenshaw form.
kak wk−1 .
k=1
101
Chapter Chapter 11. Map projections projections Now we may write the complex equation ( 11.18 11.18)) in real-valued matrix form
dx dy
where dx =
the scale ; and
dx dy
T
dz dw dz dw
Re
=
, dw =
− −
Im
T
du dv
dz dw dz dw
Im
Re
du dv
,
,
Re
dz dw
=
∂x ∂y = = K cos K cos θ, ∂u ∂w
Im
dz dw
=
∂y = ∂u
∂x = K sin K sin θ, ∂v
from which θ, the meridian convergence 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 ( 11.17)) could be reduced by using a Lambert projection for the reference latitude ϕ0 . Whether Whether the saving saving achievable achievable is worth the trouble, is a so-called Good Question.
11.7. 11.7. Curvat Curvature ure of the Earth Earth surface surface and scale Recall the equation 10.13 derived above:
− −
∂g 11 ∂g 22 + . 2 (∂x 2 ) (∂x 1 )2
R=
Because according to equation 10.12 we have R = 2K , we obtain K =
1 ∂g 11 ∂g 22 + . 2 2 (∂x 2 ) (∂x 1 )2
Let us now take a rectangular co-ordinate system in the map plane ( x, y) and transfer its coordinate 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, the metric is the unit matrix and stationary at the origin. In this case the theory in the previous previous chapter chapter applies. In the case of a conformal projection the form of the metric is gij = m
2
−
1 0 0 1
,
where m is the scale of the map projection, which thus depends on location xi . Derive K =
− 1 2
Now
∆ m−2
= =
=
− 12 ∆
m−2 .
∂ 2 ∂ 2 + m−2 (ξ, η ) = 2 2 ∂ξ ∂η ∂ ∂ ∂ ∂ 2 m−3 m 2 m−3 m = ∂ξ ∂ξ ∂η ∂η 2∆m, 2∆ m,
− ≈ −
102
∂ 2 gξξ ∂ 2 gηη + ∂η 2 ∂ξ 2
−
11.7. Curvature Curvature of the Earth surface and scale if we consider the stationarity of m 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 relationship between the scale’s second derivatives in the North-South and East-West directions: if, e.g., ∂ 2 m ∂ 2 m =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
Mercator Lambert conical Oblique stereographic ¨ger, UTM Gauss-Kr¨ Gauss-Kruger u
∂ 2 m/∂x2
∂ 2 m/∂y 2
K K K 2 0
0 0 K 2 K
In the table we used again, instead of (ξ, ( ξ, η ),(x, ),(x, y), i.e., we substituted ∂ 2 m ∂ξ 2
−→
∂ 2 m ∂x 2
and
∂ 2 m ∂η 2
−→
∂ 2 m , ∂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 coordinat ordinate, e, in other other words words,, from map West to map East. In transv transvers ersal al cylind cylindric rical al projectio projections ns 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), ¨ ger again is best for countries extending in the North-South direction whereas Gauss-Kr¨ Gauss-Kruger u (New Zealand). The conformal azimuthal projection called oblique stereographic is suitable for “square” “square” countries countries (The Netherlands Netherlands). ).
103
Chapter Chapter 11. Map projections projections
N
N
N
¨ ger (a) Gauss-Kr¨ Gauss-Kruger u
(b) Mercator, Lambert
(c) Stereographic
Figure 11.9.: The classification of map projections as “vertical”, “horizontal” and “square” pro jections
104
Chapter 12 Map projections in Finland 12.1. Traditional raditional map project projections ions In Finland, the traditional map projection has been Gauß-Kr uger u System ¨ ger with zone width 3 ◦ . System name: kkj (“National “National Map Grid Co-ordinate System”) System”) created created in 1970 (Parm, 1988). 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 o¨ (nr. 90), at latitude and longitude values from the ED50 European adjustment, and geoidal undulation from the Bomford astro-geodetic geoid (Bomfor ( Bomford, d, 1963 1963). ). ◦
◦
◦
◦
◦
◦ Map plane co-ordinates were obtained using Gauß-Kr uger u ¨ 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, (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”, System”, also “Helsinki System”) System”) co-ordinates, cf. (Ollikainen, 1993). 1993 ). Equation:
x y
=
kkj
1.00000075 0.0000 000043 4399
−
−0.00000439
x y
1.00000075
61 61..571m 95 95..693m
+
ED50 ED 50
12.2. Modern map projec projections tions The modern Finnish system is quite different:
◦ Based on the GRS80 reference ellipsoid of 1980; datum is called EUREF-FIN, created by
keeping keeping fixed four stations fixed to their ITRF96 values values at epoch 1997.0: the permanent permanent GPS stations Metsahovi, a¨hovi, Vaasa, Joensuu and Sodankyl a¨ (Ollikainen et al., 2000). 2000). Then, Then, a transformation (Boucher Boucher and Altamimi, Altamimi, 1995 1995)) was applied to obtain co-ordinates in ETRF89 ETRF89.. Thus, Thus, the datum is correct correctly ly describe described d as ETRF89 ETRF89,, but the epoch remains 1997.0, as no correction for individual station motion (mostly, glacial isostatic adjustment) was made in the transformation.
Equation:
E
X (tC ) =
I Xyy
(tC ) + Tyy +
−
0 ˙ R3 R˙ 2
R˙ 3 0 R˙ 1
−
R˙ 2 R˙ 1
−
0
I (tC ) (tC Xyy
·
− 1989 1989..0)
yy
with tC observations central epoch, yy = (19)96. (19)96. The values T96 and R˙ i,96 i,96 are tabulated in (Boucher and Altamimi, ) Tables 3 and 4.
105
Chapter Chapter 12. Map projections in Finland 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. system. This also defines defines the map sheet division. division. Howeve However, r, 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). 2003 ).
◦ For large scale maps, such as used for planning and cadastral work, the Gauß-Kr uger u ¨ger
projection continues to be used (but based on the above reference ellipsoid and datum), with a central central meridian interv interval al of only one degree: ETRS-GK ETRS-GKn n, where n designates designates the central central meridian longitude. This avoids avoids the problem of significan significantt scale distortions. distortions.
12.3. The triangulated triangulated affine transformati transformation on used in Finland 12.3.1. 12.3.1. Plane co-ordina co-ordinates tes The National Land Survey offers a facility to convert kkj co-ordinates to the new ETRS89TM35FIN system of projection co-ordinates. The method is described in the publication (Anon., 2003), 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-wis triangle-wisee 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)
(1)
(1)
where x , y 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 ∆ x, ∆y, a1 , a2 , b1 ja (1) (1) b2 . If, in the three three corners corners of the triangl triangle, e, are given 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, edges, but not differe different ntiab iable: le: the scale is discon discontin tinuou uouss across across triangle triangle edges. edges. Becaus Becausee the mapping is not conformal either, the scale will also be dependent upon the direction considered. A useful property of triangulation triangulation is, that it can be locally “patched” “patched”:: if better data is availabl availablee in the local area – a denser point set, whose co-ordinate pairs x(i) , y (i) , 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 smaller triangle, triangle, inside which the transformation 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:
x(2) y(2)
=
∆x ∆y
+
a1 a2 b1 b2
x(1) y(1)
.
Most often the co-ordinates in the (1)and (2) datums are close to each other, i.e., are small. In that case we may write the shifts δx δy
106
≡ ≡
x(2)
− x(1) = ∆x + (a (a1 − 1) x(1) + a2 y(1) , y(2) − y (1) = ∆y + b1 x(1) + (b (b2 − 1) y(1) .
∆x ∆y
T
12.3. The triangulated triangulated affine transformation transformation used in Finland Finland
7
819 230 31 2
44 85
203
301 22
441
233 503
415
440
510
5
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
≡ ∆x ∆y
, A=
a11 a12 a21 a22
≡
a1 1 a2 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 . 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)
(1)
(1)
(1)
(1)
(1)
(1)
(1)
(1)
(1)
(1)
δx A = ∆x + a11 xA + a12 yA δy A = ∆y + a21 xA + a22 yA
δx B = ∆x + a11 xB + a12 yB δy B = ∆y + a12 xB + a22 yB
δx C = ∆x + a11 xC + a12 yC δy C = ∆y + a21 xC + a22 yC
107
Chapter Chapter 12. Map projections in Finland Finland C
= ω(∆ APC ) ω(∆ ABC )
p B =
p A =
=
= ω(∆ PBC ) ω(∆ ABC )
P
= ω(∆ ABP ) ω(∆ ABC )
p C =
B
Figure 12.2.: Computing barycentric co-ordinates as the ratio of the surface areas of two triangles
or in matrix form
δx A δy A δx B δy B δx C δy C
=
(1) xA
1 0 1 0 1 0
(1) yA
0 0 0 (1) (1) 1 0 xA 0 yA (1) (1) 0 xB 0 yB 0 (1) (1) 1 0 xB 0 yB (1) (1) 0 xC 0 yC 0 (1) (1) 1 0 xC 0 yC
from which they can all be solved.
∆x ∆y a11 a21 a12 a22
,
Let us write the coordinates (x, ( 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 δx A + pB δx B + pC δx C ,
(12.1)
δy = pA δy A + pB δy B + pC δy C.
(12.2)
The set of three numbers pA , pB , pC is called the barycentric co-ordinates of point P See figure 12.2. ω (∆BC (∆BC P ) P ) They can be found as follows (geometrically pA = etc., where ω is the surface area ω (∆ABC (∆ABC )) of the triangle) using determinants:
pA =
xB xC x yB yC y 1 1 1 xA xB xC yA yB yC 1 1 1
, pB =
xC xA x yC yA y 1 1 1 xA xB xC yA yB yC 1 1 1
These equations are very suitable for coding.
108
, pC =
xA xB x yA yB y 1 1 1 xA xB xC yA yB yC 1 1 1
.
Appendix A Isometric latitude on the ellipsoid We follow the presentation from the book (Gross ( Grossman, man, 1964 1964). ). The starting formula is
ϕ
ψ (ϕ) =
ˆ 0
M ( M (ϕ ) dϕ . p (ϕ )
As a differential equation dψ =
M dϕ = p 1
−
1 e2 dϕ. e2 sin2 ϕ cos ϕ
−
The integrand is decomposed into partial fractions:
− 1
1 e2 e2 sin2 ϕ cos ϕ
−
− e2 sin2 ϕ − e2 cos2 ϕ = 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 ϕ 1
= = =
=
The integral of the first term is
ϕ
ψ=
ˆ 0
Proof by using the chain rule: dψ dϕ
= = =
1 π ϕ dϕ = lntan + . cos ϕ 4 2
d ln tan tan π4 + ϕ2 d tan π4 + ϕ2 d π4 + ϕ2 = dϕ d tan π4 + ϕ2 d π4 + ϕ2 1 1 1 = ϕ ϕ π π 2 tan 4 + 2 cos 4 + 2 2 1 1 1 = = . ϕ ϕ cos ϕ 2sin π4 + 2 cos π4 + 2 sin π2 + ϕ
This is the full solution in the case that e = 0 (solution for the sphere). In the case of the ellipsoid the second integral
ˆ
−
where we designate f
e cos ϕ 1 + e sin ϕ
dϕ =
ˆ f (ϕ) f ( f (ϕ)
dϕ = ln f ( f (ϕ) = ln (1 + e sin ϕ) ,
≡ 1 + e sin ϕ. In the same way
ˆ
−
e cos ϕ 1 e sin ϕ
−
dϕ =
− ln(1 − e sin ϕ) , 109
Appendix A. Isometric Isometric latitude latitude on the ellipsoid ellipsoid and the end result is
− −
π ϕ + + 4 2 π ϕ = ln tan + 4 2
ψ = ln t an an
110
e (ln (ln (1 + e sin ϕ) 2 e 1 + e sin ϕ . 1 e sin ϕ 2
ln(1
− e sin ϕ)) =
Appendix B Useful equations connecting the main radii of curvature When we have given the principal radii of curvature of the ellipsoid of revolution:
− − −
N ( N (ϕ) = a 1
e2 sin2 ϕ
M ( M (ϕ) = a 1
e2
1/2
−
,
e2 sin2 ϕ
1
3/2
−
,
we can calculate by brute-force derivation:
d (N ( N (ϕ)cos ϕ) = M ( M (ϕ)sin ϕ, dϕ d M ( M (ϕ) (N ( N (ϕ)sin ϕ) = + cos ϕ. dϕ 1 e2
−
−
Furthermore
− −
d 2 d M 2 (ϕ) = a2 1 e2 1 e2 sin2 ϕ dϕ dϕ 2 2 2 e M N = 3 sin2ϕ, sin2ϕ, a2 d d N 2 (ϕ)cos2 ϕ = 2N (N ( N (ϕ)cos ϕ) = dϕ dϕ = M N sin2 N sin2ϕ. ϕ.
3
−
=
−
111
Appendix C The Christoffel symbols from the metric Let us start from the definition of the metric in Chapter 9.2 9.2::
gij =
∂ x ∂ x ∂u i ∂u j
·
=
∂ x ∂ x ∂u 1 ∂u 1 ∂ x ∂ x ∂u 2 ∂u 1
· ·
∂ x ∂ x ∂u 1 ∂u 2 ∂ x ∂ x ∂u 2 ∂u 2
· ·
,
where ui = u1 , u2 is the parametrization (“co-ordinate frame”) of a curved surface in space (or more generally a sub-space). sub-space). Differenti Differentiate: ate: ∂ g jk ∂u i
= =
∂ ∂u i
∂ x ∂u j
∂ 2 x ∂u i ∂u j
· ·
Correspondingly, by interchanging indices: ∂ gki = ∂u j ∂ gij = ∂u k
∂ x = ∂u k ∂ x ∂ 2 x + ∂u k ∂u i ∂u k
∂ 2 x ∂ x j k ∂u ∂u ∂u i ∂ 2 x ∂ x ∂u k ∂u i ∂u j
·
·
+ +
·
− ∂u∂ k gij = 2
·
·
∂ 2 x ∂u i ∂u j
Let’s write the second derivatives of x on the local base (9.7 9.7), ), even if in a slightly different notation: ∂ 2 x ∂u i ∂u j
.
(C.1)
∂ 2 x ∂ x j i ∂u ∂u ∂u k ∂ 2 x ∂ x ∂u k ∂u j ∂u i
Compute equation (C.1 (C.1)) plus equation (C.2 (C.2)) minus equation (C.3 ( C.3): ): ∂ ∂ g + gki jk ∂u i ∂u j
∂ x ∂u j
∂ x · ∂u k
(C.2) (C.3)
.
∂ x ∂ x , ,n ∂u 1 ∂u 2
(C.4) as we did in equation
∂ x ≡ Γij ∂u + β ij ij n,
(C.5)
which implicitly defines the Γ symbols. Substitution into equation (C.4 ( C.4)) yields Γij
Here we recognise
· · − · · − ⇒ ⇒ −
∂ x ∂ x ∂u ∂u k
+ β ij ij n
∂ x ∂ x ∂u ∂u k
or
Γij gk = Γij =
∂ x ∂u k
=
1 2
= gk ja
∂ ∂ g + gki jk ∂u i ∂u j
n
∂ x ∂u k
∂ gij . ∂u k
= 0,
1 ∂ ∂ ∂ g + g gij jk ki 2 ∂u i ∂u j ∂u k 1 k ∂ ∂ ∂ g g + g gij , jk ki 2 ∂u i ∂u j ∂u k
equation (10.7 (10.7). ). Here, g ij g jk = δ ki , i.e., g ij is the inverse matrix of gij .
113
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 co-ordinate rectangle. Of the spatial spatial vector v is transported parallelly inside a i u -parametrized surface S , the following holds ∂ v = 0. ∂u i Let us write v on the basis of the tangent vectors: v = vi
∂ x . ∂u i
Then ∂ v 0 = j ∂u
2 ∂v i ∂ x i ∂ x +v = ∂u j ∂u i ∂u i ∂u j ∂v i ∂ x i k ∂ x + Γ v + β ij ij n, jk ∂u j ∂u i ∂u i
= =
∂v i ∂ x i k ∂ x using (C.5 (C.5). ). Here, the v derivativ derivativee consists of two parts: the “interior” “interior”part, part, j i + Γ jk v , ∂u ∂u ∂u i embedded embedded in the surface, surface, and the “exterior” “exterior” part, β ij ij n, aperpendicular to the surface. When the surface S has been given, we can only zero the internal part, i.e. ∂v i i k + Γ jk v =0 j ∂u
(D.1)
describes the parallel transport of the vector v i within the surface. Let us now consider a small rectangle ABCD, ABCD, side lengths ∆u ∆uk and ∆u ∆u (see Fig. 10.3 10.3), ), along co-ordinate curves. The sides AB and C D 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 = ∆uk = k ∂u i
i −Γkm v m ∆uk .
In the same way i ∆CD v i = +Γkm v m ∆uk .
For the side BC we obtain ∆BC v i =
∂v i ∆u = ∂u
i −Γm v m ∆u
and i ∆DA v i = +Γm v m ∆u .
115
Appendi Appendix x D. The Riemann tensor from the Christoffel symbols sum together these four terms: ∆ABCD v i = = =
− − − − − − i Γkm vm
i Γkm vm
CD
AB
∆uk
i Γm vm
DA
i Γm vm
BC
∆u =
∂ ∂ i m k i Γ v ∆ u ∆ u Γm v m ∆uk ∆u = km k ∂u ∂u i i m m ∂ Γkm ∂ Γm m i ∂v i ∂v + Γ Γ ∆u ∆uk . v km m k k ∂u ∂u ∂u ∂u
Equation (D.1 (D.1)) gives
∂v m = ∂u
−
m h Γh v ,
∂v m = ∂u k
−Γkhm vh;
substituting: ∆ABCD v
i
= =
i ∂ Γkm ∂u l
i ∂ Γkj
∂u
−
i ∂ Γlm ∂u k
∂ Γij
− ∂u k
m
v +
i + Γm Γm kj
i m Γlm Γkh
−
i − Γkm Γm j
where we have changed the names of the indices m last two terms).
i m Γkm Γlh
v h ∆u ∆uk =
v j ∆u ∆uk ,
→ j (in the first two terms) ja h → j (in the
Here we see in ready form the Riemann curvature tensor : tensor : i R jk
=
i ∂ Γkj
∂u l
∂ Γij
− ∂u k
i + Γm Γm kj
i − Γkm Γm j ,
i i aapart from the names of the indices and interchanges of type Γ jk = Γkj , just what already was given in eq. (10.9 ( 10.9). ).
116
Bibliography Anon. Anon. (2003). (2003). JHS154 JHS154.. ETRS89 ETRS89 -jarjestelm a¨rjestelma¨an a¨n liittyvat arttaprojektiot, tasokoordinaati tasokoordinaatistot stot ¨at karttaprojektiot, ja karttalehtijako karttalehtijako (Map projections, proj ections, plane co-ordinates and map sheet division in relation to the ETRS89 ETRS89 system) system).. Web site, Finnish Finnish National National Land Survey. Survey. URL URL:: http://www. jhs-suositukset.fi/suomi/jhs154 , accessed August 30, 2005. Anon. (2008). JHS153. ETRS89-j¨rjestelm¨ rjestelm¨ n mukaiset mukaiset koordinaatit koordinaatit Suomessa. Suomessa. Web site, Finnish National National Land Survey. Survey. URL: http://www.jhs-suositukset.fi/suomi/jhs153 , accessed Oct. 19, 2010. Bomf Bo mford ord,, G. (1963 (1963). ). Repor Reportt of Stud Study y Grou Group p No No.. 10. The Geoi Geoid d in Euro Europe pe and Connec Connecte ted d Countries. In Traveau Traveauxx de d e l’Ass l ’Associat ociation ion Interna Int ernatio tionale nale de G´eod´ eod´esie esi e , number 22. Bouch Boucher, er, C. and Altamimi, Altamimi, Z. Spe Specific cifications ations for refer eference ence frame frame fixing fixing in the analys analysis is of a EUREF GPS campaign . Version ersion 4, 24-10-2008 24-10-2008,, http://etrs89.ensg.ign.fr/memo-V7. pdf. Boucher, Boucher, C. and Altamimi, Z. (1995). Specifications Specifications for Reference Reference Frame Frame Fixing in the Analysis of a EUREF GPS Campaign. In VI Meeting of the TWG in Bern, March 9 and 10, Ver off. ¨ ¨ der Bayerischen Kommission f ur ¨ ¨ die Internationale Erdmessung , number 56, pages 265–267, Munchen. u ¨nchen. Gross, R. S. (2000). The excitation of the Chandler wobble. Geophys. Res. Lett., Lett., 27:2329–2332. 27:2329–2332. Grossma Grossman, n, W. (1964). (1964). Geod atische Rechnung Rechnungen en und Abbildung Abbildungen en . Verlag erlag Konrad Konrad Wittwe Wittwer, r, ¨ ¨ Stuttgart. Heikkinen, Heikkinen, M. (1981). Solving Solving the shape of the Earth by using using digital density density models. models. Report 81:2, Finnish Finnish Geodetic Institute, Institute, Helsinki. 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. http://mycoordinates.org/is-definition-of-wgs-84-correct/.
URL:
McCa Mc Cart rth hy, D. D. (199 (1996) 6).. IERS IERS Conv onvention tionss (199 (1996) 6).. (IER (IERS S Techn echnic ical al No Note te ; 21). 21). Technical repo reporrt, Central ral Bureau reau of IERS - Observato atoire de Paris. URL: http://www.iers.org/nn 11216/IERS/EN/Publications/TechnicalNotes/tn21.html. Moritz, Moritz, H. (1992). Geodetic Reference Reference System System 1980. In Tscherning Tscherning,, C. C., editor, editor, Geodesist’s Handbook 1992 , pages 187–192. International Association of Geodesy, Copenhagen. Ollikainen, M. (1993). GPS-koordinaattien muuntaminen Kartastokoordinaateiksi (Transform(Transforming GPS co-ordinates co-ordinates to the National National Grid Co-ordinate Co-ordinate System). System). Tiedote Tiedote 8, Finnish Finnish Geodetic Institute, Helsinki. Ollikainen, M., Koivula, H., and Poutanen, M. (2000). The densification of the EUREF network in Finland. Finland. Publication Publication 129, Finnish Finnish Geodetic Institute, Institute, Kirkkonum Kirkkonummi. mi. Parm, T. (1988). Kansallisen Kansallisen koordinaatti koordinaattijjarjestelm an luominen Suomessa (The creation of a ¨arjestelm¨an national national co-ordinate co-ordinate system in Finland). Finland). Maanmittaus , (1). Stevenson, P. (2008). The Truth About WGS84. URL: http://www.ashgps.com/out3/General
117
