The Mathematics of Map Projections

In projecting a picture of the world onto a planar map, there is one crucial geometric problem; the (Gauss) curvature of the R-sphere is , while that of the plane is 0. Gauss's beautiful Theorem Egregium says that two isometric surfaces have the same Gauss curvatures. Therefore, the sphere and the plane are not isometric. In other words, there is no way to map the sphere onto the plane in a (locally) one-to-one fashion so that the notion of distance on each of the surfaces corresponds. In differential geometry, distance on a surface is determined by certain tangent vectors to the surface and associated dot products. Suppose a surface in is parametrized by . The tangent vectors obtained by taking partial derivatives with respect to each of the parameters in each coordinate are and and we define the metric coefficients to be , and . If is any tangent vector, then it always has the form since the vectors and form a basis for each tangent plane. Therefore, the (square of the) length of is given by and this type of expression is often referred to as the metric of the surface. Below, we give a simple procedure to calculate the metric coefficients of a surface parametrized by and (since we will be dealing with the sphere). A small bit of area on a surface is approximated by the parallelogram spanned by the tangent vectors and . Of course, this area is given by the length of and it is a consequence of the Lagrange identity for vectors that . The procedure called area_element below calculates this quantity for a surface. The area element is exactly the quantity which is the integrand of the double integral calculating the surface area of a surface. Finally, we include a procedure to calculate the Gauss curvature of a surface directly from the metric coefficients in the case where . These procedures will allow us to analyze maps.

> with(linalg):

Warning, new definition for norm

Warning, new definition for trace

> dp := proc(X,Y)
local i;
sum(X[i]*Y[i],i=1..3)
end:

> EFG := proc(X)

> local E,F,G,Xtheta,Xphi;

> Xtheta := [diff(X[1],theta),diff(X[2],theta),diff(X[3],theta)];

> Xphi := [diff(X[1],phi),diff(X[2],phi),diff(X[3],phi)];

> E := dp(Xtheta,Xtheta);

> F := dp(Xtheta,Xphi);

> G := dp(Xphi,Xphi);

> map(simplify,[E,F,G]);

> end:

> area_element:=proc(X)
simplify(sqrt(EFG(X)[1]*EFG(X)[3]-EFG(X)[2]^2),symbolic);
end:

> GaussCurvature:=proc(X)
local E,G,sq,Ephi,Gtheta,part1,part2;
E:=EFG(X)[1]; G:=EFG(X)[3];
sq:=sqrt(E*G);
Ephi:=diff(E,phi);
Gtheta:=diff(G,theta);
part1:=diff(Ephi/sq,phi);
part2:=diff(Gtheta/sq,theta);
simplify(-1/(2*sq)*(part1+part2),symbolic);
end:

Let's look at the (unit) sphere and calculate its metric coefficients, area element, surface area and Gauss curvature. We will compare these to other surfaces below. First, let's parametrize the sphere in terms of longitude and latitude .

> sphere:=[cos(theta)*cos(phi),sin(theta)*cos(phi),sin(phi)];

The metric coefficients are

> EFG(sphere);

and the area element is

> area_element(sphere);

with total surface area as it should be.

> int(int(area_element(sphere),theta=0..2*Pi),phi=-Pi/2..Pi/2);

The Gauss curvature (note that ) is then given by

> GaussCurvature(sphere);

The first type of projection of the sphere onto a plane that mathematicians think of is stereographic projection . This projection forms a line through the North Pole of the sphere (i.e. ) and a point on the sphere and projects to the point in the plane which is intersected by the line. We can use Maple to get a formula for this procedure. Start by parametrizing the line as

> stereoline:=evalm([0,0,1]+scalarmul(sphere-[0,0,1],t));

This line hits the plane when the third coordinate becomes zero. To find the t which makes this so, we use Maple's solve command.

> T:=solve(stereoline[3]=0,t);

Finally, we substitute the T we have found into the line and obtain a formula describing where the sphere point maps.

> stereoprojection:=subs(t=T,[stereoline[1],stereoline[2],0]);

Since the North Pole itself lies outside the domain of stereographic projection, we see that the formula makes sense everywhere. Hence, this very formula provides a parametrization (albeit a very unusual one!) for the plane. We can therefore calculate this parametrization's metric coefficients.

> EFG(stereoprojection);

We can see that these coefficients are related to the sphere's metric coefficients by using some Maple algebra. Putting the minus signs in the denominator and numerator of the first and third coordinates respectively, we can simplify both to

> factor(cos(phi)^2/(2-2*sin(phi)-(1-sin(phi)^2)));

> factor((2*sin(phi)+2-(1-sin(phi)^2))/(1-sin(phi)^2)^2);

This bit of algebra shows that the metric coefficients of the plane (with the stereographic parametrization) are just multiples by of those of the sphere. Such a mapping is called a conformal map . These mappings are special because they are the curved analogues of similitudes in plane geometry. That is, they stretch tangent planes by the same amount in all directions. Just as similitudes preserve angles, so do conformal maps (and this is their key quality for map projections). So, we see that stereographic projection preserves the angles between things. But what about area: is it preserved? Unfortunately, this is definitely not the case. The area element is

