Conic Projections

In conic projections the surface of projection is a cone. The cone can be tangent to the sphere as shown below

> with(plottools):

> icecream := cone([0,0,-1],0.7,color=red), sphere([0,0,0.1],0.6,color=blue):

> plots[display](icecream, scaling=constrained, style=wireframe,orientation=[60,-100]);

Or it can be secant to the cone as shown here.

> icecream := cone([0,0,-1],0.6,1.6,color=red), sphere([0,0,0.1],0.6,color=blue):

> plots[display](icecream, scaling=constrained, style=wireframe,orientation=[60,-100]);

As might be expected there are many different types of conical projection. In the material that follows we show some of them.

Equidistant conic

The equidistant conic projection is given by the following equations

> `Equidistant conic/equations`:
Phi[0]=(Phi[1]+Phi[2])/2,
n=(cos(Phi[1])-cos(Phi[2]))/(Phi[2]-Phi[1]),
rho=r*((Phi[2]*cos(Phi[1])-Phi[1]*cos(Phi[2]))/(cos(Phi[1])-cos(Phi[2]))-phi),
rho[0]=r*((Phi[2]*cos(Phi[1])-Phi[1]*cos(Phi[2]))/(cos(Phi[1])-cos(Phi[2]))-Phi[0]),
theta=n*lambda,
x=rho*sin(theta),
y=rho[0]-rho*cos(theta);

where and are the latitudes at which the cone cuts the sphere of the earth.

The equidistant conic projection can be created using the following command

> mapcoords(`Equidistant conic`,
input = [Lambda,phi],
coords = [shift= '`if`(type(Shift,name),0,Shift)',
lambda = 'readlib(`Maps/shiftf`)'(Lambda-shift*Pi/180,Pi),
Phi_1=convert(P_1*degrees,radians),
Phi_2=convert(P_2*degrees,radians),
Phi_0=(Phi_1+Phi_2)/2,
n=(cos(Phi_1)-cos(Phi_2))/(Phi_2-Phi_1),
rho=r*((Phi_2*cos(Phi_1)-Phi_1*cos(Phi_2))/(cos(Phi_1)-cos(Phi_2))-phi),
rho_0=r*((Phi_2*cos(Phi_1)-Phi_1*cos(Phi_2))/(cos(Phi_1)-cos(Phi_2))-Phi_0),
theta=n*lambda,
x=rho*sin(theta),
y=rho_0-rho*cos(theta),
[x,y]],
params = [r,Shift,P_1,P_2],
view = [-180..180,-90..90,13,7,-180..180,-150..150]);

The coordinate system can be shown in the usual way using coordplot

> p0:=removelines(coordplot(`Equidistant conic`(1,-90,60,20),scaling=constrained)): p0;

Now, we superimpose the outline of the land/water bodies

> p1:=removelines(changecoords(world[ng,50],`Equidistant conic`(1,-90,60,20))):

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

Perspective Conic

The perspective conic is given by the following command:

> mapcoords(`Perspective conic`,
input = [Lambda,phi],
coords = [shift= '`if`(type(Shift,name),0,Shift)',
lambda = 'readlib(`Maps/shiftf`)'(Lambda-shift*Pi/180,Pi),
Phi_1=convert(P_1*degrees,radians),
Phi_2=convert(P_2*degrees,radians),
Phi_0=(Phi_1+Phi_2)/2,
n=sin(Phi_0),
rho=r*cos((Phi_2-Phi_1)/2) * (cot(Phi_0)-tan(phi-Phi_0)),
rho_0=r*cos((Phi_2-Phi_1)/2),
theta=n*lambda,
x='`if`'((phi<=(Phi_0-Pi/2) and Phi_0>0) or(phi>=(Phi_0+Pi/2) and Phi_0<0), FAIL,rho*sin(theta)),
y='`if`'((phi<=(Phi_0-Pi/2) and Phi_0>0) or(phi>=(Phi_0+Pi/2) and Phi_0<0), FAIL,rho_0-rho*cos(theta)),
[x,y]],
params = [r,Shift,P_1,P_2],
view = [-180..180,-90..90,13,7,-180..180,-150..150]);

> p0:=changecoords(world[50],`Perspective conic`(1,-90,60,30)):

> p1:=changecoords(world[ng,50],`Perspective conic`(1,-90,60,30)): p1;

> plots[display]({p0,p1},view=[-100..100,-100..80]);

Albers Equal Area Conic Projection

