Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
    
The implementation I did is for images. Would this be admissible? (I can jury-rig something that transforms a mesh, tho.) –  J. M. 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. Apr 29 '13 at 4:27
    
For more information on the cartographic applications of elliptic functions, see this. –  J. M. Apr 29 '13 at 14:03

1 Answer 1

up vote 24 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, π}, {θ, -π, π}, 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}(u\mid m)$ 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 postprocessing 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?

share|improve this answer
    
+1 Want to know something about cartographic projections long ago. –  Silvia Apr 29 '13 at 6:07
    
    
@bel, "God plays roulette with the Earth"? –  J. M. 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
2  
@rcollyer, we should then be impressed that Peirce imagined all this without having Mathematica on his hands. :) –  J. M. Apr 30 '13 at 1:32

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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