# Peirce's quincuncial projection

The Peirce quincuncial projection is the cartographic projection of a sphere onto a square.

In short, I would like to see it implemented in Mathematica.

Here is my code:

Clear[InversePeirceQuincunicalMapping];
InversePeirceQuincunicalMapping[x_Real, y_Real] :=
Module[{cn, m = 1/2, th, phi},
cn = JacobiCN[x 2 EllipticK[m] + I 2 y EllipticK[1 - m], m];
{th, phi} = {2 ArcTan[Abs[cn]], Arg[cn]};
{th, phi}]

zf[{th_, phi_}] := Cos[th]
xf[{th_, phi_}] := Sin[th] Cos[phi]
yf[{th_, phi_}] := Sin[th] Sin[phi]


The contours that correspond to the equation $\theta=\frac{\pi}{2}$, besides the expected tilted square, also includes the diagonal and anti-diagonal. This is not an artifact and can be proven. These lines make me suspect that I am not doing the projection the right way. I would be glad to see a prototype implementation of it in Mathematica.

Thank you.

• The implementation I did is for images. Would this be admissible? (I can jury-rig something that transforms a mesh, tho.) – J. M.'s technical difficulties Apr 29 '13 at 4:23
• Yes, completely. Additionally, there is a fun way to use this projection to build spherical panoramas, which is what I am ultimately after. arXiv:1011.3189 and wikipedia provide examples. – Sasha Apr 29 '13 at 4:26
• I haven't tried turning my implementation into something for producing panoramas, but okay, I'll put in my take... – J. M.'s technical difficulties Apr 29 '13 at 4:27
• For more information on the cartographic applications of elliptic functions, see this. – J. M.'s technical difficulties Apr 29 '13 at 14:03

In demonstrating how the quincuncial projection works, consider first the following complex mapping:

With[{ω = N[EllipticK[1/2], 20]},
ParametricPlot[{Re[InverseJacobiCN[Tan[φ/2] Exp[I θ], 1/2]],
Im[InverseJacobiCN[Tan[φ/2] Exp[I θ], 1/2]]},
{φ, 0, π}, {θ, -π, π}, Mesh -> Automatic,
MeshStyle -> {Orange, Green}, PlotPoints -> 75,
PlotRange -> {{0, 2 ω}, {-ω, ω}}]] The components should be easily recognizable (e.g. the argument of the elliptic integral is a stereographic projection/Riemann sphere representation of a complex number).

The only fixed thing in applying the projection is the use of either the cosine amplitude function $\mathrm{cn}\left(u\mid \frac12\right)$ or its inverse; the choice of which stereographic projection to do is dependent on what coordinate convention you take for spherical coordinates.

For instance, if your convention is longitude/latitude (as with the output of CountryData[]), here's how one might implement the projection:

