Azimuthal Projections

Azimuthal (or zenithal) projections are projections onto a plane that is tangent to some reference point on the globe. If the reference point is one of the poles the projections are polar azimuthal (zenithal) views. If the reference point lies on the equator the projections are termed normal. For all other reference points the projections are oblique . There are five azimuthal projections considered below: equal area, gnomonic, equidistant, stereographic, and orthographic.

A generalized set of equations may be written for the family of azimuthal projections (Snyder, 1993):

> `Azimuthal/x` := rho* (cos(phi)*sin(lambda))/sin(z);

> `Azimuthal/y` := rho* (cos(Phi)*sin(phi)- sin(Phi)*cos(phi)*cos(lambda))/sin(z);

> `Azimuthal/z`:=z=arccos(sin(Phi)*sin(phi) + cos(phi)*cos(Phi)*cos(lambda));

where is the longitude of the point of tangency and z is the great circle distance from the center.

One of the important azimuthal projections is the equal area projection developed by J.H. Lambert, for which projection is defined as

> `Azimuthal/equal area/rho` := rho = 2 * r*sin(z/2);

The gnomonic projection is a perspective view of the globe as seen by an observer at its center. For this projection is defined by

> `Azimuthal/gnomonic/rho` := rho = r*tan(z);

The equidistant projection has parallels that are the same distance apart and for which is defined by

> `Azimuthal/equidistant/rho` := rho = r*z;

The stereographic projection is a perspective view of the globe from a point on the surface that lies on the far side from the reference point. For example, if the reference point is the North Pole then the viewpoint is the South Pole. For this projection is defined by

> `Azimuthal/stereographic/rho` := rho = 2*r*tan(z/2);

The orthographic projection is a perspective view from infinitely far above the surface of the globe. is defined by

> `Azimuthal/orthographic/rho` := rho = r*sin(z);

It is interesting to compare the different projections by plotting these functions

> eqnlist := map(rhs,subs(r=1,[seq(`Azimuthal/`.pname.`/rho`,pname=[`equal area`,gnomonic,equidistant,stereographic,orthographic])]));

> plot(eqnlist,z=0..Pi/2-0.0001,view=[0..Pi/2,0..Pi],color=[red,blue,green,black,yellow],thickness=2);

The red line is the equal area projection, the blue is the gnomonic, the green is for the equidistant, black represents the stereographic and yellow the orthographic. This illustration tells us that the gnomonic projection cannot be used to project very large areas (e.g. entire hemispheres). We may also need to exercise some caution with the stereographic projection as well.

We make the projections as follows. Note how the computations are optimized to eliminate repeated trigonometric calculations involving the same angles.

> readlib(optimize):

> for pname in [`equal area`,gnomonic,equidistant,stereographic,orthographic] do
optcalc:=[optimize([`Azimuthal/z`,`Azimuthal/`.pname.`/rho`,x=`Azimuthal/x`,y=`Azimuthal/y`])];
newcalc:=[op(optcalc[1..nops(optcalc)-2]),map(rhs,optcalc[nops(optcalc)-1..nops(optcalc)])];
mapcoords(`Azimuthal `.pname,
input = [alpha, phi],
coords = [Phi = '`if`(type(_Phi,name), 0.0, _Phi/180*Pi)',
Lambda = '`if`(type(_Lambda,name), 0.0, _Lambda/180*Pi)',
lambda = 'readlib(`Maps/shiftf`)'(alpha - Lambda, Pi),
op(newcalc)],
params = [r, _Phi,_Lambda],
view = [-180..180,-90..90,13,7,-120..120,-120..120]);
od:

The code for the equal area projection is shown below for illustrative purposes.

> print(`Maps/Azimuthal equal area`);

The parameters are, respectively, the latitude and longitude (in degrees) of the point to appear at the centre of the projection. The default value of both angles is zero.

Below we show the equal area projection from a point of view above the North Pole. Note that the equations given above involve a division by zero if we specify the reference latitude as 90 degrees. We can get around this using a traperror command as we have done elsewhere but in this case we avoid the problem by using a latitude that is just slightly less than 90 degrees. The resulting projection is not affected by this simple minded trick. Note, further, that the grid is constructed separately in order to avoid the appearance of some spurious lines that would otherwise appear.

> p1:=changecoords(world[ng,50], `Azimuthal equal area`(1,89.999)):

> p2 := removelines(coordplot(`Azimuthal equal area`(1,89.999),scaling=constrained)):

> plots[display]({p1,p2});

Note how the entire surface of the earth fits into a circle. Viewed from the North Pole as done here we see that the outer "circle" is, in fact, the outline of Antarctica. Below is the equal area projection from the equator centered on North and South America.

> p1:=changecoords(world[50], `Azimuthal equal area`(1,0,-90)):

> plots[display]({p1});

Here is the azimuthal equidistant projection from latitude 0 degrees, longitude 0 degrees. The conversion of the grid leads to the appearance of a few straight lines that should not appear but which have been retained as they don't significantly detract from the appearance of the projection.

> changecoords(world[50], `Azimuthal equidistant`);

The stereographic projection from the same reference point.

> splot:=changecoords(world[50], `Azimuthal stereographic`): splot;

We see here what was hinted at by our earlier comparison of the projections, some points are so far away that portion of interest is greatly compressed. Use of the view option allows us to to display a projection that we find quite appealing.

> plots[display](splot,view=[-200..200,-200..200]);

We now look at the world using the orthographic projection viewed from a point very close to the North pole.

> changecoords(world[50], `Azimuthal orthographic`(1,89.99,0));

