Fix and additional options for Peirce Quincuncial projections

docs/plot/plotdefs.json

@@ -1053,11 +1053,22 @@
     {
         "filename": "peirce_q.png",
         "latmax": 90,
-        "latmin": 0,
+        "latmin": -90,
         "lonmax": 180,
         "lonmin": -180,
         "name": "peirce_q",
-        "projstring": "+proj=peirce_q",
+        "projstring": "+proj=peirce_q +lon_0=25",
+        "res": "low",
+        "type": "poly"
+    },
+    {
+        "filename": "peirce_q_d.png",
+        "latmax": 90,
+        "latmin": -90,
+        "lonmax": 180,
+        "lonmin": -180,
+        "name": "peirce_q_d",
+        "projstring": "+proj=peirce_q_d +lon_0=25",
         "res": "low",
         "type": "poly"
     },

docs/source/operations/projections/peirce_q.rst

@@ -1,9 +1,22 @@
 .. _peirce_q:
 
 ********************************************************************************
-Peirce Quincuncial
+Peirce Quincuncial (Square)
 ********************************************************************************
 
+The Peirce Quincuncial projection is a conformal map projection
+that transforms the circle of the northern hemisphere into a square,
+and the southern hemisphere split into four triangles arranged
+around the square to form a quincunx. The resulting projection
+is a regular diamond shape that is usually rotated 45 degrees to form a square.
+The resulting tile can be infinitely tessellated. Though this implementation
+defaults to a central meridian of 0, it is more common to use a central
+meridian of around 25 to optimise the distortions. Peirce's original
+published map from 1879 used a central meridian of approx -70.
+
+Use Peirce Quincuncial (Diamond) to skip the rotation of the diamond to
+a square.
+
 +---------------------+----------------------------------------------------------+
 | **Classification**  | Miscellaneous                                            |
 +---------------------+----------------------------------------------------------+
@@ -24,9 +37,9 @@ Peirce Quincuncial
 .. figure:: ./images/peirce_q.png
    :width: 500 px
    :align: center
-   :alt:   Peirce Quincuncial
+   :alt:   Peirce Quincuncial (Square)
 
-   proj-string: ``+proj=peirce_q``
+   proj-string: ``+proj=peirce_q +lon_0=25``
 
 Parameters
 ################################################################################

docs/source/operations/projections/peirce_q_d.rst

@@ -0,0 +1,54 @@
+.. _peirce_q_d:
+
+********************************************************************************
+Peirce Quincuncial (Diamond)
+********************************************************************************
+
+The Peirce Quincuncial projection is a conformal map projection
+that transforms the circle of the northern hemisphere into a square,
+and the southern hemisphere split into four triangles arranged
+around the square to form a quincunx. The resulting projection
+is a regular diamond shape that is usually rotated 45 degrees to form a square.
+The resulting tile can be infinitely tessellated. Though this implementation
+defaults to a central meridian of 0, it is more common to use a central
+meridian of around 25 to optimise the distortions. Peirce's original
+published map from 1879 used a central meridian of approx -70.
+
+Use Peirce Quincuncial (Square) to rotate diamond to a square.
+
++---------------------+----------------------------------------------------------+
+| **Classification**  | Miscellaneous                                            |
++---------------------+----------------------------------------------------------+
+| **Available forms** | Forward spherical projection                             |
++---------------------+----------------------------------------------------------+
+| **Defined area**    | Global                                                   |
++---------------------+----------------------------------------------------------+
+| **Alias**           | peirce_q_d                                                 |
++---------------------+----------------------------------------------------------+
+| **Domain**          | 2D                                                       |
++---------------------+----------------------------------------------------------+
+| **Input type**      | Geodetic coordinates                                     |
++---------------------+----------------------------------------------------------+
+| **Output type**     | Projected coordinates                                    |
++---------------------+----------------------------------------------------------+
+
+
+.. figure:: ./images/peirce_q_d.png
+   :width: 500 px
+   :align: center
+   :alt:   Peirce Quincuncial (Diamond)
+
+   proj-string: ``+proj=peirce_q_d +lon_0=25``
+
+Parameters
+################################################################################
+
+.. note:: All parameters are optional.
+
+.. include:: ../options/lon_0.rst
+
+.. include:: ../options/R.rst
+
+.. include:: ../options/x_0.rst
+
+.. include:: ../options/y_0.rst

src/pj_list.h

@@ -122,7 +122,8 @@ PROJ_HEAD(ortel, "Ortelius Oval")
 PROJ_HEAD(ortho, "Orthographic")
 PROJ_HEAD(pconic, "Perspective Conic")
 PROJ_HEAD(patterson, "Patterson Cylindrical")
-PROJ_HEAD(peirce_q, "Peirce Quincuncial")
+PROJ_HEAD(peirce_q, "Peirce Quincuncial (Square)")
+PROJ_HEAD(peirce_q_d, "Peirce Quincuncial (Diamond)")
 PROJ_HEAD(pipeline, "Transformation pipeline manager")
 PROJ_HEAD(poly, "Polyconic (American)")
 PROJ_HEAD(pop, "Retrieve coordinate value from pipeline stack")

