Space Oblique Mercator Projection 

Mathematical Development 

GEOLOGICAL SURVEY BULLETIN 1518 '• 

I 

: ' 

Space Oblique Mercator Projection 

Mathematical Development 

By JOHN P. SNYDER 

GEOLOGICAL SURVEY BULLETIN 

Refined and improved equations that can be applied 

to Landsat and other near-polar-orbiting 

satellite parameters to describe a continuous 

true-to-scale groundtrack and that will contribute 

to the automated production of image-base maps 

1 5 1 8 
UNITED STATES DEPARTMENT OF THE INTERIOR 

JAMES G. WATT, Secretary 

GEOLOGICAL SURVEY 

Doyle G. Frederick, Acting DirectM 

Library of Congress catalog-card No. 81-607029 

For sale by Superintendent of Documents, U.S. Government Printing Office 

Washington, D.C. 20402 
FOREWORD 

On July 25, 1972, the National Aeronautics and Space Administration launched an Earth-sensing satellite named ERTS-1. This satellite, which was renamed Landsat-1 early in 1975, circled the Earth 

in a near-polar orbit and scanned continuous 185-km-wide swaths of 

the surface. Landsat-1 was followed by Landsats-2 and -3 and thus 

established a pattern by which the Earth is to be viewed from space. 

Conventional map projections are based on a static Earth, but Landsat imagery involves the parameter of time and the need for a dynamic 

map projection for its cartographic portrayal. 

During 1973, U.S. Geological Survey personnel conceived a new dynamic map projection which would, in fact, display Landsat imagery 

with a minimum of scale distortion. Moreover, the map projection was 

continuous so that the 185-km swath would involve no zone boundaries 

regardless of its length. This was accomplished by defining the groundtrack of the satellite as the centerline of the map projection. The 

projection was made conformal and named the Space Oblique Mercator (SOM). The SOM was relatively easy to define in geometric 

terms, and NASA began casting Landsat imagery on it early in 1975. 

The mathematical definition of the SOMis, however, a complex matter 

which took several years to resolve. In October 1976, at a symposium 

on "The Changing World of Geodetic Science" sponsored by the Ohio 

State University, I expressed concern that the equations for the SOM 

had not been developed in spite of considerable effort on the part of the 

Geological Survey. I described the task as "a challenge to the academic 

community." John P. Snyder, a chemical engineer, was among the 

audience. Map projections had been Snyder's hobby for some 35 years, 

and since he had found few that he could not mathematically define, 

the SOM presented a challenge. Moreover, Snyder had recently bought 

his first programable pocket calculator and was eager to put it to a real 

test. By May 1977 Snyder had contacted me on this matter, and I was 

more than willing to provide him all possible assistance, even though 

we had just established a contract with John L. Junkins of the University of Virginia to perform the same task. 

By January 1978 both Junkins and Snyder, acting independently, 

had developed mathematical solutions to the SOM. Junkins' solution is 

intended to be more general, for application to non-circular orbits, and 

is quite complex. His map projection is not conformal, although it 

is believed this aspect can be rectified. Snyder, on the other hand, 

III 
IV 

developed a relatively simple set of equations which applies to circular 

orbits (such as Landsat) and provides a rigorous near-conformal map 

projection. 

Snyder's development of an SOM solution is remarkable because he 

had little formal training in the sciences related to map projections, he 

developed and tested his equations using only a pocket calculator, and 

he worked alone and without compensation. He has now developed 

more complex equations for non-circular orbits. 

Snyder's work has been recognized by the U.S. Geological Survey 

through the John Wesley Powell Award which was presented to him 

in 1978. As the award indicates, Snyder "provided the sought-after 

link by which Earth surface data obtained from orbiting spacecraft 

can now be transformed to any of the common map projections. 

This is an essential step in automated mapping systems which we see 

developing." 

Reston, Va. 

July 15, 1981 

Alden P. Colvocoresses 
CONTENTS 

Page 

Forevvord__________________________________________________________________ iii 

Introduction_______________________________________________________________ 1 

Space Oblique Mercator projection for the sphere_____________________________ 3 

Fourier series__________________________________________________________ 15 

Inverse equations __ ___ _ _ _ __ _ ___ ___ _ _ _ ___ ___ ___ __ _ ___ ___ ___ _ _ _ ______ ___ _ 18 

Distortion analysis _______ --------_---------____________________________ 21 

Formulas for the satellite groundtrack___________________________________ 29 

Iteration procedures____________________________________________________ 30 

Radius of curvature analysis-------------------------------------------- 30 

Space Oblique Mercator projection for the ellipsoid--------------------------- 35 

Elements of the ellipse ------------------------------------------------- 36 

Orbital plane on the ellipsoid ------------------------------------------- 39 

Plotted satellite groundtrack on the projection--------------------------- 42 

Scan line projections____________________________________________________ 46 

Transformation of geodetic coordinates---------------------------------- 48 

Inverse equations __ _ _ _ __ _ ______ __ _ ______ ___ __ _ __ ____ ____ _ _ ___ ___ __ _ _ _ _ _ 50 

Correction for vertical scanning----------------------------------------- 55 

Summary of closed equations ------------------------------------------- 58 

Fourier series__________________________________________________________ 61 

Distortion analysis_____________________________________________________ 67 

Incorporation of Landsat rovv and path numbers______________________________ 73 

Summary of simplified equations for the ellipsoidal SOM for Landsat---------- 74 

Modification of ellipsoidal SOM equations for noncircular orbits--------------- 76 

Space Oblique Mercator projection for the general elliptical orbit-------------- 84 

Inverse equations __ ___ _ _ _ __ _______ ______ ______ ___ _______ _ _ ___ ___ _ _ _ __ _ _ 102 

Summary of equations__________________________________________________ 104 

References cited----------------------------------------------------------- 108 

ILLUSTRATIONS 

Page 

FIGURE 1. Orthographic view of the Earth----------------------------------- 5 

2. Element (AB) of an Earth-stationary groundtrack__________________ 8 

3. Space Oblique Mercator projection applied to the sphere 

for tvvo orbits ------------------------------------------------- 14 

4. Exaggerated elliptical cross section of the Earth ------------------- 37 

5. Equatorial vievv of ellipsoid shovving intersection with 

orbital plane ____ __ _ ___ ___ _ _ _ __ _ ___ ___ _ _ _ _________ _ _ _ ___ ___ ___ _ 39 

6. Plane of satellite orbit ------------------------------------------- 40 

7. Octant of ellipsoid with geographic latitude 4> and satellite-apparent 

longitude "-t plotted-------------------------------------------- 41 

8. Octant of ellipsoid vvith orbit------------------------------------- 48 

v 
VI CONTENTS 

Page 

FIGURE 9. Quadrant of plane of satellite orbit for ellipsoid____________________ 49 

10. Correction for vertical scanning on an ellipsoid_____________________ 56 

11. Enlargement of second quadrant of Space Oblique Mercator 

projection applied to the ellipsoid------------------------------- 62 

12. Elements of groundtrack curvature________________________________ 86 

13. Orthographic projection of the ellipsoid____________________________ 89 

14. Groundtrack and extended sample scanlines for Space Oblique 

Mercator projection, elliptical orbit______________________________ 101 

TABLES 

TABLE 1. Scale factors and angular deviation 1 o from Landsat 

groundtrack given by preliminary formulas for the 

Page 

spherical Space Oblique Mercator projection--------------------- 11 

2. Scale factors and angular deviation away from Landsat 

groundtrack given by final formulas for the spherical 

Space Oblique Mercator projection ------------------------------ 26 

3. Location of pole of curvature for Landsat groundtrack 

relative to sphere---------------------------------------------- 34 

4. Scale factors and angular deviation away from geocentric Landsat 

groundtrack for the ellipsoidal Space Oblique Mercator 

projection_____________________________________________________ 72 
SPACE OBLIQUE MERCATOR PROJECTION 

MATHEMATICAL DEVELOPMENT 

By JOHN P. SNYDER 

INTRODUCTION 

In 1974 A. P. Colvocoresses announced a new map projection which 

would permit continuous mapping of satellite imagery, especially 

Landsat. He named it the Space Oblique Mercator (SOM) projection 

(Colvocoresses, 1974). 

Until the SOM, no map projection had been devised which showed 

the ground track of an Earth-orbiting satellite continuously and true to 

scale. In defining the SOM, it was preferred that the areas scanned by 

the satellite be mapped conformally with a minimum of scale error. 

Since the satellite track proceeds essentially along an oblique circumference of the Earth, the oblique Mercator was selected as the closest 

member of the map projection family. 

The mapping problem was simplified by the relatively narrow swath 

scanned by the Landsat (formerly ERTS) satellite at any given time, 

officially 100 nautical miles (about 185 km), but the mathematical 

formulas had not been developed when Colvocoresses presented the 

projection. He offered a geometric analogy using an oscillating cylinder (Colvocoresses, 1975), and made numerous appeals, but adequate 

formulas were not developed until1977. 

In May 1977 a rigorous derivation was initiated voluntarily by the 

writer which was to lead to a satisfactory solution of the SOM 

projection problem by January 1978. The active interest and encouragement shown by A. P. Colvotoresses throughout the work were extremely important factors in its completion. 

Coincidentally, in March 1977 a government contract had been issued to John L. Junkins, then of the University of Virginia, to carry out 

the derivations. His work (Junkins and Turner, 1977) proceeded independently and almost simultaneously with that of the writer. Junkins 

took a more general, theoretical approach-his formulas were much 

more complex, but were designed to be more universal, permitting 

application to any noncircular orbit and finite scanning time. Moreover, Junkins' equations, in present form, do not result in a conformal 

projection. 

1 
2 SPACE OBLIQUE MERCATOR PROJECTION 

During the present work the Landsat orbit was stipulated to be 

circular. Actually it was planned to have a nominal eccentricity of 

0.0001. Its corrected orbit in 1974 had a reported eccentricity of 0.0020, 

which is one-fortieth the eccentricity of the Earth ellipsoid (Colvocoresses, 1974). The Landsat orbit model has also been stipulated to 

have a constant altitude, but this is inconsistent with the circular orbit 

and an ellipsoidal Earth, and would in fact require an orbit with two 

perigees and two apogees, an unlikely and certainly non-Keplerian 

situation; therefore, it was agreed that the circular-orbit stipulation 

be followed. 

It should be stressed that the writer does not claim to have found a 

perfectly conformal projection. It is only a very close approximation, 

becoming closer as the groundtrack is approached. The groundtrack 

itself is perfectly true to scale using either set of formulas. Within 1 o 

of the groundtrack (somewhat farther than the 0.88° limit of the 

Landsat scan), the scale factors using the final formulas below for the 

sphere are within six-millionths of what they should be for a truly 

conformal projection. (For the ellipsoid they are within thirtymillionths of being correct for conformality.) They are thus well within USGS mapping requirements. 

The following derivations are arranged to provide four successively 

more complicated levels of equations, including with each level only 

those features needed. Formulas for the more complex levels may be 

simplified fairly readily to the formulas for any of the more elementary levels: 

1. The first level assumes the Earth to be a sphere and the orbit to 

be circular. This level is recommended only for understanding the 

concepts of the Space Oblique Mercator, since it leads to errors of over 

one-half percent in tropical areas. Formulas are summarized on p. 3. 

2. The second level assumes the Earth to be an ellipsoid and the 

orbit to be circular. This is recommended for normal use with Landsat 

and normal mapping from satellites. Formulas are summarized on p. 

58-61. 

3. The third level assumes the Earth to be an ellipsoid and the orbit 

to be elliptical, but with an eccentricity of 0.05 or less. It is useful for 

more rigorous plotting of Landsat and similar data, although the satellites normally vary slightly from a perfectly elliptical planar orbit, 

giving these formulas more theoretical than practical accuracy. Formulas are summarized on p. 81-82. 

4. The final level assumes the Earth to be an ellipsoid and the orbit 

to be elliptical with any eccentricity. As the eccentricity increases, the 

number of terms for practical computation increases substantially, 

but the formulas provide a correctly curved groundtrack which is true 
PROJECTION FOR THE SPHERE 3 

to scale. Mapping may be accomplished as accurately as the mapping 

for levels 2 and 3 above. Formulas are summarized on p. 104-107. 

In each case, the ellipsoid is assumed to be centered at the center of 

the Earth, and an elliptical orbit is assumed to follow Kepler's laws for 

two bodies. Scanning is assumed to be instantaneous. 

SPACE OBLIQUE MERCATOR PROJECTION FOR THE SPHERE 

The final forward equations for the Space Oblique Mercator projection for the sphere are as follows: 

X - JA' H - s 2 d" s I t ( 11. + lL A.') (1) R - o (1 + S2)112 1\ - (1 + S2)1/2 n an 147r r2'f' 

Y _ JA' s 1 R - (H + 1) 0 (1 + S2)112 dX 

+ (1 + 1

8 2)112 In tan ( 'A1r + V2 4>') (2) 

where S = (P21 P1) sin i cos X' 

H = 1 - (P2/P1) cos i 

tan X' = cos i tan At + sin i tan 4>/ cos At 

sin 4>' = cos i sin 4> - sin i cos 4> sin At 

At = X + (P2/ P1) X'. 

4> = geodetic (or geographic) latitude. 

(3) 

(4) 

(5) 

(6) 

(7) 

X = geodetic (or geographic) longitude, measured east from 

intersection of satellite orbit with Earth's Equator, as satellite crosses in northerly direction (the ascending node). 

P2 = time required for revolution of satellite (103.267 min. for 

Landsat). 

P1 = length of Earth's rotation with respect to the precessed 

ascending node. For Landsat, the satellite orbit is Sunsynchronous; that is, it is always the same with respect to 

the Sun, equating P1 to the solar day (1,440 min.). (Landsat 

is orbited to complete exactly 251 revolutions in 18 days.) 

i = angle of inclination between the plane of the Earth's 

Equator and the plane of the satellite orbit, counterclockwise from the Earth's Equator to the orbi~al plane at 

the ascending node (99.092° for Landsat). 

R = radius of Earth, or of globe at scale of map. 

x, y = rectangular coordinates on the map. 

The significance of S, H, X', 4>', At, as well as calculating techniques 

for the integrals of equations (1) and (2), iteration and quadrant adjust-
4 SPACE OBLIQUE MERCATOR PROJECTION 

ment for equation (5), inverse equations, and equations for scale factors and distortion will all be developed below. Inverse equations are 

numbered (80), (87), (19), (92), and (93). Fourier equivalents are equations (73) through (77), (86), and (90). 

The regular Mercator projection for a stationary sphere is frequently expressed with the following formulas (Thomas, 1952, p. 2): 

X = R~ (8) 

y = R In tan ( ~ 1r + 1h <J>) (9) 

where ~ is measured in this case east from the origin of the x and y 

coordinates, and ~ and <J> are given in radians, as are all angles 

throughout this paper, except where the degree symbol (0 ) or word is 

used. Either radians or degrees may be assumed, of course, if used with 

a trigonometric function such as cos <J>. 

The regular Mercator projection consists of meridians projected 

onto a cylinder tangent to the globe at the equator. The parallels 

cannot be projected directly from the globe, but must be spaced for 

conformality according to equation (9). For the transverse Mercator 

projection, the cylinder is rotated so that it is tangent along a meridian 

from pole to pole. For the oblique Mercator projection, the cylinder is 

rotated so that it is tangent along a great circle which crosses the 

equator at an angle. Direct equations may be provided which show 

little resemblance to equations (8) and (9) (Thomas, 1952, p. 6). It is 

preferable, however, for the sake of subsequent derivations to use 

intermediate transformation equations, substituting ~' and </>' from 

equations below in place of ~ and <J> of (8) and (9). 

Referring to figure 1, an orthographic view of the globe centered at 

the crossing of the transformed equator with the actual equator, the 

equator of the globe CBDA is AB, the regular north pole is C, and the 

regular south pole is D. The transformed equator (inclined i to the 

regular equator) is NG; the transformed poles are Hand J. Angle FOB 

is geographic latitude</>, and EF is the projection of the parallel. Arc 

CMT is a meridian at longitude ~ with respect to central meridian CD, 

so that angle O'C'M' (polar view) or spherical angle OCM (equatorial 

view) equals ~. Point M on the transformed globe has transformed 

latitude</>', represented by KL, and transformed longitude ~', represented by arc HMU, with spherical angle OHM equal to~'. 

In the spherical triangle HCM, side HC is an arc of length i, side CM 