A passing glance at the above might suggest a succesful projection. More careful inspection of this image, however, reveals that all is not well. We can see the tip of South America occupying the same space as parts of the USA and Canada; New Zealand is unreasonably close to the Siberian peninsula. We are, in fact, seeing southern hemisphere through the northern. The projection is quite correct from the mathematical point of view; in practice, however, the earth is not transparent. What we need to do is exclude the land/water bodies that are not visible from infinitely far above the reference point. For an orthographic projection of the northern (southern) hemisphere we can do this very simply be redefining our coordinate system so that only points with a positive (negative) latitude are considered. The following, for example, would suffice for a view of the northern hemisphere.

> mapcoords(Az,
input = [alpha, phi],
coords = [q='`if`(phi<0, 0, RETURN([FAIL,FAIL]))',
Phi = '`if`(type(_Phi,name), 0.0, _Phi/180*Pi)',
Lambda = '`if`(type(_Lambda,name), 0.0, _Lambda/180*Pi)',
lambda = 'readlib(`Maps/shiftf`)'(alpha - Lambda, Pi),
`Azimuthal/z`,
`Azimuthal/`.pname.`/rho`,
[r*`Azimuthal/x`, r*`Azimuthal/y`]],
params = [r, _Phi,_Lambda],
view = [-180..180,-90..90,13,7,-120..120,-120..120]);

Simply change the < to > for views of the southern hemisphere. Only slightly more complicated constructions will suffice for equatorial views. Sadly, such simple modifications will not work for oblique projections (from points that are not above the poles or equator) and we need something more rigorous.

All points south of the equator are invisible to an observer stationed far above the North pole. Conversely, all points north of the equator are invisible to an observer located far below the south pole. For an observer stationed high over the Atlantic Ocean just south of the North African bulge (0 degrees latitude, 0 degrees longitude) then all points within +/- 90 degrees east and west are visible. It should be obvious that any point that is more than 90 degrees removed from the reference point directly below the observer cannot be seen by that observer. It is important to realize that the 90 degrees must be measured along the great circle that joins the point in question to the reference point. Thus, in order to check the angular separation between two points on the surface of the reference sphere we need to compute the great circle distance between them. This is a relatively simple problem in spherical or differential geometry with the following result.

> cos(theta) = collect(cos(Lambda)*cos(Phi)*cos(lambda)*cos(phi)+sin(Lambda)*cos(Phi)*cos(phi)*sin(lambda)+sin(Phi)*sin(phi),[cos(Phi),cos(phi)]);

where is the great circle distance (in angular units), are the longitude and latitude of the reference point and are the longitude and latitude of the point of interest. We encapsulate this result in a procedure for repeated use in what follows.

> greatcircleangle:=proc(Lambda,Phi,lambda,phi)
arccos((cos(Lambda)*cos(lambda)+sin(Lambda)*sin(lambda))*cos(phi)*cos(Phi)+sin(Phi)*sin(phi));
end;

The next task is to design a procedure to test whether or not a given point is more or less than 90 degrees ( radians) away from the reference point.

> greatcircletest := proc(A::list,Lambda,Phi)
local gca;
gca := evalf(greatcircleangle(Lambda*Pi/180,Phi*Pi/180,A[1]*Pi/180,A[2]*Pi/180));
if type(A,list(numeric)) and gca > evalf(Pi/2)+0.0000001 then
RETURN([FAIL,FAIL])
else
RETURN(A)
fi
end;

Note that the procedure takes as arguments the coordinates of the point in question in the form of a list of two numbers (this is the form used everywhere else in this material on map projections), and the longitude and latitude (both in degrees) of the reference point. If the point is more than radians from the reference point then the procedure returns [FAIL, FAIL], otherwise it returns the coordinates of the point unchanged. Maple's plot routines will ignore points with FAIL coordinates and lines that include points with such coordinates somewhere in the sequence will be broken (which is exactly the desired behavior here). We now need a procedure that applies greatcircletest to all of the coordinates of our world maps:

> Hidehemisphere := proc (areamap::PLOT,Lambda,Phi)
local newmap, allcurves, numcurves, datalist, numpoints, i, newpoints;
newmap := areamap;
allcurves := select(has,[op(areamap)],CURVES);
numcurves := nops(allcurves);
for i to numcurves do
for datalist in [op(allcurves[i])] do
if type(datalist,list) then
newpoints:=map(greatcircletest,datalist,Lambda,Phi);
newmap := subs(datalist = newpoints,newmap) ;
fi;
od;
od;
RETURN(newmap)
end;

Let us test our procedure by developing projections from a point somewhere near the Black Sea. The first step is to recreate the basic map showing only those points of interest.

> wplot := Hidehemisphere(world[50],45,45):

The simple equirectangular projection of this map is shown below.

> wplot;

We can now display this map in any coordinate system.

> `Maps/changecoords`(wplot,`Azimuthal orthographic`(1,45,45));

> `Maps/changecoords`(wplot,`Azimuthal equidistant`(1,45,45));

Finally, we come to the gnomonic projection:

> gplot:=changecoords(wplot,`Azimuthal gnomonic`(1,45,45)): gplot;

The above shows the very large distortion that results from using the gnomonic projection over too great an area. This is another of those projections where we must restrict the view window.

> plots[display](gplot,view=[-200..200,-200..200]);

Note the excessive distortion as we move further away from the reference point. Let us try again and restrict the view window still more.

> plots[display](gplot,view=[-80..80,-80..80]);

As a second example of application of hiding parts of the world map we present azimuthal projections from above the North Pole.

> nplot := Hidehemisphere(world[50],0,90):

> nplot;

> nplot1:=changecoords(nplot,`Azimuthal stereographic`(1,89.99,0)):

> plots[display](nplot1,view=[-120..120,-120..120]);

> changecoords(nplot,`Azimuthal equal area`(1,89.99,0));

We leave it as an exercise for readers to create south polar projections and normal equatorial projections of eastern and western hemispheres.