src/projections/adams.cpp

@@ -6,6 +6,11 @@
 * PROJ by Kristian Evers. Original code found in file src/proj_guyou.c, see
 * https://github.com/rouault/libproj4/blob/master/libproject-1.01/src/proj_guyou.c
 * for reference.
+* Corrections added for Peirce Quincuncial projection by Toby C. Wilkinson to
+* correctly flip out southern hemisphere into the four triangles of Peirce's
+* quincunx. The fix inspired by a similar rotate and translate solution applied
+* by Jonathan Feinberg for cartopy, see
+* https://github.com/jonathf/cartopy/blob/8172cac7fc45cafc86573d408ce85b74258a9c28/lib/cartopy/peircequincuncial.py
 *
 * Copyright (c) 2005, 2006, 2009 Gerald I. Evenden
 * Copyright (c) 2020 Kristian Evers
@@ -33,7 +38,8 @@
 #include "proj_internal.h"
 
 PROJ_HEAD(guyou, "Guyou") "\n\tMisc Sph No inv";
-PROJ_HEAD(peirce_q, "Peirce Quincuncial") "\n\tMisc Sph No inv";
+PROJ_HEAD(peirce_q, "Peirce Quincuncial (Square)") "\n\tMisc Sph No inv";
+PROJ_HEAD(peirce_q_d, "Peirce Quincuncial (Diamond)") "\n\tMisc Sph No inv";
 PROJ_HEAD(adams_hemi, "Adams Hemisphere in a Square") "\n\tMisc Sph No inv";
 PROJ_HEAD(adams_ws1, "Adams World in a Square I") "\n\tMisc Sph No inv";
 PROJ_HEAD(adams_ws2, "Adams World in a Square II") "\n\tMisc Sph No inv";
@@ -43,6 +49,7 @@ namespace { // anonymous namespace
 enum projection_type {
     GUYOU,
     PEIRCE_Q,
+    PEIRCE_Q_D,
     ADAMS_HEMI,
     ADAMS_WS1,
     ADAMS_WS2,
@@ -117,11 +124,9 @@ static PJ_XY adams_forward(PJ_LP lp, PJ *P) {
             sn = lp.phi < 0.;
         }
         break;
+    case PEIRCE_Q_D:
     case PEIRCE_Q: {
-            if( lp.phi < -TOL ) {
-                proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
-                return proj_coord_error().xy;
-            }
+      /* Note that the original Peirce model used a central meridian of around -70, but the default for proj is, atypically, +lon0=0 */
             const double sl = sin(lp.lam);
             const double cl = cos(lp.lam);
             const double cp = cos(lp.phi);
@@ -173,7 +178,21 @@ static PJ_XY adams_forward(PJ_LP lp, PJ *P) {
     xy.x = ell_int_5(m);
     xy.y = ell_int_5(n);
 
-    if (Q->mode == ADAMS_HEMI || Q->mode == ADAMS_WS2) { /* rotate by 45deg. */
+    if (Q->mode == PEIRCE_Q || Q->mode == PEIRCE_Q_D) {
+      /* For Quincuncial projections, spin out southern hemisphere to triangular segments of quincunx */
+      if (true && (lp.phi < 0.)) {
+        /* Constant complete elliptic integral of the first kind with m=0.5 https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.ellipk.html as shift distance */
+        double shd = 1.8540746773013719 * 2;
+
+        if (lp.lam < ( -0.75 * M_PI )) xy.y = shd - xy.y; /* top left segment, shift up and reflect y */
+        if ( (lp.lam < (-0.25 * M_PI)) && (lp.lam >= ( -0.75 * M_PI ))) xy.x = - shd - xy.x; /* left segment, shift left and reflect x */
+        if ( (lp.lam < (0.25 * M_PI)) && (lp.lam >= ( -0.25 * M_PI ))) xy.y = - shd - xy.y; /* bottom segment, shift down and reflect y */
+        if ( (lp.lam < (0.75 * M_PI)) && (lp.lam >= ( 0.25 * M_PI ))) xy.x = shd - xy.x; /* right segment, shift right and reflect x */
+        if (lp.lam >= (0.75 * M_PI)) xy.y = shd - xy.y; /* top right segment, shift up and reflect y */
+      }
+    }
+
+    if (Q->mode == ADAMS_HEMI || Q->mode == ADAMS_WS2 || Q->mode == PEIRCE_Q ) { /* rotate by 45deg. */
         const double temp = xy.x;
         xy.x = RSQRT2 * (xy.x - xy.y);
         xy.y = RSQRT2 * (temp + xy.y);
@@ -230,6 +249,10 @@ PJ *PROJECTION(peirce_q) {
     return setup(P, PEIRCE_Q);
 }
 
+PJ *PROJECTION(peirce_q_d) {
+    return setup(P, PEIRCE_Q_D);
+}
+
 PJ *PROJECTION(adams_hemi) {
     return setup(P, ADAMS_HEMI);
 }