# 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.) Apr 29, 2013 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. Apr 29, 2013 at 4:26
• I haven't tried turning my implementation into something for producing panoramas, but okay, I'll put in my take... Apr 29, 2013 at 4:27
• For more information on the cartographic applications of elliptic functions, see this. Apr 29, 2013 at 14:03

(I had been meaning to write a blog entry about this myself, but since this question has come up, I suppose I'll just write about it here instead...)

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[#[[1]] + I #[[2]], 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[#[[1]] + I #[[2]]]}, {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. Apr 29, 2013 at 6:07
• eugeneparnell.com/wp-content/uploads/2011/07/… Apr 29, 2013 at 12:32
• @bel, "God plays roulette with the Earth"? Apr 29, 2013 at 12:38
• +1 for linking to a paper from 1879. I'd give more, but there are limits to this sort of thing. :) Apr 29, 2013 at 19:53
• @rcollyer, we should then be impressed that Peirce imagined all this without having Mathematica on his hands. :) Apr 30, 2013 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"]


• That's neat! Since when was it incorporated? May 11, 2016 at 0:14
• @J.M. Very recently I think. May 11, 2016 at 6:59
• @VitaliyKaurov Wondering if you can extract the projection as a function for use in ImageTransformation?
– M.R.
Oct 5, 2020 at 20:19