> mapcoords(`Albers equal area conic`,
input = [Lambda,phi],
coords = [shift= '`if`(type(Shift,name),0,Shift)',
lambda = 'readlib(`Maps/shiftf`)'(Lambda-shift*Pi/180,Pi),
Phi_1=convert(P_1*degrees,radians),
Phi_2=convert(P_2*degrees,radians),
Phi_0=(Phi_1+Phi_2)/2,
n=(1+sin(Phi_1))/2,
rho=r*sqrt(2*(1-sin(phi))/n),
rho_0=r*sqrt(2*(1-sin(Phi_0))/n),
theta=n*lambda,
x=rho*sin(theta),
y=rho_0-rho*cos(theta),
[x,y]],
params = [r,Shift,P_1,P_2],
view = [-180..180,-90..90,13,7,-180..180,-150..150]);

> p0:=coordplot(`Albers equal area conic`(1,-90,60,30),scaling=constrained): p0;

> p1:=changecoords(world[ng,50],`Albers equal area conic`(1,-90,60,30)): p1;

> plots[display]({p0,p1},view=[-100..100,-100..80]);

Polyconic Projection

> mapcoords(Polyconic,
input = [Lambda,phi],
coords = [shift= '`if`(type(Shift,name),0,Shift)',
lambda = 'readlib(`Maps/shiftf`)'(Lambda-shift*Pi/180,Pi),
x='`if`'(phi=0,r*lambda,r*cot(phi)*sin(lambda*sin(phi))),
y='`if`'(phi=0,0,r*(phi+cot(phi)*(1-cos(lambda*sin(phi))))),
[x,y]],
params = [r,Shift],
view = [-180..180,-90..90,13,7,-180..180,-150..150]);

> p0:=coordplot(Polyconic(1,0),scaling=constrained):p0;

> p1:=changecoords(world[ng,50],Polyconic(1,-90)): p1;

> p1:=removelines(p1):

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

Below we use the polyconic projection for a map of North America

> plots[display]({p0,p1},view=[-50..50,0..90]);

Rectangular Polyconic Projection

> mapcoords(`Rectangular polyconic`,
input = [Lambda,phi],
coords = [n= '`if`(type(_n,name),0,_n)',
lambda = 'readlib(`Maps/shiftf`)'(Lambda-n*Pi/180,Pi),
theta = 2*arctan(lambda*sin(phi)/2),
x='`if`'(phi=0,r*lambda,r*cot(phi)*sin(theta)),
y='`if`'(phi=0,0,r*phi+r*cot(phi)*(1-cos(theta))),
[x,y]],
params = [r,_n],
view = [-180..180,-90..90,13,7,-180..180,-150..150]);

> p0:=coordplot(`Rectangular polyconic`(1,0),scaling=constrained): p0;

Now combine the coordinate system with the world:

> p2:=changecoords(world[ng,50],`Rectangular polyconic`(1,90)):
p2:=removelines(p2):

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

A projection of Asia is shown below.

> plots[display]({p2,p0},view=[-50..50,0..90]);

Bonne Projection

> mapcoords(Bonne,
input = [Lambda,phi],
coords = [shift= '`if`(type(Shift,name),0,Shift)',
lambda = 'readlib(`Maps/shiftf`)'(Lambda-shift*Pi/180,Pi),
A=cot(Phi)+Phi-phi,
B=lambda*cos(phi)/A,
x=r*A*sin(B),
y=r*(cot(Phi)-A*cos(B)),
[x,y]],
params = [r,Shift,Phi],
view = [-180..180,-90..90,13,7,-180..180,-150..150]);

> removelines(changecoords(world[50],Bonne(1,-90)));

Werner Projection

The Werner projection is a special case of the Bonne

> mapcoords(Werner,
input = [Lambda,phi],
coords = [shift= '`if`(type(Shift,name),0,Shift)',
lambda = 'readlib(`Maps/shiftf`)'(Lambda-shift*Pi/180,Pi),
Phi=Pi/2,
A=cot(Phi)+Phi-phi,
B1=lambda*cos(phi),
x='`if`'(phi=evalf(Pi/2),0,r*A*sin(B1/A)),
y='`if`'(phi=evalf(Pi/2),0,r*(cot(Phi)-A*cos(B1/A))),
[x,y]],
params = [r,Shift],
view = [-180..180,-90..90,13,7,-180..180,-150..150]);

Make the coordinate system

> coordplot(Werner,view=[-150..150,-180..80]):

> wp0:=removelines(%): wp0;

Reflect the above image around the center line and display with the aboveto make a more complete coordinate system

> wp1:=plottools [reflect](wp0,[[0,-180],[0,180]]):

> wp2:=plots[display]({wp0,wp1}): wp2;

Now add the world.

> wp:=removelines(changecoords(world[ng,50],Werner(1,-90))):

> plots[display]({wp,wp2});