> area_element(stereoprojection);

Below, we give a procedure for stereographically projecting curves on the sphere to curves in the plane, and vice versa! From this it can be seen that stereographic projection is probably not a preferred way to create a map.

In the first part of this article, we used a cylindrical projection (often called Lambert's projection). This consisted of projecting outward from the z-axis horizontally from the sphere to a cylinder tangent to the sphere along the equator. The projecting line has a form similar to the stereographic line above. Note that is the third coordinate of the point on the sphere that is being projected. For this projection, we require the value of the parameter t which puts the first two coordinates on a unit circle. Again the solve command is useful.

> cylline:=evalm([0,0,sin(phi)]+scalarmul(sphere-[0,0,sin(phi)],t));

> S:=solve(cylline[1]^2+cylline[2]^2=1,t);

The cylindrical projection is then given by substituting the value S into the line projection.

> cylprojection:=subs(t=1/cos(phi),[cylline[1],cylline[2],cylline[3]]);

This projection is very simple, but what are its properties? First, we see that the metric coefficients do not bear an easy relationship to those of the sphere (except that they are a scrambled version). The cylindrical projection is not conformal.

> EFG(cylprojection);

However, when we compute the area element of the cylinder with this parametrization, we find that it is exactly the same as that of the sphere!

> area_element(cylprojection);

> int(int(area_element(cylprojection),theta=0..2*Pi),phi=-Pi/2..Pi/2);

Therefore, the cylindrical projection is area preserving and this is what recommends it as a map projection. We note that, however,

> GaussCurvature(cylprojection);

and this means that essential information must be lost in projecting from the sphere to the cylinder (as we see from the metric coefficients).

Put simply, there is no way to make a perfect map! This is the content of the Theorem Egregium in this situation. So, map makers must always decide what aspect of the world they want to preserve. For example, they can preserve angles by conformal mappings and they can preserve areas by area preserving ones --- but it is a fact (an exercise for the reader) that they can never do both simultaneously, for this would entail an isometry, contradicting the Theorem Egregium . There are yet other properties which recommend types of maps. In particular, it is often desired to have longitudes as vertical lines on a map and latitudes as horizontal lines. The Mercator projection has this property as well as being conformal. To see that Mercator is conformal, we will show that it is the composition of stereographic projection with the complex logarithm. Of course, we have seen that stereographic projection is conformal and it is a basic fact that complex analytic functions, such as the logarithm, are conformal (where there derivatives are non-zero). Hence, the composition is conformal as well. Recall that, for complex, the log is defined by where we note that and (essentially). With this in mind, apply the complex log to stereographic projection.

> expand(evalc(ln(stereoprojection[1]+I*stereoprojection[2])));

This expression simplifies to and the second term to (but it is not easy to get Maple to do this in a straightforward way). Thus, we end up with a parametrization of the plane and this is the Mercator projection. It is easy to check that Mercator is not area preserving (as it cannot be since it is conformal). Also, the parametrization says that a latitude = fixed goes to a horizontal line (since the second coordinate is fixed) and a longitude = fixed goes to a vertical line. This has the following consequence which proved vital for navigation in years past: if a curve on the sphere keeps a constant direction (as can be measured by the angle with meridians say), then the curve goes to a straight line under the Mercator projection. To see this, note that Mercator is conformal, so the image of the curve makes a constant angle with the images/mapmath of meridians - vertical lines. Since these lines are all parallel, there is only one choice for the image curve: it must be a straight line! Conversely, under the inverse of the Mercator projection, straight lines go to curves having constant direction, the so-called loxodromes . Below we give a procedure for taking curves on the sphere into planar curves via stereographic projection. The last two arguments allow the user to restrict the viewing area. For this procedure and the Mercator procedure below, we have reverted to using Cartesian coordinates since the curves are inputted parametrically in three Cartesian coordinates. In Cartesian coordinates, stereographic projection is given by and the Mercator projection is given by where is a point on the sphere. The second coordinate of the Mercator projection uses the fact that .

> Stereo:=proc(alpha,a,b,view1,view2)
local sqleng,curv,sph,projcurv;
sqleng:=nrm(alpha)^2;
curv:=[alpha[1]/(1-alpha[3]),alpha[2]/(1-alpha[3]),t=a..b];
projcurv:=plot(curv,color=black,thickness=2,numpoints=200);
plots[display]({projcurv},scaling=constrained,view=[-view1..view1,-view2..view2]);
end:

Here is an interesting curve to try. Take the intersection of a cylinder with the unit sphere. As the reader can show, a parametrization for the intersection curve (known as Viviani's curve) is and shown below:

> sphereplot:=plot3d(sphere,theta=0..2*Pi,phi=-Pi/2..Pi/2,color=cyan):
cylplot:=plot3d([cos(theta)/2+1/2,sin(theta)/2,v],theta=0..2*Pi,v=-1..1,color=orange):
viviani:=plots[spacecurve]([cos(t)^2,sin(t)*cos(t),sin(t)],t=0..2*Pi,color=black,thickness=2):
plots[display]({sphereplot,cylplot,viviani},scaling=constrained,style=wireframe,orientation=[0,70]);

Stereographic projection then gives

> Stereo([cos(t)^2,sin(t)*cos(t),sin(t)],-Pi,Pi,2,4);

We can also check that latitude circles map to concentric circles about the origin.

> Stereo([cos(t)*cos(Pi/4),sin(t)*cos(Pi/4),sin(Pi/4)],0,2*Pi,4,4);

The inverse of stereographic projection is found by parametrizing a line from any point on the plane to the North Pole of the unit sphere and then determining the value of the parameter which gives a point on the sphere by the condition that the sum of the squares of the coordinates be 1. The reader should carry out this procedure and derive the formula below (i.e. curv).

> InvStereo:=proc(alpha,a,b)
local sqleng,curv,sph,plcurv,sphcurv;
sqleng:=alpha[1]^2+alpha[2]^2;
curv:=[2*alpha[1]/(1+sqleng),2*alpha[2]/(1+sqleng),(sqleng-1)/(1+sqleng)];
sphcurv:=plots[spacecurve](curv,t=a..b,color=black,thickness=2,numpoints=200);
sph:=plot3d([cos(theta)*cos(phi),sin(theta)*cos(phi),sin(phi)],theta=0..2*Pi,phi=-Pi/2..Pi/2):
plots[display]({sph,sphcurv},scaling=constrained,shading=zhue,style=wireframe);
end:

We now give examples of the inverse of stereographic projection applied to a line and to the standard astroid.

> InvStereo([t+1,2*t+3,0],-8,8);

> InvStereo([cos(t)^3,sin(t)^3,0],-8,8);

Of course, now we can write a procedure to Mercator project curves and apply it to curves such as Viviani's curve:

> Merc:=proc(alpha,a,b,view1,view2)
local curv,projcurv;
curv:=[arctan(alpha[2]/alpha[1]),1/2*ln((1+alpha[3])/(1-alpha[3])),t=a..b];
projcurv:=plot(curv,color=black,thickness=2);
plots[display]({projcurv},scaling=constrained,view=[-view1..view1,-view2..view2]);
end:

> Merc([cos(t)^2,sin(t)*cos(t),sin(t)],-Pi,Pi,10,10);

as well as to meridians and parallels (where we can check our results):

> Merc([cos(t)*cos(Pi/4),sin(t)*cos(Pi/4),sin(Pi/4)],0,2*Pi,2,2);

> Merc([cos(t)*cos(Pi/4),cos(t)*sin(Pi/4),sin(t)],-Pi/2,Pi/2,2,2);

It is more fun to create curves on the sphere by using the inverse of Mercator projection. Since Mercator projection is accomplished by first doing stereographic projection and then taking the complex logarithm, to achieve its inverse we must undo these operations in the reverse order.That is, we first do the inverse of the complex log, the complex exponential , to the input curve. This is called invlogalpha in the procedure which follows (but to have the correspondences latitude=horizontal, longitude=vertical, we have interchanged and ). We then apply the inverse of stereographic projection to this invlogalpha to obtain the inverse of Mercator projection.

> InvMerc:=proc(alpha,a,b)
local sqleng,curv,sph,plcurv,sphcurv,invlogalpha;
invlogalpha:=[exp(alpha[2])*cos(alpha[1]),exp(alpha[2])*sin(alpha[1]),0];
sqleng:=(exp(alpha[2])*cos(alpha[1]))^2+(exp(alpha[2])*sin(alpha[1]))^2;
curv:=[2*invlogalpha[1]/(1+sqleng),2*invlogalpha[2]/(1+sqleng),(sqleng-1)/(1+sqleng)];
sphcurv:=plots[spacecurve](curv,t=a..b,color=black,thickness=2,numpoints=200);
sph:=plot3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],u=0..2*Pi,v=-Pi/2..Pi/2):
plots[display]({sph,sphcurv},scaling=constrained,shading=zhue,style=wireframe);
end:

Now we can see a loxodrome by taking the image of a straight line under inverse Mercator.

> InvMerc([3*t+3,t+1,0],-5,5);

We can also see that vertical lines go to longitudes and horizontal lines to latitude circles.

> InvMerc([3,t,0],-5,5); InvMerc([t,0.5,0],-5,5);

Parabolas may be pulled back to the sphere as well as other familiar curves.

> InvMerc([t,t^2,0],-5,5); InvMerc([t^2,t,0],-5,5);

> InvMerc([t+1,sin(t),0],-5,5);

> InvMerc([sin(t),cos(t),0],-5,5);

> InvMerc([cos(t)^3,sin(t)^3,0],-5,5);

>