world = {{Switch[CountryData[#, "Continent"],
"Asia", Yellow, "Oceania", Green, "Europe", Pink,
"NorthAmerica", Red, "SouthAmerica", Orange, "Africa", GrayLevel[1/10],
"Antarctica", GrayLevel[9/10], _, Blue],
CountryData[#, {"FullPolygon", "Equirectangular"}]} & /@
Append[CountryData[], "Antarctica"]} /.
{θ_?NumericQ, φ_?NumericQ} :>
Through[{Re, Im}[InverseJacobiCN[Cot[φ °/2 + π/4] Exp[I θ °], 1/2]]];

With[{ω = N[EllipticK[1/2], 20]},
tile = Image[Graphics[Prepend[world,
{ColorData["Legacy", "PowderBlue"], Rectangle[{0, -ω}, {2 ω, ω}]}],
ImageResolution -> 300]] The only snag in this is that some post-processing is necessary if one wants to remove the polygons that are turned into "slivers" by the transformation, as can be seen when one tries to tile the image given above:

Graphics[{Texture[tile],
Polygon[{{0, 0}, {1, 0}, {1, 1}, {0, 1}},
VertexTextureCoordinates -> {{1, 0}, {1, 1}, {0, 1}, {0, 0}}],
Polygon[{{1, 0}, {2, 0}, {2, 1}, {1, 1}},
VertexTextureCoordinates -> {{0, 1}, {0, 0}, {1, 0}, {1, 1}}],
Polygon[{{0, 1}, {1, 1}, {1, 2}, {0, 2}},
VertexTextureCoordinates -> {{0, 1}, {0, 0}, {1, 0}, {1, 1}}],
Polygon[{{1, 1}, {2, 1}, {2, 2}, {1, 2}},
VertexTextureCoordinates -> {{1, 0}, {1, 1}, {0, 1}, {0, 0}}]}] (You can do the sliver removal yourself, if you want it.)

If, like me, you prefer the longitude/colatitude convention, the stereographic projection proceeds a bit differently. For this example, I'll transform an image instead of transforming polygons. ImageTransformation[] does a nice job for this route:

earthGrid = Import["http://i.stack.imgur.com/Zzox0.png"];
With[{ω = N[EllipticK[1/2], 20]},
eg = ImageTransformation[earthGrid,
With[{w = JacobiCN[#[] + I #[], 1/2]}, {Arg[w], 2 ArcCot[Abs[w]]}] &,
DataRange -> {{-π, π}, {0, π}}, Padding -> 1., PlotRange -> {{0, 2 ω}, {-ω, ω}}]] Using code similar to the one given above, we can see the tiling for this as well: For image transformation purposes, however, I have found the execution of the built-in JacobiCN[] for complex arguments to be a bit slow, so I wrote my own implementation of a function that can replace JacobiCN[z, 1/2]:

SetAttributes[cnhalf, Listable];
cnhalf[z_?NumericQ] := Block[{nz = N[z], k, zs},
k = Max[0, Ceiling[Log2[4 Abs[nz]]]];
zs = (nz 2^-k)^2;
Nest[With[{cs = #^2}, -(((cs + 2) cs - 1)/((cs - 2) cs - 1))] &,
(1 - zs/4 (1 + zs/30 (1 + zs/8)))/(1 + zs/4 (1 - zs/30 (1 - zs/8))), k]]


which works quite nicely. (Exercise: try to recognize the formulae I used.)

Using, for instance, the ETOPO1 global relief,

etopo1 = Import["http://www.ngdc.noaa.gov/mgg/image/color_etopo1_ice_low.jpg"];


we finally present the other way to demonstrate the quincuncial projection:

With[{ω = N[EllipticK[1/2], 20]},
etp = ImageTransformation[etopo1,
With[{w = cnhalf[#[] + I #[]]}, {Arg[w], 2 ArcCot[Abs[w]]}] &,
DataRange -> {{-π, π}, {0, π}}, Padding -> 1.,
PlotRange -> {{0, 2 ω}, {-ω, ω}}]];

Graphics[{Texture[etp],
Polygon[{{0, 0}, {1, 0}, {1, 1}, {0, 1}},
VertexTextureCoordinates -> {{1, 1}, {0, 1}, {0, 0}, {1, 0}}],
Polygon[{{1, 0}, {3/2, 0}, {3/2, 1/2}, {1, 1/2}},
VertexTextureCoordinates -> {{0, 0}, {1/2, 0}, {1/2, 1/2}, {0, 1/2}}],
Polygon[{{1/2, -1/2}, {1, -1/2}, {1, 0}, {1/2, 0}},
VertexTextureCoordinates -> {{1/2, 1/2}, {1, 1/2}, {1, 1}, {1/2, 1}}],
Polygon[{{1, -1/2}, {3/2, -1/2}, {3/2, 0}, {1, 0}},
VertexTextureCoordinates -> {{1, 1/2}, {1/2, 1/2}, {1/2, 0}, {1, 0}}]}] (I had previously posted this image in chat. Now you know where it came from. ;))

For an out-of-this-world bonus, here's another image that I have quincuncially projected: Can you guess the original equirectangular image I made this from?

• +1 Want to know something about cartographic projections long ago. – Silvia Apr 29 '13 at 6:07
• eugeneparnell.com/wp-content/uploads/2011/07/… – Dr. belisarius Apr 29 '13 at 12:32
• @bel, "God plays roulette with the Earth"? – J. M.'s technical difficulties Apr 29 '13 at 12:38
• +1 for linking to a paper from 1879. I'd give more, but there are limits to this sort of thing. :) – rcollyer Apr 29 '13 at 19:53
• @rcollyer, we should then be impressed that Peirce imagined all this without having Mathematica on his hands. :) – J. M.'s technical difficulties Apr 30 '13 at 1:32

Just an update. This is now built in in WL:

GeoProjectionData["PeirceQuincuncial"]


{"PeirceQuincuncial", {"Centering" -> {90, -90}, "GridOrigin" -> {0, 0}, "ReferenceModel" -> 1}}

GeoGraphics["World",
GeoProjection -> "PeirceQuincuncial",
GeoBackground -> "ReliefMap"] 