Map Distortion

It is not possible to develop a map projection that does not lead to some disortion of the patterns that appear on the surface of the earth. The extent to which different projections distort the original images/distortion can be studied with Maple. The fundamental idea behind an invesigation of map distortion is to plot a simple shape (X's and O's, for example) and see how that shape changes when plotted in a different coordinate system.

A Simple, but Incorrect, Way to Look at Map Distortion

The simplest way to look at map distortion is plot simple shapes in rectangular coordinates. Here we plot a number of circles (and semi-circles) on the centers of a 30 degree equirectangular graticule.

> p1:=coordplot(Equirectangular,scaling=constrained):

> c1:=seq(seq(plottools[circle]([k*30,j*30], 10, color=red),k=-5..5),j=-2..2):

> s1:=seq(plottools[arc]([-180,j*30], 10, Pi/2..-Pi/2,color=red),j=-2..2):
s2:=seq(plottools[arc]([180,j*30], 10, 3*Pi/2..Pi/2,color=red),j=-2..2):

> ti:=plots[display]({c1,p1,s1,s2}): ti;

To see how this image changes we simply replot using a different coordinate system. Some examples follow.

> changecoords(ti,Sinusoidal,view=[-180..180,-90..90]);

> changecoords(ti,Robinson);

> changecoords(ti,`Eckert I`);

> changecoords(ti,Cassini);

> removelines(%);

The above method of looking at map distortion is not the best approach because the original image (the equirectangular projection) already distorts the original outlines on the Earth's surface.

Tissot Method

A better approach, developed by Tissot, is to project the image of circles drawn on the surface of a sphere. In what follows we show (admittedly without much explanation) one way to draw circles on the surface of a sphere.

> with(linalg):

Warning, new definition for norm

Warning, new definition for trace

> T[1] := matrix(3,3,[[sin(theta)*cos(phi), sin(theta)*sin(phi), cos(theta)], [cos(theta)*cos(phi), cos(theta)*sin(phi), -sin(theta)],
[-sin(phi), cos(phi),0]]):

> T_1[1] := matrix(3,3,[[sin(theta)*cos(phi), cos(theta)*cos(phi), -sin(phi)], [sin(theta)*sin(phi), cos(theta)*sin(phi), cos(phi)],
[cos(theta), -sin(theta),0]]):

> Tz := matrix(3,3,[[cos(lambda), -sin(lambda), 0],
[sin(lambda), cos(lambda), 0],[0, 0,1]]);

> Tz2 := matrix(3,3,[[cos(phi), -sin(phi), 0],
[sin(phi), cos(phi), 0],[0, 0,1]]);

> Ty := matrix(3,3,[[cos(theta), 0, -sin(theta)],
[0,1,0], [sin(theta),0, cos(theta)]]);

> Tx := matrix(3,3,[[1,0,0],[0,cos(phi), -sin(phi)],
[0,sin(phi), cos(phi)]]);

> mm := vector(3,[0,r,R]);

> temp := R = sqrt(1 - r^2):temp;

> tt1 := linalg[multiply](Tz,mm): convert(tt1, list);

First, the sphere.

> ss := plottools[sphere]([0,0,0], 0.99):

then a circle around the north pole

> pp := plots[spacecurve]( subs(temp,r = 0.1, convert(tt1, list)) , lambda = -Pi..Pi, color = black, thickness = 3, axes = boxed):

> display({pp,ss},scaling=constrained);

> tt2 := multiply(Ty,tt1): convert(tt2, list);

> tt3 := multiply(Tz2,tt2): convert(tt3, list);

Here we draw two more circles on the sphere.

> with(plots,spacecurve):

> pp := spacecurve( subs(temp,r = 0.1, theta = Pi/4, phi = Pi/4, {convert(tt1, list), convert(tt2, list), convert(tt3, list)}) , lambda = -Pi..Pi, color = black, thickness = 3, axes = boxed,orientation=[-140,60]):
display(pp,ss);

Finally, we show how to draw circles all over the sphere.

> dd := evalf(seq(op({seq(subs(temp,r = 0.2, theta = 2*Pi/12*ii, phi = Pi/6*jj, convert(tt3, list)), ii = [2,3,4])}), jj = 1..12)): dd:

> ee := evalf(seq(op({seq(subs(temp,r = 0.2, theta = 2*Pi/12*ii, phi = Pi/6*jj*2, convert(tt3, list)), ii = [1,5])}), jj = 1..6)): ee:

> gg := evalf(seq(op({seq(subs(temp,r = 0.2, theta = 2*Pi/12*ii, phi = Pi/6*jj*2, convert(tt3, list)), ii = [0,6])}), jj = 1)): gg:

> pp := spacecurve({dd,ee,gg}, lambda = -Pi..Pi, color = black, thickness = 2, axes = boxed):display(pp,ss);

We now wish to project the surface of the above sphere using one of the map coordinate systems discussed elsewhere. We have automated the creation of distortion images/distortion in a procedure included in the Maps package. The procedure must be accessed using the readlib command.

> readlib(`Maps/distortion`):

It is now a trivial exercise to create distortion images/distortion of any other projection. In what follows we show some representative examples. Note that in an equal area projection the area of all complete circles is the same, although the shape suffers somewhat. The degree of distortion introduced by the Mercator projection is evident from the illustrations below.

> mp:=[`Cylindrical Equal Area`, `Eckert IV`, `Lambert Conical Equal Area`, Mercator, Mollweide, Sinusoidal, `Transverse Mercator`, `interrupted Boggs`];

> for m in mp do
`Maps/distortion`(m)
od;

>