Miscellaneous Projections - Peirce Quincuncial
Notice how the south pole appears in all four corners of a square. We begin our discussion of this interesting projection by defining the following integral:
> `Pierce/integral` := unapply(int(1/sqrt(1-1/2*sin(t)^2),t=0..v),v);
Pierce's projection, also known as the quincuncial has x and y coordinates defined by
> x = `Pierce/integral`(v);
> y=`Pierce/integral`(u);
where
> sin(u) = (sqrt(2)* cos((a+b)/2));
> sin(v) = (sqrt(2)* sin((a-b)/2));
with parameters a and b given by
> cos(a) = (cos(phi)*cos(Pi/4 + lambda));
> cos(b) = (cos(phi)*cos(Pi/4 - lambda));
The evaluation of these coordinates appears to be somewhat involved. Fortunately Maple can provide some help with the integration step. We may evaluate the integral numerically by enclosing the function inside a call to evalf . Let us evaluate the function for some numerical values of the argument.
> for s in [1/2,1,3/2,2,3] do `Pierce/integral`(s) od;
From these results we see that the result of the integral is a combination of complete and incomplete elliptic functions depending on the magnitude of the argument. Further experimentation will show that the switch from the single elliptic function to the combination of two occurs at . Thus, we can combine the above results into a single procedure:
> PierceElliptic := proc(xx) option remember; if evalf(abs(xx) > Pi/2) then evalf(signum(xx)*(2*EllipticK(1/2*2^(1/2))-EllipticF(sin(abs(xx)),1/2*2^(1/2)))) else evalf(EllipticF(sin(abs(xx)),1/2*2^(1/2))); fi; end;
We now proceed to define the coordinate system as a series of equations as shown below. Note the use of the signum function to assign the appropriate sign to the coordinates and the enclosure of the entire algorithm within quote marks to delay evaluation of the functions.
> `Pierce/equations`:= 'a = arccos(cos(phi)*cos(Pi/4 + lambda)), b = arccos(cos(phi)*cos(Pi/4 - lambda)), u = arcsin(sqrt(2)* cos((a+b)/2)), v = arcsin(sqrt(2)* sin((a-b)/2)), x = 'signum'(Pi/2-abs(lambda))*'PierceElliptic'(u), y = 'signum'(lambda)*'PierceElliptic'(v)';
> mapcoords(Pierce, input = [lambda,phi], coords = [`Pierce/equations`, [x,y]], params = [r], view = [-180..180,-90..90,13,7,-90..90,-180..180]):
We will not use coordplot to isplay the coordinate system; rather, we shall make a grid:
> tt1 := graticule([seq(ii*10, ii = -18..18)],[-90,-60,-40,-20,-5,5,20,40,60,90]):
and convert it to the Pierce coordinate system.
> pp3 := display(`Maps/changecoords`(tt1,Pierce)): pp3;
Now we can test these new coordinates on the whole world.
> changecoords(world[50],Pierce);
We see here a problem encountered with some other projections, the northern and southern hemispheres overlap. If we redefine the coordinate system as follows:
> mapcoords(`Pierce northern`, input = [lambda,phi], coords = [temp = '`if`(phi < 0, RETURN([FAIL,FAIL]), 0)', `Pierce/equations`,[x,y]], params = [r], view = [-180..180,-90..90,13,7,-90..90,-180..180]):
we may obtain an attractive view of just that one hemisphere.
> northhemi:=changecoords(world[50],`Pierce northern`): northhemi;
Simply repeating the above projection for the southern hemisphere will give an image analogous to that shown above centered on the south pole.
> mapcoords(`Pierce southern`, input = [lambda,phi], coords = [temp = '`if`(phi > 0, RETURN([FAIL,FAIL]), 0)', `Pierce/equations`,[x,y]], params = [r], view = [-180..180,-90..90,13,7,-90..90,-180..180]):
> southhemi:=changecoords(world[50],`Pierce southern`): southhemi;
The next question is how to attach the southern hemisphere to the northern. Pierce's solution was to cut the southern hemisphere into four along the diagonals shown below.
> diagonals := plot([x,-x],x=-105..105,color=red,thickness=2):
> plots[display]({southhemi,diagonals});
Each triangle is folded outwards along it's outer edge using the following procedure that tests into which quadrant the point just computed falls and then modifies the coordinates appropriately.
> `Pierce/foldout` := proc(X) local x, y, border; border := 1.854046393; x := X[1]; y := X[2]; if not type(x,numeric) or not type(y,numeric) then RETURN([FAIL,FAIL]) fi; if y > x and y > -x then RETURN([x,2*border-y]) elif y > x and y < -x then RETURN([-2*border-x,y]) elif y < x and y > -x then RETURN([2*border-x,y]) elif y < x and y < -x then RETURN([x,-2*border-y]) else RETURN([FAIL,FAIL]) fi; end:
> mapcoords(`Pierce southern`, input = [lambda,phi], coords = [temp = '`if`(phi > 0, RETURN([FAIL,FAIL]), 0)', `Pierce/equations`, t2 = '`Pierce/foldout`'([x,y]), [t2[1],t2[2]]], params = [r], view = [-180..180,-90..90,13,7,-90..90,-180..180]):
Warning: coordinates already exists, system redefined.
Warning: default information already exists, redefining.
Now the southern hemisphere appears as shown below.
We can patch this onto the Northen hemisphere to obtain:
> plots[display]({northhemi,southhemi});
Here we see that part of South America cut off from the rest. The solution to this unsatisfying image is to shift the origin 20 degrees east. We can accomplish all of this in one single projection as follows:
> mapcoords(Pierce, input = [Lambda,phi], coords = [lambda = 'readlib(`Maps/shiftf`)'(Lambda-20,180), `if`(phi > 0,`Maps/Pierce northern`(lambda,phi,r), `Maps/Pierce southern`(lambda,phi,r))], params = [r], view = [-180..180,-90..90,13,7,-207..207,-207..207]):
The coordinate system appears as shown below
> piercecoords:=coordplot(Pierce,scaling=constrained): piercecoords;
> piercegrid:=changecoords(tt1,Pierce): piercegrid;
and the entire world:
> pierceworld := removelines(changecoords(world[ng,50],Pierce)): pierceworld;
> plots[display]({pierceworld,piercecoords});
One of the interesting characteristics of the Pierce projection is that it can be tiled indefinitely.
> maxval :=2*106.2283620;
> shiftvalues := [[maxval,maxval],[maxval,-maxval],[-maxval,maxval],[-maxval,-maxval],[0,2*maxval],[0,-2*maxval],[2*maxval,0],[-2*maxval,0]];
> plots[display]({pierceworld,seq(plottools[translate](pierceworld,op(v)),v=shiftvalues)});
>