has a length (90° - </>),and side HM has a length (90° - </>').Angle C 

or (HCM) equals (90° + ~),and angle H or (CHM) is (90° - ~').From 

the Law of Sines, as applied to this triangle, 

sin C/sin (HM) = sinH/sin (CM) (10) 
A' E' 

PROJECTION FOR THE SPHERE 

C' 

I 

I 

I 

POLAR VIEW 1 

\ 

........ ~ 

"' \ 

\ 

\ 

I 

c 

\ I 

EQUATORIAL VIEW 

FIGURE !.-Orthographic view of the Earth. 

5 

F' B' 
6 SPACE OBLIQUE MERCATOR PROJECTION 

or sin (90° + A)/sin (90° - 4>') 

= sin (90° - A')/sin (90° - 4>) (11) 

or cos A/cos 4>' = cos A' /cos cJ> (12) 

or cos A' = cos 4> cos AI cos 4>' (13) 

From the Law of Cosines, 

cos (HM) = cos (HC) cos (CM) + sin (HC) sin (CM) cos C (14) 

or cos (90° - 4>') = cos i cos (90° - cJ>) + sin i sin (90° - cJ>) 

cos (90° + A) (15) 

or sin 4>' = cos i sin cJ> - sin i cos cJ> sin A (16) 

Using the Law of Cosines for a different side, 

cos (CM) = cos (HC) cos (HM) + sin (HC) sin (HM) cosH (17) 

or cos (90° - cJ>) = cos i cos (90° - 4>') + sin i sin (90° - 4>') 

or 

or 

cos (90° - A') (18) 

sin cJ> = cos i sin 4>' + sin i cos 4>' sin A' 

sin A' = (sin cJ> - cos i sin 4>')/sin i cos 4>' 

Dividing (20) by (13), then substituting from (16), 

tan A' = (sin cJ> - cos i sin 4>')/sin i cos cJ> cos A 

= [sin cJ> - cos i (cos i sin cJ> - sin i cos cJ> sin A)]/ 

sin i cos cJ> cos A 

= [sin cJ> (1 - cos2i) + sin i cos i cos cJ> sin A]/ 

sin i cos cJ> cos A 

(19) 

(20) 

= sin i tan cJ>/cos A + cos i tan A (21) 

Equations (16) and (21) provide the transformation equations in a 

form suitable for subsequent derivations. This transformed longitude 

A' is the angular distance along the transformed equator, measured 

north from the same intersection of equators from which A is measured. The transformed latitude 4>' is the angular distance from the 

transformed equator, measured positively in the same direction as i, 

namely counterclockwise from the transformed equator, viewed from 

the ascending node. 

If a satellite is following a uniform circular orbit inclined i to the 

Earth's Equator and the Earth is assumed to be a stationary sphere, 

the groundtrack will follow the transformed equator, and A' will be 

directly proportional to the time elapsed since the satellite crossed the 

plane of the Earth's Equator in a northerly direction; that is, the 

ascending node. (This ascending node for Landsat is on the dark side 

of the Earth; imaging occurs as the satellite is proceeding southward 
PROJECTION FOR THE SPHERE 7 

on the opposite, lighted face.) The angular distance along the scan line 

will be </>', measured from the ground track. Scanning by the satellite is 

assumed to be instantaneous throughout these derivations. 

For a stationary spherical Earth with a revolving satellite, these 

formulas are exact. Let us think, however, of the satellite orbit as fixed 

in space, and the satellite revolving and the Earth rotating with respect to the orbit. At the time the satellite has reached X', starting 

from zero, the parallels will not have changed latitude, but the meridians will have rotated so that, from the viewpoint of the satellite in 

space, the actual longitude X will appear along the scanning line at the 

point where some other longitude At would have appeared if the Earth 

had remained stationary. The "satellite-apparent" longitude At is thus 

the geodetic longitude increased by the angle of Earth rotation at the 

time the satellite has reached X', or, as listed in the earlier summary 

(p. 3), 

(7) 

where P2 and P1 are the relative times of orbit and Earth rotation as 

described after the original listing of equation (7). 

The term At should then be substituted for X in transformation 

equations (21) and (16) to account for Earth rotation, to give equations 

(5) and (6) listed earlier: 

tan X' = cos i tan At + sin i tan <t>l cos At 

sin </>' = cos i sin </> - sin i cos </> sin At 

(5) 

(6) 

Since At is a function of X', finding X' for a given</> and X involves trial 

and error (iteration), and X' must be calculated before<!>' is determined 

so that the proper value of 'At may be inserted into equation (6), but 

equations (5), (6), and (7) are exact for the transformation in the case 

of the rotating sphere. Iteration techniques are discussed on p. 30. 

It is not appropriate, however, to substitute X' and</>' from equations 

(5) and (6) in place of the </> and X of equations (8) and (9) for the 

rotating sphere and revolving satellite. This substitution is only suitable for the stationary sphere. If this were done for the rotating 

sphere, the satellite groundtrack would be shown as a straight line and 

the scan lines as a series of parallel lines perpendicular to the groundtrack. Such a groundtrack would have neither conformality nor true 

scale. 

Due to Earth rotation, the direction of the Landsat ground track on 

the Earth is about 86° from that of the scanlines at the Earth's Equator, even though the orbital track is perpendicular to the scanlines in 

space. The groundtrack becomes perpendicular to the scanlines when 

the satellite is closest to the poles (at 180° - 99.092° or about 81 o 

North or South latitude). 
8 SPACE OBLIQUE MERCATOR PROJECTION 

To derive the curved satellite groundtrack, the instantaneous slope 

of the ground track relative to the parallels of latitude on the sphere is 

determined first: 

For the groundtrack, ¢'of equation (6) must be zero, or 

cos i sin c/>0 = sin i cos cf>o sin ~t (22) 

sin ~t = tan cf>o/tan i (23) 

where cf>o is the latitude cf> at the groundtrack. Differentiating, 

cos ~td~t = (sec2cf>o/tan i) dcf>o (24) 

Let i' be the inclination of the groundtrack to a given parallel of 

latitude if the Earth were stationary. From figure 2, showing an element AB of this Earth-stationary groundtrack near latitude cf>o and 

longitude ~' BE = d¢0 and AE = cos cf>od~ if the sphere is given unit 

radius. For the stationary situation~ = ~'so AE = cos cf>od~. Then 

tan i' = d¢0 /cos cf>o d~t (25) 

y~ 

cflo 

FIGURE 2.-Element (AB) of an Earth-stationary groundtrack near latitude 4>o and longitude )., 
PROJECTION FOR THE SPHERE 

Substituting from equation (24), 

tan i' =cos At tan i cos2¢o dAti cos cPo dAt 

=tan i cos cPo cos At 

9 

(26) 

Working further with figure 2, the following elements may be defined: 

AB=dA', element of satellite motion along groundtrack, 

disregarding Earth rotation, during a given time 

interval dt. 

BD = (P21 P1) cos¢0 d A', element of distance resulting from 

change of longitude due to Earth rotation, along a 

given parallel of latitude during dt. 

AD=projection of line AB onto map to account for change 

of longitude from B to D during dt. 

AC=dx l elements of rectangular coordinates for groundCD=dy) track projection AD on SOM map. 

Then AC+CB=AB 

or 

Then 

dx+CB=dA' 

CD=BD sin i' 

dy=(P2/P1) cos c/>osin i'dA' 

CB=BD cos i' 

dx =dA'[1-(P2/Pl) cos cPo cos i'] 

To eliminate i', rearrange equations (23) and (26) as follows: 

tan cPo =tan i sin At 

COS c/>0 =1/(1 +tan2¢0 )

112 

=cos i/(1-sin2i cos2At)112 

tan i' =tan i cos i cos A/ (1-sin2i cos2At )

112 

=sin i COS Atl(1-sin2i cos2At)112 

sin i'=tan i'/(1 +tan2i')112 

(27) 

(28) 

(29) 

(30) 

(31) 

=sin i cos At (32) 

COS i'=(1-sin2i COS2At)112 (33) 

Substituting from (33) into (28), 

dx=dA'[1-(P2/Pl) cos cPo (1-sin2 i cos2At)112] (34) 

Substituting from (30), 

dx =dA'[1-(P2/P1) cos i] (35) 

Integrating, 

(36) 
10 SPACE OBLIQUE MERCATOR PROJECTION 

Letting (4) 

and reentering the radius R of the sphere, 

x=RH"A' (37) 

Substituting from (32) and (30) into (27), 

dy =[(P2/P1) cos i sin i cos "A/(1-sin2i cos2"At)112] d"A' (38) 

Substituting for q, from equation (23) into (5), 

tan "A'= cos i tan At +sin i sin At tan i I cos At 

=tan At· (cos i +sin2i/cos i) 

=tan At ( cos2i + sin2i )/cos i 

=tan "A/cos i (39) 

tan At =tan "A' cos i (40) 

cos At =1/(1 +tan2"A'cos2i)112 (41) 

Substituting from (41) into (38), 

dy =[(P2/P1) cos i sin i/(1/cos2"At -sin2 i)112] d"A' 

=[(P2/P1) cos i sin i/(1 +tan2"A'cos2i -sin2i)112 ] d"A' 

=[(P2/P1) cos i sin i/(cos2i sec2"A')112] d"A' 

=(P2/P1) sin i cos "A'd"A' (42) 

Comparing equation (3), dy=Sd"A' 

Integrating and including R, 

y=R (P2/P1) sin i sin "A' (43) 

Equations (37) and (43) provide coordinates for the groundtrack with 

satellite scanlines running perpendicular to the X axis. The transformed parallels of latitude¢' are at a distance yin equation (9) from 

the groundtrack, so equation (9) must be added to equation (43) for the 

full equation for y. Thus, 

x=RH"A' (37) 

y=R [(P2/P1) sin i sin "A'+ln tan (J:41r+lh¢')] (44) 

where symbols are as given after equations (1) through (7), and His 

found from equation ( 4). · 

Equations (37) and (44) provide a groundtrack which is sinusoidal, 

true to scale, and with conformality (letting ¢'=0). It is inclined about 

4° to the X axis at the crossing of the equator, as the following calculation shows. Dividing equation (42) by (35), 

dyldx=(P2/P1) sin i cos "A'/[1-(P2/Pl) cos i]. (45) 

If "A'=O, at the equatorial crossing, 

dy/dx =(P2/P1) sin i/[1-(P2/P1) cos i ]. 
PROJECTION FOR THE SPHERE 11 

If P2/P1 =18/251, and i =99.092°, as indicated following equations (1) 

through (7), 

dy I dx = 0.07002, 

of which the arctangent is 4.01 o, the angle of slope to the X axis. 

Time is proportional to x, but it is not quite proportional to distance 

along the ground track and should not be because of the effect of Earth 

rotation. 

This change from a straight to a curved transformed equator (or 

groundtrack) is not enough to place scale errors within the satellite 

swath at less than the desired limit of about one part in 10,000 for 

mapping accuracy, even though the groundtrack itself is satisfactory. 

The normal scale factor for a tangent conformal cylindrical projection, 

whether regular, oblique, or transverse Mercator, is sec cJ>', or 1.000152 

at 1 o from the transformed equator (somewhat farther than the 0.83° 

limit of the Landsat scan). Dividing this factor out of the values of h 

and k in table 1 for various points along one quadrant of the projection 

based on formulas (37) and ( 44), and using a scale factor derivation 

following the same pattern as the differentiation discussed later, it is 

found that residual errors are over one part in 1,000 at the polar 

approaches when ¢' = ± 1 o. Since the parallels of latitude run parallel 

to the ground track at the polar approach, it may be noted that k is too 

small at¢' =1 o (away from the north pole) at 'A' =90°, and too large 

at ¢' = -1 o (toward the north pole), and by about the same amount: 

0.998930/1.000152 = 0.998778, while 1.001375/1.000152 = 1.001223, each 

TABLE 1.-Scale factors and angular deviation r l/J' from the Landsat groundtrack 

given by preliminary formulas for the spherical Space Oblique Mercator projection 

[Seanlines drawn parallel. From equations 37 and 44) 

X' 

10 -10 

h k w h k 

oo 1.000154 1.000151 0.0006° 1.000154 1.000151 

15° .999847 1.000142 .0181 1.000462 1.000159 

30° .999566 1.000128 .0350 1.000747 1.000169 

45° .999339 1.000102 .0495 1.000983 1.000187 

600 .999212 1.000035 .0606 1.001131 1.000233 

75° .999350 .999776 .0676 1.001059 1.000428 

900 1.000152 .998930 .0700 1.000152 1.001375 

Note: h =scale factor along meridian of longitude 

k =scale factor along parallel of latitude 

w=maximum angular deviation. 

X'=angular distance along groundtrack from ascending node. 

l/J'=angular distance from groundtrack (positive in direction away from 

north pole). 

On groundtrack l/J'=O, h=k=l.O and w=0.0°. 

w 

0.0006° 

.0181 

.0350 

.0495 

.0607 

.0677 

.0700 
12 SPACE OBLIQUE MERCATOR PROJECTION 

quotient being about 0.001223 from unity. The maximum angular deformation is consequently almost identical each place: 0.0700°. At the 

equatorial crossing (A.' =0°), this discrepancy disappears, so that k is 

less than h by 0.000003, whether Q>' is 1 o or -1 o. 

Therefore, after much conceptualizing, the retrospectively simple 

answer became apparent: the ground track should be bent more sharply on the projection in the polar areas but not in the equatorial areas, 

and the scanlines should continue to intersect the groundtrack at the 

same angles as before to prevent distortion along the ground track. The 

scanlines would thus be skewed with respect to the Y axis. 

By observing the magnitude of the residual scale error at various 

values of A.', it was decided to double the slope of the groundtrack and 

to rotate the scanlines so that they still intersected the ground track at 

the proper angle. This greatly reduced the residual errors, but it was 

found that doubling the angle of the slope instead of doubling the slope 

itself resulted in both simpler formulas and a slightly improved scale 

factor. This form had scale-error residuals (after dividing out sec A.') 

of up to 0.000020. The last slight modification brought the maximum 

residual down to about 0.000006, without complicating the formulas. 

The approach was strictly empirical, with the preceding formulas and 

the desire to retain the exact scale and conformality of the groundtrack already achieved. 

Following is the derivation of equations for x and y retaining the 

proper groundtrack, but altering its slope to the final form. 

By substituting equations (3) and (4) into equation (45), 

dyldx=S/H (46) 

for the ground track of equations (37) and ( 44). Calling the angle of this 

slope()= arctan (S/ H), this angle is to be increased by angle 01 =arctan 

S. Temporarily calling the revised coordinates x' andy', to contrast with 

the x andy of (37) and (44), the new angle of slope, 

Since 

then 

()'=arctan (dy' ldx') 

=0+01. 

tan ()' = (tan ()+tan 01)/ (1-tan () tan Ot) 

dy'ldx'=(SIH +S)/(1-82/H) 

=(S+HS)I(H -82) 

To maintain the same scale along the ground track, 

(dx )2+(dy )2=(dx')2+(dy')2 

Rearranging, 

[1 +(dyldx )

2] (dx )

2=[1 +(dy' ldx')2] (dx')2 

(dx')2= { [1 +(dyldx)2]/[1 +(dy~/dx') ]} (dx)2 

= t [1+(S/H)2 ]/[1+((S+HS)/(H-S2))2]} (dx)2 

(47) 
PROJECTION FOR THE SPHERE 

This simplifies to 

dx'=[(H -82)/H(l +8:2

)

112 ]dx 

Substituting from equation (35) into (48), 

Integrating, 

dx'=[(H -82) H/H (1+82}

112]dA' 

=[(H -82}/(1 +82)

1l

2]dA' 

f>.' x'= Jn [(H -82

)/(1 +8:~} 112 ]dA' 

Substituting equation ~49) into (47), 

dy'=[(8 + H8)/(H -82)]dx' 

=[8 (1 +H) (H -82}/(H -82) (1 +82)

112]dA' 

=[8 (H + 1)/(1 +82}

112 ]dA' 

13 

(48) 

(49) 

(50} 

(51) 

Integrating, 

y'=(H +1) I:[S/(l+S2

)

112]dA' (52) 

For points away from the satellite path, the second term of equation 

( 44) must be used, but corrected for inclination at angle 01 with respect 

to vertical. Since tan 01 = 8, as above, 

sin 81 =tan 81/(l +tan201}

112=S/(l +82)112 

cos 81 =1/(1 +tan201}

112=1/(1 +82)112 

Thus, from equation (50) must be subtracted sin 01 times ln tan 

( 1A 1r+ 1/2<1>'), and to equation (52) must be added cos 01 times ln tan 

( 1A r + '12 </>'). Then the final equations for the new x and y, eliminating 

the prime marks and including the radius R of the sphere, become 

x/R= I~'~ dA'- kIn tan ('Ar+'t.!t/>') (1) 

y/R=(H+l) I~' y&s.dA'+ Vf~s ln tan ('A,.-+'f.!t/>') (2) 

The resulting scale factors, after dividing by sec </>', as calculated 

below, are within 0.000006 of correct for a conformal projection, in a 

zone 1 o on either side of the satellite groundtrack. This is still not 

perfectly conformal, but for the sphere the discrepancy within the 

required scanning range is negligible. 1~he X axis is not now proportional to time, but }..' remains so. For Landsat, the poles are near to, but 

not on, the new X axis. The scanlines vary between ± 4.05° from vertical, at counterclockwise angle arctan 8, and ± 4.01 o (at clockwise angle 

arctan (8/ H)) from normal to the satellite ground track, which now is 

inclined 8.06° to the X axis at the Equator. The 4.01 o is as calculated 

following equation (45), while 4.05° is the arctan of 8=(18/251)Xsin 

99.092°=0.07081. Figure 3 shows a 30° graticule extended to most of 
---SCAN LINES (EXTENDED 30°) 

-GROUNDTRACK 

y 

I 

FIGURE 3.-Space Oblique Mercator projection applied to the sphere for two orbits. Heavy line indicates satellite groundtrack; dashed lines, scan; 

orbit inclination, 99°; period, 103 minutes. 

..... 

~ 

00 

~ 0 

~ 

0 

1:1:1 

t"" 

~ tr:l 

~ tr:l 

~ 

0 

~ ~ 

"'tt 

~ 

~ tr:l 

0 

~ 

0 

z 
PROJECTION FOR THE SPHERE 15 

the globe for two orbits. The progressive rotation of meridians may be 

observed. 

Numerical integration is now required, but this can be reduced to a 

one-time calculation of constants for a Fourier series for any given 

satellite orbit. 

FOURIER SERIES 

The Fourier series for a repeating function may be found thus: 

2~x 4~x 2~nux F(x)=lhao+~ cosT+~ cosT+ .. . +am cos~+ ... 

+" . 2~x+ +b . ·2~nx+ V1 Sln T . . . n Sln -rr- ... 

where 2 IC+T 

ao=T c F (x) dx 

2 Jc+T 2~nux am =T c F (x) cos ~ dx 

2 Jc+T . 2~nx bn =T c F (x) s1n -;y- dx 

and F(x+T)=F(x}, with C and T constants. 

(53) 

(54) 

(55) 

(56) 

(57) 

Taking the first term of equation (1), 

x.IR= J:'[(H -S2)/(l+S2

)

112]d>-.' (58) 

Let 

(59) 

The repeating cycle for F(~'), for use in equation (57}, is T=2~. Let 

C=O. Letting lhao=B, and because of symmetry integrating only to 

~12 and multiplying by 4, equation (54) becomes 

J

r/2 

B=(2/~) 0 Fx(~') d~' (60) 

The term F(x) in equations (55) through (57) should in theory be 

[F(~')- B] to remove the noncyclical portion of x. After first publication of the formulas (Snyder, 1978b), it was realized that B in these 

cases cancels out and may be omitted from equations for am, etc. 

Substituting in equations (53) through (56), 

, _ m 2~m~' " . 2~n~' (61) F;x;(~ )-B+ ~ am cos -2-+ ~ bn Sln -2- m-1 ~ n-1 ~ 

where 

(62) 
16 SPACE OBLIQUE MERCATOR PROJECTION 

1 J211' bn =; 0 

Fxf>..') sin n'A.'d'A.' (63) 

In calculating these constants, it is found that bn, and Gt, as, as, etc., 

are all zero. Therefore, eliminating the unnecessary terms from (61), 

m 

Fx('A')=B+ 2; am cos m'A.' (64) m-2 

where m is even. Equation (58) requires the integration of Fx ('A'), or of 

(64). Since J cos m'A.' d'A.'= ~sin m'A.' (64a) 

we may combine equations (58), (59), (64) and (64a) for the following: 

X m a 

---! = B 'A.'+ 2; __??! sin m 'A.' R m-2 m (65) 

where am is found from equation (62), m is an even number, and B is 

found from (60). It is found that three terms of the series give 10-place 

accuracy for Landsat and other orbits tested. As a multiplier of B, 'A.' 

must be in radians, or the calculated value of B must be multiplied by 

1r /180 to use 'A' in degrees. 

For similar handling of the first term of (2), inserting (H + 1) into 

the integral (it is a constant for the sphere, but a variable for the 

ellipsoid): 

let 

y.fR = J)<H +l)S/(l+S2)

112]dX' 

F11 ("A')=(H + 1)8/(1 +82

)

112 

As in equations (53) and (61), 

m n 

(66) 

(67) 

F11 ('A.')= lhao+ 2; am cos m'A.'+ 2; bn sin n'A.' (68) m-1 n-1 

where (69) 

(70) 

(71) 

This time, it is found that ao, bn and~' a.,~, etc., are all zero. With 

integration as before, combining equations (66), (67), (68) and (64a), 

Y m a 

---! = 2; __??! sin m'A.' R m-1 m (72) 
PROJECTION FOR THE SPHERE 17 

where am is found from equation (70) this time, and m is an odd 

number. It is found that three terms give 10-place accuracy for this 

orbit. 

To recapitulate, the Fourier equations as applied to equations (1) 

and (2) give the following (incorporating (59), (60), (62), (65), (67), (70) 

and (72), and consolidating constants): 

~=BA.'+A2 sin 2A.'+A4 sin 4A.'+ ... 

- k ln tan ( 1A1r+ %~') (73) 

~ = C1 sin A.'+ Ca sin 3A.' +. . . 

+ v'1 ~S2 ln tan ( 14 1r+ %~') (74) 

where 

(75) 

(76) 

(77) 

It is recommended that Simpson's rule be used for calculating any of 

the numerical integrations encountered in the equations in this study. 

In the general forms, y = Jb F(x )dx, the interval from x =a to x =b a 

must be divided into an even number n of equal intervals, each ~x 

wide, and the value of F(x) is calculated for each value of x. Then 

~ y =3 [F(xa) +4F(xa +~X)+ 2F(xa + 2~x) +4F(Xa + 3~x) 

+2F(xa+4~x )+ ... +4F(xb- ~X)+ F(xb)] (78) 

This is increasingly accurate as ~xis reduced, but is far more accurate than trapezoidal integration. In the case of the Landsat orbit and 

probably most others, a sufficiently small interval for sphere or ellipsoid is 9° in A.' for 10-place accuracy in determining all Fourier coefficients, while 15° is a satisfactory interval for 7- to 8-place accuracy, 

which is sufficient for most purposes. 

Another interesting aspect of the Fourier constants for the spherical or ellipsoidal SOM formulas is that while the 21r or 360° calculation is necessary to determine whether the constant is zero or not, 

once this is determined the needed constant may be calculated by 
18 SPACE OBLIQUE MERCATOR PROJECTION 

integrating from oo to goo or 1r/2, instead of to 360°, and then multiplying the answer by 4. This results from symmetry of the quadrants. 

Thus the computer time for calculating the Fourier constants for the 

SOM needs to be only about 1/36, using a /i'A.' of go and one quadrant, 

instead of say 1 o intervals for the entire circumference. 

For the Landsat on the sphere, using P2/P1=18/251 and i=gg,Qg2o, 

the Fourier constants are found to be as follows (showing decimal 

places in excess of those needed): 

B = 1.0075654142 

0.017585333g 

A2=- 0.001881g808 

A,= 0.0000006868 

A6=- 0.0000000004 

Ct= 0.14215g7410 

Ca=- 0.00002g6182 

C5= 0.0000000167 

for >..' in radians 

for >..' in degrees 

INVERSE EQUATIONS 

Equations (1) through (7) give x and y in terms of cf>' and >..', and c/>' 

and >..' in terms of cf> and >... It is also desirable to have the inverse. The 

equations for cf>' and >..' in terms of x and y will be derived first. The 

equations are easily inverted, although Fourier series equivalents are 

the most practical means of calculating with the resulting forms. 

Multiplying equation (2) by S, we obtain 

- I~ s I s 11" ~ ySIR-S(H+l) 0 

yl+S2 d>.. + v1+82 ln tan (4+2) (79) 

Adding this equation to equation (1), the last terms, and thus c/>', cancel 

out; 

(80) 

Since cf>' is eliminated, >..' may be found in terms of x and y from this 

equation, but both numerical integration and iteration are involved, 

resulting in excessive computer time. Thus equation (80) is an almost 

academic equation of little practical use. By changing it to Fourier 

series, however, repeated numerical integration is eliminated, and 

while trial and error (or iteration) is required for the remaining equation, convergence is rapid. To convert to the Fourier series, substitute 

the first portions of equations (73) and (74) in place of the respective 

integrals of (80), multiplying the series of (74) by S: 

X ~yS = B >..' + A2 sin 2 >..' +A, sin 4 )..' + ... + s ( cl sin >..' 

+C3 sin 3 >..' +C5 sin 5>..' + ... ) (81) 
PROJECTION FOR THE SPHERE 

Since S = D cos 'A' , where D = (P2/ P1) sin i, from equation (3), 

x~yS =B 'A' +A2 sin 2 'A' +A4 sin 4 'A'+ ... 

19 

+D cos 'A' (C1 sin 'A' +C3 sin 3 'A' +Cs sin 5 'A'+ ... ) (82) 

Taking a standard trigonometric identity, 

sin a +sin b =2 sin V2(a+b) cos 1

h(a-b) (83) 

let a=(m+n)'A' and b=(m-n)'A' 

Substituting, sin (m +n)'A'+sin (m -n)'A'=2 sin 1h(2m)'A' cos lh(2n)'A' 

or sin m'A' cos n 'A'=lhsin (m+n)'A'+V2sin (m-n)'A' (84) 

Applying (84) to (82), 

where 

x~yS =B'A'+A2 sin 2'A'+A4 sin 4'A'+As sin 6'A'+ ... 

+ V2DC1 sin 2'A'+ V2DC3 (sin 4'A'+sin 2'A') 

+ lh DC5 (sin 6'A' +sin 4'A') + ... 

= B 'A'+ E2 sin 2'A' + E4 sin 4'A' + E6 sin 6'A' +. . . (85) 

E2=A2+ lhD(C1+Cs) 

E4=A4+ lhD(Ca+Cs) 

Es=As+ lhD(C5+C1) 

For the same values of P2/ Ph i, and the Fourier constants as before, 

D = 0.0708121348 

E2= 0.003150288 

E4 = -0.0000003619 

Es= 0.0000000002 (negligible). 

Solution of (85) may then be accomplished by transposing as follows: 

B'A'=(x/R)+(y/R) D cos 'A'-E2 sin 2'A'-E4 sin 4'A'-... (86) 

After substituting the known values of x, y, and R, almost any value 

of 'A' may be tried in the right side of the equation, followed by solving 

for 'A' on the left side and using the latter 'A' for the next trial, etc., until 

there is no significant change between successive trial 'A"s. 

To find 4>' in terms of x andy, it is most convenient to find 'A' first 

from equation (86) and then to substitute it into a rearranged equation 

(2): 

In tan (14r+'hq,')=v/l+S2 [~-(H+l)J."' vf+s. dX'] (87) 

While no iteration is involved in (87), it is preferable to convert to 

Fourier series to eliminate the numerical integration, first converting 

yl +S2 to combine with the integral. Using binomial expansion, 
20 SPACE OBLIQUE MERCATOR PROJECTION 

(1 + 8 2)

112 =[1 + (D cos X')2]112 (88) 

=1 + lhD2 cos2 X'- ~D 4 cos4 X'+M6D6 cos6 X'+ ... 

=1 + lhD2(1h) (1 +cos 2X')- ~D (~) (3+4 cos 2X' 

+cos 4X')+M6D6(l32) (10+15 cos 2X'+6 cos 4X' 

+cos 6X')+ ... 

=(1 + D 4/4-3D'/64+5D6/256)+(D2/4-D 4/16 

+ 15D6/512) cos 2X' +(-D 4/64+3D6/256) cos 4X' 

+ (D 6 /512) cos 6X' + ... 

=Go+G2 cos 2X'+G, cos 4X'+G& cos 6X'+... (89) 

where G0 =1 + D 2/4-3D4/64+(5D6

/256) 

G2=D2/4-D 4/16+(15D6/512) 

G4=-D 4/64+(3D6/256) 

G&=(D6

/512) 

Substituting from equations (88), (89), and (74) into (87), and dropping G6, which is only 0.0000000002: ·--='""""-.,...----,. 

In tan ( 1A 1r + 1h <P') = (y / R) y'1 + D2 cos2 X'-(Go+ G2 cos 2X' 

+G4 cos 4X'+ ... ).(C1 sin X' 

+C3 sin 3X'+C5 sin 5X'+ ... ) 

Using equation (84) to combine the last part, 

In tan (lAr+lh<P')=(y/R) y'1+D2 cos2 X'-[C1G0 sin X' 

where 

+CaGo sin 3X'+C5G0 sin 5X'+ lhC1G2 sin 3X' 

+ 1h C1 G2 sin (-X')+ 1h CaG2 sin 5X' 

+ 1h C3G2 sin X'+ 1h C5G2 sin 7X' 

+ lhC5G2 sin 3X' + lhC1G4 sin 5X' 

+ lhC1G4 sin ( -3X')+ lhCaG4 sin 7X' 

+ lhCaG4 sin (-X')+ lhC5G4 sin 9X' 

+ lhC5G4 sin X'+ ... ] 

=(y/R) y'1+D2 cos2 X'-Lt sin X' 

-La sin 3X'-Lr, sin 5X'-... 

Lt = C1 Go -lh C1G2 + 1h CaG2-( 1h CaG4) + ( 1h C5G4) 

La=CaGo+ lhC1G2+(1hC5G2)-lhC1G4 

Lr,=(C5Go+ lhCaG2+ lhC1G4) 

(90) 

The following values of constants G. and L. for Landsat have been 

calculated. The parenthetical products in the formulas for G. and 4., 

and thus G6 and 4,, may be ignored for 7 -place accuracy: 

Go= 1.001252413 

G2 = 0.0012520218 

G4 = - 0.0000003914 
PROJECTION FOR THE SPHERE 

G6= 0.0000000002 (negligible) 

Lt = 0.1422487716 

La= 0.0000593661 

L5=- 0.0000000296 (negligible). 

21 

The inverses of equations (5), (6), and (7) are also desirable so 

that cJ> and X may be calculated for a given 4>' and X'. 

Equation (19) already provides cJ> in terms of 4>' and X': 

sin cJ> =cos i sin 4>' +sin i cos 4>' sin X' 

Rearranging equation (13), 

cos X= cos 4>' cos X' I cos cJ> 

Rearranging equation (16), 

sin X=(cos i sin cJ>-sin 4>')/sin i cos cJ> 

Dividing (91a) by (91), then substituting from (19), 

tan X=(cos i sin c~>-sin 4>')/sin i cos 4>' cos X' 

=[cos i (cos i sin cJ>'+sin i cos 4>' sin X') 

-sin 4>']/ sin i cos 4>' cos X' 

=cos i tan X'-sin i tan 4>' /cos X' 

(19) 

(91) 

(91a) 

(91b) 

This is the inverse for the stationary globe. As in equation (5), At must 

be substituted for X, so that 

and 

tan At =tan X' cos i -tan 4>' sin i/cos X', 

X=Xt-(P2/P1) X' 

(92) 

(93) 

provides the geodetic longitude from a transposition of (7). Equations 

(92) and (93) require adjustment to place X in the proper quadrant. (See 

step 1 following equation (119).) 

Equations (19), (92), and (93) are the formulas for cJ> and X in terms 

of 4>' and X'. 

DISTORTION ANALYSIS 

To indicate the accuracy of map projection formulas it is important 

that scale factors and angular distortion be determined for representative points throughout the map. This is equally needed on the SOM 

formulas, where conformality is sought but not entirely achieved. The 

fundamental formulas for analysis of scale factors on the sphere involve comparing differential lengths between points on the map with 

the corresponding distance on the globe. 

The distance on the map for a differential change of latitude is: 

(iJs I iJcJ> )m =[(iJx/ iJcl> )

2+ (iJy I iJcJ> )2]112 (94) 
22 SPACE OBLIQUE MERCATOR PROJECTION 

and for a differential change of longitude, 

(asloA)m =[(ax!oA)2+(ay!aA)2]112 

The corresponding distance on the globe, 

(as/ac/>)g =R 

(asloA)g =R cos cJ> 

(95) 

(96) 

(97) 

Dividing (94) by (96) gives the scale factor h along a meridian (constant longitude): 

h = [ (ax 1 ocf> )

2 + ( ay 1 ocf> )2]112 1 R (98) 

The corresponding scale factor k along a parallel (constant latitude) is 

found by dividing (95) by (97): 

k =[(axlaA)2+(aylaA)2]112/R cos cJ> (99) 

A more general scale fa~tor m is found by combining (98) and (99), 

but using finite increments, 

_1 [ (dx)2+(dy)2 ]1/2 

m -R (de/>)2+(cos c/> dA)2 (100) 

This gives the scale factor in any direction, allowing variation in both 

c/> and A. It is useful if small finite increments are given for a map 

projection that is complex to differentiate (especially if applied with 

proper changes to the ellipsoid, as described later). If desired in the 

differential form, the d's can be changed to a's and then divided 

through by oe/> or OA to give calculable expressions. 

In addition to h and k, a third function, the maximum angular 

deformation w first described by Tissot in 1881, should be included. To 

obtain a composite equation for the sphere, we begin with certain 

formulas in Maling (1973), corrected for two errors. 

sin 1

hw=(a -b )/(a +b) (101) 

cos 8'=c:: ~ + :: :~) 1 h k cos cP 

a +b =(h2+k2+2hk sin 8')112 

a -b =(h2+k2-2hk sin 8')112 

(102) 

(103) 

(104) 

where 8' is the angle at which the meridian and parallel intersect at a 

given latitude and longitude, and h and k are determined from equations (98) and (99). Solving for 2hk sin 8' in equation (104) by substituting differentials, letting R =1, since it cancels out, 
PROJECTION FOR THE SPHERE 

2hk sin ()' =2hk V1-cos2 ()' 

[ C

ay ay ax ax)2

/ ]

112 

=2hk 1- oct> oA. +oct> aA. h2k2 cos2 c/> 

=co! c/> {[ c ~:Y +( ~~YJ [ c ~~y +( ~~YJ -[ay ay ax ax]2

} 

112 

ocp oA. + ocp oA. 

This simplifies to 

. '- 2 cay ox ox oy) 2hk sin()-cos c/> ocp oA.- ocp oA. 

Substituting into equation (104), 

(a-b Y=C~:Y+C~:Y+ co;2 c/> [(~~Y+C~~) J 

2 cay ax ax oy)-[ax + 1 cay)J2 

-cos c/> oct> oA.- oct> oA. - oct> cos c/> oA. 

+ [ay 1 cax)J2 oct> cos c/> oA. 

Similarly, equation (103) becomes 

(a+bf=[~: -co! cJ> (~)J +[~~+co! cJ> (~~)J 

23 

(105) 

(106) 

Substituting from (105) and (106) into (101) and multiplying through 

by cos2 c/>, 

. 1 -[(cos c/> ~+~Y +(cos c/> ~-~X]112 

Sln Y2w- ( OX oy)2 ( oy ox)2 (107) 

cos c/> oct>- iJA. + cos c/> oct>+ iJA. 

The tangent form is slightly simpler, using the sine-to-tangent 

identity: 

tan lhw=sin 1hw/(1-sin2 lf2w)112 

[( ax ay)2 ( ay ax)2]1

1

2 

= lh cos c/> a;p+-ax + cos c/> a;p--ax 

(

ay ax ax ay) (

108) 

cos c/> ocp oA.- oct> oA. 

With equations (98), (99), and (108), a fairly complete picture of the 

distortion at any point may be obtained. For a truly conformal 

projection, h, k, (and m) are equal, and w is zero. 

To evaluate these for the final spherical SOM formulas (1) to (7), and 

inverses (19), (73) through (77), (80), (86), (87), (92), and (93), differ-
24 SPACE OBLIQUE MERCATOR PROJECTION 

entiation is carried out so that h, k, and w may be calculated preferably 

for any values of¢' and 'A', since these variables enable us to pinpoint 

the relation of a point to the ground track of zero distortion. Distortion 

as a function of cJ> and 'A may be determined with little additional 

derivation, but it is omitted here as superfluous. 

Differentiating equations (3) and (7), 

dS= -(P2/P1) sin i sin 'A' d'A.'= -Stan 'A.'d'A.' (109) 

d'At =d'A. +(P2/P1) d'A.' (110) 

From equation (5), 

sec2'A.'d'A.'=cos i sec2'Atd'At +(sin i sec2¢/cos'At) d¢ 

+sin i tan cJ> sec At tan'Atd'At (111) 

Substituting (110) into (111), 

sec2'A.'d'A.=cos i sec2'At [d'A+(P~P ) d'A.']+(sin i sec2¢/cos'At)d¢ 

+sin i tan cJ> sec'Attan'At [d'A+(P2/P1) d'A.'] 

Separating by differentials, 

[sec2'A.'-(P2/P1) (cos i/cos2'At +sin i tan cJ> sin 'At/cos2'At)l d'A.' 

= (cos i I cos2'At +sin i tan cJ> sin At I cos2'At) d 'A 

+(sin i/cos2 cJ> cos 'At) d¢ (112) 

Taking partials, 

iJ'A.' sin i 

iJ¢ cos2¢ [cos ~!-pp2 ~ (cos i +sin i tan cJ> sin At)] 

cos 1\ 1 cos 1\t 

(113) 

iJ'A.' cos i +sin i tan cJ> sin At 

iJ'A. cos

2 ~-p 'At P2 ( • + . . t . '\ ) COS 't SID 't an cP SID At 

cos 1\ 1 

(114) 

Differentiating equation (6), and then substituting from (110), 

cos ¢'d¢'=cos i cos cJ> d¢-sin i cos cJ> cos 'Atd'At +sin i sin cJ> sin 'Atdc/> 

=(-sin i cos cJ> cos At) [d'A+(P2/P1) d'A.']+(cos i cos cJ> 

+sin i sin cJ> sin 'At) d¢ (115) 

Taking partials, 

~~=(cos ¢/cos¢') [cos i +sin i tan cJ> sin At 

p2 • • '\ iJ'A.'] - p 1 

SID 't COS At iJcJ> (116) 

iJ¢' _ sin i cos cJ> cos At ( 1 + P2 iJ'A.') 

iJ'A. - cos ¢' P1 iJ'A. (117) 
PROJECTION FOR THE SPHERE 25 

Differentiating equation (1), 

H-82 

[ -S2 

dx/R (l +S2)112 dX'-ln tan ( ~1r+ lh¢') (l +S2)a12 

+ (l +12)112] dS (l +~2)112 sec cp'dcp' 

Taking partials, 

ax _ 1 [ 2 ax' , aq,' s 

Raq, - (1 +S2)172 (H -s ) aq, -s sec q, aq, + 1 +S2 

iJX'J In tan ( ~ 1r + 1h ¢') tan X' iJcp (118) 

The partial ax/ RiJX is the same equation with X in place of q, (at four 

places). 

Differentiating equation (2), 

dy/R = {i~st)})2 dX' (l +~2)a/2ln tan (~1r+ 1h¢') dS+ (l +12}112 

sec ¢' dcp' 

iJy - 1 [ iJX' I iJcp' Raq,- (l+S2)1;2 S(H+l) aq,+sec¢ aq, 

+ 1 !~2 In tan ( ~1r+ lh¢') tan X'~~] (119) 

The partial iJy/Rax is the same equation with X in place of q, (at four 

places). 

These equations permit calculation of distortion for the spherical form 

of the SOM. Assembling the steps, 

1. For a given¢' and X', calculate¢, At, and X from equations (19), (92), 

and (93), making sure that the At calculated from the arctan in equation (92) is placed in the proper quadrant. To do this, subtract from the 

At initially computed (falling between -90° and +90°) this factor: 90° 

times (1 ± 1, taking the opposite sign of cos X') times ± 1, taking the 

sign of the initial At (assume minus if Xt is zero). If cos X' is zero, add 

lo-s degree to X' for the calculation. Thus the subtrahend is 0, ± 180°, 

± 180°, 0, 0, etc., in successive quadrants of the orbit, depending on the 

sign of the initial At. 

2. Calculate S and H from equations (3) and ( 4). 

3. Calculate ( iJX' I iJcp ), (ax' I iJX), ( iJ¢' I aq, ), and ( aq,' I iJX) from equations 

(113), (114), (116), and (117), respectively. 

4. Calculate (axliJcp), (ax/aX), (iJy/iJcp), and (iJy/aX) from equations 

(118) and (119). 

5. Calculate h, k and w from equations (98), (99), and (108). 
26 SPACE OBLIQUE MERCATOR PROJECTION 

Table 2 shows these values for every 5° of a Landsat quadrant at 1 o 

from the ground track, along the ground track, and each 15° of a Landsat quadrant at greater distances from the groundtrack. 

TABLE 2.-Scale facWr's and angular deviati<m away from Landsat groundtrack 

given lYy final farmulas fqr the spherical Space Oblique Mercator projecti<m 

[See notes at end of table] 

At 1° from ground track 

"A' f/)'=10 f/)'=-10 

h k w m •. h k w m •. 

0° 1.000154 1.000151 0.0006° 1.000152 1.000154 1.000151 0.0006° 1.000152 

5 1.000153 1.000151 0.0006 1.000151 1.000154 1.000151 0.0006 1.000152 

10 1.000153 1.000151 0.0006 1.000151 1.000155 1.000151 0.0006 1.000153 

15 1.000153 1.000151 0.0006 1.000150 1.000155 1.000151 0.0006 1.000153 

20 1.000152 1.000151 0.0006 1.000150 1.000156 1.000151 0.0006 1.000154 

25 1.000152 1.000151 0.0006 1.000150 1.000156 1.000151 0.0006 1.000154 

30 1.000152 1.000151 0.0005 1.000149 1.000156 1.000151 0.0005 1.000154 

35 1.000152 1.000150 0.0005 1.000149 1.000156 1.000151 0.0005 1.000154 

40 1.000152 1.000150 0.0005 1.000150 1.000156 1.000151 0.0005 1.000154 

45 1.000152 1.000150 0.0004 1.000150 1.000156 1.000151 0.0005 1.000154 

50 1.000152 1.000150 0.0004 1.000150 1.000156 1.000151 0.0004 1.000154 

55 1.000152 1.000150 0.0004 1.000150 1.000155 1.000151 0.0004 1.000154 

60 1.000153 1.000151 0.0003 1.000151 1.000155 1.000151 0.0003 1.000154 

65 1.000153 1.000151 0.0003 1.000151 1.000155 1.000151 0.0003 1.000153 

70 1.000153 1.000151 0.0002 1.000152 1.000154 1.000151 0.0002 1.000153 

75 1.000153 1.000151 0.0002 1.000152 1.000154 1.000151 0.0002 1.000153 

80 1.000153 1.000151 0.0001 1.000152 1.000153 1.000152 0.0001 1.000153 

85 1.000153 1.000152 0.0001 1.000152 1.000153 1.000152 0.0001 1.000152 

90 1.000152 1.000151 0.0001 1.000152 1.000152 1.000152 0.0000 1.000152 

Beyond 1° from groundtrack 

"A' f/)' h k w "A' f/)' h k w 

oo 15° 1.03567 1.03489 0.142° 450 15° 1.03582 1.03455 0.099° 

10 1.01560 1.01526 0.062 10 1.01563 1.01515 0.044 

5 1.00386 1.00378 0.015 5 1.00386 1.00376 0.011 

-5 1.00386 1.00378 0.015 -5 1.00386 1.00379 0.011 

-10 1.01560 1.01526 0.062 -10 1.01553 1.01539 0.044 

-15 1.03567 1.0348~ 0.142 -15 1.03540 1.03535 0.102 

15 15 1.03574 1.03474 0.136 60 15 1.03575 1.03459 0.071 

10 1.01561 1.01521 0.060 10 1.01562 1.01517 0.031 

5 1.00386 1.00377 0.015 5 1.00386 1.00377 0.008 

-5 1.00387 1.00378 0.015 -5 1.00386 1.00380 0.008 

-10 1.01558 1.01530 0.060 -10 1.01550 1.01542 0.032 

-15 1.03558 1.03505 0.137 -15 1.03531 1.03546 0.073 

30 15 1.03579 1.03461 0.122 75 15 1.03551 1.03483 0.038 

10 1.01562 1.01518 0.053 10 1.01554 1.01525 0.016 

5 1.00386 1.00377 0.013 5 1.00385 1.00378 0.004 

-5 1.00387 1.00379 0.013 -5 1.00385 1.00380 0.004 

-10 1.01556 1.01535 0.054 -10 1.01548 1.01544 0.017 

-15 1.03549 1.03521 0.124 -15 1.03526 1.03552 0.039 
PROJECTION FOR THE SPHERE 27 

Table 2.-Continued 

Beyond 1° from ground track 

X' if>' h k w X' 'if>' h k w 

goo 15° 1.03528 1.03505 0.012° goo -50 1.00382 1.00383 0.000° 

10 1.01543 1.01536 0.004 -10 1.01543 1.01549 0.004 

5 1.00382 1.00381 0.000 -15 1.03528 1.03551 0.013 

Notes: X' = angular position along groundtrack, from ascending node. 

if>' = angular distance away from groundtrack, positive in direction away 

from north pole. 

h = scale factor along meridian of longitude. 

k = scale factor along parallel of latitude. 

w = maximum angular deformation. 

m~, = scale factor along line of constant if>'. 

m,., = scale factor along line of constant X'. 

= sec if>', or 1.000152 at if>'=1 °. 

If if>' = 0°, h, k, and m~,=l.O, while w=O. 

It is also of interest to calculate the scale factor on the SOM at constant cl>' or at constant "A', especially since the equations simplify somewhat in the case of the sphere. 

At constant cJ>', we return to equation (100), letting R =1, and using 

differentials with respect to "A': 

[ 

( ax)2 cay)2 ]1/2 a>..' + a>..' 

m= ac~> 2 2 a>.. 2 

(a>..') +cos cJ> (a>..') 

(120) 

For the numerator of this, we may refer to equations (118) and (119), 

but replacing cJ> in each with "A'. Also let 

V =[S/(1 +S2)] tan "A' In tan ( 1A1r+ 1/2c1>')+ H (121) 

Then from equations (118) and (119), with R = 1, and since ac~>' is zero, 

ax/a>..'=(1+S2)-112 (V-S2

) (122) 

ay/a>..'=(1+S2)-112 S(V+1) (123) 

Combining for the numerator of equation (120), 

(ax/a>..')2+(ay/a>..')2 = (1+S2)-1 [(V-S2

)

2+S2 (V+1)2] 

= (1 +S2)-1 [y2-2VS2+S4+S2V2+2VS2+S2] 

= (1 +S2)-1 [y2 (1 +S2)+S2 (1 +S2)] 

= V2+S2 (124) 

For the denominator of equation (120), it is not correct to use reciprocals of previously derived partial derivatives, since cl>' is not held 

constant in those derivations. Instead, it is necessary to differentiate 

inverse equations (19) and (92). For (92), 
28 SPACE OBLIQUE MERCATOR PROJECTION 

sec ~td~t = sec ~'cos i d~'-(sec cj>' sin i/cos ~')del>' 

-tan 4>' sin i sec ~' tan ~' d ~~ 

Substituting from (110) for d~tJ 

sec ~t [ d ~ + (P 21 P 1) d ~'] = sec ~' cos i d ~~- ( sec2cj>' sin i I cos~') del>' 

-tan 4>' sin i sec~' tan~' d ~~ 

Separating differentials and rearranging, 

d ~/ cos ~t = [cos i I cos ~'- tan 4>' sin i sin ~'I cos ~' 

-(P2/P1)/cos ~t] d~'-(sin i/cos ~' cos2c/>') dq,' 

Then, using partial differentials, 

a~_ cos ~t ( 2. t , . . . '\') P2 ~,- cos2~' cost- an cJ> sin t sin 1\ - p 1 

;~, = -sin i cos ~/ cos ~' cos2cj>' 

Differentiating equation (19), 

cos cJ> dcj>=cos 4>' cos i dcJ>'+cos ~~cos 4>' sin i d~' 

-sin ~'sin 4>' sin i del>' 

Rearranging and using partials, 

(125) 

(126) 

!f,=cos 4>' cos ~'sin i/cos cJ> (127) 

;:,=(cos 4>' cos i -sin 4>' sin i sin ~')/cos cJ> (128) 

While the denominator of equation (120) cannot be readily consolidated, it is found from equations (127) and (125): 

( Oc/>)2

+ 2.-~. ca~)2 = (cos cl>' cos ~'sin i)2+ 2A, [cos ~t ( • o~' cos 'I' o~' cos cJ> cos 'I' cos2~1 cos '1, 

- tan cl>' sin i sin ~')-~:r (129) 

Equations (124) and (129) may be used in equation (120) to calculate m 

for lines of constant cj>', after calculating steps 1 and 2 listed after 

equation (119). Table 2 shows the results for Landsat, using the SOM 

spherical formulas. 

At constant ~', substituting 4>' in place of ~~ in equation (120), 

[ 

( ox )2 ( lJy )2 ]1/2 _ ocJ>' + ocJ>' 

m,..,- ocJ> 2 2 a~ 2 

(ac~>') +cos cJ> (ac~>') 
PROJECTION FOR THE SPHERE 29 

For the numerator, replacing <P with <!>' in equations (118) and (119), 

and since d'A' is zero, and again letting R = 1, 

(ax/8<J>')2+(ay 18</>')2 =(1 +S2)-1 {S2 sec2<f>' +sec2<J>') 

=(1 +S2)-1 (1 +S2

) sec2

<f>' 

=sec2<t>' 

For the denominator, equations (128) and (126) may be used. Actually 

the denominator equals 1, although it takes considerable work on the 

formulas to simplify them. It is more easily verified by making several 

representative calculations on a programable calculator. Thus, 

m,.,=sec <!>' 

in accordance with the actual design of the scanlines, which are lines 

of constant 'A', since the incorporation of equation (9) into (44) and then 

into (1) and (2) changed the direction but not the scale along any given 

scan line. 

FORMULAS FOR THE SATELLITE GROUNDTRACK 

If </>'=0°, previous equations may be readily reduced for calculations 

along the groundtrack. The symbols <Po and 'Ao will be used for values 

of </> and 'A, respectively, along the groundtrack. 

Equations (1) and (2) reduce to the portion determined as equations 

(50) and (52), namely, 

X [' H-82 R = 0 y1 +S2 d'A.' (130) 

~=<H +1) J.' ..;1! 8 , d>..' (131) 

with H and S unchanged from equations (3) and (4). 

Equations (19), (92), and (93) become 

sin <Po= sin 'A' sin i 

'Ao = 'Ato-(P2/ P1) 'A' 

(132) 

(93) 

where tan At0 =tan 'A' cos i (133) 

The latter three equations may be easily inverted to give 'A' for a given 

<Po or 'Ao: 

sin 'A'=sin </>0 /sin i (134) 

tan 'A'=tan 'At/cos i (135) 

where At 0 = 'Ao + (P21 P1) 'A' (136) 

Equations (135) and (136) require combined iteration, since 'A' cannot 

be isolated, just as in equations (5) and (7). This is discussed below. To 

find <Po for a given 'Ao, use equations (135) and (136), then (132). To find 

'Ao for a given </>0 , use equations {134), (133), and (93) in order. 
30 SPACE OBLIQUE MERCATOR PROJECTION 

If )..' should be desired for a given x or y along the groundtrack, 

iteration is also required, although the flatness of the groundtrack 

curve and its periodicity make the calculation from y undesirable. For 

)..' from x, transpose the series part of equation (73): 

B>..'=(x01R)-A2 sin 2)..'-A4 sin 4>..'-... (137) 

Iteration may proceed in the manner described after equation (86). 

ITERATION PROCEDURES 

There are two types of iterations in the equations so far: 

1. Equations (86) and (137) follow the very simple procedure given just 

after equation (86). 

2. Equations (5) with (7) and (135) with (136) may be rapidly converged with the procedure below. Note that unless a procedure such as 

this is used, attempts to calculate )..' automatically will result in quadrant problems and discontinuities near )..' of 90°, 270°, etc., in the 

arctan calculation. 

The equations converge rapidly if, after selecting the desired q, and 

>.., the )..' of the nearest polar approach, )..' P, is used as the first trial )..' 

on the right side of the equation, solving for the corresponding )..' on 

the left side, adding a factor (see below), and using that result as the 

next trial >..'. The )..' P value may be calculated as )..' P =goo X ( 4N + 2 ± 1 ), 

where N is the number of orbits completed at the last ascending node 

before the satellite passes the nearest pole, and the ± is minus in the 

northern hemisphere and plus in the southern (either for the equator). 

(There is actually substantial overlap, so that the minus may be extended almost, but not quite, to the south pole, and the plus almost to 

the north pole; the arbitrary switch of signs at the equator is convenient and covers all situations.) Since the computer normally calculates the arctan as an angle between -90° and 90°, it is necessary to 

add the proper factor. Each)..' given on the left side by the computer 

must be increased by )..'p minus the following: goo times sin )..'P times 

±1 (the ±taking the sign of cos~' where~ =>..+(P2/P1) )..'p). If cos 

)..'Pis zero, the final)..' is )..'p· Thus )..'P is goo, goo, 270°, 270°, 450°, etc., 

and the adder to arctangent is 0°, 180°, 180°, 360°, 360°, etc., for each 

successive quadrant beginning at the origin (q,, >.., q,', and )..'=0). These 

quadrants automatically change along the equator, rather than along 

the scanline which crosses the equator at the node. 

RADIUS OF CURVATURE ANALYSIS 

Junkins and Turner (1g77, p. 18) applied two basic constraints to 

their development of the SOM formulas. First, they maintained the 

groundtrack perfectly true to scale and conformal, as did the present 

writer. Secondly, they made the curvature of the groundtrack on the 
PROJECTION FOR THE SPHERE 31 

map the same as its curvature on the Earth. The writer did not attempt 

to do this until much later and in fact did not think of this constraint 

until Junkins mentioned it. The writer's ground track very closely approaches this concept of curvature, judging from the very small differences in coordinates between the two derivations, but it is as a result 

of empirical reduction of scale error, as previously described. 

Junkins and Turner's report shows larger errors than the calculations based on the present formulas, even in the ellipsoidal forms 

(below), but this appears to be due to errors in some of the concepts for 

developing the formulas for distance from the groundtrack (Junkins 

and Turner, 1977, p. 115-118). The writer's own formulas embodying 

the curvature restraint are derived at the end of this work. Because 

they are much lengthier than the formulas using an empirical curvature, even when simplified for a circular orbit, they are recommended 

only for a non-circular orbit. The increased scale accuracy for a circular orbit appears to be negligible. 

Following is derived the equation for the curvature of the groundtrack on the sphere. It will be rederived for the ellipsoid later and 

incorporated into mapping formulas. The standard formula for radius 

of curvature in plane coordinates is 

rc=[1 +(dyldx)''f12/(d 2yldx2

) (138) 

To convert this to spherical coordinates, it was simplest to use an 

oblique orthographic projection, the azimuthal with the simplest formulas, tangent to the globe at the point where the groundtrack curvature is to be calculated. To follow conventional practice, x and y are 

used in equations (138) through (145), but they are not the SOM coordinates derived earlier. These formulas for a sphere of radius 1.0 are 

as follows: 

x =cos cf> sin A (139) 

y =cos a sin cp-sin a cos cf> cos A (140) 

where a is the latitude at the tangent point, and A is zero. 

Since 

dy/dx =f'(dy/dcp)l(dx/dcp) 

differentiating each with respect to cf>, 

Dividing, 

dx/ d cf> = -sin cf> sin A+ cos cf> cos A d AId cf> 

dy/dcf> = cos a cos cJ>+sin a sin cf> cos A 

+sin a cos cf> sin A dAidcJ> 

(141) 

(142) 

dy _ cos a cos cJ>+sin a sin cf> cos A+ sin a cos cf> sin A (dAidcJ>) 

dx - -sin cf> sin A+cos cf> cos A (dAidcJ>) 

cos a sec A+sin a tan cJ>+sin a tan A (dAidcJ>) 

-tan cf> tan A+(dAidcJ>) (143) 
32 SPACE OBLIQUE MERCATOR PROJECTION 

If a=cJ> and ~=0 at the point of tangency, where the slope on the map 

is the same as the slope on the sphere itself, 

dy cos cJ> + sin2 cJ>/ cos cJ> 

dx d~/dcJ> 

=sec cJ>/(d~/dcJ>) (144) 

For the second derivative, we differentiate equation (143) with respect 

to c/>, and divide by (dx/dcJ>) from equation (141): 

d2y d (dy/dx )/dx dx2 de/> de/> 

[-tan cJ> tan ~+(d~/dcJ>)] [cos a sec~ tan~ (d~/dcJ>) 

+sin a sec2 ¢+sin a sec2 ~ (d~/d¢) +sin a tan ~ (d ~/dcJ> )] 

-[cos a sec ~+sin a tan ¢+sin a tan~ (d~/dcJ>)] 

[ -sec2 c/> tan ~-tan p sec2 ~ (d~/dp)+(d ~/dp )] 

[-tan cJ> tan ~+(d~/d¢)]2 [-sin cJ> sin ~+cos c/> cos~ (d~/dcJ>)] 

If a=cJ> and ~=0, 

(d~/dcJ>) [sin cJ>/cos2 ¢+sin cJ> (d~/d¢)2]-[cos ¢+ 

d2y sin2 ¢/cos p] [-tan p (d~/dcJ>)+(d ~/dp )] 

dx 2 (d~ldc/>)2 cos c/> (d~/dcJ>) 

Multiplying numerator and denominator by cos2¢/cos2c/>, 

(d~/dcJ>) (sin ¢+sin cJ> cos ¢+sin cJ>)-(d ~/dcJ> ) 

(cos3p+sin2p cos p)+sin p cos2p (d~/d¢)3 

COS3c/> (d~/d¢)3 

2 sin c/> (d~/dcJ>)-cos cJ> (d ~/d¢ )+sin cJ> cos2¢ (d~/d¢)3 

COS3c/> (d~/dcJ>)3 

(145) 

Substituting equations (145) and (144) into (138), and inserting subzeros, since only the groundtrack is being considered: 

[1 +sec2cJ>ol(dAoldc/>o)?f12 COS3c/>0 (dAoldc/>0 )

3 

rc 2 sin c/>0 (dAofdc/>0 )-cos c/>0 (d ~jdc/> )+sin c/>o COS2c/>o (dAofdc/>0 )

3 

(146) 

For differentials in (146), we go to groundtrack formulas. Combining 

(133) and (93), 

Ao=arctan (cos i tan ~')-(P /Pl)~' (147) 

Differentiating, 

dAold~'=cos i sec ~' /(1 +cos2 i tan ~')- P2/P1 (148) 

Differentiating equation (134), 

cos ~'(d~'/dcJ>o)=cos c/>0 /sin i (149) 
PROJECTION FOR THE SPHERE 33 

Multiplying (148) by transposed (149), 

(dXoldcJ>o)=[cos i sec2X'I(1+cos2 i tan2X')-P2IP1] cos cJ>jsin i cos X' 

=[cos il(cos2

X'+cos2 i sin2

X')-P21P1] cos ¢jsin i cos X' 

=[cos il(1-sin2 i sin2X')-P21Pd cos ¢jsin i cos X' 

Substituting from equation (134) for X', 

(dXo/dcJ>o)=(cos ilcos2 cl>o-P21P1) cos ¢/(sin2 i-sin2 c/>0 )

112 (150) 

Differentiating this, 

(cos ilcos2cJ>o-P21P1) [(sin2 i-sin2

¢0 )

112 (-sin cl>o)-

tL ( • 2 • • 2 ) -1/2 ( 2 • )] d2Xold¢ 2= 12COSpoSin t -sin .Po . ~- Sin c/>0 cos </>0 

0 Sin2 t - Sin2 cl>o 

+ ( . 2 ~os 2 <I> ) 172 (2 cos i sec cl>o sec cl>o tan cl>o) sin t -sin o 

1 {[ . (. 2• • 2 ) 

( . 2 · . 2A.. )3;2 -sin cl>o Sin t -sin cl>o sin t -sin .,o 

+sin c/>0 cos2 cl>o] (cos i I cos2c/>o-P 2l P 1) 

+2 cos i sin c/>0 sec2¢o (sin2 i -sin2 c/>0 )} 

sin c/>0 cos i [ 2 .

1 2 (P IP ) . 

( • 2 · • 2 A. )3/2 COS t COS c/>0 - 2 1 COS t Sin t -sin .,0 

+ 2 sin2 i I cos2 cl>o- 2 tan2 cl>o] 

sin cl>o cos i [2 2 .

1 2 (PIP) "] ( • 2 · • 2 A. )3/2 -COS t COS c/>0 - 2 1 COS t Sin t -sin .,0 

(151) 

Equations (150) and (151) may be substituted in equation (146) to give 

the radius of curvature rc in terms of groundtrack latitude c/>0 • Without 

giving the intermediate steps, the combination can be somewhat simplified to the following: 

r =cos2 Po [1-2 (PiP1) cos i+(cos cl>o PiP1)Z'f12 

c sin cl>o {(cos i I cos2 cl>o-P 21 P 1) [2( cos2 <l>o- cos2 i) 

+(cos i -cos2 ¢~ P21P1)2]-cos i (2-cos2 ilcos2 <l>o 

-cos i P2IP1)} (152) 

Actually, rc=tan Zo, (153) 

where Z0 is the polar distance of the small circle with radius of curvature rCJ· that is, the pole of the groundtrack curve is Zo away from it. It 

is also desirable to know the latitude and longitude of this pole. This 

involves determining the direction of the groundtrack at a given point 

and moving a distance Zo perpendicular to it along a great circle path. 

For the direction, or azimuth, the direction Zg east of north of the 

groundtrack is determined from the fundamental relationship, 

tan Zg =cos cl>o dXoldcJ>o 
34 SPACE OBLIQUE MERCATOR PROJECTION 

The azimuth of the momentary pole of the circle of curvature from the 

groundtrack is the perpendicular Z, east of north: 

Z=arctan [cos cl>o (d~/dcp )]+90° 

From equation (150), 

Z=arctan [(cos i/cos2 cp0 -P2/P1) cos2 cp0 /(sin2 i-sin2 c/>0 )

112]+90° 

=arctan {[cos i -(P2/P1) cos2 cl>oY(cos2 cl>o -cos2 i)112 +90° (154) 

If ( c/>t, X1) are coordinates of the pole, the standard formula for the 

great circle distance of the groundtrack is written 

cos zo=sin cl>o sin cp1+cos cl>o cos c/>1 cos (~-X ) (155) 

If a spherical triangle is drawn with the geographic pole, the pole of 

the circle of curvature, and the point on the groundtrack as vertices, 

the Law of Sines provides an azimuth formula as follows: 

sin (~-X )/sin zo=sin Z/cos c/>1 (156) 

Since c/>o, ~' Z, and Z0 are known or determined from equations previously given, we have two equations (155) and (156) in two unknowns. 

To eliminate (~- X1), rearrange (155) and (156) as follows: 

cos (~_X ) cos Z0 - sin c/>0 sin c/>1 (

157) 1 COS c/>0 COS c/>1 

sin (~-X )=sin Z sin zjcos c/>1 (158) 

Squaring each and adding the two equations gives 1.0 for the left side: 

1 = (cos Z0 - sin c/>0 sin c/>1)

2 I cos2 c/>0 cos2 c1>t + sin2 Z sin2 zj cos2 c/>1 

Solving for c/>1 leads to the following: 

sin c/>1 =sin c/>0 COS Z0 - COS c/>0 sin Z0 COS Z (159) 

Equation (158) may be used to find X11 after determining c1>t from (159). 

The (~- X1) must fall between 90° and 180° for a first-quadrant portion of the ground track. Table 3 shows the interesting pattern of these 

points for the Landsat groundtrack, with i=99.092° and 

P2/ P1 = 18/251. Each quadrant is symmetrical. 

TABLE 3.-Locaticm of pole of curvature for Landsat groundtrack relative to sphere 

[i==99.092°; P2"P1=18/251. Each quadrant symmetrical] 

tPo Ao A' zo z tPt AI 

oo 0.0° 0.0° -90.0° 76.90276° 13.09724° 90.0° 

10 - 2.34332 10.12857 -88.60379 76.82351 13.21777 88.57343 

20 - 4.79247 20.26557 -87.25038 76.56102 13.56506 87.23070 

30 - 7.48298 30.42181 -85.98117 76.03133 14.09788 86.04627 

40 -10.62969 40.61451 -84.83482 75.05302 14.75274 85.07670 

50 -14.64319 50.87686 -83.84607 73.23903 15.45126 84.35388 

60 -20.48721 61.28790 -83.04467 69.64980 16.10931 83.88063 
tPo 

PROJECTION FOR THE ELLIPSOID 

-31.25479° 

-71.32738 

-96.45418 

X' 

72.11023° 

85.81878 

90.0 

* <f>o=180°-i=80.908°. 

TABLE 3.-CONTINUED 

-82.45460° 

-82.09342 

-82.07238 

z 

61.25029° 

24.20186 

.00000 

35 

16.64701° 83.62960° 

16.99861 83.54691 

17.01962 83.54582 

Note: <f>o=geodetic latitude along groundtrack. 

>.c,=geodetic longitude along groundtrack, relative to longitude of ascending node. 

'A'= transformed longitude along ground track. 

Zc,=polar distance of small circle with radius of curvature of groundtrack. 

Z =azimuth of radius of curvature of groundtrack, east of north. 

4>1 =geodetic latitude of pole of radius of curvature of ground track. 

'A1 =geodetic longitude of pole of radius of curvature of groundtrack, relative to 

longitude of ascending node. 

The Space Oblique Mercator projection may be more correctly considered equivalent to an oblique Lambert Conformal Conic projection 

with one constantly changing standard transformed parallel of distance Z0 from a constantly shifting transformed pole with coordinates 

(c/>11 ).1). This concept results in a substantial reduction (by a factor of 

10 to 500 in various tests) in non-conformality for the sphere. With the 

shifting of the transformed pole, true conformality is still not 

achieved. Since the scanlines are skewed with respect to the groundtrack, they are not radii of the conic projection, nor are they truly 

straight using this concept. The formulas involve a lengthy iteration to 

find the radius of curvature and related parameters for each latitude 

and longitude, in addition to using Fourier constants for calculations 

of x and y. The considerably expanded complexity merely for the 

sphere does not appear justified, since the empirical approach provides 

such an accurate projection. Therefore, the rest of the formulas for the 

conic interpretation of the SOM for the sphere are not given in this 

work, and the concept has not been applied to the ellipsoid. 

SPACE OBLIQUE MERCATOR PROJECTION 

FOR THE ELLIPSOID 

Accurate though the above formulas may be for the sphere, there are 

errors of over one-half of a percent in using these (or any accurate 

spherical formulas) in place of formulas based on the ellipsoid, especially in tropical latitudes (Snyder, 1978a). In maps of very large 

areas, errors fundamental to plane projection far outweigh the effect 

of the ellipsoid, but for topographic mapping of small areas or strips, 

the use of the ellipsoid is essential for high-quality mapping. Therefore, as indicated earlier, the development of SOM formulas for the 
36 SPACE OBLIQUE MERCATOR PROJECTION 

sphere served largely to avoid even more complex derivations and 

calculations in determining the feasibility of various concepts. 

At first, the writer assumed that Hotine's classic work on the ellipsoidal oblique Mercator would be the logical link in changing from 

sphere to ellipsoid (Hotine, 1946, 1947). Indeed, John B. Rowland (written commun., 1977) of the Geological Survey has applied it in five 

stationary zones to approximate each north-to-south pass of the satellite, with consequent discontinuities. Hotine pointed out that an 

ellipsoidal version of the oblique Mercator cannot be exactly developed 

(both the regular and transverse Mercator projections may be derived 

exactly for the ellipsoid), but he prepared an approximation by using 

an "aposphere" tangent to the ellipsoid at a chosen point. The scale of 

his central line is true only at the point of tangency. Within several 

degrees of the tangent point, the central line is essentially true to 

scale, but it gradually deviates. The projection is exactly conformal, 

however. For normal usage of the ellipsoidal oblique Mercator, this is 

completely satisfactory. In the case of a revolving satellite, the increasing error of the central line soon becomes unacceptable. The 

curve presented by the plotting of the groundtrack presents an additional serious problem in attempting to adapt the Hotine to the SOM 

for the ellipsoid. 

Thus it was soon evident that it was better to disregard the Hotine 

derivations in developing the ellipsoidal SOM formulas, and to work 

instead from the basic geometry of the ellipse and ellipsoid, with 

assistance from Thomas (1952). It will be assumed that the center of 

the ellipsoid coincides with the center of the Earth's mass, slightly 

incorrect for most official ellipsoids as used with local datums, but 

accurate for satellite-determined ellipsoids such as GRS 1980. 

ELEMENTS OF THE ELLIPSE 

In the ellipse of figure 4, a is the major semi-axis AC or CE, b is the 

minor semi-axis BC or PC, and the eccentricity e, 

e =(1-b2/a2

)

112 (160) 

The geocentric latitude for point L is c/>g, with the north pole at P and 

the equator along AE. 

The geographic latitude cl>n is the inclination to the equator of a 

perpendicular to the surface at point L. Radii shown to point L are Pg 

or CL, the radius to the center; Pm, the radius of curvature of meridian 

PE at point L; and p, the radius of curvature of the surface in a plane 

perpendicular to the plane of the meridian PE, and also perpendicular 

to a plane tangent to the Earth at L. Both Pm and Pp lie along LD, 

although they are unequal. In much of the literature, Pm and Pp are 

given the symbols R and N, respectively, but the latter are not used 

here, to avoid confusion with other symbols. 
PROJECTION FOR THE ELLIPSOID 37 

p 

B 

FIGURE 4.-Exaggerated elliptical cross section of the Earth. 

The basic equation for the ellipse (using conventional x andy coordinates, which are not those for the SOM map, in equations (161) 

through (167)) 

is differentiated: 

Since, from (161), 

then 

(2xla 2)dx +(2ylb2)dy =0 

dy/dx= -xb2/ya 2. 

Y =b(a2-x2)112fa 

dy/dx = -xb/a(a2-x2)

112 

(161) 

(162) 

This is the slope of the ellipse, and is therefore the negative reciprocal 

of the slope of the perpendicular DL, or 

dyldx= -1/tan cl>n 

Therefore, 

tan cl>n =a(a 2-x2

)

1121bx 

For geocentric coordinates: 

tan cPu =ylx =b(a 2-x2

)

112/ax 

(163) 
38 SPACE OBLIQUE MERCATOR PROJECTION 

Substituting from (163) and (160), 

tan cPu =b2 tan cJ>,Ja 2 

=(1-e2) tan cPn (164) 

Solving for x from equation (163), 

b2x2 tan2 cPn =a2(a2-x2) 

x2(b2 tan2 cPn +a2)=a 4 

x =a/(b2 tan2 cJ>,./a 2+1)112 

=a/[(1-e2) tan2 cl>n + 1]112 

=a cos cJ>,./(1-e2 sin2 c/>,.)112 (165) 

Solving for y from equations (162) and (165), 

y =b [a 2-a 2 cos2 cJ>,./(1-e2 sin2 c/>,.)]112/ a 

=b [1-cos2 cJ>,./(1-e2 sin2 c/>,.)]112 

=b sin c/>,. [(1-e2)/(1-e2 sin2 c/>,.)]112 

=a(1-e2) sin cJ>,./(1-e 2 sin2 c/>,.)112 (166) 

For geocentric radius Pu, using the Pythagorean theorem and equations (165) and (166), 

Pu =(x2+y2)112 

=a , /cos2 cPn +(1-e2

)

2 sin2 c/>,. 

V 1-e2 sin2 cPn 

=a , /1-e2 (2-e2

) sin2 c/>,. 

V 1-e2 sin2 cl>n (

167) 

It is also useful to have Pu in terms of c/>(/' 

Multiplying numerator and denominator of (167) by sec c/>,., 

=a , /sec2 c/>,. -e 2 (2-e2

) tan2 cl>n 

Pu V sec2 c/>,. - e 2 tan2 cPn 

=a , /1 +tan2 c/>,. (1-e2)

2 

V l +tan2 c/>,. (1-e2) 

Substituting from (164), 

_ y 1+tan2 cPa 

Pu -a1 +tan2 c!>u1(1-e 2) 

=a[(1-e2)/(1-e2 cos2 cl>u)]112 (168) 

For radii Pm and Pp the derivations are given in several standard works 

and will be omitted here (see Thomas, 1952, p. 59, eqn. 162, 154). 

Pm =a(1-e2)/(1-e2 sin2 c/>,.)312 (169) 

Pp =a/(1-e2 sin2 c/>,.)112 (170) 
PROJECTION FOR THE ELLIPSOID 39 

ORBITAL PLANE ON THE ELLIPSOID 

Referring to figure 5, the satellite orbit cuts the Earth with an 

inclination i (99.092° for Landsat) at the equatorial crossing. The cross 

section in the orbital plane is also an ellipse having the same major 

semi-axis a, but a different minor semi-axis r (fig. 6). If the satellite 

orbit is assumed to be circular, the geocentric angle X', as measured 

from the equatorial crossing, will be changed at a uniform rate with 

respect to time. Therefore, the elements of the orbital plane are as 

follows: 

POLE 

FIGURE 5.-Equatorial view of ellipsoid showing intersection with orbital plane. 

The minor semi-axis r in figure 6 or figure 5 (where r = AB or AC) 

is the same as pg of figure 4, if </Jg is made equal to i Therefore, using 

equation (168), 

(171) 

Again using equation (168), but calculating the radius Po in figure 6 to 

the surface of the Earth in terms of X' and the ellipse of dimensions a 

and r, we may use (1-r2/a 2

) in place of e 2 (see equation (160) ): 
40 

Then 

SPACE OBLIQUE MERCATOR PROJECTION 

Po=a(r/a)/[1-(1-r 2/a 2

) cos2 X']112 

h'=po sin X' 

=a(r/a) sin X'/[1-(1-r 2/a 2

) cos2 X']112 

FIGURE 6.-Plane of satellite orbit. 

(172) 

(172a) 

(173) 

Relating h' to rectangular coordinates x' and z' in the plane of figure 

5, following the usual convention for the vertical axis of the ellipsoid, 

for points on the satellite groundtrack BC at a given time angle X', 

x'=h' cos i (173a) 

=a(r/a) sin X' cos i/[1-(1-r/a2

) cos2 X']112 (174) 

z'=h' sin i 

=a(r/a) sin X' cos i/[1-(1-r2/a2)cos X']112 (175) 

Taking the same octant (upper left) as that used in figure 5, let us 

plot geodetic (or geographic) latitude cf> (the same as cf>.n) and satelliteapparent longitude Xt, using equations (165) and (166) for the outer 

meridian, but multiplying x' by sinXt to give coordinates for any intermediate point (fig. 7). Longitude Xt is measured from the central meridian, and is calculated from geodetic longitude from equation (7), 

just as for the sphere. 

x' =a cos c1> sin X/ (1-e2 sin2¢ )

112 (176) 
PROJECTION FOR THE ELLIPSOID 41 

POLE 

FIGURE 7.-0ctant of ellipsoid with geographic latitude 4> and satellite-apparent longitude X, plotted. 

z'=a(1-e2) sin cf>/(1-e 2 sin2¢)112 

(Note that x' is negative, but so is At) 

(177) 

By equating (174) to (176) and (175) to (177) and solving simultaneously, cf> and X may be determined at various places along the 

groundtrack. This will be done later. To determine the true distance 

along the satellite groundtrack on the ellipsoid as a function of X', but 

for a stationary Earth (to be corrected shortly for rotation), we return 

to figure 6 and apply equation (169) using (1-r2

/a2

) for e2 and cf>'n as 

a fictitious perpendicular equivalent of X', just as cl>n relates to cf>u· Let 

us call the elemental distance along the groundtrack ds. 

Then ds/dcf>'n =pm adapted to the orbital plane, or 

dsldct>'n =a[1-(1-r2/a2)]/[1-(1-r2/a 2

) sin2 cf>'n]312 

=a(r2/a 2)/[(1-t(1-r2/a 2

) sin2 c/>' nf12 (178) 

Applying (164) to figures 4 and 6, 

tan X'=(r2/a2

) tan cf>'n (179) 

Equation (179) may be rearranged to give 

sin cf>'n=sin X'/[1-(1-r4fa 4

) cos2 X']112 (180) 

Differentiating (179), 

sec2 X' dX'=(r 2/a2

) sec2 cf>'n dct>'n 

dct>'nldX'=cos2 ct>'nl(r2/a2

) cos2

X' 
42 SPACE OBLIQUE MERCATOR PROJECTION 

Substituting from (180) for c/J'n, 

dct>'nld~'=(r /a )/[1-(1-r /a ) cos ~'] 

Multiplying this equation by (178), and substituting from (180) into 

(178), 

[ (1-r2

/a2

) sin ~' ]a12 , 

dsld~'=a(r /a )1 1 1_(1-r4/a4) cos ~, [1-(1-r4fa4

) cos ~] 

(181) 

PLOTTED SATELUTE GROUNDTRACK ON THE PROJECTION 

Referring again to figure 2, but defining the elements with respect to 

the ellipsoid, 

AB = ds, element of satellite motion during a given time interval 

proportional to angular satellite motion d~'. Earth rotation is 

disregarded along line AB. 

BD=(P~ ) a cos cf>od~' 1(1-e sin2c/>0 )

112, element of distance resulting from change of longitude along a given parallel of latitude 

cf>o during d~', using equation (170) multiplied by cos cf>o to give 

distance along the parallel, and by (P2/P1) d~' to give distance 

during satellite motion d~'. P2 and P1 are defined after equations (1) through (7). 

AD=projection of line AB onto map to account for change of longitude from point B to D during time of satellite travel through 

d~'. 

AC=dx l elements of rectangular coordinates for 

CD=dy groundtrack projection AD on the SOM map. 

Following a derivation similar to that for equations (27) through (45), 

but applicable to the ellipsoid, 

AC+CB=AB 

dx+CB=ds 

CD=BD sin i' 

where i' is, as for figure 2, the inclination of the groundtrack to a given 

parallel of latitude if the Earth were stationary. 

dy=(P2IP1) a cos cf>o sin i' d~'/(1-e2 sin2 c/>0 )

112• (182) 

CB=BD cos i' 

dx=ds-(P2/P1) a cos cf>o cos i' d~'/(1-e2 sin2 c/>0 )

112• (183) 

To determine i', using equations (169) and (170), 

tan i'=BE/ AE 

=pm dc/>of Pp COS c/Jo d~t 

=(1-e2

) dct>of(1-e2 sin2 cf>o) cos c/Jo d~t· (184) 
PROJECTION FOR THE ELLIPSOID 43 

Equating (175) to (177), 

z'la=(rla) sin A' sin i/[1-(1-r2/a2) cos2 A']112, and also 

=(1-e2) sin Q>0 /(1-e2 sin2 l/>0 )

112. (185) 

Equating (174) to (176), 

x'la=(rla) sin A' cos i/[1-(1-r2/a2) cos2 A1112, and also 

= cos l/>o sin Atl(1-e2 sin2 l/>0 )

112. (186) 

Rearranging (185) and also (186), several common terms may be moved 

to the left side of each: 

(ria) sin A'(1-e2 sin2 Q>0 )

112/[1-(1-r2/a2) cos2 A']112 

=(1-e2) sin Q>/sin i 

and also =cos l/>o sin At! cos i 

Equating the two right sides of (186a), 

(186a) 

sin At =tan l/>0 (1-e2)/tan i (187) 

Differentiating (187), 

cos At dAt =(1-e2

) sec2 l/>o dQ>/tan i 

and substituting into (184), 

tan i'=cos At tan i cos Q>/(1-e2 sin2 l/>0 ) (188) 

Converting to the sine function, 

sin i' = cos At tan i cos l/>ol[ (1- e2 sin2 l/>0 )

2 + cos2

At tan2 i cos l/>0]

112 

Substituting for At from (187), after a few steps we obtain 

. . ysin2 i - cos2 i tan2 Q> (1- e2)

2 

sin t'= · o 1-e4 sin2 l/>o cos2 i (189) 

Substituting into equation (182), 

dy =a P2 , /sin2 i -sin2 Q>o[1-e2 cos2 i (2-e2)] (

190) 

dA' P1 V (1-e2 sin2 l/>o)(1-e 4 sin2 l/>o cos2 i) 

Reverting to (185) to solve for l/>o in terms of A', first let 

yA =(ria) sin A'/[1-(1-r2/a2) cos2 A']112 (191) 

Then, from (185), 

(1- e2) sin l/>o -

sin i(1-e2 sin2 l/>0 ) 112 yA 

or sin l/>o= yA sin i/[(1-e2)

2+ Ae2 sin2 i]l12 (192) 

Substituting from (171) to simplify A, 

A 1-e2 . 2A'1[1 (1 1-e2 ) 2A'] 1-e2 cos2 i sin - 1-e2 cos2 i cos 

=(1-e2) sin2A'I[1-e2(1-sin2 i sin2 A')] (192a) 
44 SPACE OBLIQUE MERCATOR PROJECTION 

Substituting into (192), this reduces to 

sin <J>o=sin i sin X' /[(1-e2)

2+e2 sin2 i sin2

X'(2-e 2)]112 (193) 

The tangent form is simpler: 

tan <J>o=sin i sin X' /(1-e2)(1-sin2 i sin2X')112 (194) 

Using an auxiliary angle, such that 

sin 1f=sin i sin X', (194a) 

this simplifies still further: 

tan <J>o=sin 1f/(1-e2)(1-sin2 1f)112=tan 1f/(1-e2

) (194b) 

Comparing this to (164), it is apparent that 1f is the same as </>g, the 

geocentric latitude. Therefore, 

tan <Po= tan </> g/ (1- e2

) (195) 

where 

sin <J>g =sin i sin X' (196) 

Substituting (193) into (190) and carrying out several laborious steps 

of algebra and trigonometry, a harmonious exact relationship may be 

developed, using certain operators which are constants for a given 

orbit: 

~~,=~~sin i cos X'[(1 + T sin2X')/(1 + W sin2 X')(1 +Q sin2 X')]112, 

(197) 

where T=e2 sin2 i (2-e2)/(1-e2

)

2 (198) 

Q=e2 sin2 i/(1-e2

) 

W =[(1-e2 cos2 i)/(1-e2)]2-1 

(199) 

(200) 

This is the ellipsoidal version of equation ( 42) for the sphere, and may 

be given the symbol S, corresponding to equation (3). That is, 

S=dyladX' (201) 

Similarly equation (189) may be converted to the cosine form, 

., cos i (1-e2 sin2 </>0 ) 

cos ~ cos </>o(1-e4 sin2 <Po cos2 i)172 (

202) 

and substituted into equation (183) to give 

dx _ ds P2 cos i(1-e2 sin2 </>0 )

112 

dX'-dX'-a P1 (1-e4 sin2 <Po cos2 i)112 (

203) 

Substituting from equation (193) for sin <Po and (181) for ds/dX', equation (203) also may be harmonically simplified after many laborious 

steps to 

dx _ , I 1 + Q sin2 X' [ 1 + W sin2 X' P2 • ] 

ad X'- V 't + W sin2 X' (1 +Q sin2 X')2 P 1 cos~ (204) 
PROJECTION FOR THE ELLIPSOID 45 

where W and Q are as in equations (199) and (200). This corresponds 

to equation (38) for the sphere, and may be given the symbol H, corresponding to equation (4), or 

H=dxlad"A'. (205) 

Unlike the spherical H, which is a constant for a given satellite 

orbit, the ellipsoidal H varies slightly with position along the groundtrack. This gives an ellipsoidal SOM ground track corresponding to the 

spherical groundtrack plotting of equations (39) and (43). Like them, 

equations (197) and (204) give a conformal groundtrack which is true 

to scale, but this is again not sufficient to reduce the scale error to the 

proper range away from the ground track. The slope must be changed, 

rotating scanlines at the same time, as for the sphere, but the optimum 

correction is slightly more complicated than before. Although the calculation of scale errors requires formulas for cJ>' and "A' in terms of cJ> 

and "A, which are derived subsequently, the rotation formulas are included at this point for convenience. 

Combining equations (201) and (205), equation (46) is again obtained: 

dyldx=SIH (46) 

After changing the slope with the same formula used for the sphere, 

and calculating scale factors at constant cJ>' and "A' along the sides of 

quadrilaterals 0.01 o on a side, the slope was further revised so that the 

new angle ()' of the slope, temporarily using x' and y' again for revised 

coordinates, is found as follows: 

()'=arctan (dy' /dx') 

=8+81 

where ()=arctan (SIH) or the angle of the original slope 

81 =arctan (S/ J) (206) 

and J =(1-e 2

)

3 (207) 

J is an empirical factor giving almost correct scale factors for the 

ellipsoid, but reducing as desired to J = 1 for the sphere, since e is then 

zero. 

Since 

then 

tan ()'=(tan 8+tan 81)/(1-tan ()tan 81) 

dy' /dx'=(SIH +SI J)/(1-82/HJ) 

=(JS + HS )I(HJ -82) (208) 

Further derivation is pursued exactly as in the case of equations (47) 

through (52), but using equation (208) instead of (47). The following 

equivalent to (48) is obtained: 

dx'=[(HJ -S2)/H(J2+S2

)

112] dx (209) 

Substituting from (205), 
46 SPACE OBLIQUE MERCATOR PROJECTION 

(210) 

Integrating, 

x'=a J.[(HJ -S')I(J'+S')11"J d>.' (211) 

Substituting from (210) into (208), and reducing, 

dy'=a[S(H +J)I(J2+S2

)

112] d'A' (212) 

or y'=a J.[s(H +J)I(J'+S2

)

11'] d>.' (213) 

Equations (211) and (213), if primes are dropped from x' and y', apply 

only to the satellite groundtrack. The factors for the scanlines are 

derived below. 

SCANLINE PROJECTIONS 

The normal (or perpendicular) radius of curvature Pn at a given 

point and in a given direction on the ellipsoid, using the nomenclature 

of equations (169) and (170), is as follows (Thomas, 1952, p. 60, eqn. 

157): 

Pn = Pm Ppl (Pm sin2 'Y + Pp cos2 'Y) 

where 'Y is the azimuth east or west from north. 

(214) 

The azimuth of the scanlines is perpendicular to the satellite 

groundtrack in space (but not perpendicular on the rotating globe or 

projected map). The space-fixed azimuth of the ground track north 

from east is i' (equations (182) ff). The azimuth of the scanline west 

from north (space-fixed or Earth-fixed) is also i'. Therefore, 

Pn =pm Ppl(pm sin2 i'+pp COS2 i') (215) 

The scale in the scanline direction is based on this radius. The differential of the map distance S, along the scanline, as a function of the 

angular distance cf>' from the satellite groundtrack, positive (as for the 

sphere) in the direction,of increasing i, is 

dS,=pn dcf>'/cos c/>' (216) 

for a conformal cylindrical projection. Actually, Pn should vary with 

angle c/>', but for simplicity it is taken as the radius at the satellite 

groundtrack, in the proper direction, however. This error practically 

cancels out (less than one-millionth) in final equations. 

Substituting into (216) from (215), 

dS =-1- Pm Pv de/>' 

' cos c/>' Pm (1-cos2 i')+pp cos2 i' 

=-1- Pm Pp de/>' 

COS c/>' Pm -(pm -pp) cos2 i' 
PROJECTION FOR THE ELLIPSOID 47 

Substituting for Pm, Pp, and cos i' from (169), (170), and (202), and 

combining, 

d~ a 1 

dct>'= cos ct>' (

1- 2 . 2 .A.. ) 112 [

1 + e

2 

cos

2 i (1-e2 sin2 c/>0 )

2 J 

e Sin 'Po (1-e 2) (1-e 4 sin2 cl>o cos2 i) 

Since all these values are constants for a given scanline, except for ct>', 

a In tan (hr + ~ct>') 

Ss [ 2 2 · (1 2 • 2 )2 ] (1- 2 · 2 cJ> )112 1 + e cos t - e Sin cl>o 

e Sin o (1-e2) (1-e4 sin2 cl>o cos2 i) 

(217) 

Another symbol is used to simplify this and future expressions. Let 

F=(1-e2 sin2 cJ> }112 1 + e cos t -~ sin cl>o . [ 

2 2 . (1 2 . )2 ] 

o (1-e2) (1-e4 sin2 cl>o cos2t) 

By substituting for sin cl>o from equation (193), and using several steps 

to simplify, it is found that 

_ y.1 +Q sin2 "A' [ U(1 +Q sin2 "A')2 ] 

F -1 + T sin2 "A' 1 + (1 + W sin2 "A') (1 + T sin2 "A') (218) 

where Q, T, and W are as in equations (198) - (200) and 

U =e 2 cos2 i/(1-e2) (219) 

Thus equation (217) becomes 

Ss=(a/F) In tan (i1r+~ct>') (220) 

Geometrically, (a/F) is Pn or the radius of curvature of the ellipsoid 

at the satellite groundtrack in the direction of the scan line. 

The scan lines must be rotated 81 counterclockwise from vertical, 

where 

81 =arctan (S/ J). (206) 

As in the discussion following equation (52), sin 8t times 88 must be 

subtracted from equation (211), and cos 81 times Ss must be added to 

equation (213). 

sin 81 =tan 8/ (1 + tan2 81)

112 

=(S/ J)/(1 +S2fJ2)112 

= S/ (J2 + s 2)112 

cos 81 =1/(1 +tan2 81)

112 

=J /(J2+82)1;2 

Dropping the pri.mes on x' andy' for the final rectangular coordinates, 

f.~: HJ-82 S 

x/a = 0 (J2+S2)1;2 d"A' F (J2+S2)172ln tan (l1r+lct>') (221) 

(222) 
48 SPACE OBLIQUE MERCATOR PROJECTION 

TRANSFORMATION OF GEODETIC COORDINATES 

In order to determine the geodetic latitude and longitude at a given 

cp' and 'A', figure 8 may be used. This figure shows an octant of the 

equatorial aspect of an orthographic view of the ellipsoid, with cp and 

'At as in figure 7, but with the groundtrack, a transformed latitude cp' 

and transformed longitude ).'. A given scanline is a portion of the 

transformed meridian for 'A'. 

I 

I 

z,l I 

I 

I 

+ 

FIGURE 8.-0ctant of ellipsoid with orbit. 

The scanline is perpendicular to the groundtrack in the plane tangent to the ellipsoid, but in figure 8 it is only exactly perpendicular at 

the Equator. At the closest polar approach it is slightly skewed in 

order to remain tangent to the surface of the ellipsoid. At intermediate 

latitudes, its skewness will be less. For the point (cp, 'At), equations for 

x' and z' have already been given: 

x' =a cos cp sin 'A/(1-e 2 sin2 cp)112 

z' =a (1-e 2

) sin cp/(1-e 2 sin2 cp)112 

(176) 

(177) 

The equation for y', in a direction perpendicular to the plane of figure 

8, and increasing toward the reader, is found by multiplying x in 

equation (165) by cos 'At, or 

y'=a cos cp cos 'Atl(1-e 2 sin2 c/>)112 (223) 

Since y' is the same length when viewed in the plane of the orbit (fig. 

9), the distance f from the ascending node to the projection of (cp, 'At) 

in figure 8 onto the groundtrack is 

f=y' tan ).' (224) 
PROJECTION FOR THE ELLIPSOID 49 

FIGURE 9.-Quadrant of plane of satellite orbit for ellipsoid. 

The distance g as projected onto figure 8, from the discussion for 

equations (216) through (220) is nearly, but not exactly, 

g =(a/F) sin cp' (225) 

with cp' positive to the left of the groundtrack in this view. Subsequent 

scale analysis shows that this approximation gives scale factors along 

the scanlines that are in error by no more than one part in a million 

at the ground track. Greater errors away from the ground track are due 

to other approximations as well. 

By simple trigonometric relationships in figure 8, remembering that x' 

is shown as a negative value, 

x'=f cos i -g sin i 

z'=f sin i +g cos i 

(226) 

(227) 

Solving for f and g by multiplying (226) by sin i, (227) by cos i, and 

subtracting for g, then multiplying (226) by cos i, (227) by sin i, and 

adding for f, and substituting from (224) and (225) for f and g, 

(a/F) sin cp'=z' cos i-x' sin i (228) 

y' tan X'=x' cos i+z' sin i (229) 

Substituting for x' and z' from (176) and (177) into (228), and solving 

for cp': 

sin cp' F [(1 2) • . • . • '" ] (1 2 • 2 A..)ll2 -e s1n cJ> cos -r, -cos cJ> s1n -r, s1n "t 

-e s1n '~' 

(230) 
50 SPACE OBLIQUE MERCATOR PROJECTION 

Substituting into (229), and solving for )/: 

tan A'=tan At cos i +(1-e2) tan q, sin i/cos At (235) 

with equation (7) as before to relate At to A and A'. For iteration, follow 

type (2) in the section "Iteration procedures." 

For formulas for A' along the groundtrack in terms of t1>o or ~, with 

tl>' of course equal to zero, we may easily invert (195) and (196), 

sin A'= sin tl>ulsin i (236) 

where tan tl>u=(1-e2) tan tl>o (237) 

to get A' in terms of t/>0 • For A' in terms of ~, we may equate (176) with 

(173a), substituting (172a) for h': 

a cos t/> sin At • , • 

(

1 2 · 2 )1;2 Po Sin A cos ~ (238) -e s1n t/> 

Since in figure 9, when measured to B at the groundtrack, 

y'=AB cos A' 

=p0 COS A1 (238a) 

Equating (223) with (238a), 

a cos t/> cos At , (1 2 . 2 4>}1;2 Po cos A (239) -e s1n 

Dividing (238) by (239), and transposing, we get a simple rigorous 

equation which is identical with (135), the equivalent formula for the 

sphere: 

where 

as before. 

tan A'=tan At/cos i, 

At 0 =Ao+(P2/P1) A' 

INVERSE EQUATIONS 

(135) 

(136) 

Equations (230) and (235) may be inverted as follows to find t/> and 

A from 4>' and A': Equation (235) is transposed and rearranged to solve 

for tan q,: 

(1-e2

) tan tJ> sin i/cos At =tan A' -tan At cos i 

tan tJ>=cos At (tan A'-cos i tan At)/(1-e2

) sin i 

=(tan A' cos At -cos i sin At)l(1-e2) sin i 

(240) 

This is a satisfactory form in which to leave this equation (except for 

the case cos A'=O, discussed later) but it requires the calculation of At 

from t/>' and A' first. To invert equation (230) for this purpose requires 

many more steps. Transposing (230), we isolate sin At: 
PROJECTION FOR THE ELLIPSOID 

cos cf> sin i sin At =(1-e2) sin cf> cos i -sin cf>' (1-e2 sin2 cf>)~IF 

sin At=(1-e 2) tan ct>ltan i-(sin ct>'IF sin i) (1-e2 sin2 cf>)~lcos cf> 

=(1-e2) tan ct>ltan i-(sin cf>'IF sin i) (sec2 cf>-e 2 tan2 cf>)~ 

=(1-e2) tan ct>ltan i -(sin cf>' IF sin i) [1 +(1-e 2) tan2 cf>]~ 

Substituting from (240) for tan c/>, 

tan A' cos At -cos i sin At sin At sin i tan i 

sin cf>' [

1 + (tan A' cos At -cos i sin At )

2 ]~ 

F sin i (1-e 2) sin2 i 

Multiplying through by sin2 i, 

sin2 i sin At =tan A' cos At cos i - cos2 i sin At 

sin cf>' [ • 2 . +(tan A' cos At -cos i sin At) ]~ -----p- s1n ~ (

1-e2) 

51 

Combining the left and the second right terms, then dividing by cos At, 

t '\ -t v . sin cf>' [ sin2 i +(tan A'-cos i tan At) ]~ an 1\t- an 1\ cos t--F ~ (

1 2) cos 1\t -e 

Transposing and squaring, 

sin2 cf>' tan2 At -2 tan At tan A' cos i +tan2 A' cos2 i =~ 

[ . 2 . (1 + t 2 '\ ) + (tan A'-cos i tan At )

2

] stn ~ an 1\t (

1-e2) 

Expanding and collecting, terms 

2 [ sin2 cf>' ( • 2 . cos

2 i )J ( sin2 cf>' ) tan At 1-F2 stn t+ 1_e 2 -tanAt 1 F2(1-e2) 

sin2 cf>' ( tan2 A') (2 tan A' cos i)-~ sin2 i + 1_e2 +tan2 A' cos2 i =0 

Solving for tan At as a quadratic equation, 

tan At=[-B ± (B2 -4AC )~]12A, (241) 

where 

( sin2 f/J' ) , · B=- 1 F 2 (

1-e2) (2 tan A cost) 

A _ 1 stn '+' • 2 • 2 A..' ( . + cos 2 t ') 

- -~ s1n t 1_e2 

C=--- sin2i+-- +tan2A'cos2i sin2 cf>' ( tan2 A') 

F 2 1-e 2 

To simplify these terms before combining them, first let 

K =sin c/>' I F. (242) 

Then the term - B becomes, 

- B =2[1-K 21(1-e2)] tan A' cos i (243) 
52 SPACE OBLIQUE MERCATOR PROJECTION 

while 

2A =2 [1-K2 ( 1-cos2 i+ ~~:~)] 

= 2 [ 1_ K 2 ( 1 + cos2 i -c~s~:te 2 cos2 i)] 

= 2 [ 1-K 2 ( 1 + e ~ ~:: i)] 

=2[1-K2 (1+U)] (244) 

using the term U from equation (219). 

The radical (B2-4AC )lis more complicated: 

B2 -4AC=4[1-K 2/(1-e 2)f tan2 X' cos2 i 

-4[1-K 2 (1 + U)] [-K 2 

( sin2 i + t:~ e~')+tan2 X' cos2 i] 

!(B2 -4AC )=tan2 X' cos2 i -2K2 tan2 X' cos2 i/(1-e 2

) 

+ K 4 tan2 X' cos2 i/(1-e 2

)

2+ K 2 sin2 i 

+ K 2 tan2 X' /(1-e 2

)-tan2 X' cos2 i 

- K 4 (1 + U) [sin2 i +tan2 X' /(1-e 2)] 

+ K 2 (1 + U) tan2 X' cos2 i 

Cancelling the first and sixth terms on the right, multiplying by 

(1-e 2)/ K 2

, and subsequently resubstituting its trigonometric form 

for U, 

(B2 -4AC) (1-e 2)/4K2= -2tan2 X' cos2 i + K 2 tan2 X' cos2 i/(1-e 2

) 

+sin2 i (1-e 2)+tan2 X'-K 2 (1 + U) 

[sin2 i (1-e 2)+tan2 X']+(1 + U) tan2 X' 

cos2 i (1-e 2

) 

= tan2 X' [-2 cos2 i + 1 + cos2 i (1-e 2

) 

+e 2 cos4 i]+sin2 i (1-e 2)+ K 2 

[tan2 >: cos2 i/(1-e 2)-sin2 i (1-e 2)-tan2 X' 

-e 2 sin2 i cos2 i -e 2 cos2 i tan2 X' /(1-e 2)] 

=tan2 X' [ -2(1-sin2 i)+ 1 +(1-sin2 i) 

(1-e 2)+e 2 (1-sin2 i)2]+sin2 i (1-e 2

) 

+ K 2 { tan2 X' [cos2 i/(1-e 2)-1-e 2 cos2 i/ 

(1-e 2)]-sin2 i (1-e 2)-e 2 sin2 i (1-sin2 iH 

=tan2 X' sin2 i (1-e 2+e 2 sin2 i)+ K 2 

[tan2 X' (cos2 i -1)-sin2 i (1-e 2 sin2 i)] 

+sin2 i (1-e 2

) 

=tan2 X' sin2 i (1-e 2 cos2 i)-K 2 sin2 i 

(tan2 X'+ 1-e2 sin2 i)+sin2 i (1-e 2

) 

Multiplying by cos2 X' /sin2 i, 
PROJECTION FOR THE ELLIPSOID 

(B2 -4AC) (1-e 2

) cos2 A.'/4K2 sin2 i =sin2 A.' (1-e 2 cos2 i)-K 2 

(sin2 A.' +cos2 A.' -e 2 sin2 i 

cos2 A.') +cos2 A.' (1-e 2

) 

=sin2 A.' -e 2 sin2 A.' (1-sin2 i) 

+(1-sin2 A.') (1-e 2)-K 2 

[1-e 2 sin2 i (1-sin2 A.')] 

=1-e 2+e 2 sin2 i sin2 A.'-K 2 

53 

(1-e 2 sin2 i +e 2 sin2 i sin2 A.') 

=(1-e 2)+Q sin2 A.'(1-e 2)-K 2 

[1-e 2+e 2 cos2 i +Q sin2 A.' 

(1-e 2)], 

making use of Q from equation (199). We may again divide out (1-e 2): 

(B 2 -4AC) cos2 A.'/4K2 sin2 i =(1 +Q sin2 A.') (1-K 2)-K 2U. 

The radical of equation (241) is then 

(B 2 -4AC )~=(2K sin i/cos A.') [(1 +Q sin2 A.') (1-K 2)- U]~ (245) 

Combining equations (241) through (245) and dividing out the 2 (calculating experience shows that the ± needs to be a - ): 

[1-K 2/(1-e 2)] tan A.' cos i -(K sin i/cos A.') 

[(1 +Q sin2A.') (1-K 2)-K2U]i tan At (246) 

This is the final formula for At in terms of A.' and cp', where K is found 

from (242). 

Then (93) 

If cos A.'= 0 in equation (240) the equation is indeterminate, since cos 

At is also zero. From equation (235), if A.t is 90°, A.' is -90°, and vice 

versa. This applies to each previous or subsequent cycle of 360° in A.'. 

Thus, in (230), if cos A.'=O, sin A.t = ± 1, taking minus the sign of sin A.', 

or, for cos A.'=O, 

sin cp'=[F /(1-e 2 sin2 cp)l] [(1-e 2

) sin cJ> cos i ±cos cJ> sin i] 

taking the sign of sin A.'. To find cJ> in terms of cp', i, and F, first let (sin 

cp' IF)= K. Transposing and squaring, 

K 2 - K 2e2 sin2 cp=(1-e 2

)

2 sin2 cJ> cos2 i ±2 sin cJ> cos cJ> sin i 

cos i (1-e 2

)+cos2 cJ> sin2 i 

=(1-e 2

)

2 sin2 cJ> cos2 i +sin2 i -sin2 cJ> sin2 i 

±2 sin cJ> (1-sin2 cp)i sin i cos i (1-e 2

) 

sin2 cJ> [(1-e 2

)

2 cos2 i -sin2 i + K 2

e2]+[sin2 i-K2]= 

Let 

and 

±2 sin cJ> (1-sin2 cp)i sin i cos i (1-e 2

) 

A =(1-e 2

)

2 cos2 i -sin2 i + K 2

e 2 

B=sin2 i-K2 

(246a) 
54 SPACE OBLIQUE MERCATOR PROJECTION 

Substituting and squaring, 

A 2 sin4 cf>+2AB sin2 cf>+ B 2=4 sin2 cf> (1-sin2 cf>) sin2 i cos2 i (l-e2

)

2 

[A 2+4 (l-e 2

)

2 sin2 i cos2 i] sin4 cf> 

+[2AB -4 (l-e 2

)

2 sin2 i cos2 i] sin2 cf>+ B 2

=0 

Solving for sin2 cf> with the binomial theorem, 

sin2 cf>= {-2AB+4 (l-e 2

)

2 sin2 i cos2 i ±[4A 2

B 2 -16AB 

(l-e 2

)

2 sin2 i cos2 i + 16(1-e 2

)

4 sin4 i cos4 i -4A 2

B 2 

-4B2 (1-e 2

)

2 sin2 i cos2 i]I}/[2A 2+8 (l-e 2

)

2 sin2 i cos2 i] 

= {-AB +2 (l-e2)

2 sin2 i cos2 i ±2 (l-e 2

) sin i cos i 

[(l-e 2

)

2 sin2 i cos2 i-B (A+ B )]l} I 

[A 2+4 (l-e 2

)

2 sin2 i cos2 i] (246b) 

Resubstituting for A and B in the radical, 

(l-e 2

)

2 sin2 i cos2 i-B (A+ B )=(l-e 2

)

2 sin2 i cos2 i -(l-e 2

)

2 

cos2 i sin2 i +(l-e 2

)

2 K 2 cos2 i 

- K 4 (l-e 2

)+ K 2 (l-e 2

) sin2 i 

=K2 (1-e 2) [(l-e 2)cos2 i +sin2 i-K2] 

=K2 (l-e2) (l-e2 cos2 i-K 2

) 

Letting A'=(l-e 2

) sin i cos i (246c) 

and returning to equation (246b ), the formula for cf> in terms of c/J' when 

cos >..' = 0 is as follows: 

sin c/J= ± {[-A (sin2 i-K2

) 

+2A' [A ±K (1-e 2)1 (1-e 2 cos2 i-K 2)1]]!(A 2+4A '

2

) }I (246d) 

Experience with the formula indicates that each of the two "±" signs 

should take the sign of sin >..". 

In the unlikely event that the inclination i of the orbit is 0°-an 

equatorial orbit-equation (240) will be indeterminate. In that case, 

equation (280) may be rewritten fori =0 as follows: 

sin c/J' (l 2 ~ 2 c/>)112 (l-e 2

) sin c/J (247) -e stn 

By considering the geometry, it becomes evident that cf> and c/J' are 

independent of>.. in this case. Solving (247) for c/J, 

sin cf> sin c/J' 

(l-e 2 sin2 c/>)112 F (l-e 2

) 

Substituting K from equation (242), 

sin c/J K 

(l-e 2 sin2 c/>)172 1-e2 

(1-e 2

)

2 sin2 cf>=K2 (l-e 2 sin2 cf>) 

sin2 cf> [(l-e 2)

2+e 2

K2]=K2 

sin cf>=K/[(l-e 2

)

2+e 2

K2]1 (248) 
PROJECTION FOR THE ELLIPSOID 55 

only if i =0. Under these conditions, equation (246) may still be used, 

but it simplifies merely to ~t = ~~, after substituting equation (219) for 

u. 

In a previous paper (Snyder, 1978b) inverse formulas different from 

(240) and (246) are shown as equations (34) through (36) and (38). 

These are partially empirical equations, not exact inverses of present 

equations (230) and (235) although they are satisfactory within the 

Landsat scanning range. Equations (240) and (246), derived after the 

earlier publication, are exact inverses (even though the projection 

itself involves some empirical assumptions) and are slightly simpler to 

compute. Therefore, they replace the earlier formulas, for which the 

derivation is omitted here. 

Inverse equations for </>1 and~~ in terms of x andy are derived from 

ellipsoidal formulas (221) and (222) in a manner completely analogous 

to and as simply as the spherical equivalents, although terms F and J 

are in one set and not the other. Comparable to steps for equations (80) 

and (87), equation (222) is multiplied by (S/ J) and added to (221), 

cancelling out the right term of each: 

x+(S/J)y [' HJ-82 I S [' S(H+J) 1 

a 

0 

(J2+s2)u2 d~ +J 

0 

(J2+S2)112 d~ (249) 

This allows solution for ~~ in terms of x and y, but like equation (80), 

it is impractical as it stands, because it requires both numerical integration and iteration unless Fourier series are used, in which case the 

iteration converges rapidly, as for (80). These series are developed 

below. For </>1 in terms of x andy, equation (222) is transposed: 

1 t (

1 +lA..I) F (J2+S2)~ [U -I.~' S (H +J) d'~] n an 411" 2"' J a (J2+S2)112 1\ 

0 

(250) 

Like equation (87), this first requires calculation of~~ from (249), or its 

Fourier equivalent, but another Fourier series is recommended to 

eliminate numerical integration in (250), even though iteration is not 

required. 

CORRECTION FOR VERTICAL SCANNING 

The above equations were derived for Landsat scanning in which the 

sensor would follow a groundtrack directly between the satellite and 

the center of the ellipsoid, or geocentric scanning. Actually, the sensor 

follows a groundtrack which is vertical, or normal to the plane at the 

surface of the ellipsoid. This is a small but significant correction. (In 

the case of the sphere, geocentric and vertical ground tracks are identical.) As with other groundtrack calculations, the positions for the 

vertical groundtrack may be rigorously determined. 

In figure 10, the satellite is at S. (SA is not a projection of the orbital 

plane, but merely a vector in the orbital plane at some position on the 
56 SPACE OBLIQUE MERCATOR PROJECTION 

NORTH 

POLE 

PLANE TANGENT :----.8 

AT H 

FIGURE 10.-Correction for vertical scanning on an ellipsoid. 

orbit crossing meridian CB.) The radius of the circular orbit is a constant Ro =SA. The satellite scans vertically with true ground track at 

G, at which geodetic latitude is (/>0 • The geocentric ground track is at J, 

at geodetic latitude ¢, and geocentric latitude Q>0 • Geometrically, 

SD=SF+GE=Ro sin ¢u (251) 

FH=FG+EA=R0 cos¢0 

SF /FG =tan 4>o 

(252) 

(253) 

From equations (165) and (166), for rectangular coordinates of point G, 

EA=x =a C08f¢c/(1-e 2 sin2 (/>0 )1 (254) 

GE=y =a (1-e 2

) sin (/>0 /(1-e 2 sin2 4>o)i (255) 

Combining (251) through (255), 

tan ¢o=[Ro sin 4>u-a (1-e 2

) sin ¢j(1-e 2 sin2 (/>0 )i]l 

[Ro cos 4>u -a cos (/>0 /(1-e 2 sin2 (/>0 )i] 

Transposing and multiplying by cos (/>0 , 

Rosin 4>o cos ¢u -a sin 4>o cos ¢o/(1-e 2 sin2 4>o)i=Ro cos 4>o sin ¢u 

-a sin (/>0 cos (/>0 /(1-e 2 sin2 (/>0 )i+ae 2 sin (/>0 cos ¢j(1-e 2 sin2 4>o)i 
PROJECTION FOR THE ELLIPSOID 57 

Combining terms, 

sin ( 4>o- 4>u) = ae 2 sin 4>o cos 4>ol Ro ( 1 - e 2 sin 2 4>o)~ 

. [ ae

2 sin 24>o ] or 4>u = 4>o- arcsin 2R (

1- 2 . 2 ,.. )112 o e Sin '+'o 

(256) 

Since we have previously obtained (236), 

sin X'=sin 4>ulsin i, (236) 

we can now calculate X' for a given 4>o along a vertical groundtrack, 

using (236) and (256). For the inverse, 4>o in terms of X', it was found 

that the direct inverse of (236) and (256), which requires iteration, 

converges rapidly enough to be more satisfactory than an approximation. Thus: 

_ · · , . . . [ ae 

2 sin 24>o ] 4>o- arcsin (sin X sin t) +arcsin 2R (

1 _ 2 . 2 4> )112 o e sin o (257) 

For X' in terms of Xo and vice versa, there is no change from the 

geocentric formulas to correct for vertical scanning, since the groundtrack merely changes latitude along the same meridian when it drops 

from a geocentric to vertical position under the satellite. Thus (135) 

and (136) apply, with iteration required in the direction shown, and no 

iteration for the inverses: 

where 

Inversely, 

tan X'=tan At/cos i 

At 0 =Xo+(P2/P1) X' 

(135) 

(136) 

(93) 

tan At 0 =tan X' cos i (133) 

The best means of correcting the foregoing formulas for the geocentric groundtrack with respect to the vertical groundtrack was a matter 

of some concern. It was apparent from scale factor calculations that 

the earlier equations are quite accurate (though not perfectly, of 

course), as far as plotting meridians and parallels onto the grid of the 

map projection, and they are based on a groundtrack which cuts the 

ellipsoid with a plane, while the "plane" is "warped" for the vertical 

groundtrack. After various approaches were tried, it was decided to 

retain the geocentric formulas exactly as derived, but to change 4>' and 

X' in them to 4>" and X", defined as pseudotransformed latitude and 

longitude, respectively, relative to a geocentric groundtrack. Transformed coordinates 4>' and X' relative to the vertical groundtrack are 

then converted to the equivalent 4>" and X" as follows. 

Equations (257), (93), and (133) are used to calculate 4>o and Xo for a 

given X'. Then equations (235) and (7) are used to calculate A;; (instead 
58 SPACE OBLIQUE MERCATOR PROJECTION 

of ~') for the same f/Jo and ~' where ~ is the adjusted position of ~~ in 

a direction parallel to the groundtrack: 

tan ~=tan ~to cos i +(1-e2

) tan f/Jo sin i/cos ~to (258) 

~t 0 =~o+(P2/P1) ~" (259) 

Equation (230) is similarly used to find f/J~ (instead of f/J') for this f/Jo 

and ~' where f/J~ is the angular displacement of the true groundtrack 

from the geocentric groundtrack: 

sin f/J~ (

1 2 ~ 2 f/J )112 [(1-e2

) sin f/Jo cos i 

-e sin o 

-cos f/Jo sin i sin ~tJ (260) 

Thus f/J" and~" for the groundtrack can be determined as functions of 

~'. Then f/J" for any given angular displacement f/J' in the satellite 

scanning is 

f/J"=f/J' +f/J~ (261) 

and ~" for a given time-angle ~' in the satellite orbit is 

~"=~'+(~-~'). (262) 

These corrections from f/J' to f/J" and ~' to ~" are small and cyclical. 

Therefore, it is better to eliminate the repetitive lengthy calculations 

for each value by using Fourier series as derived below. 

SUMMARY OF CLOSED EQUATIONS 

At this point it is well to bring together the various closed equations 

derived above for the ellipsoidal SOM before giving the Fourier equivalents, and at the same time changing f/J' and ~' to f/J" and ~"so that they 

can be used with the vertical groundtrack corrections. 

1. For x andy in.terms of f/J" and~" (from equations (221) and (222) ): 

a X -r· HJ -82 dV s 1 t (1 +l..L") (263) - 0 (J2 + s2)112 1\ F (J2+ s2)112 n an 411" 2fP 

y -I.~· S (H +J) d " J 1 (1 +1 ") (264) a- 0 

(J2+S2)112 ~ + F(J2+S2)112 n tan 411" 2<P 

where, from equations (201) and (197), 

S=(P2/P1) sin i cos~" [(1 + T sin2 ~")/(1 + W sin2 ~") (1 +Q sin2 ~")]l 

(265) 

from equations (205) and (204) 

-V1+Q sin2 ~" [ 1+ Wsin2 ~" ] H (P2/P1) cos i - 1 + W sin2 ~" (1 +Q sin2 ~")2 (266) 

from equation (218), 

r:-~-:---:..--

- y'1 +Q sin2 ~" [ U(1 +Q sin2 ~")2 ] 

F -1 + T sin2 ~" 1 + (1 + W sin2 ~") (1 + T sin2 ~") (267) 
PROJECTION FOR THE ELLIPSOID 59 

and from equations (207), (198) through (200), and (219), 

J =(1-e 2)

3 (207) 

Q =e 2 sin2 i/(1-e 2

) (199) 

T=e 2 sin2 i (2-e 2)/(1-e 2)

2 (198) 

W=[(1-e 2 cos2 i)/(1-e 2)]2-1 (200) 

U =e 2 cos2 i/(1-e 2

) (219) 

2. For c/> and X in terms of c/>" and X", from equations (240), (242), 

(246), (248), and (93), 

X= At-(P21 P1) X" (268) 

where 

(1-R) tan X" cos i 1-e2 

K si~,,i [(1 +Q sin2 X") (1-K 2)-K 2U]l cos (269) 

K=sin cf>"/F (270) 

tan cf>=(tan X" cos At -cos i sin At)/(1-e 2

) sin i (271) 

If cos X"=O, (271) is indeterminate, but 

sin cf>= ±{ [-A (sin2 i-K 2

)+2A' [A' 

±K (1-e 2)i (1-e 2 cos2 i-K 2)i]]I(A 2+4A'2)}l (246d) 

where the two ± signs take the sign of sin X", 

A =(1-e 2

)

2 cos2 i -sin2 i + K 2

e 2 (246a) 

A'=(1-e 2

) sin i cos i (246c) 

If i =0, (271) is indeterminate, but 

sin cf>=K/[(1-e 2)

2+e2

K2]l (248) 

See step 1 following equation (119) for quadrant adjustment. 

3. For c/>" and X" in terms of c/>' and X', using equations (261) and (262), 

cf>"=cf>'+cl>~ (261) 

X"=X'+(~-X') 

with equations (258) through (260), 

tan ~=tanAt0 cos i +(1-e 2

) tan c/>0 sin i/cos At 0 

At 0 = Xo + (P2/ P1) X~ 

(262) 

(258) 

(259) 

sin c/>~ (

1 2 ~ 2 c/> )112 [(1-e 2

) sin <Po cos i -cos <Po sin i sin AtJ (260) -e s1n o 
60 SPACE OBLIQUE MERCATOR PROJECTION 

and equations (257), (93), and (133), 

_ . . , . . . [ ae 2 sin 2</>0 ] </>0 -arcsin (sin X SID i) +arcsin 2Ro (

1-ez sin2 </>o)ll2 

Xo=Xto -(P1iP1) X' 

tan Xto =tan X' cos i 

(257) 

(93) 

(133) 

Equations (259) and (257) require iteration, and the use of Fourier 

equivalents for the combination of these equations is especially recommended. 

4. For </>' and X' in terms of </>" and X". While the inverse of the above 

equations could be listed, the coefficients of the Fourier series are so 

small that they can be used for forward or inverse merely by changing 

signs, without significant difference. (See below under the heading 

"Fourier series.") 

5. For </>" and X" in terms of </> and X, from equations (235), (7), and 

(230), 

tan X"=tan Xt cos i +(1-e2) tan</> sin i/cos Xt 

Xt=X+(P2/P1) X" 

following type 2 under the section "Iteration procedures." 

(272) 

(273) 

sin</>" (l-e2 ~n2 <1>)112 [(1-e2

) sin</> cos i -cos</> sin i sin Xt] (274) 

6. For </>" and X" in terms of x and y from equations (249) and (250), 

x+(SIJ)y r· HJ-82 " S r· S(H+J) " 

a 

0 

(J2+S2)112 dX +J 

0 

(J2+S2)112 dX (275) 

In tan (br+W) F(J•js•) 112 ~-{<~~!i~J.d>:] (276) 

As stated before, the solution of (275) for X" is possible but totally 

impractical unless a Fourier series is substituted. This series is equation (298) below. 

7. For the satellite groundtrack, 

a. For X' in terms of </>0 , using (236) and (256), 

sin X'= sin <f>ul sin i 

where 

. [ ae 

2 sin 2</>0 ] <l>u = <l>o- arcsin 2Ro (1-e 2 sin2 </>o)l/2 

(236) 

(256) 

b. For X' in terms of Xo, using (135) and (136) with iteration, 

tan X'= tan Xw/ cos i (135) 

Xto =Xo+(P1iP1) X' (136) 
PROJECTION FOR THE ELLIPSOID 61 

c. For cl>o in terms of 'A', using (257), 

. . , . . . [ ae 2 sin 2cp0 ] 

cp0 =arcsin (sin 'A sin ~)+arcsin 2Ro (

1-e2 sin2 c/>o)l/2 (257) 

d. For "Ao in terms of 'A', using (93) and (133), 

"Ao="At 0 -(P2/P1) 'A' (93) 

tan At 0 =tan 'A' cos i (133) 

See step 1 following equation (119) for quadrant adjustments. 

Symbols c/>, 'A, P2, Ph and i are explained after equation (7). In addition, 

a= major semi-axis of ellipsoid (see fig. 4). This is 6378206.4 m for 

the Clarke 1866 ellipsoid. 

e 2=square of eccentricity of ellipsoid (see equation (160) ). The 

square is 0.00676866 for the Clarke 1866 ellipsoid. 

Ro=radius of circular satellite orbit. This is 7,294.69 km for Landsat. 

Other symbols such as cp', 'A', cJ>" and 'A" are intermediate quantities 

calculated from formulas above. 

Iteration procedures given for the spherical SOM, after equation 

(137), apply identically to the corresponding formulas for the ellipsoidal SOM. The only additional ellipsoidal formula requiring iteration is (257), but the iteration procedure for it is the same as that 

following equation (86): trying any value of cl>o in the right side, solving 

for c/>0 on the left side, and using that as the next trial cJ>o on the right 

side, until negligible change takes place. 

Figure 11 shows an enlargement of the second quadrant of the 

projection after the ascending node. This graticule is 10°, calculated 

for the ellipsoid, but almost imperceptibly different from the sphere at 

this scale. 

FOURIER SERIES 

The principle of the Fourier series calculations is given in the section 

under the spherical SOM b~ginning with equation (53). For corresponding ellipsoidal SOM formulas, the Fourier equivalents are very 

similar, except that the constants are somewhat different. Furthermore, since ellipsoidal equations are more complicated than the corresponding spherical equations, there is a benefit in converting some of 

the ellipsoidal terms to series, whereas the corresponding spherical 

terms should be left in the closed form. In any case, the repeating cycle 

is 360° or 21r. 

Beginning with the integral term of (263), for x in terms of 'A", the 

Fourier equivalent is like the first part of equation (73), coupled with 

equations (75) and (76): 

Xala=B"A"+A2 sin 2"A"+A4 sin 4"A"+ ... (277) 
SCAN LINES (EXTENDED 15° FROM GROUNDTRACK) 

~~~SCANNING RANGE 

- GROUNDTRACK 

FIGURE 11.-Enlargement of second quadrant of Space Oblique Mercator projection applied to the ellipsoid. Heavy line indicates satellite groundtrack; dashed scan lines extended to 4>' of 15°. 

oN 

C/) 

"t' 

>- 0 

t'l1 

0 

t:x:l 

t"" 

'8 

t'l1 

~ 

t'l1 

~ 

0 

>-

~ 

0 

~ 

"t' 

2 

t'l1 

0 

~ 

0 

z 
PROJECTION FOR THE ELLIPSOID 63 

where 

211fi2HJ-S2 B =; 0 

(J2+S2)112 dX" (278) 

1 127( HJ-S2 An= 1rn 0 

(J2+S2)112 cos nX" dX" (279) 

As with the spherical cases, the B in the Snyder (1978b) form of 

equation (279) may be omitted. 

Since the coefficient of the In tan term in (263) involves calculating 

S and F from the lengthy equations (265) and (267) for each point, 

although numerical integration is not involved, a Fourier series is 

desirable here. Using the same principles, 

F(J2fS2)112 ~cos X"+ba cos 3X"+b5 cos 5X"+ ... , (280) 

where 

b _1127( s '"d'" n -; F (J2+S2)112 cos n" " 0 

(281) 

Combining (277) and (280) with (263), the Fourier equivalent of the 

entire (263) is 

xla=BX"+A2 sin 2X"+A4 sin 4X"+ ... 

-(~cos X"+ba cos 3X" 

+b5 cos 5X"+ ... ) In tan (br+~<t>"), (282) 

where Bn, An, and bn are found from (278), (279), and (281), respectively. 

Using the constants P2/P1=18/251, i=99.092°, and the Clarke 1866 

ellipsoid, 

B = 1.005798138 for X" in radians 

= 0.0175544891 for X" in degrees 

A2= -0.0010979201 

A4 = -0.0000012928 

A6= -0.0000000021 (negligible) 

~ = 0.0721167880 

ba= -0.0000471816 

bs= -0.0000001313 

~= -0.0000000003 (negligible). 

Similarly, for equation (264), 

y/a =C1 sin X"+C3 sin 3X"+ ... 

+Ouo+o2 cos 2X"+g4 cos 4X"+ ... ) 

In tan (11r+l<t>"), (283) 
64 SPACE OBLIQUE MERCATOR PROJECTION 

where 

_ 1 J.2,.. S(H+J) " " C"- 1rn 0 (J2+S2)112 cos nA dA 

1 J.21r J on=; 0 

F(J2+S2)u2 cosnA"dA" 

With the above constants, 

c1 = 0.1434409899 

Ca = 0.0000285091 

Cs=- 0.0000000011 (negligible) 

uo= 2.000384416 

g2=- 0.0029599346 

04 = - 0.0000032363 

g&= -0.0000000052 (negligible). 

(284) 

(285) 

In equations (270), (260), and (274), there is some benefit in replacing 

F with the following: 

where 

Calculating, 

F=ldo+lh. cos 2A"+d4 cos 4A"+ ... 

1[11" dn =;. 0 

F cos nA" dA" 

do= 1.99702154 76 

lh,=0.0016545831 

d4 =0.0000047966 

d.;= 0.0000000142 (generally negligible). 

For equation (261), 

(286) 

(287) 

cl>" = cl>' + jl sin A'+ ja sin 3A' + js sin 5A' +. . . (288) 

where 

j. =~ r. q,: sin n>.' d>.' (289) 

jl = 0.0085556722 for c/>" and c/>' in degrees 

j3 = 0.0008178430 II 

j5 = - 0.0000026304 II 

j7=- 0.0000000094 11 (negligible). 

For equation (262), 

A"= A' +rl'l2 sin 2A' +m4 sin 4A' +'mt; sin 6A' +. . . (290) 
PROJECTION FOR THE ELLIPSOID 65 

where 

m,. =~ r. (X:- X') sin n X' d X' (291) 

'fn2= -0.0238400531 for X" and X' in degrees 

m4 = 0.0001060576 11 

'1ni; = 0.0000001932 II 

ms= 0,0000000007 II (negligible), 

For the inverses of (288) and (290), the signs are reversed: 

cJ>'=cJ>"-jl sin X"-ja sin 3X"-js sin 5X"-... 

X'= X"-'fn2 sin 2X"-m4 sin 4X"-'1rl1i sin 6X"-. . . 

(291a) 

(291b) 

While these are not true inverses, they are within 0.000003° and 

0.000009°, respectively, of the true inverses of (288) and (290). 

For equation (275), three Fourier series may be combined into two. The 

two integrals are equivalent to equations (277) and the first series of 

equation (283), whileS/ J (in two places) may be converted with the 

equations 

S/ J = Dt cos X"+ D3 cos 3X" + D5 cos 5X" +. . . (292) 

where 

D =__!__ [11" S cos nX" dX" n J7r 0 

Substituting the three series into (275), 

x/a+(y/a) (D1 cos X"+D3 cos 3X"+Ds cos 5X"+ ... )= 

BX"+A2 sin 2X"+A4 sin 4X"+As sin 6X"+ ... 

+(D1 cos X"+ D8 cos 3X"+ D5 cos 5X"+ ... ) 

(C1 sin X"+Ca sin 3X"+Cs sin 5X"+ ... ) 

(293) 

Multiplying the last two parenthetical expressions together and applying equation (84), a simplification similar to that of equation (85) 

may be made to give a right expression as follows: 

x/a+(y/a) (D1 cos X"+D3 cos 3X"+D5 cos 5X"+ ... )= 

BX"+ E2 sin 2X"+ E4 sin 4X"+ Es sin 6X"+. . . (294) 

where 

E2=A2+l(C1 D1 +Ca D1-C1 Da+Ca Da*-Ca Ds*) (295) 

E4=A4+~(Ca D1+C1 Da+CsDl-Cl Ds) (296) 

Es=As+l(CsDl+Cl Ds+CaDa) (297) 

(Terms with asterisks are less than 10-10.) 

Transposing equation (294) gives the final equivalent to (275), 

BX"=x/a+(y/a) (D1 cos X"+D3 cos 3X"+D5 cos 5X"+ ... ) 

- E2 sin 2X"-E4 sin 4X"-E6 sin 6X"-. . . (298) 
66 SPACE OBLIQUE MERCATOR PROJECTION 

Starting with a trial "A" of (x/ aB) on the right side, "A" on the left is 

calculated and substituted for the previous trial value on the right 

until change in "A" is minimal. 

With the same basic constants, Dn and En are as follows: 

D1 = 0.0722098655 

D3=0.0000597921 

D6 = 0.00000007 42 

E2= 0.0040777 482 

E4 = 0.0000040248 

Es=0.0000000040 (negligible). 

For equation (276), the first series of (283) can again replace the integral, but the expression ahead of the bracket can be replaced as follows: 

F (J

2

;S2

)l ~Go +G2 cos 2"A"+G4 cos 4>-."+Gs cos 6>.."+. . . (299) 

where 

G _1 ['If F(J2

+S2)i '-"d"" n -; J cos n 1\ 1\ 

0 

(300) 

Substituting into (276), 

In tan (l1r+lct>")=(!Go+G2 cos 2"A"+G4 cos 4>-."+Gs cos 6>.."+ ... ) 

[(y/a)-(C1 sin "A"+C3 sin 3"A"+C6 sin 5>.."+ ... )] 

Multiplying the parenthetical expressions with G's and C's and using 

equation (84) to combine the terms, the final Fourier equivalent of 

(276) may be written 

In tan (br+~¢")=(~Go+G2 cos 2"A"+G4 cos 4"A"+G6 cos 6>.."+ ... ) 

(y/ a)-~ sin "A"-La sin 3>.."-L5 sin 5>.." + .. . 

where Gn is as in (300) and 

~ =!(Cl Go-Cl G2+Ca G2-Ca G4+C6 G4*-Cs Gs*) 

La=i(Ca Go+C/ G2-Cl G4-Ca Gs*+C6 G2*) 

Ls=i(C6 Go+C1 G4-Cl Gs+Ca G2) 

.Lr=l(C7 Go*+C1 Gs+Ca G4+C6 G2*) 

(Terms with asterisks are less than 10-10.) 

For Landsat constants, 

Go= 1.999624414 

G2=0.0029588261 

G4 =0.0000076132 

Gs = 0.0000000212 

(301) 
PROJECTION FOR THE ELLIPSOID 

G8=0.0000000001 (negligible) 

£,. =0.1432018864 

La=0.0002401662 

L5 = 0.0000005882 

Lr=0.0000000016 (negligible). 

67 

Although integrals are shown from 0 to 211", the comments following 

equation (78) merit a reminder. For a circular orbit, all these Fourier 

coefficients which are not zero may be determined by integrating from 

0 to 1r/2 and multiplying by 4, due to quadrant symmetry. 

DISTORTION ANALYSIS 

The ellipsoidal distortion analysis is also analogous to that for the 

spherical SOM, except that the formulas are, again, more complicated. 

Paralleling the discussion beginning with equation (94), but applying 

the principles to the ellipsoid, the distance on the map for a differential change of latitude is still 

(as/acJ>)m =[(axlacJ>)2+(ay/acJ>)2Ji (94) 

and for a differential change of longitude is still 

(as/ aA)m =[(ax/ aA)2+ (ay I aA)2]i (95) 

The corresponding distances on the ellipsoid, though, are Pm and Pp 

cos cJ> as determined from equations (169) and (170): 

(as/acJ>)e=a(1-e2)/(1-e2 sin2 cJ>)I 

(aslaA)e=a cos cJ>/(1-e2 sin2 cJ>)i 

(302) 

(303) 

Dividing (94) by (302) gives h, the scale factor along a meridian of the 

ellipsoid: 

h =[(ax/acJ>)2+(ay/acJ>)2Ji (1-e2 sin2 cJ>)I/a(1-e 2) (304) 

Dividing (95) by (303) gives k, the scale factor along a parallel: 

k =[(ax!aA)2+(aylaA)2]i (1-e2 sin2 cJ>)i!a cos cJ>. (305) 

Combining (304) and (305) for a general scale factor m, using finite 

increments, 

m -[ (~x)2+(~y)2 ]i - [a(1-e2) ~ cJ>f+ (a cos cJ> ~ A)2 

(1-e 2 sin2 c/>)3 1-e2 sin2 cJ> 

(306) 

The discussion following equation (100) applies equally well here, 

through equation (104), with 0' defined the same way, but with h and 

k determined from (304) and (305). It is not useful to try to combine 

equations to obtain an equation for w similar to (107) or (108); such a 

combination is about as complicated as calculating w from (101) 

through (104) and (304) through (305), so the existing equations should 

be used instead. 
68 SPACE OBLIQUE MERCATOR PROJECTION 

To apply these ellipsoidal scale factors and maximum angular deformation to the SOM formulas, a sequence of differentiations is carried 

out. They are straightforward, subject to the inadvertent errors which 

can and did plague so many of the writer's initial derivations of foregoing equations. These have also been calculator tested and are believed to be free of errors. 

In differentiating equation (265), the parenthetical expressions may 

be treated separately: 

d (1 + T sin2 ).")/d)."= 2T sin )."cos )." = T sin 2)." 

For a calculating shorthand, let 

(1 + T sin2 ).") = T' (307) 

Then 

dT'Id>-."=Tsin 2>-." (308) 

Likewise, the shorthand may be used with expressions (1 + W sin2 >-.") 

and (1 +Q sin2 >-."). Then (265) may be rewritten 

S=(P2/P1) sin i cos>-." (T'/W'Q')i (309) 

dS/d>-."= -(P2/P1) sin i sin>-." (T'IW'Q')I+(P~Pt) sin i cos>-." 

[T sin 2>-." /2(T'W'Q')I-l(T' /Q')I W sin 2>-." /( W')l 

-l(T' IW')I Q sin 2>-."/(Q')~ 

=S[sin >-."cos>-." (T/T'- W/W'-Q/Q')-tan >-."] (310) 

Rewriting (267), then differentiating, 

F=(Q'/T')~ [1 +U(Q')2/W'T'] (311) 

dF/d>-."=[1 + U(Q')2/W'T'] [Q sin 2>-."/2(T'Q')I 

- ( Q')l T sin 2>-." /2( T')~ 

+(Q'IT')~[2UQ'Q sin 2>-."IW'T' 

- U ( Q')2 W sin 2>-." I ( W')2 T' 

-U(Q')2 Tsin 2>-."IW'(T')2] 

After some intermediate steps this simplifies to 

dF / d >-."=sin 2>-." t F ( Q/Q'-T /T')/2 

+[F-(Q'IT')l] (2Q/Q'- W/W'-T/T'H (312) 

Differentiating (272) and (273), 

sec2 )."d>-."=cos i sec2 >-.td>-.t +(1-e2) (sin i sec2 q,/cos >-.t) dq, 

+(1-e2) sin i tan q, sec ).t tan >-.td>-.t (313) 

(314) 
PROJECTION FOR THE ELLIPSOID 69 

Combining (313) and (314), and separating differentials, 

sec2 'A"-(P2/P1) [cos i/cos2 At 

+(1-e2) sin i tan cp sin 'A/cos2 'At] d'A." 

=[cos i/cos2 'At +(1-e2) sin i tan cp sin 'At/cos2 'At] d'A. 

+(1 '-e2) (sin i/cos2 cp cos 'At) dcp (315) 

The partial derivatives are then 

a'A."/acp=(1-e2) sin i/cos2 cp [cos 'A/cos2 'A"-(P2/P1) [cos i 

+(1-e2) sin i tan cp sin 'At]/ cos 'At] (316) 

a'A"/a'A.=[cos i +(1-e2) sin i tan cp sin 'At]/[cos2 'Atlcos2 X" 

-(P2/P1) (cos i +(1-e2) sin i tan c/J sin >-t)] (317) 

Differentiating (27 4), 

cos cp"dcp"=[(1-e2) cos i sin cp-sin i cos cp sin >-t] dF /(1-e2 sin2 cp)~ 

+ F { (1-e2 sin2 cp)~ [(1-e2) cos i cos cp dcp+sin i sin cp 

sin 'At dcp-sin i cos cp cos >-t d Xt]-[(1-e2) cos i sin c/J 

-sin i cos c/J sin >-tj (1-e2 sin2 c/J)-~ 

( -2e2 sin c/J cos cp dcp)}/(1-e 2 sin2 cp) 

Separating differentials, and then giving partial derivatives as before, 

~~ =B A(a>."/acp)+(1-e 2

) cos i +sin i tan c/J sin >-t 

+e 2 sin cp [(1-e2) cos i sin cp-sin i cos cp sin 'At] (318) 

~f =B [A(a>."/a>.)-sin i cos >-t] (319) 

where 

A =(1/F) [(1-e2) cos i tan cp-sin i sin 'At] dF/dX" 

-(P2/Pl) sin i cos 'At (320) 

B=F cos cp/cos c/J" (1-e2 sin2 cp)~ (321) 

Differentiating (263), 

dx HJ-8 2 S 

a (J2+S2)112 dX" F (J2+S2)172 sec cp"dcp"-ln tan (br+lc/J") 

{ l[(J2+S2)~-iS(J2+S2)-~2SldS+ S (-F-2)dF} F J2+S2 -J (J2+S2)112 

(J2+

1

S 2)u2 [(HJ -S2

) dX"-~ sec cp"dcp"-(F (J2~S2)1;2 dS 

-:2 dF) ln tan (i1r+lc/J")] (322) 
70 SPACE OBLIQUE MERCATOR PROJECTION 

Taking the partial derivative with respect to c/>, using intermediate 

differentials derived earlier for Sand F, then rearranging, 

~ F(J•:s•)m }[F(HJ-S') 

( J 2 dS SdF) , ifA" - J2+S2dA"-FdA" In tan (hr+l</> )] ac~> 

s ac~>"} cos c~>" ac~> (323) 

The formula for (iJx/aA) is the same as (323), but with A substituted for 

cf> at all three places: 

:~ F(J2~S2)112 { [F(HJ -82

) 

( J 2 dS SdF) , aA" - J2+S2 dA"-FdA" In tan (hr+lc!> )] aA 

s ac~>"} cos c~>" aA (324) 

Obtaining the same sort of relationships for y, equation (264) is differentiated: 

dy 

a 

S(H+J) d'"+ J "d " (J2+ 8 2)112 1\ F (J2+ s2)1/2 sec cf> cf> 

+In tan (i1r+lc1>") · { ~ [(J2~:2)s12] dS 

+ (J2.:S2)112 (-F -

2

) dF} 

(J2+1S2)112[s (H +J) dA"+~ sec q,"dq," 

-(F(Jf+S2) dS+ ;2dF )In tan (hr+lc!>")] (325) 

(326) 

As with x, the formula for (ay/aA) is the same as (326), but with A 

substituted for q, at all three places. 

~ F(J2~S2)112 {[FS(H+J)-(J2!s2~+~:~) 

, ] aA" J aq," } In tan (l1r+lc~> ) aA +cos cf>" aA (327) 
PROJECTION FOR THE ELLIPSOID 71 

The equations for distortion analysis of the ellipsoidal SOM have 

now been completed. There appears to be no advantage in further 

combining the equations, but each portion may be calculated separately as part of a multipart program and then combined. Specifically, 

paralleling the analysis for the spherical SOM, 

1. For a given ¢" and X", calculate ¢, 'At and X from equations (268) 

through (271), placing 'At in the proper quadrant by the procedure 

detailed in step 1 following (119). 

2. Calculate J, S, H, and F from equations (265) through (267), (198) 

through (200), (207), and (219). 

3. Calculate (dS/d'A") from (310) and (dF/d'A") from (312), using the 

notation of (307), etc. 

4. Calculate (aX"/a¢), (aX"/aX), (a¢"/a¢), and (a¢"/ax) from (316) 

through (321). 

5. Calculate (ax/a¢), (ax/aX), (ay/a¢), and (ay/aX) from (323), (324), 

(326), and (327). 

6. Calculate h, k, and w from equations (304), (305), and (101) through 

(104). 

While these calculations are obviously lengthy, they are important to make, but need to be made only once for a given satellite 

orbit (like Fourier constants) to establish the degree of accuracy. They 

do not involve computer time during mapping operations. 

Table 4 lists h, k, and w for each 15° of a quadrant of the ellipsoidal 

SOM for Landsat. The values are first given within 1 o of the groundtrack to show normal ranges of errors, and then given up to 15° from 

the groundtrack (with fewer decimals) to show ranges of errors when 

the SOM is used for other satellites. While other orbital characteristics would affect these values slightly, they have been found to be 

representative of errors to be expected. In the 1 o range, sin iw is also 

given, since it is the variatiol\ which occurs from the normal conformal 

scale factor at that position, due to the lack of perfect conformality of 

the projection. For example, at '/..."=15° and ¢"=1 °, the scale varies up 

to 0.000019 more or less from the normal scale factor of about 1.000152, 

because of the maximum angular deformation of 0.0022°. 

The nonsymmetry of the errors with respect to both X" and ¢" is 

interesting, due largely to the manner in which the radius of curvature 

changes away from the groundtrack. The slight but negligible variation of the ground track from true conformality and scale is due to the 

approximation used for equation (225) and deviates from unity in the 

crosstrack direction only. 
72 SPACE OBLIQUE MERCATOR PROJECTION 

TABLE 4.-Scale factors and angular deviation away from geocentric Landsat graundtrack 

for the eUipsoidal Space Oblique Mercator projection 

At 1° from groundtrack 

'A" 4>" h k w sin lh w 

oo 10 1.000154 1.000151 0.0006° 0.000005 

0 1.000000 1.000000 .0000 .000000 

-1 1.000154 1.000151 .0006 .000005 

15 1 1.000161 1.000151 .0022 .000019 

0 1.000000 1.000000 .0001 .000000 

-1 1.000147 1.000151 .0011 .000010 

30 1 1.000167 1.000150 .0033 .000029 

0 1.000000 1.000000 .0001 .000001 

-1 1.000142 1.000150 .0025 .000021 

45 1 1.000172 1.000150 .0036 .000031 

0 .999999 1.000000 .0001 .000001 

-1 1.000138 1.000150 .0031 .000027 

60 1 1.000174 1.000150 .0031 .000027 

0 .999999 1.000000 .0002 .000001 

-1 1.000136 1.000150 .0028 .000025 

75 1 1.000174 1.000152 .0019 .000016 

0 .999999 1.000000 .0001 .000000 

-1 1.000135 1.000150 .0019 .000016 

goo 1 1. 000170 1.000156 .0008 .000007 

0 .999999 1.000000 .0000 .000000 

-1 1.000133 1.000151 .0010 .000009 

Beyond 1° from ground track 

'A' 4>" h k w 

oo 15° 1.03565 1.03490 .136° 

10 1.01559 1.01526 .060 

5 1.00386 1.00378 .015 

-5 1.00386 1.00378 .015 

-10 1.01559 1.01526 .060 

-15 1.03565 1.03490 .136 

15 15 1.03589 1.03472 .157 

10 1.01570 1.01521 .074 

5 1.00390 1.00377 .022 

-5 1.00383 1.00378 .006 

-10 1.01551 1.01530 .041 

-15 1.03549 1.03503 .107 

30 15 1.03616 1.03453 .162 

10 1.01583 1.01514 .080 

5 1.00395 1.00376 .027 

-5 1.00381 1.00378 .003 

-10 1.01548 1.01530 .024 

-15 1.03543 1.03510 .076 

45 15 1.03636 1.03440 .147 

10 1.01594 1.01509 .075 

5 1.00398 1.00375 .027 

-5 1.00380 1.00377 .007 

-10 1.01549 1.01529 .014 
INCORPORATION OF LANDSAT ROW AND PATH NUMBERS 73 

TABLE 4.-Continued 

Beyond 1° from ground track 

)>.' r/J" h k w 

-15° 1.03546 1.03511 .052° 

15 1.03635 1.03446 .112 

10 1.01597 1.01509 .059 

5 1.00400 1.00374 .022 

-5 1.00381 1.00376 .009 

-10 1.01552 1.01526 .015 

-15 1.03553 1.03510 .041 

75 15 1.03596 1.03489 .061 

10 1.01583 1.01526 .032 

5 1.00398 1.00378 .012 

-5 1.00381 1.00376 .008 

-10 1.01556 1.01525 .018 

-15 1.03554 1.03514 .040 

90 15 1.03553 1.03533 .011 

10 1.01561 1.01549 .007 

5 1.00391 1.00385 .004 

-5 1.00373 1.00384 .006 

-10 1.01524 1.01557 .019 

-15 1.03449 1.03571 .040 

Notes: >."=angular position along geocentric groundtrack, from ascending node. 

r/J"=angular distance away from geocentric groundtrack, positive in direction 

away from north pole. 

w=maximum angular deformation 

sin 1f2 w=maximum variation of scale factors from true conformal values. 

h =scale factor along meridian of longitude 

k =scale factor along parallel of latitude 

INCORPORATION OF LANDSAT ROW AND PATH NUMBERS 

Rather than using the number of orbits N discussed under the second 

type of equation in the section "Iteration procedures," the user of 

Landsat information is more likely to refer to the row and path numbers of a given Landsat orbit. A given Landsat orbit has the same path 

number from one north polar approach to the next north polar approach. Therefore, there are 251 path numbers before the orbits exactly repeat. The paths are numbered not in the order of the occurrence of the orbit, but consecutively for each path encountered as one 

follows a given parallel of latitude westward around the Earth from a 

specified starting point. Thus path 2 is 1/251 of the circumference, or 

360°/251=1.434° of longitude west of path 1. Chronologically, each 

successive orbital journey of a Landsat satellite has a path number 18 

greater than the previous one, subtracting 251 if the path number 

would exceed this limit. (With each satellite revolution, the Earth 

turns 18/251 of a complete rotation, or 18 path numbers.) 
74 SPACE OBLIQUE MERCATOR PROJECTION 

Since the path number changes near the north polar approach, the 

path cycle approximately covers a change in X' of 90° to 450° if X' is 

related to the preceding ascending node. The row numbering is proportional to the time along the orbit and therefore, proportional to X'. 

Since 248 rows are used, each corresponds to 360°/248 or 

45°/31 = 1.452° of X'. Numbering of rows does not begin until a X' of 

94.355°, so that a X' of 180° (the descending node) corresponds to row 

60. It should be noted that the ascending node on a given path occurs 

at row 184 (60+248/2), when X'=360°. The ascending node at a X' of0° 

occurs at a fictitious row number 184-248 or -64. 

Relating the above information to formulas for the satellite groundtrack, or to arithmetic considerations: 

For X' in terms of r, the number of the row, 

X'=90°+(45°/31)(r+2) (328) 

and conversely, 

r = (31/ 45° )X'-64. (329) 

For cf>o or Xo along the ground track in terms of the row r, use equations 

(257), (93), and (133), adjusted for quadrant as described in step 1 

following equation (119). 

To convert from longitude Xa east ( +) or west (-) of Greenwich to the 

longitude X relative to the ascending node, row "-64" of a given Landsat path p, use the following: 

X=Xa-128.87°+(360°/251)p. (330) 

Conversely, 

Xa =X+ 128.87°-(360° /251)p. (331) 

This is the value of X which must be used in the various SOM equations, 

not a value based on the later ascending node. 

SUMMARY OF SIMPLIFIED EQUATIONS 

FOR THE ELUPSOIDAL SOM FOR LANDSAT 

While some simplifications and empirical treatment were used in 

the development of the foregoing equations, several factors make the 

equations appear even more complicated than necessary when compared with equations for other projections, such as the ellipsoidal 

transverse and oblique Mercator projections. By deleting the following 

functions, the equations may be somewhat simplified to a set of equations for converting from latitude and longitude to rectangular coordinates and vice versa, especially for Landsat: 

1. Remove all equations for finding positions along the groundtrack. 
SUMMARY OF ELLIPSOIDAL SOM FOR LANDSAT 75 

2. Remove all Fourier constants which are not to eliminate numerical integration. Show calculation only for integration of one 

quadrant. 

3. Eliminate the factor F, since it almost completely cancels out in 

the region of Landsat sensing, and has a very slight effect in 

broader areas. 

The formulas can then be summarized as follows: 

Constants to be calculated for the given satellite orbit: 

B = (2hr) [

12

[ H J-S2)/ (.I'+ S')l'ylx• 

A. =(4hrn) ['l<H J -S2)/(.f+S2)l]cos nX'dX' 

where n =2 and 4. 

J.'!r/2 1 

Cn=(4hrn) 

0 

[S(H +J)I(J2+S2)2] cos nX"dX" 

where n = 1 and 3. 

(278) 

(279a) 

(284a) 

Simpson's rule (equation (78)) is recommended for integration, using 

go steps for each of the above. For S, H, J, Q, T, and W, equations 

(265), (266), (207), (199), (198), and (200) may be used unchanged. They 

are shown together following (265). For Landsat, i =99.092° and 

P2/P1 =18/251. 

Forward equations (geodetic to rectangular coordinates): 

Using equation (330) to change AG to X, 

At= X+ (P21 P1)X" (273) 

X"=arctan [cos i tan Xt+(1-e2

) sin i tan cJ>/cos At] (272a) 

iterating together, using type 2 under "Iteration procedures." 

cJ>"=arcsin{[(1-e2

) cos i sin cJ>-sin i cos cJ> sin Xt]l 

(1-e2 sin2 cJ>)~} (274a) 

using the final At above. 

x!a=BX"+A2 sin 2X"+A4 sin 4 X"-[S/(J +S )~] 

In tan (n/4+cJ>"/2) (282a) 

y/a =C1 sin X"+C3 sin 3X"+[J /(J2+S )~] In tan (n/4+cJ>"/2) (283a) 

where S is found, as in the integration, from equation (265). 

Inverse equations (rectangular to geodetic coordinates): 

X"=[(xla)+(SIJ)(yla)-A2 sin 2 X"-A4 sin 4X"-(S/J)(C sin X" 

+ C3 sin 3X")]/B (275a) 

with S found from equation (265), and iteration as described following 

equation (298). 
76 SPACE OBLIQUE MERCATOR PROJECTION 

4J"=2 arctan {exp((1 +S2/J2)i(y/a-C1 sin 'A." 

-C3 sin 3'A.")]} -1r/2 (276a) 

'A.t =arctan {[(1-sin2(/J" /(1-e2)) tan 'A." cos i -(sin 4J" sin i/cos'A.") 

((1 +Q sin2 'A.") cos2 4J"- U sin24J")1]![1-sin24J"(1 + U)]} (269a) 

adjusting 'A.t for quadrant as described in step 1 following equation 

(119). 

U =e2 cos2 i/(1-e2

) 

'A.= 'A.t-(P2/ P1)'A." 

Conversion of 'A. to 'A.a is carried out according to equation (331). 

(219) 

(268) 

(/J=arctan [(tan 'A." cos 'A.t -cos i sin 'A.t)/(1-e2

) sin i] (271) 

using the quadrant-adjusted value of 'A.t obtained above. 

MODIFICATION OF ELLIPSOIDAL SOM EQUATIONS 

FOR NONCIRCULAR ORBITS 

For a noncircular orbit, 'A.' is not proportional to time, but if the orbit 

is elliptical, 'A.' may be calculated by Keplerian formulas from the mean 

longitude in orbit of the satellite. The factor (P2/P1)'A.' of equation (7) 

is changed, and this change necessitates changes in the derivations of 

the earlier formulas. The changes are actually not very difficult for 

obtaining a correct groundtrack, but the proper curve to give this 

groundtrack on the plotted map becomes much more complicated, except for small orbital eccentricities of 0.05 or less. For larger eccentricities, see the next section, although several formulas in this section 

apply to an orbit of any eccentricity and are referred to in the next 

section. 

In terms of the discussion following equation (28) and leading to equation (7), the "satellite-apparent" longitude 'At will be displaced from 

the actual longitude 'A. by an angle proportional to time, called the 

mean longitude of the satellite along its orbit. If we call the mean 

longitude L for the satellite with respect to the perigee, and if the 

perigee of the orbit is at a geocentric angle w with respect to the 

ascending node of the orbit, the difference between 'A. and 'At becomes 

(P2/ P1)(L + w ), or 

(332) 

To relate 'A.' to L, standard Keplerian equations are used (Moulton, 

1914): 

tan lh('A.'-w)=(tan lhE')y(1 +e')/(1-e') (333) 

where 

E-e' sin E'=L, (334) 
SOM EQUATIONS FOR NONCIRCULAR ORBITS 77 

e' is the eccentricity of the satellite orbit, 

E' is an intermediate angle, or the "eccentric anomaly", 

and w is the longitude of the perigee relative to the ascending node. 

E' may be found from L by iteration, first rearranging (334) to 

E'=L+e' sinE' (335) 

assuming a value forE', inserting it in the right side, solving for the 

left side, and substituting the new value of E' into the right side, until 

there is sufficient convergence. Convergence is very rapid when e' is 

small. This is the same technique used in much of the iteration in 

preceding formulas. 

The radius vector of the elliptical orbit, R0 , is found from the following 

equation: 

Ro=a'(1-e' cos E') (336) 

where a' is the mean distance of the satellite from the center of the 

Earth, in the same units as a. 

Equations for the satellite ground track as plotted on the globe thus 

vary from those given for the ellipsoid, as follows: 

Equation (236)-no change. 

(256)-no change, except that Ro is found from (336). 

(135)-no change. 

(136) becomes At 0 = Ao + (P2/ P1) (L + w ). (337) 

(257)-no change, except that Ro is found from (336). 

(93) becomes Ao = At 0 - (P2/ P1) (L + w ). (338) 

(133)-no change. 

In deriving the general equations for the ellipsoidal SOM, the relationship between A', L, and At for the ellipsoid first affects the definition of BD following equation (181). Now, 

BD=(P2/P1) a cos cl>o dL/(1-e 2 sin2 c/> )~, 

so that equation (182) becomes 

dy =(P2/P1) a cos cl>o sin i' dL/(1-e 2 sin2 cp )~, 

and (183) becomes 

dx =ds -(P2/P1) a cos cl>o cos i' dL/(1-e 2 sin2 c/> )~. 

The next change is at equation (190), where A' is merely replaced with 

L Equation (197) becomes 

dy/adL =(P2/ P1) sin i cos A'[(1 + T sin2 A')/ 

(1 + W sin2 A') (1 +Q sin2 A')]~. 
78 SPACE OBLIQUE MERCATOR PROJECTION 

It is more convenient to change this to a function of d X': 

dy/adX'=(P21P1) sin i cos X1(1 + T sin2 X')/(1 + W sin2 X') 

(1 +Q sin2 X'):f dL/dX'=S (339) 

Equation (203) becomes 

dx _ ds P2 cos i (1-e2 sin2 c/>0 )

112 

adL - dL -a P1 (1- e 4 sin2 cf>o cos2 i)112 

_ ds dX' P2 cos i (1-e2 sin2 c/>0 )

112 

- dX' dL -a P1 (1-e 4 sin2 cf>o cos2 i)112 

The only change this introduces into (204) is the insertion of dX' /dL: 

dx _ y'1 +Q sin2 X' [ 1 + W sin2 X' dX' P2 . ] 

adL -1 + W sin2 X' (1 + Q sin2 X')2 dL - P1 cos "' 

As in the case of equation (339), it is preferable to multiply by dL/dX': 

dx _ y'1+Qsin2X' [ 1+Wsin2X' P2 dL ·]-

adX'- 1+Wsin2X' (1+Qsin2X')2 P1dX'cost -H <

340) 

The value of dL/ dX' is derived from equations (333) and (334): 

Differentiating (334), 

dL=(1-e'cosE')dE' (341) 

Differentiating (333), 

lsec2 MX' -w)dX'=lsec2lE' y'(1 +e')/(1-e') dE' .~--::~----:---;;- dE' /dX'=sec2l(X' -w) cos2lE' y'(1-e')/(1 +e') (341a) 

Substituting into (341), 

dL/dX'=(1-e' cos E') sec2l(X'-w) cos2lE'y'(1-e')/(1 +e') (342) 

From (333), however, 

tan2 ~(x'=-w)=(tan lE') [(1 +e')/(1-e')] 

sec2l(X' -w)=1 +(tan2lE') [(1 +e')/(1-e')] 

Substituting into (342l 

dL/ dX' =(1-e' cos E') [ cos2lE' y'(1-e')/(1 +e') 

+sin2lE' y'(1 +e')/(1-e')] 

=1(1_ , E') [<1 +cos E') (1-e')+(1-cos E') (1 +e')] 2 e cos (

1 _ e'2)u2 

_t( , ') (2-2e' cos E'] - 2 1-e cos E (

1-e'2}112 

=(1-e' cos E') /(1-e' )~ 

where E' is found from X' in equation (333). 

(343) 
SOM EQUATIONS FOR NONCIRCULAR ORBITS 79 

As with both spherical and ellipsoidal equation derivations, the 

values of x and y calculated from equations (340), (339), and (343) will 

give the groundtrack true to scale, but the groundtrack must be bent 

to a different shape to give sufficiently accurate mapping of the region 

near the groundtrack. For small orbital eccentricities the bending of 

the groundtrack can follow exactly the same approach as that used in 

ellipsoidal equations (206a) through (213). The eccentricity of the orbit 

does not affect this, except for its incorporation into the new formulas 

for H (340) and S (339); J does not change. Therefore, equations (211) 

and (213) remain the same: 

x'=a I,'' [(HJ -S2)/(J2+S')l] dA' (211) 

y'=a f' [S (H +J)/(J'+S')l] dA' (213) 

An additional complication is introduced by the elliptical orbit, however. With x' andy' starting at (0,0) when X' is zero, y' will not be zero 

when X' is 360°, since the quadrants are not symmetrical. For convenience, it is preferred that y return to zero at each 360° repetition. 

This is accomplished by rotating the X' and Y' axes slightly so that y 

becomes zero at each 360° interval. Calling the rotational angle 02, 

x =x' cos 02+y' sin 02 

y =y' cos 02-x' sin 02 

(344) 

(345) 

using standard formulas. The tangent of 02 is the ratio of y' at 360° to 

x' at 360°, or 

where 

B1 =_!_J.

2

11" [(HJ -S )/(J +S )~] dX' 21r 0 

B,= ~r.· [S (H +J)/(J'+S')l] dA' 

(346) 

(347) 

dividing each equation by 21r to provide convenient numbers usable 

later, since 21r cancels out in their use for 02• 

Then 

and 

sin 02=tan 02/(1 +tan2 02)i 

=B2/(B12+ B22)i 
80 SPACE OBLIQUE MERCATOR PROJECTION 

Assigning symbols H1 to cos 82 and & to sin 82, and deleting a from 

(211) and (213) for simplicity, equations (344) and (345) become 

where 

x/a=x'H1+y'S1 (348) 

y/a=y'H1-x'S1 (349) 

H1 =B/(B12

+ B22)

112 

& =B2/(B12+ B22)1t2 

(350) 

(351) 

B1 and B2 are found from (346) and (347), and x' andy' are found from 

(211) and (213) omitting the a. 

In developing positions of points along the scanlines, the equations 

for calculating F (214) through (220) are unchanged by orbital eccentricity. Rotation of scanlines in the x' andy' frame occurs as before, 

but the above rotation of X' and Y' axes to X and Y must include the 

scanlines as well. Therefore, equations (221) and (222) undergo only a 

slight change to the following in place of equations (211) and (213): 

J.x' HJ-82 S 

x'= 0 (J2+S2)112 dX' F(J2+S2)u2ln tan (br+lq,') (352) 

'-J.>.' S (H +J) d'\'+ J 1 t (1-+1..1..') (353) Y - 0 

(J2+S2)112 1\ F (J2+S2)112 n an ;pr 2fP 

Subsequent derivations of equations for transformation of geodetic 

latitude and longitude into c/J' and X' are unchanged from the ellipsoidal 

form except for the substitution of (P2/P1) (L +w) in place of (P2/P1)X' 

wherever the latter appears. This occurs in equations (7) and (136), and 

wherever the like formula is referred to for calculation of Xt, etc. The 

iteration required for solving equations (235) and (7), or (135) and 

(136), remains as before using successive substitution, except that 

additional steps are involved. 

Procedure (3) under the section "Iteration procedures" is followed, 

but the first trial X'(X'p) is used in place of X' in equation (333), which 

is solved for E', for insertion into (334), which is solved for L This 

value of Lis then used in (332), solving for Xt, and Xt is in turn used in 

(235) or (135), from which a new value of X' is found for use as the trial 

X' in equation (333). This is repeated until sufficient convergence occurs. 

For inverse equations for the noncircular SOM, the foregoing discussion coupled with approaches taken for circular inverse equations, 

makes the modifications fairly obvious: 

For c/J and X in terms of Q>' and X', only equation (93) is changed, to 

X=Xt-(P2/Pl) (L+w) (354) 

the inverse of (332). 
SOM EQUATIONS FOR NONCIRCULAR ORBITS 81 

For ~' and X' in terms of x and y, the axis rotation is the principal 

change. The inverse of (344) and (345) involves standard formulas, so 

that 

x'=x cos 82-y sin 82 

y'=y cos 82+x sin 82 

which are converted to our terminology as the inverse of (348) and 

(349): 

x'=(x/a) H1-(y/a) S1 

y'=(y/a) H1 +(x/a) S1 

(355) 

(356) 

Equations (249) and (250) remain the same, except that x' and y' are 

substituted for (x/ a) and (y/ a), respectively: 

I '-J.'A' HJ -82 I s f.'A' s (H +J) I 

X +(S/ J)y - 0 (J2+S2)1/2 dX +J 0 (J2+S2)112 dX (357) 

1 1 I F (J2+S2

)

112 [ I f.'A' s (H +J) '] In tan (11r+2Q>) J y -

0 

(J2+S2)112 dX (358) 

Solution using Fourier series is almost mandatory, as it was for the 

circular orbits. 

To correct for vertical rather than geocentric scanning, the same approach is taken as with circular orbits, but with the revised formulas 

above. As before, this changes single primes to double primes in most 

of the revised formulas, but a slight modification is shown below in the 

formula summary. 

In summary, for equations giving a true-to-scale groundtrack for any 

Keplerian elliptical orbit, and sufficiently accurate plotting of points 

near the groundtrack (see later discussion) for an orbital eccentricity 

of about 0.05 or below, the following formulas may be used: 

1. For x and y in terms of ~" and X" use equations (348), (349), (350), 

(351), and the following modifications of (333), (339), (340), (343), (346), 

(347), (352), and (353): 

B1 =(1/2,.) f.l<HJ -S2)/(.J'+S')!JdA" (359) 

B2=(1/21r) f.21r[S(H +J)/(J2+S )~]dX" (360) 

0 

S = (P2/ P1)L' sin i cos X'1(1 + T sin2 X")/ 

(1 + W sin2 X")(1 +Q sin2 X")J 

_ y'1 + Q sin2 X" [ 1 + W sin2 X" , . ] H-1 + W sin2 X" (1 +Q sin2 X")2 (P2/P1)L cost 

L'=(dL/dX")=(1-e' cos E') /(1-e' )~ 

(361) 

(362) 

(363) 
82 SPACE OBLIQUE MERCATOR PROJECTION 

tan lhE'=tan lh(X"-w) [(1-e')/(1 +e')]~ (364) 

Equations (352) and (353) are altered only by replacing X' and 4>' with 

X" and 4>", respectively. 

Equations (267), (207), (199), (198), (200), and (219) shown in sequence 

in the summary of circular-orbit equations for F, J, W, Q, T, and U, 

remain unchanged. 

2. For 4> and X in terms of 4>" and X", use equations (269), (270), (271), 

(248), and (354) without change, where L is found from equation (334) 

and E' from (364). 

3. For 4>" and X" in terms of 4>' and X', use equations (261), (262), (258), 

(260), (257), and (133) in the sequence shown in the circular-orbit 

summary, but equation (259) is replaced by (337), where L is calculated 

from (334) and (364) from X", and equation (93) is replaced by 

Xo=Xt -(P2/P1)(L +w) 

where L is calculated from (334) and (333) from X'. The radius vector 

Ro is obtained from equation (336). Fourier series are strongly recommended here, as with the circular-orbit formulas, except that equations (289) and (291) must be supplemented with integrations involving cos n X', as in equations (53) through (56). 

4. For 4>' and X' in terms of 4>" and X", the same Fourier series as that 

used for 4>" and X" in terms of 4>' and X' may be used with reversal of 

signs, since the correction is so small. 

5. For 4>" and X" in terms of 4> and X, use equations (272), (332), and 

(274), finding L from equation (334) and E' from (364). 

6. For 4>" and X" in terms of x and y, use equations (355), (356), (350), 

(351), and (359) through (364) as they stand. Equations (357) and (358) 

are altered only by replacing X' and 4>' with X" and 4>", respectively. 

7. For the satellite groundtrack, 

a. For X' in terms of 4J0 , use (236) and (256), but determine Ro from 

(336), and E' from (364), using X' in place of X". Iteration by 

successive substitution is required for the group of equations. 

b. For X' in terms of Ac,, use (135), (337), (334), and (364) with 

iteration and quadrant adjustment essentially as followed for 

circular orbits. 

c. For 4>o in terms of X', use (257), determining Ro from (336) and 

(364). 

d. For Xo in terms of X', use (338) and (133), finding L from (334) 

and (364). For (b), (c), and (d), X' should replace X" in (364). 

A check of scale factors for various values of X", taken at constant 4>" 

(parallel to the groundtrack) and constant X" (along a scanline), shows 

correct scale along the groundtrack and along the scanline at the 

groundtrack, regardless of orbital eccentricity. The scale along the 
SOM EQUATIONS FOR NONCIRCULAR ORBITS 83 

scanlines is also the same as it is for a circular orbit with the same 

inclination, with a slight variation from sec cj>" due to ellipsoidal approximations. The scale parallel to the ground track, but 1 o away from 

it, deviates from sec c/>" by a quantity increasing with orbital eccentricity, when using the formulas in this section. These formulas are 

sufficient only if the eccentricity is as low as 0.05. For greater eccentricities the next section should be consulted. 

In any case, as stated before, the groundtrack is shown correctly with 

respect to geodetic latitude and longitude and to scale. The fact that it 

is not curved quite correctly on the map affects the scale away from it, 

but not along it. 

The rapidly converging and minimal number of Fourier constants 

sufficient for circular orbits is not the case with noncircular orbits. In 

the case of a fictitious orbit with i and P2/ P1 the same as those for 

Landsat, and an orbital eccentricity of 0.1 and w of 20°, it was found 

that 13 coefficients were needed for the x' integral and 12 for the y' 

integral, instead of 3 and 2 respectively for the corresponding functions of the circular orbit, to provide 6- to 7-place accuracy in computing x and y. Assuming that similar numbers of coefficients are 

needed for the non-integral functions, it would be more efficient to 

calculate the latter directly without Fourier coefficients and to reserve 

the Fourier for the integrals only. 

The approach in determining the Fourier coefficients is similar, of 

course, to preceding derivations, except that more terms are involved. 

Taking the term under the integral in equation (352), using X" in place 

of X', 

F("") HJ-S2 

1 + ~ ""+ ~ b . "" 1\ (J2+S2) 112 2ao n~1 an cos n 1\ n~ 1 n sm n 1\ 

where 

a =! f

2

1r F (X")d X" 0 7r Jn 0 

a =!J.21r F (X") cos n X"dX" n 7r 

0 

b, =~f. F (X') sin n X'd X' 

Integrating, 

I n n 

F (X")dX"=~oX"+ n~l (anln) sin nX"- n~/bnln) cos nX"+C 

where C is the constant of integration, which will be determined to 

make the integral zero when X" is zero, by making it the sum of the 
84 SPACE OBLIQUE MERCATOR PROJECTION 

coefficients of cos n>-.". Replacing ~ac, with Bh (anln) with Am and 

(bnln) with A'n, 

i

x• HJ-S2 n 

0 

(J2+S2)172 d>-."=B1>-."+ n~ An sin n>-." 

n n 

- ~ A'n cos n>-."+ ~ A'n (365) n-1 n=1 

where 

1 J.2,.. HJ-S2 

B1 = 27r 

0 

(J2+ s2)112 d>-." (359) 

1 r.,.. HJ-82 

An= 1rn 0 (J2+S2)112 cos n>-."d>-." (366) 

A I - 1 r.,.. HJ -S2 . ""d"" (367) n-1rn 0 (J2+S2)172 Sln n" " 

It will be noted that B1 is the same term as that used in the rotation 

of axes earlier. An analogous approach leads to Fourier constants for 

the integral term for y' in equation (353): 

where 

i

x· S(H +J) d""-B ""+ ~ C • "" 0 (J2+S2)112 " - 2" n~1 n Sln n" 

n n 

- ~ C'n cos n>-."+ ~ C'n (368) n-1 n=1 

SPACE OBUQUE MERCATOR PROJECTION 

FOR THE GENERAL ELLIPTICAL ORBIT 

(360) 

(369) 

(370) 

As stated in the preceding section, the semi-empirical formulas used 

for the curvature of the groundtrack with a circular or near-circular 

orbit are not satisfactory if the eccentricity of the orbit exceeds about 

0.05. While the formulas relating q, and >-.to cp" and >-.", as derived above 

for the elliptical orbit, may be used here without change, it is important to incorporate the constraint pointed out by Junkins, namely 

that the groundtrack on the map should have the same radius of 
PROJECTION FOR THE GENERAL ELLIPTICAL ORBIT 85 

curvature as the ground track on the surface of the ellipsoid as viewed 

vertically above the surface. 

Since tests using subsequent formulas show that the use of this 

constraint does not appreciably improve the scale accuracy of the SOM 

for the circular orbit over the formulas derived already in this work, 

the formulas based on radius of curvature are proposed only for the 

general elliptical orbit, to convert c/>" and A" to x andy and the inverse. 

They may be fairly readily simplified for the circular orbit and sphere 

as desired. 

For at least two reasons the radius of curvature restraint (combined 

with a groundtrack true to scale) does not alone produce a perfectly 

conformal projection: (1) The position c/>" along the scanline is calculated by somewhat empirical formulas for the ellipsoid independent of 

the groundtrack curvature, and (2) the scanlines are assumed to be 

straight, but this is not exactly correct, even for the sphere. 

In figure 12, line ABCE is a satellite ground track as plotted on a flat 

map with origin at A, the Y axis arbitrarily tangent to the groundtrack at A, and the X axis perpendicular. Element BC, of length ds, has 

a radius of curvature BF=CF=rco Angle d8, CFB, is 

d8=ds/rc (371) 

Angle 8 is the cumulative or integral of d8's from A to B: 

8= f.' d8 

= J.8 (ds/rc) 

0 

To relate 8 to A", the expression is multiplied by dA"IdA": 

8= J.x· (1/rc) (ds/ dA")dA" 

0 

(372) 

Element BC has a slope such that angle BCD equals 8, and since BD 

and DC are elemental changes of x andy, respectively, in moving from 

B to C, 

and 

dx=BD=BC sin 8 

=ds sin 8 

=(ds/dA") sin 8 dA" 

dy=DC=BC cos 8 

=ds cos 8 

=(ds/dA") cos 8 dA" 

(373) 

(374) 

Since ds =[(dx )

2+ (dy )2]~, for which the X and Y axes may be in any 

perpendicular orientation, we may use equations (339) and (340) or 

earlier derivations to show that 

(ds/dX")=a (H +S )~ (375) 
86 SPACE OBLIQUE MERCATOR PROJECTION 

Y axis 

FIGURE 12.-Elements of groundtrack curvature. 

where H and S are found for the non-circular orbit from (361) and 

(362) in terms of >..". 

Combining (372) through (375), and integrating, the basic equations 

for the groundtrack are found: 

x/ a = s.•· (H2+ 82)! sin 8 d >." (376) 

y/a = s.•· (H2+S2)l cos 8 d>." (377) 

where 8=a s.•· [(H2+S2)h,]d>." (378) 

For reasons which will develop later, the X and Y axes will not be 

rotated from the orientation shown in figure 12, except for convenience 

in certain cases. 
PROJECTION FOR THE GENERAL ELLIPTICAL ORBIT 87 

Formulas for H and S have been developed, but the radius of curvature rc is yet to be derived for the general case. The practical evaluation 

of Fourier constants for the double integral also must be delineated. 

Before this is done, the adjustment for position along the scan lines 

will be derived. 

The slope of the ground track in figure 12 at B is dy/ dx, and from 

(373) and (374), 

Thus 

dy/ dx =cos 0/sin 0 

=1/tan 0 

LCBD=arctan (1/tan O) 

Angle JBC, between scanline JBH and the groundtrack, must be the 

same as that derived for the circular orbit, first for the sphere and 

then for the ellipsoid, except for using equations (361) and (362) for H 

and S. From equation (46) preceding either (47) or (206), since the slope 

of the groundtrack with vertical scanlines is (S/H), 

LJBC=arctan (H/S) 

For the slope of the scanline with respect to the X axis of figure 12, we 

take the tangent of LJBD, which is the sum of LJBC and LCBD, or 

(dy/dx)sz=tan LJBD=tan (LJBC+LCBD) 

=(H/8+1/tan 0)/(1-H/S tan 0) 

=(H sin O+S cos 0)/(S sin 0-H cos 0) (379) 

This function is applied to equations (376) and (377) by adding 

cos LJBD times In tan ('n/4+cJ>"/2) to (376) and adding sin LJBD times 

In tan (7r/4+cJ>"/2) to (377). To obtain cosLJBD, 

cos LJBD=1/(tan2 LJBD+ 1)~ 

=(S sin 0-H cos 0)/[(S sin 0-H cos 0)2 

+(H sin O+S cos 0)2]~ 

=(S sin 0-H cos 0)/(8 2 sin 2 0-2SH sin 0 cos 0 

+ H 2 cos 2 O+ H 2 sin2 0+2HS sin 0 cos 0+82 cos 2 O)~ 

=(S sin 0-H cos O)/(H +S )~ (380) 

sin LJBD =tan LJBD cos LJBD 

=(H sin O+S cosO) (S sin 0-H cos 0)/ 

(S sin 0-H cos O) (H2 + )~ 

=(H sin O+S cos 0)/(H +S )~ 

Equations (376) and (377) then become, for the complete projection, 

f.]+.• 

x/a = 

0 

(H +S )~ sin 0 dX" 

-[(H cos O-S sin O)/(H +S )~] In tan (7r/4+cJ>"/2) (381) 
88 SPACE OBLIQUE MERCATOR PROJECTION 

yla= [" (H'+S2)l cos 8 dll' 

+[(H sin 8+8 cos 8)/(H +S )~] In tan (1r/4+<t>"/2) (382) 

These equations, looking simpler than they are in practical use, do not 

simplify when applied to the circular orbit for sphere or ellipsoid. 

Instead, the calculations of H, S, 8, </>",and~" from <1> and~ (as well as 

the Fourier equivalents) simplify for the latter cases. 

Following is the derivation for rc, the radius of curvature used to 

find 8 in equation (378). This will follow the principles used to obtain 

equation (152) for the radius of curvature for the sphere and a circular 

orbit. First the equations for rectangular coordinates of the ellipsoidal 

orthographic projection will be obtained. 

In figure 13, an octant of the ellipsoid is shown, with semi-major axis 

AD=a. The rectangular coordinates of random point K, with ellipsoidal coordinates <1> and ~' and of the proposed center G of the tangent 

plane GL, at latitude <Po and longitude 0, are obtained relative to the X', 

Y', and Z' axes from earlier formulas, reversing the functions of ~ 

since ~=0 on the X' axis rather than the Y' axis. For point K, 

KJ =x'=a cos <1> cos ~/(1-e 2 sin2 <t>)i (176a) 

SK=y'=a cos <1> sin V(1-e 2 sin2 </>)~ (223a) 

JD=z'=a (1-e 2

) sin <f>/(1-e 2 sin2 </>)~ (177) 

For point G, these formulas reduce to 

x~ =a cos <Pol (1-e2 sin2 </> )~ 

y~=O 

z~=a(1-e ) sin <f>0 /(1-e2 sin2 </>0 )1 

(383) 

(384) 

(385) 

On the tangent plane, with the Y axis pointing north and the X axis 

as shown in figure 13, the rectangular coordinates are as follows: 

the x coordinate is the same as y', or 

x =a cos <1> sin ~/(1-e2 sin2 <t>)l (386) 

while y equals GL in the equatorial view. To find this, 

tanLKGM = KM/GM 

=(z'-z~)l(x~-x') 

LLGM=90°-</>o 

GL=GK cosLLGK 

=[(z'-Zo') +(x~-x')2]~ cos (LLGM -LKGM) 

=[ (z'-z~2 + (x~- x')2]l![1 + tan2(LLGM-LKG M) ]l 

(387) 
PROJECTION FOR THE GENERAL ELLIPTICAL ORBIT 89 

X' axis 1 

I 

I 

I 

I 

cpo 

NORTH 

POLE 

Y' axis 

Z'axis 

NORTH 

E POLE 

~~--:---:------lj 

EQUATORIAL VIEW 

FIGURE 13.-0rthographic projection of the ellipsoid. 
90 SPACE OBLIQUE MERCATOR PROJECTION 

For the trigonometric term, 

tan (LLGM-LKGM)=(tanLLGM -tanLKGM)/ 

(1 +tanLLGM tanLKGM) 

=[1/tan cl>o -(z'-z:JI(x'o -x')V 

[1 +(z'-z:J/(x~-x')tan cl>o] 

=[(x~-x')-(z'-z:J tan cl>oV 

[(x~-x') tan cl>o+(z'-z:J] (388) 

Substituting into (387) and multiplying through by the denominator of 

(388), 

GL=y=[(z'-z:J +(x~-x')2J [(x~-x')tan cl>o+(z'-z:JV 

i { [(x~- x')tan cl>o + (z'- z:Jf+[ (x~- x') 

-(z'-z:J tan cl>of~l 

Expanding the new denominator, canceling like terms and combining, 

y=[(z'-Zo') +(x~-x')1 [(x~-x') tan cl>o+(z'-z:JV 

sec cpo[(x~-x') +(z'-z:J1 

=(x'o-x') sin cl>o+(z'-z'o) cos cl>o 

Substituting from (176a), (177), (383), and (385), and canceling one pair 

of terms, 

y =[a/(1-e2 sin2 q,)lJ(1-e 2

) cos cl>o sin cp-sin cl>o cos cp cos ~] 

+ae2 sin c/>0 cos cpof(1-e2 sin2 cl>o)l (389) 

To satisfy the formula for radius of curvature, 

rc=[1 +(dy/dx )2]812/(d2y/dx2

) (138) 

equations (386) and (389) are differentiated with respect to q,, and the 

latter derivative is then divided by the former for dy I dx ( cl>o is a 

constant): 

dx/dcJ>=a[ -(1-e2 sin2 q,)l sin cJ> sin ~-(cos cJ> sin ~)(~)(1-e 2 sin2 cJ>)-l 

( -2e2 sin cJ> cos cJ>)]/(1-e2 sin2 cJ>)+a cos cJ> cos Md~/dcp)/ 

(1-e2 sin2 q,)l 

=a[(1-e2 sin2 f) cos cJ> cos Md~/dcJ>)-(1-e ) sin cJ> sin ~V 

(1-e2 sin2 cJ>)! (390) 

dy/dcJ>=[a/(1-e2 sin2 cJ>)l] [(1-e2

) cos cl>o cos c~>+ sin cl>o sin cJ> cos ~ 

+ sin cl>o cos cJ> sin ~(d~/dcJ>)]+a[(1-e ) cos cl>o sin cJ> 

- sin cl>o cos cJ> cos ~]( -~)(1-e2 sin2 q,)-1( -2e2 sin cJ> cos c/>) 

dyldx =(dy/dcJ>)/(dx/dcJ>) 

= [ (1-e2

) cos cPo cos cJ> + sin cl>o sin cJ> cos ~ + sin cPo cos cJ> sin ~( d ~/ d cJ>)] 

(1-e2 sin2 cJ>)+[(1-e2

) cos c/>0 sin c~>-sin cl>o cos cJ> cos ~] 

e2 sin cJ> cos cJ> /[(1-e2 sin2 cp) cos cJ> cos ~(d~/dq,) 

-(1-e2

) sin cJ> sin ~] 
PROJECTION FOR THE GENERAL ELLIPTICAL ORBIT 91 

= [(1-e2

) cos cl>o cos cJ>+ sin cl>o sin cp cos "A+ sin cl>o cos cp sin "A(d"Aidcp) 

-e2(1-e 2)cos cl>o cos cp sin2 cp-e2 sin cl>o sin3 cp cos "A 

-e2 sin cl>o cos cp sin2 cp sin "A(d"A/dcp)+e 2(1-e 2

) sin2 cp cos cp cos cl>o 

-e2 sin cl>o cos2 cp sin cp cos "A]/[(1-e2 sin2 cp) cos cp cos "A(dAidcp) 

-(1-e2) sin cp sin "A] 

= [(1-e2)(cos cl>o cos cJ>+ sin cl>o sin cp cos "A)+(1-e 2 sin2 cp) 

sin cl>o cos cp sin "A(d"A/dcp)]/[(1-e2 sin2 cp) cos cp cos A(dA!dcp) 

-(1-e2

) sin cp sin "A] (391) 

At the point of tangency, where the radius of curvature is being determined, cJ>o=cJ>, and "A=O. Equation (391) simplifies to 

(dyldx) 0=(1-e2)/(1-e2 sin2 cp) cos cp(d"A/dcp) (392) 

For (d2yldx2), equation (391) rather than (392) must be differentiated 

with respect to cJ> and divided by (dx/dcp) from (390): 

d2yldx2=[d(dyldx)ldcp]/(dxldcp) 

= {[(l-e2 sin2 cp) cos cp cos "A (d"A/dcp)-(1-e 2

) sin cJ> sin "A] 

t (1-e2)[ -cos cl>o sin cp + sin cl>o cos cp cos "A 

-sin cl>o sin cp sin "A (dAidcp)] 

-[sin cl>o sin cp sin A(d"A/dcp)-sin cl>o cos cp cos "A (d"Aidcp) 2 

-sin cl>o cos cp sin "A (d2Aidcp2)](1-e2 sin2 cp) 

+sin cl>o cos cp sin "A (d"Aidcp)( -2e2 sin cp cos cp)} 

-[(1-e2)(cos cl>o cos cp + sin cl>o sin cp cos "A)+(1-e2 sin2 cp) 

sin cl>o cos cp sin A(d"Aidcp)] (1-e2 sin2 cp) 

[-sin cp cos "A(dAidcp) 

-cos cJ> sin A (d"Aidcp)2+ cos cp cos A (d2"A/dcp2)] 

+ cos cJ> cos "A (d"Aidcp)( -2e2 sin cJ> cos cp) 

-(1-e2)[cos cp sin "A+ sin cp cos "A (d"A/dcp)]}}{l-e 2 sin2 cp)i/ 

a[(l-e2 sin2 cp) cos cp cos A(d"A/dcp)-(1-e2) sin cp sin "A]3 

Since this equation does not need to be further differentiated, the value 

at tangent point cp=cJ>o and "A=O may be determined now to reduce 

handling of the equation. 

(d2yldx2

) 0= {{1-e2 sin2 cp) cos cp (d"Aidcp)[(1-e 2)( -cos cp sin cp 

+ sin cp cos cp) 

+sin cJ> cos cJ> (d"Aidcp)2(1-e 2 sin2 cp)] 

-(1-e2)(cos2 cJ> + sin2 cp)[(l-e2 sin2 cp)( -sin cJ> (d"Aidcp) 

+cos cJ> (d2"Aidcp2))+ cos cJ> (d"A/dcp)( -2e2 sin cJ> cos cp) 

-(1-e2

) sin cJ> (d"A/dcp)]}(1-e 2 sin2 cp)l/ 

a[1-e2 sin2 cp) cos cJ> (d"A/dcp)"f 

After intermediate steps, this "simplifies" to 

(d2yldx2

) 0= {(1-e2

)[2 sin cp (d"Aidcp)-cos cp(d2"Aidcp2) 

-e2 sin cp(3 sin2 cp(dAidcp)- sin cp cos cp(d2"Aidcp2

) 

-(d"Aidcp))]+(d"Aidcp)3 (l-e2 sin2 cp)2 sin cJ> cos2 cp}/ 

a(1-e2 sin2 cp)i cos3 cp(d"A/dcp)3 (393) 
92 SPACE OBLIQUE MERCATOR PROJECTION 

Before substituting (392) and (393) in (138), equation (138) for rc should 

be inverted to give 11 rc, which shall be called (-Ri ). This results in 

finite values throughout, instead of infinity at the points of inflection 

of the groundtrack at the equatorial crossings or nodes. It is also 

advisable to multiply (d'Aid(j>) and (d2

'Aid4>2), derived subsequently, by 

cos 'A" and cos3 'A", respectively, and to multiply the numerator and 

denominator of the equation for Ri by cos3 'A", thus avoiding indeterminate values or dividing by zero. The sign will also be reversed to give 

initial positive values of~ consistent with figure 12. Subscript zero 

should be applied to 4> and 'A to denote groundtrack conditions. As a 

further simplification, let a =1. 

Therefore, let 

Substituting, 

G =cos 'A"(d'Ac,ld4>o) 

M =cos3 'A"(d2Ac,/d(j>0

2

) 

C =l-e2 sin2 4>o 

K=l-e2 

(394) 

(395) 

(396) 

(397) 

Ri= -lire= -(d2yldx2)of[l +(dyldx)!.]l 

= - { K [2G sin 4>o cos2 'A"-M cos 4>o- e2 sin 4>o 

(3G sin2 4>o cos2 'A"-M sin 4>o cos 4>o 

- G cos2 'A")]+ G3C2 sin 4>o cos2 4>o } I 

G3cl cos3 4>o (1 + K2 cos2 'A" /C2G2 cos2 4>o)i 

= -{ K[G sin 4>o cos2 'A"(2-3e2 sin2 4>o +e2

) 

- M cos 4>o (l-e2 sin2 4>o)]+G3C2 sin 4>o cos2 4>o} I 

( G2 C cos2 4>o + K2 cos2 'A" IC )i 

= - { K[ G sin 4>o cos2 'A"(3C-K)-MC cos 4>o] 

+G3C2 sin 4>o cos2 4>ol/(G2C cos2 4>o 

+ K2 cos2 'A" IC)i (398) 

This formula for the reciprocal of the radius of curvature applies to 

any path on the surface of the ellipsoid, if G and M can be determined 

for the specific path desired. It is now necessary to derive (d'Ac,/d4>o) 

and (d "Aald4>~) for the elliptical satellite orbit. For simplicity, since the 

geocentric ground track rather than the vertical groundtrack has been 

used for determining 4>", 'A", H, and S, the same will be used here. The 

slight correction for the vertical groundtrack will be made as before, 

as described in the summary of equations for a slightly elliptical orbit 

in the preceding section. 

For the geocentric groundtrack, equations for 4>o in terms of 'A" are 

unchanged from (195) and (196), using 'A" in place of 'A': 

tan 4>o=tan (j>0 /(l-e2

) 

sin 4>0 =sin i sin 'A" 

(195) 

(196) 
PROJECTION FOR THE GENERAL ELLIPTICAL ORBIT 93 

Equations for longitude for geocentric or vertical groundtrack were 

derived in the preceding section, again using >-.": 

'Ao='Ato-(Pz/Pl)(L+w) (338) 

tan At 0 =tan 'A" cos i (133) 

L = E'-e' sinE' (334) 

tan~ E'= tan~ ('A"-w)[(1-e')/(1 +e')]112 (364) 

L'=(dL/d>-.")=(1-e' cos E')2/(1-e'2

)

112 (363) 

Combining (338) and (133) and differentiating, 

>-."=arctan (tan 'A" cos i)-(P2/P1)(L+w) 

(d'Ao/d'A")=cos i sec2 'A"/(1 +cos2 i tan2 'A")-(P2/PI)(dL/d'A") (399) 

Rearranging (196) and differentiating, 

sin >-."=sin </>glsin i (236) 

cos 'A"(d'A"/d<Po)=(cos </>g/sin i)(d</>gld<Po) (400) 

Differentiating (195), 

sec2 <Pod <Po =[sec2 </>gl(1-e2

)]d</>g 

(d</>gld<Po)=(1-e2)(tan2 <Po+1)/(tan2 </>g+1) 

= (1-e2)[tan2 </>g/(1-e 2

)

2+ 1]/(tan2 </>g + 1) 

Simplifying and assigning the symbol D, 

D=(d</>gld<Po)=[1-e2(2-e2

) cos2 </>g]/(1-e 2

) (401) 

Rearranging (400) and multiplying (d'A/d'A") from (399) by (d'A"Id<t>o) 

from (400), 

(d'Ao/d<t>o)=[cos i/cos2 >-."(1 + cos2 i tan2 'A")-(P2/P1)(dL/d'A")] 

(cos </>glsin i cos >-.")(d</>gld</>0 ) 

The first denominator may be changed to (cos2 'A"+cos2 i sin2 >-."), 

(1-sin2 'A"+cos2 i sin2 >-."), and (1-sin2 i sin2 >-."). Substituting from 

(196), 

(d'A/d<Po)=[cos i/cos2 </>g-(P2/P1)(dL/d>-.")] (cos<j>g/ 

sin i cos 'A")(d</>gld</>0 ) (402) 

To find G of equation (394), 

G =cos 'A"(d'Aold<Po) 

=[cos i/cos2 </>g -(P2/P1)(dL/d'A")](cos </>glsin i) 

(d</>gfd</>0 ) (403) 

where the differentials are found from (363) and (401) above. 
94 SPACE OBLIQUE MERCATOR PROJECTION 

For (d2"Ao/d¢0

2

) equation (402) is differentiated: 

(d2"Aoldcp0

2

)=[cos i/cos2 cf>u -(P2/P1)(dLidA")]{sin i cos A" 

[-sin cf>u(dct>uldc/>0 )

2

+cos cf>u(d 2

ct>uld¢0

2

)] 

+cos cl>u(dct>uldcp0 ) sin i sin A"(dA"Idct>oHI 

sin2 i cos2 A"+[cos cf>u(dct>uldc/>0 )/sin i cos A"] 

[2 cos i sin cl>u(dct>uldc/>0 )1cos3 cf>u 

- (P2/ P1)(d2LI dAY2)(dA" I dct>o)] 

Substituting from (400) for the last differential, and combining some 

other terms, 

(d2"Ao/dcp0

2

)=[cos i/cos2 cf>u -(P2/P1)(dLidA")] {-sin cl>o 

(dct>uldct>o)2

+cos cf>u(d2

¢uldc/>0

2

) 

+sin cl>u cos2 cf>u (dct>uldcp0 )

2

/sin2 i cos2 A"}/ 

sin i cos A"+2 sin cf>u(dcf>uldc/>0 )

2/tan i cos A" cos2 cf>u 

-(P2/Pl)(d2

LidAY2) cos2 cf>u(dct>uldct>o)2

/sin2 i cos2 A" 

To obtain M in equation (395), 

M =cos3 A"(d "Ao/d¢~) 

M ={[cos i/cos2 cf>u -(P2/P1)(dL!dA")] [cos cf>u cos2 A"(d ct>uldct>~) 

+sin cf>u(dct>uldct>o)2(cos2 cl>ulsin2 i-cos2 A")] 

+2 sin cf>u cos i cos2 A"(dcf>uld¢0 )

2

/cos2 cf>u 

-(P2/Pl)(d2LidAY2) cos2 cl>u cos A" (dct>uldc/>0 )

2

/sin i }/sin i 

The expression (cos2 cf>u/sin2 i -cos2 A") may be reduced by using (196): 

cos2¢ulsin2 i -cos2 A"=(cos2 cl>u -sin2 i +sin2 A" sin2 i)/sin2 i 

=( cos2 cf>u -sin2 i + sin2 cf>u)lsin2 i 

=1/tan2 i 

Thus 

M={[cos i/cos2 cf>u-(P2/P1)(dLidA")] [cos cf>u cos2 A"(d ¢uldct>~) 

+sin cf>u(dcf>uldct>o)2/tan2 i]+2 sin cf>u cos i cos2 A"(dct>uld¢0 )

2/ 

cos2 cf>u -(P~P )(d LidAY2) cos2 cf>u cos A"(dct>uld¢0 )

2

/sin i }/sin i 

(404) 

Solution of ( 404) requires two more second derivatives. First differentiating (dct>uld¢0 ) in equation (401) and assigning the symbol D', this 

is simply 

D'=(d ¢uldct>~)=[21(1-e )]e 2 (2-e2)sin cl>u cos cf>u(dct>uldct>o) (405) 

Differentiating (363), 

(d2LidAY2)=[2/(1-e'2)f](1-e' cos E') e' sinE' (dE' IdA") (406) 

Repeating equation (341a) in the previous section, 

(dE' ldA")=sec2 ~ (A"-w) cos2 ~ E' [(1-e')/(1 +e')]i (341a) 

From (333), 

tan2 l (A"-w)=tan2 l E' (1 +e')/(1-e') 
PROJECTION FOR THE GENERAL ELLIPTICAL ORBIT 95 

Thus sec2 l (A"-w)=1 +tan2 l E'(1 +e')/(1-e'2

) 

Substituting this in (341a) and the latter in turn in (406), and calling 

this L", 

L"=(d2LidA'12)=[2/(1-e 12)i] (1-e' cos E') e' sinE' cos2 ~ E' 

[(1-e')/(1 +e')]i [1 +tan2 iE'(1 +e')/(1-e')] 

=[2/(1-e 12 )~] (1-e' cos E') e' sinE' { cos2 l E' [(1-e')/ 

(1 +e')]i+sin2 ~ E' [(1 +e')/(1-e')]i} 

=[1/(1-e 12)l] (1-e' cos E') e' sinE' [(1 + cos E')(1-e') 

+(1-cos E')(1 +e')]/(1-e12)i 

=2e' sinE' (1-e' cos E')2/(1-e'2

) (407) 

Now the fundamental() can be found with equation (378) by numerical 

integration: HandS are found in terms of A" from (361) and (362) and 

( -1/rc) or Ri from (398). To find (398) in terms of A", equations (396), 

(397), (196), (195), (401), (405), (364), (363), (407), (403), and (404) must 

be used, as well as (198), (199), and (200). These formulas will be 

assembled later, and they present no serious computing problem except that they are lengthy. 

Fourier constants are again almost indispensable for practical computation. The usual procedure is described with equations (53) and 

those following. Simpson's rule is recommended for numerical integration. Instead of an integrating interval of 9°, as used for the circular 

orbit, an interval of 3° is recommended for the noncircular orbit, and 

integration must be carried out for the full 360° cycle due to lack of 

symmetry of the groundtrack for each quadrant. 

Three types of constants must be determined, and none are zero for the 

general case, except that convergence leads ultimately to negligible 

coefficients: 

where 

(J=gA"+ ~dn sin nA" + ~ fn(1-cos nA") (408) 

g=(l/2.-) [R,(H'+S')idx• 

d,. =(1/...n) fR;(.H'+S')! cos nX'dX' 

f. =(lhrn) r·R,(H'+S")I sin nX'dX' 

(409) 

(410) 

(411) 

For circular orbits g and all dn, coefficients are zero, as are coefficients 

In for even values of n. For elliptical orbits with perigee w coinciding 

with either ascending or descending node, g is zero. The number of 

terms required for convergence increases with eccentricity: for a circular orbit fi, fa, and is are sufficient for six-place accuracy; for a hypo-
96 SPACE OBLIQUE MERCATOR PROJECTION 

thetical eccentricity of 0.2, about 18 coefficients (g, dn, and fn) are 

required for the same accuracy, and for an eccentricity of 0.7, about 40. 

Once the coefficients are determined, the second numerical integration 

may occur in order to evaluate the Fourier coefficients for the integrals 

of (376) and (377), or (381) and (382). If the coefficient g from (409) is 

zero, the coefficients are determined in the same way as the coefficients 

for fJ or for circular orbits in earlier sections of this study. For the 

elliptical orbit, it is highly unlikely that the perigee will coincide with 

a node. Therefore, in general, the simple forms of (376) and (377) must 

be expanded as follows: 

A requirement of Fourier coefficient determinations is the repeatability of the function with every cycle, as in equation (57). If g is not 

zero, fJ does not return to its original value at each 360° cycle. It is 

necessary, then, to separate the repeatable portion of ( 408) from the 

nonrepeating portion. This may be done by letting 

()='Y+g'A" (412) 

where 'Y=~ dn sin n'A"+~fn (1-cos n'A") (413) 

The integral for (xl a) from (376), for the ground track only, may thus 

be written 

(xlah= J.'"(H+S ~ sin (-y+gX")dX" 

f.,.. = 

0 

(JI2+S )~ (sin 'Y cos g'A" + cos 'Y sin g'A")d'A" 

f.,.. = 

0 

[(JI2+S2)i sin 'Y] cos g'A"d'A" 

+ J.'"[(If2+S2)! cos 'Y] sin gX" d>." (414) 

The bracketed expressions repeat every 360°, and Fourier series may 

be determined for each. (A similar equation for y/a will follow later.) 

Let 

and 

where 

(JI2+S2)l sin 'Y =b + ~ h~ sin n'A"+ ~ hn cos n'A" (415) 

(JI2+S2)l cos 'Y =b' + ~ Cn sin n'A"+ ~ c~ cos n'A" 

[

11' 1 

b = (1/2?r) (JI2 + S2)1 sin 'Y d 'A" 

0 

b'=(l/2r) f.<H"+S")i cos 'Y dX' 

hn =(1/?r) i 211'(JI2+S2)l sin 'Y cos n'A"d'A" 

0 

(416) 

(417) 

(418) 

(419) 
PROJECTION FOR THE GENERAL ELLIPTICAL ORBIT 97 

h', =(1/,·) ['<H"+S')l sin 'Y sin nA"dA" (420) 

c~=(1hr) ['<H"+S')! cos 'Y cos nA"dA" (421) 

c,. =(11.-) ['<H"+S')! cos 'Y sin nA"dA" (422) 

Note that the left sides of (415) and (416) are not integrals, while() of 

equations (378) and ( 408) involves an integral. Consequently the form 

of these coefficients is somewhat different from that of the constants 

in (409) through (411)., They may be determined, however, with the 

same application of Simpson's rule as that for all the other Fourier 

coefficients developed in this bulletin. 

Now (415) and (416) may be substituted into (414): 

(xlah = r· [b cos g'A"+ ~ h~ sin n'A" cos g'A" 

0 

f.~· 

+ ~ hn cos n'A" cos g'A"]d'A"+ [b' sin g'A" 

0 

+~en sin n'A" sin g'A"+ ~ e~ cos n'A" sin g'A"]d'A" 

This equation may be integrated using standard trigonometric forms: 

(xla)1={(blg) sin g'A"+~ h~ [-cos (n-g)'A"/2(n-g) 

-cos (n +g)'A"/2(n +g)]+~ hn[sin(n -g)'A"/2(n -g) 

+sin (n+g)'A"/2(n+g)]-(b'lg) cos g'A" 

+~ en[sin (n-g)'A"/2(n-g)-sin (n+g)'A"/2(n+g)] 

+~ e~cos (n-g)'A"/2(n-g)-cos (n+g)'A"/2(n+g)]}x· 0 

= {(1/g)(b sin g'A"-b' cos g'A") 

- ~ [( -hn +en)(sin n'A" cos g'A"+cos n'A" sin g'A") 

+(h~+e~)(cos n'A" cos g'A"-sin n'A" sin g'A")]/2(n +g) 

- ~ [( -hn -en)(sin n'A" cos g'A"-cos n'A" sin g'A") 

+(h~-e~)(cos n'A" cos g'A"+sin n'A" sin g'A")]/2(n -g)} x· 

= { (1/g)(b sin g'A"-b' cos g'A")- ~ [( -hn +en) cos g'A"/ 

2(n +g)-(h~+e~) sin g'A"/2(n +g) 

+(hn +en) cos g'A"/2(n -g)-(h~-e~) sin g'A"I 

2(n -g)] sin n'A"+ ~ [( -hn +en) sin g'A"/2(n +g) 

-(h~+e'n) cos g'A"/2(n +g)+(hn +en) sin g'A"I 

2(n-g)-(h~-e~) cos g'A"/2(n-g)]1x· 0 

0 
98 SPACE OBLIQUE MERCATOR PROJECTION 

=(1/g)(b sin g)."-b' cos g).")+ ~{[(nc~-gh~)l 

(n2 -g2)] sin g"A"+[(nhn+ucn)l(n2-g2)] cos g)."}sin n)." 

- ~ l[(ncn+ghn)l(n2 -g2)] sin g"A"+[(nh~-gc~)l 

(n2 -g2)] cos g)."} cos n"A"+b' /g 

+ ~ (nh~-gc~)l(n -g ) (423) 

The last two terms are the negative value of the integral when )."=0. 

It is convenient to determine new constants to replace the coefficients 

of (423): 

B=b/g (424) 

B'=b'/g (425) 

An =(nhn +ucn)l(n2 -g2

) (426) 

A~=(nc~-gh~)l(n -g ) (427) 

Cn =(ncn +ghn)l(n2 -g2

) (428) 

c~ = ( n~-gc~)/ ( n 2- g2

) (429) 

Cx=B'+~ C~ (430) 

These symbols may be substituted into (423), but since (423) is the 

Fourier equivalent of (376), the final term of (381) may also be attached for the complete equation for (xla), the equivalent of (381): 

x/a =B sin g)."-B' cos g"A"+ ~ (A~ sin g)."+ An cos g).") sin n"A" 

- ~ (Cn sin g"A"+C~ cos g"A") cos n'A" + Cx 

-[(H cos 8-S sin 8)/(H2+S2)l] ln tan (1r/4+cJ>"/2) (431) 

Using symbols which will also be found useful in later formulas, let 

H'=H cos 8- S sin 8 (432) 

A"= ~(An sin n)."-C~ cos n"A") (433) 

C"=~(A~ sin n)."-Cn cos n'A") (434) 

Then (431) may be rearranged and rewritten as follows: 

xla=(B+C") sin g)."-(B'-A") cos g"A"+Cx 

-[H'/(H2+S2)6] ln tan (1r/4+cJ>"/2) (435) 

In computer programing, care should be taken, of course, to minimize 

the number of calculations of repeated trigonometric terms. 

The Fourier equivalent of equation (382) for (y I a) is derived in an 

analogous manner, so most of the steps are omitted here. From (377) 

and (412), 
PROJECTION FOR THE GENERAL ELLIPTICAL ORBIT 99 

(y/a),= J,' .. <H'+S2)l cos b+uX')dX' 

r>-· = Jo [(H +S )~ cos 'Y] cos gX" dX" 

r>-· - Jo [(W+S )~ sin 'Y] sin gX" dX" 

Substituting (415) and (416), 

r>-· (y/a)1= Jn [b' cos gX"+~cn sin nX" cos gX"+~c~ cos nX" cos gX"]dX" 

0 

r>-· - Jo b sin gX"+ ~~sin nX" sin gX"+ ~hn cos nX" sin gX"]dX" 

With integration as accomplished for (x/a)h combining of terms, and 

final substitution into (382), 

yla=(B'-A") sin gX"+(B+C") cos gX"-Cy 

+[S' /(W+S )~] ln tan (·n/4+cp"/2) (436) 

where coefficients have been previously defined, except for Cy and S': 

Cy=B- ~Cn (437) 

S'=H sin 8+S cos 8 (438) 

Equations ( 435) and ( 436) have the interesting variation from equations for the circular orbit in that, aside from being much more complicated, the central line about which the sinusoidal-like groundtrack 

oscillates is no longer a straight line, but is now a circular arc. Its 

equation is found from (435) and (436) by letting nX" be any multiple 

of 360°, so that sin nX"=O and cos nX"=1, while c/>"=0: 

(xla)o=Cy singX"+Cx (1-cosgX") 

(yla )o=Cx sin gX"-Cy (1-cos gX") 

Transposing constants, 

(xla )o-Cx =Cy sin gX"-Cx cos gX" 

(yla )o+Cy =Cx sin gX"+Cy cos gX" 

Multiplying the first by Cx, the second by Cy and subtracting, then 

multiplying the first by Cy, the second by Cx and adding, the two new 

equations may be squared and added to eliminate the angle. Simplifying the result, 

[(x/a )0 -Cx:f+[(y Ia )0 +Cy f=Cx 2+Cy 2 (439) 

the equation of a circle with radius (Cx +Cy )~ and center at 

(xla=Cx), and (yla=-Cy). Generally, if Fourier constant g, equation (409), is positive, this "axis" curves to the right as the ground track 
100 SPACE OBLIQUE MERCATOR PROJECTION 

proceeds upward (y increasing) from the origin, and left if g is negative. In examples checked, a positive g occurs if w (longitude of perigee) 

is between 180° and 360°, and a negative g, if w is in the first or second 

quadrants. Empirically, g is within 1 percent less than [-(P2/ P1) e' sin 

w sin i]. The "axis" curves approximately 360° x g for each orbital 

revolution. An example of elements of the projection for the general 

elliptical orbit is shown in figure 14. 

If g is zero, for the circular orbit or a perigee coinciding with a node 

(as stated earlier), equations (435) and (436) are indeterminate, since 

there is division by zero in (424) and (425). This affects the terms B sin 

g>..", B' cos g>..", and Cz of equation (435) and the corresponding terms 

of (436). To resolve this, (424) and (425) may be used to rewrite these 

terms as follows: 

B sin g>.."=b (sin g>..")/g 

= b >.." as g approaches zero. 

-B' cos g>.."+Cz=-(b'/g) cos g>.."+(b'/g)+"i:,C~ 

=b'(1-cos g>..")lg+"i:,C~ 

To evaluate (1-cos g>..")lg as g approaches zero, the expression 

may be differentiated with respect to >..": 

d(1-cos g>..")/gd>.."=g(sin g>..")/g=sin g>.."=O, when g is zero. 

Thus (435) is written as follows, if g =0. 

xla=b>.."+"i:,C~+A"-[H'I(W+S )l] In tan (1r/4+¢"/2) (440) 

Similarly, (436) is written, for g =0, 

y/a =b'>.."+ "i:,C" +C"+[S' /(W+S2)l] In tan (1r/4+¢" /2) (441) 

Since, when g =0, the ground track oscillates about a straight line, it is 

preferable to rotate the axes so that the central line is the X axis, or 

y =0 at each 360° of>..". This may be done to obtain revised coordinates 

x' and y' as follows (derivation not shown, but conventional): 

x'=(b'y +bx)/(b12+b2)l (442) 

y'=(by -b'x)l(b12+b2)l (443) 

For the circular orbit, h:t and c,. are zero, and equations (442) and (443) 

may be combined with (440) and (441) to give the following equations. 

These correspond to (263), (264), (282a), and (283a) except that the 

latter are based on empirical curvature, while the following are based 

on calculated radius of curvature: 

x'la=B">.."+"i:,A; sin n>.." 

+[(b'S'-bH')/B"(.H2+S2)1] In tan (1r/4+¢"/2) (444) 

y'la="i:,C; sin n>.." 

+[(bS'+b'H')IB"(H2+S2)i] In tan (1r/4+¢"/2) (445) 
PROJECTION FOR THE GENERAL ELLIPTICAL ORBIT 101 

Y axis 

""--

CENTRAL "AXIS")\ o 

~~ 

~ 

GROUNDTRACK \ ;'\ 

SCAN 

FIGURE 14.-Groundtrack and extended sample scanlines for Space Oblique Mercator 

projection, elliptical orbit. Eccentricity is 0.5, perigee at north polar approach, 

(i=80°, P2=275 min., P1=1440 min.). Three orbits are shown. 
102 

where 

SPACE OBLIQUE MERCATOR PROJECTION 

B"=(b +b' )~ 

A:=(bhn +b'c~)/Bn for even n only 

C;:=(bc~-b'hn)!Bn for odd n only 

(446) 

(447) 

(448) 

and other symbols are as calculated above. Note that 'Y=() if g=O. 

When Landsat constants are applied to these Fourier coefficient calculations (P2/ P1 = 18/251; i = 99.092°, Clarke 1866 ellipsoid), the following values are found: 

h = 8.131909918° 

13=0.0007134417° 

fr, = 0.0000101583° 

J7 = 0.000000017 4 ° 

B"=1.005806754 for X" in radians 

A2"= -0.0010940568 

A4" = -0.0000015160 

As"= -0.0000000027 

Cl"=0.1433202229 

C{ = 0.0000327398 

C5" = 0.0000000046 

b = 0.1422865425 

b'=0.9956916017 

As check points in calculating ]4 for the above, if A."=45°, 

Ri =0.0992990361, and if A."=90°, Ri =0.1407675647. 

INVERSE EQUATIONS 

Necessary forward equations have now been derived, including several 

in the preceding section. For inverse equations, the new equations 

necessary are those to determine c/J" and X" from x and y, or the 

inverses of (435) and (436), and the special inverses of (440), (441), 

(444), and (445). 

For the former, following a procedure analogous to the earlier derivations for the circular orbit, c/J" in (435) and (436) may be eliminated 

by multiplying both sides of ( 435) by S' and both sides of ( 436) by H', 

and then adding the two equations. The terms containing c/J" cancel, 

and 

S'xla+H'yla=(S'B+S'C"+H'B'-H'A") sin gA." 

-(S'B'-S'A "-H'B-H'C") cos gA."+S'Cx-H'C11 

Since g is generally small, sin g A." will vary almost proportionately 

with X" for several orbits past the origin. It will not become insensitive 

to change of A." until g A." approaches 90°. Therefore, it is appropriate 

(and satisfactory from subsequent calculator tests) to separate out the 
PROJECTION FOR THE GENERAL ELLIPTICAL ORBIT 103 

terms containing sin gX" for a calculation of X" with few iterations. 

Conversely, cos gX" changes very little in the first few orbits and may 

be left where it is. Transposing and dividing by the coefficient of sin 

gX", 

sin gX"=[(S'B'-H'B) cos gX"+S'(xla-Cx-A" cos gX") 

+ H'(yla +Cy -C" cos gX")]I[S'(B +C")+ H'(B'- A")] (449) 

This iterates satisfactorily if (y I a) is chosen as the first trial X" in 

( 433), ( 434), and the right side of ( 449), then sin g X" is found, from 

which a new X" may be determined for resubstitution, until the new X" 

varies from the preceding by less than a chosen convergence value. 

While the iterations are few, the calculation is lengthy because (433) 

and ( 434) have to be evaluated for each n at each new trial value of X". 

After X" is determined, the final values of X", H, S, H', A", and C" may 

be substituted in a transposition of equation (435) to find cp": 

ln tan {-n-/4+cp"/2)=[(~+S )~/H'] [Cx -x/a +(B+C") sin gX" 

- (B'-A") cos g X"] ( 450) 

Here the equation for x was transposed rather than the equation for 

y (used for the circular orbits) because the axes (figure 14) for the 

elliptical orbit are almost perpendicular for the first few orbits to 

those for the circular orbits, and x is therefore more sensitive to 

changes in cf>" than y is. 

For inverse equations in the range gX"=90° and greater, equation 

( 449) must be inverted to solve for cos g X", but correcting for quadrant 

of the circular central "axis." This has not been done in view of the 

unlikely need for the equations. 

Inversion of equations (440) and (441) follows a similar pattern, more 

like the earlier inversions for the circular orbit. Again multiplying 

(440) by S', (441) by H', and adding, but now merely isolating the linear 

functions of X", 

X"=[S'(xla- };C~- A")+ H'(yla- };Cn -C")]I(S'b+ H'b') (451) 

in which X" is again found by successive substitution, beginning with 

a trial (X"=yla ). Inverting (440) to find cf>" from the final X" and related 

parameters, 

ln tan (-n-14+cp"/2)=[(H +S )~/H'](bX"+ };C~+ A"-xla) (452) 

Similarly, equations (444) and (445), for a circular orbit, lead to the 

following, except that cf>" should be found from inversion of ( 445) because of rotation of the axes: 

X"=[(bS'+b'H')(x'/a- };A: sin nX")+(b'S'-bH')(y'/a 

-};C: sin nX")]IB"(bS'+b'H') (453) 
104 SPACE OBLIQUE MERCATOR PROJECTION 

solved by the same type of iteration, but starting with a trial 

("A"=x'la), and 

In tan (1ri4+<J>"/2)=[B"(.H2+S2)1!(bS'+b'H')] (y'la 

- 'l;C; sin n "A") (454) 

SUMMARY OF EQUATIONS 

Since the formulas are scattered throughout the last few pages, they 

are assembled· here in the order of calculation. First, for constants 

applying to the entire orbit, the following are calculated just once for 

all the integration: 

Q =e 2 sin2 i/(1-e 2

) (199) 

T=e 2 sin2 i (2-e 2)/(1-e 2

)

2 

W=[(1-e2 cos2 i)/(1-e 2>:f-1 

K=1-e2 

(198) 

(200) 

(397) 

The following are calculated for each 3° of "A", from 0° to 360°, to 

determine Ri, H, and S for numerical integration in subsequent formulas: 

<l>u =arcsin (sin i sin "A") (196) 

<l>o=arctan (tan <l>uiK) (195) 

C=1-e 2 sin2 <l>o (396) 

D=[1-e2 (2-e 2

) cos2 <l>u]IK (401) 

D'=2De2 (2-e 2

) sin <l>u cos <l>uiK (405) 

E'=2arctan {tan l("A"-w) [(1-e')/(1 +e')]i} (364) 

L'=(1-e' cos E')2/(1-e12)l (363) 

L"=2e' sinE' (1-e' cos E')2/(1-e12) (407) 

G =D [cos i/cos2 <l>u-L' (P2/P1)] cos <l>ulsin i (403) 

M = {[cos i/cos2 <l>u-L' (P2/ P1)] [D' cos <l>u cos2 )." 

+ D 2 sin <l>ultan2 i]+2D2 sin <l>u cos i cos2 "A" /cos2 <l>u 

- L"D2 cos2 <l>u cos)." (PiP1)/sin i l /sin i (404) 

Ri=- K[G sin <l>ocos2 "A" (3C-K)-MC cos <l>o] 

+G3C2 sin <l>o cos2 <l>o /(G2C cos2 <l>o+K2 cos2 )."/C)I (398) 

H =[(1 +Q sin2 "A")/(1 + W sin2 "A")]l [(1 + W sin2 "A")/ 

(1 +Q sin2 "A")2 -(P2/P1) L' cos i] (362) 

S=(P2/P1)L' sin i cos "A"[(1 +T sin2 "A")/ 

(1 + W sin2 "A") (1 +Q sin2 "A")]i (361) 
PROJECTION FOR THE GENERAL ELLIPTICAL ORBIT 105 

Numerical integrations: 

g=(1/2r) r. R, (H'+S2)ld>..' 

d,. =(1/n.,.) r. R, (H +S~! cos n A'dX' 

f. =(1/n.,.) r. R, (H'+S"P sin n X'd>." 

'Y="l;fin sin n :\"+~In (1-cos n :\") 

b=(112?r) 12.,.. (H2+S2)i sin 'Y d:\" 

0 

b'=(112?r) 12.,.. (H2+S2)i cos 'Y d:\" 

0 

h. =(1hr) r. (H'+S2)l sin 'Y cos n X'd X' 

h',=(1/,..) [' (H'+S~l sin 'Y sin n X'd X' 

c'.=(l/.,.) r. (H'+S'~ cos 'Y cos n X'd>..' 

c,. =(1/.,.) [' (H'+S2)l cos 'Y sin n A'dX' 

Convenient terms: 

An =(nhn +gcn)l(n2 -g2

) 

A'n =(nc'n-gh'n)l(n2 -g2

) 

Cn =(ncn +ghn)l(n2 -g2

) 

C'n =(nh'n -gc'n)l(n 2 -g2

) 

If g is not zero (orbit elliptical and perigee not at a node): 

B=b/g 

B'=b'/g 

Cx=B'+"l;C'n 

11 =B-~Cn 

(409) 

(410) 

(411) 

(413) 

(417) 

(418) 

(419) 

(420) 

(421) 

(422) 

(426) 

(427) 

(428) 

(429) 

(424) 

(425) 

(430) 

(437) 

If g =0 (perigee at a node or circular orbit), b, b', "l;C'n, and "l;Cn are 

used separately as constants. See equations (444) - (448) for further 

reduction to circular orbits. 

All the above produce constants for the projection. For an individual 

point, to find cJ>" and:\" in terms of cJ> and:\: 
106 SPACE OBLIQUE MERCATOR PROJECTION 

E'=2 arctan tan ~(~"-w) [(1-e')/(1 +e')Jl 

L=E'-e' sinE' 

(364) 

(334) 

~t=~+(P2IP ) (L+w) (332) 

~"=arctan [tan ~t cos i + K tan cJ> sin i/cos ~t] (272) 

The last four equations above are iterated together, following type (2) 

under "Iteration procedures." 

cJ>"=arcsin { (K sin cJ> cos i -cos cJ> sin i sin ~t)l 

(1-e 2 sin2 cJ>)~} (455) 

This equation is like (274), but F is omitted as an unnecessary complication, not affecting accuracy significantly. To find x andy in terms 

of cJ>" and ~", 

8=g~"+ !,fln sin n~"+ !,fn (1-cos n~") (408) 

A"=!,(An sin n~"-C'n cos n~") (433) 

C"=!,(A'n sin n~"-Cn cos n~") (434) 

L'=(1-e' cos E')2/(1-e'2)1 (363) 

H =[(1 +Q sin2 ~")/(1 + W sin2 ~")]I [(1 + W sin2 ~")/ 

(1 +Q sin2 ~") -(P /P )L' cos i] (362) 

S =(P2/P1)L' sin i cos~" [(1 + T sin2 ~")/ 

(1 + W sin2 ~") (1 +Q sin2 ~")]~ (361) 

H'=H cos 8-S sin 8 (432) 

S'=H sin 8+8 cos 8 (438) 

If g is not zero, 

x/a=(B+C") sin g~"-(B'-A") cosg~"+Cx 

-[H'I(H +S )~]1n tan (7r/4+cJ>"/2) (435) 

y /a= (B'-A") sin g ~" + (B + C ") cos g ~"- Cy 

+[S' /(H2+S2)i] In tan (7r/4+cJ>"/2) (436) 

If g=O, 

x/a=b~"+!,C'n+A"-[H'I(H +S )~] In ta~ (7r/4+cJ>"/2) (440) 

y/ a =b'~"+ !,Cn +C"+[S' /(H2+S2)1] In tan (7r/4+cJ>" /2) (441) 

For the inverse equations, to find cJ>" and~" in terms of x andy: If 

g is not zero, limited to several orbits from the origin, 

sin g~"=[(S'B'- H'B) cos g~"+S' (x/a -Cx-A" cos g~") 

+ H' (y/a +Cy -C" cos g~")]I[S'(B +C") 

+H'(B'-A")] (449) 
PROJECTION FOR THE GENERAL ELLIPTICAL ORBIT 107 

using iteration in which (y/a) is the first trial "A", finding H 1

, 81, A", 

and C" for each trial "A" from the formulas just above. 

ln tan (·n-/4+ct>"/2)=[(H +S )~/H ] [Or -xla +(B +C") sin g"A" 

-(B1

- A") cos g"A"] (450) 

If g =0, without limit, 

"A"=[S1 (x/a -~C n -A")+ HI (y/a- ~Cn -C")]I(S1b+ H 1b1

) (451) 

iterated like (449) above. 

ln tan (·n-/4+ct>"/2)=[(H +S )~/H ] (b"A"+~C n +A"-x/a) 

For cJ> and "A in terms of cf>" and "A" (ignoring Fin (269)): 

E 1=2 arctan{ tan ~("A"-w) [(1-e 1)/(1 +e )]~} 

L=E1-e 1 sinE' 

U=e 2 cos2 ilK 

tan "At= {[1-sin2 cf>"/(1-e 2)] tan "A" cos i 

-(sin cf>" sin i/cos "A") [(1 +Q sin2 "A") (1-sin2 cf>") 

- U sin2 cf>"]~} /[1-sin2 cf>" (1 + U)] 

If cos "A" is not zero, 

tan cf>=(tan "A" cos At -cos i sin "At)IK sin i 

If cos "A" is zero, 

A =K2 cos2 i -sin2 i +e 2 sin2 cf>" 

A 1=K sin i cos i 

sin cJ> = ± { -A ( cos2 cf>"- cos2 i) + 2A 1 [A 1 

±sin cf>" [K (cos2 cj>"-e 2 cos2 i)]i]}i/(A 2+4A12)i 

taking the sign of sin "A" in both places. 

(452) 

(364) 

(334) 

(219) 

(269) 

(271) 

(246a) 

(246c) 

(246d) 

The relationship between c/>", "A". c/>1

, and "A', and the formulas for cf>o 

and "Ao are identical to those described in the previous section for a 

slightly elliptical orbit. 

To evaluate scale error, equations may be differentiated as they were 

for the ellipsoid (equations (302) through (327)), but it is simpler to 

evaluate scale factors as was done for the published paper (Snyder, 

1978b ). To do this, a quadrilateral 0.01 o on a side is assumed: two sides 

0.01 o of cf>" and two sides 0.01 o of "A", located at any point for which the 

scale is to be determined. For example, the quadrilateral may extend 

from a "A" of 45° to a "A" of 45.01 o, and from a cf>" of 1 o to a cf>" of 1.01 o. 

For the scale factor at constant "A", values of x, y, cf> and "A are calculated 

for ("A", c/>")=(45°, 1 °) and (45°, 1.01 °). The numbers are inserted into 

equation (306), using the differences in place of the ~'s, and the arithmetical average of the two cj>'s for cf>. 
108 SPACE OBLIQUE MERCATOR PROJECTION 

Likewise, for the scale factor at constant</>", values of x, y, </>and X are 

calculated for (X", </>")=(45°, 1 °) and (45.01 °, 1 °), and inserted into 

(306). This does not provide w, the maximum angular deformation, but 

the diagonal scale factors may also be evaluated for (X", 4>"}=(45°, 1 °) 

and (45.01 °, 1.01 °), as well as for (X", </>")=(45.01 °, 1 °) and (45°, 1.01 °), 

to give distortion information which the sides of the quadrilateral do 

not. This procedure may be followed for elliptical or circular orbits, 

and for the ellipsoid or sphere, and was used to verify the above 

formulas for the elliptical orbit. 

REFERENCES CITED 

Colvocoresses, A.P., 1974, Space Oblique Mercator: Photogrammetric Engineering, v. 40, 

no. 8, p. 921-926. 

--1975, Evaluation of the cartographic application of ERTS-1 imagery: American 

Cartographer, v. 2, no. 1, p. 5-18. 

Hotine, Martin, 1946, The orthomorphic projection of the spheroitt: Empire Survey 

Review, v. VIII, no. 62, p. 300-311. 

--1947, The orthomorphic projection of the spheroid: Empire Survey Review, v. IX, 

no.63,p.25-35;v.IX,no.64,p.52-70;v.IX,no.65,112-123;v.IX,no.66,p.157-166. 

Junkins, J. L., and Turner, J.D., 1977, Formulation of a Space Oblique Mercator map 

projection: Final rept. UV A/525023/ESS77 /105, Submitted to Computer Science 

Laboratory, U.S. Army Engineer Topographic Laboratories. 

Maling, D. H., 1973, Coordinate systems and map projections: London, George Philip and 

Son, Ltd., p. 58-67. Specific formulas are (4.21), corrected for sign, (4.33), corrected 

for function, (4.51), and (4.59). 

Moulton, F. R., 1914, An introduction to celestial mechanics, 2d ed.: New York, MacMillan, p. 159, 165. 

Short, N. M., Lowman, P. D., Jr., Freden, S.C., and Finch, W. A., Jr., 1976, Mission to 

Earth: Landsat views the world: National Aeronautics and Space Administration, 

p. 437. 

Snyder, J. P., 1978a, Equidistant conic map projections: Assoc. Am. Geog. Annals, v. 68, 

no. 3, p. 373-383. 

--1978b, The Space Oblique Mercator projection: Photogrammetric Engineering and 

Remote Sensing, v. 44, no. 5, p. 585-596. 

Thomas, P. D., 1952, Conformal projections in geodesy and cartography: U.S. Coast and 

Geodetic Surv. Spec. Pub. 251. 

-k U.S. GOVERNMENT PRINTING OFFICE: 1981-341-614/223