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.

contour of inverse Peirce mapping

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. is computer-less 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. is computer-less Apr 29 '13 at 4:27
  • For more information on the cartographic applications of elliptic functions, see this. – J. M. is computer-less Apr 29 '13 at 14:03
up vote 44 down vote accepted

(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 ω}, {-ω, ω}}]]

Peirce's projection

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 ω, ω}]}],
                       ImagePadding -> None, PlotRangePadding -> None],
              ImageResolution -> 300]]

quincuncially-projected world map

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}}]}]

A Peirce tiling

(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 ω}, {-ω, ω}}]]

quincuncially-projected image

Using code similar to the one given above, we can see the tiling for this as well: Peirce tiled image

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}}]}]

Another Peirce projection

(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:

http://www.jpl.nasa.gov/spaceimages/details.php?id=PIA15482

Can you guess the original equirectangular image I made this from?

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"]

enter image description here

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.