Fix and additional options for Peirce Quincuncial projections (#2978)

This fixes the current forward implementation of Peirce Quincuncial proj to correctly flip/reflect out the southern hemisphere to four triangles, and rotate entire result to a square or diamond. (It there resolves the issues identified with pull request #2230 , where southern hemisphere was wrongly projected over northern, and reverses the restriction to northern hemisphere introduced there). It also adds additional lateral projection of the hemispheres.

- This PR adds an optional parameter `+type` which allows selection of projection. The `+type=square` and `+type=diamond` types match in principle ESRI's twin implementations of square and diamond PQ projs. The **default** if not specified is `+type=diamond`.

- The previous behaviour restricted to the northern hemisphere can be reproduced using the `+type=nhemisphere`, though this is an edge case only.

- An additional `+type=horizontal` and `+type=vertical` rectangular lateral versions have been added that place each hemisphere side-by-side. This is primarily to allow creation of projections such as Greiger Triptychial, which also require the additional optional params `scrollx` or `scrolly` in order to shift parts of the projection from one side of the map to the other. 

- Additional documentation has been added to proj description, including quoting the usual meridian used in common usage of projection, and images showing the different types.

docs/plot/plotdefs.json

@@ -1051,13 +1051,46 @@
         "type": "poly"
     },
     {
-        "filename": "peirce_q.png",
+        "filename": "peirce_q_square.png",
         "latmax": 90,
-        "latmin": 0,
+        "latmin": -90,
+        "lonmax": 180,
+        "lonmin": -180,
+        "name": "peirce_q_square",
+        "projstring": "+proj=peirce_q +lon_0=25 +type=square",
+        "res": "low",
+        "type": "poly"
+    },
+    {
+        "filename": "peirce_q_diamond.png",
+        "latmax": 90,
+        "latmin": -90,
         "lonmax": 180,
         "lonmin": -180,
         "name": "peirce_q",
-        "projstring": "+proj=peirce_q",
+        "projstring": "+proj=peirce_q +lon_0=25 +type=diamond",
+        "res": "low",
+        "type": "poly"
+    },
+    {
+        "filename": "peirce_q_horizontal.png",
+        "latmax": 90,
+        "latmin": -90,
+        "lonmax": 180,
+        "lonmin": -180,
+        "name": "peirce_q_horizontal",
+        "projstring": "+proj=peirce_q +lon_0=25 +type=horizontal",
+        "res": "low",
+        "type": "line"
+    },
+    {
+        "filename": "grieger_triptychial.png",
+        "latmax": 90,
+        "latmin": -90,
+        "lonmax": 180,
+        "lonmin": -180,
+        "name": "grieger_triptychial",
+        "projstring": "+proj=pipeline +step +proj=ob_tran +o_proj=peirce_q +o_lat_p=-45 +o_lon_p=45 +type=horizontal +scrollx=-0.25 +step +proj=affine +s11=-1 +s12=0 +s21=0 +s22=-1",
         "res": "low",
         "type": "poly"
     },

docs/source/operations/projections/peirce_q.rst

@@ -4,6 +4,23 @@
 Peirce Quincuncial
 ********************************************************************************
 
+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 or can be rotated 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.
+The diamond and square versions can be produced using the
+``+type=diamond`` and ``+type=square`` options respectively.
+This implementation includes an alternative lateral projection
+which places hemispheres side-by-side (``+type=horizontal`` or
+``+type=vertical``). Combined with a general oblique transformation,
+this can be used to produced a Grieger Triptychial projection
+(see example below).
+
 +---------------------+----------------------------------------------------------+
 | **Classification**  | Miscellaneous                                            |
 +---------------------+----------------------------------------------------------+
@@ -21,12 +38,33 @@ Peirce Quincuncial
 +---------------------+----------------------------------------------------------+
 
 
-.. figure:: ./images/peirce_q.png
+.. figure:: ./images/peirce_q_square.png
+   :width: 500 px
+   :align: center
+   :alt:   Peirce Quincuncial (Square)
+
+   proj-string: ``+proj=peirce_q +lon_0=25 +type=square``
+
+.. figure:: ./images/peirce_q_diamond.png
+   :width: 500 px
+   :align: center
+   :alt:   Peirce Quincuncial (Diamond)
+
+   proj-string: ``+proj=peirce_q +lon_0=25 +type=diamond``
+
+.. figure:: ./images/peirce_q_horizontal.png
+   :width: 500 px
+   :align: center
+   :alt:   Peirce Quincuncial (Horizontal)
+
+   proj-string: ``+proj=peirce_q +lon_0=25 +type=horizontal``
+
+.. figure:: ./images/grieger_triptychial.png
    :width: 500 px
    :align: center
-   :alt:   Peirce Quincuncial
+   :alt:   Grieger Triptychial
 
-   proj-string: ``+proj=peirce_q``
+   proj-string: ``+proj=pipeline +step +proj=ob_tran +o_proj=peirce_q +o_lat_p=-45 +o_lon_p=45 +o_type=horizontal +o_scrollx=-0.25 +step +proj=affine +s11=-1 +s12=0 +s21=0 +s22=-1``
 
 Parameters
 ################################################################################
@@ -35,6 +73,41 @@ Parameters
 
 .. include:: ../options/lon_0.rst
 
+.. option:: +type=square/diamond/horizontal/vertical/nhemisphere/shemisphere
+
+    .. versionadded:: 8.2.1
+
+    *Defaults to diamond.*
+
+    Indicates the type of transformation applied to the southern hemisphere:
+    ``square`` and ``diamond`` represent the traditional quincuncial form suggested
+    by Peirce with the southern hemisphere divided into 4 triangles and reflected
+    outward from the northern hemisphere. The ``square`` type is rotated by 45
+    degrees to produce the conventional square presentation. The origin lies at
+    the centre of the square or diamond.
+
+    By contrast, the ``horizontal`` and ``vertical`` forms reflect the southern
+    hemisphere laterally across the x or y axis respectively to produce a rectangular
+    form. The origin lies at the centre of the rectangle.
+
+    The other two types, ``nhemisphere`` and ``shemisphere``, discard latitudes of less
+    than 0 or more than 0, respectively, to allow single hemispheres to be selected.
+    The origin lies at the centre of the square or diamond.
+
+.. option:: +scrollx=<value>
+
+    For ``horizontal`` type allows a scalar circular scroll of resulting x coordinates
+    to shift sections of the projection to the other horizontal side of the map.
+
+    *Defaults to 0.0. Must be a scale between -1.0 and 1.0.*
+
+.. option:: +scrolly=<value>
+
+    For ``vertical`` type allows a scalar circular scroll of resulting y coordinates
+    to shift sections of the projection to the other vertical side of the map.
+
+    *Defaults to 0.0. Must be a scale between -1.0 and 1.0.*
+
 .. include:: ../options/R.rst
 
 .. include:: ../options/x_0.rst

src/projections/adams.cpp

@@ -6,9 +6,21 @@
 * 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.
+* Fix for Peirce Quincuncial projection to diamond or square 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
+* Added original code for horizontal and vertical arrangement of hemispheres by Toby
+* C. Wilkinson to allow creation of lateral quincuncial projections, such as Grieger's
+* Triptychial, see, e.g.:
+* - Grieger, B. (2020). “Optimized global map projections for specific applications:
+* the triptychial projection and the Spilhaus projection”. EGU2020-9885.
+* https://doi.org/10.5194/egusphere-egu2020-9885
 *
 * Copyright (c) 2005, 2006, 2009 Gerald I. Evenden
 * Copyright (c) 2020 Kristian Evers
+* Copyright (c) 2021 Toby C Wilkinson
 *
 * Related material
 * ----------------
@@ -48,8 +60,20 @@ enum projection_type {
     ADAMS_WS2,
 };
 
+enum peirce_type {
+  PEIRCE_Q_SQUARE,
+  PEIRCE_Q_DIAMOND,
+  PEIRCE_Q_NHEMISPHERE,
+  PEIRCE_Q_SHEMISPHERE,
+  PEIRCE_Q_HORIZONTAL,
+  PEIRCE_Q_VERTICAL,
+};
+
 struct pj_opaque {
     projection_type mode;
+    peirce_type pqtype;
+    double scrollx = 0.0;
+    double scrolly = 0.0;
 };
 
 } // anonymous namespace
@@ -118,9 +142,18 @@ static PJ_XY adams_forward(PJ_LP lp, PJ *P) {
         }
         break;
     case PEIRCE_Q: {
-            if( lp.phi < -TOL ) {
+            /* lam0 - note that the original Peirce model used a central meridian of around -70deg, but the default within proj is +lon0=0 */
+            if (Q->pqtype == PEIRCE_Q_NHEMISPHERE) {
+              if( lp.phi < -TOL ) {
+                proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
+                return proj_coord_error().xy;
+              }
+            }
+            if (Q->pqtype == PEIRCE_Q_SHEMISPHERE) {
+              if( lp.phi > -TOL ) {
                 proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
                 return proj_coord_error().xy;
+              }
             }
             const double sl = sin(lp.lam);
             const double cl = cos(lp.lam);
@@ -173,7 +206,71 @@ 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) {
+      /* Constant complete elliptic integral of the first kind with m=0.5, calculated using https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.ellipk.html . Used as basic as scaled shift distance */
+      constexpr double shd = 1.8540746773013719 * 2;
+
+      /* For square and diamond Quincuncial projections, spin out southern hemisphere to triangular segments of quincunx (before rotation for square)*/
+      if( Q->pqtype == PEIRCE_Q_SQUARE || ( Q->pqtype == PEIRCE_Q_DIAMOND )) {
+        if (lp.phi < 0.)  { /* fold out segments */
+          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 */
+        }
+      }
+
+      /* For square types rotate xy by 45 deg */
+      if( Q->pqtype == PEIRCE_Q_SQUARE ) {
+            const double temp = xy.x;
+            xy.x = RSQRT2 * (xy.x - xy.y);
+            xy.y = RSQRT2 * (temp + xy.y);
+      }
+
+      /* For rectangle Quincuncial projs, spin out southern hemisphere to east (horizontal) or north (vertical) after rotation */
+      if( Q->pqtype == PEIRCE_Q_HORIZONTAL ) {
+        if (lp.phi < 0.)  {
+          xy.x = shd - xy.x; /* reflect x to east */
+        }
+        xy.x = xy.x - (shd / 2); /* shift everything so origin is in middle of two hemispheres */
+      }
+      if( Q->pqtype == PEIRCE_Q_VERTICAL ) {
+        if (lp.phi < 0.)  {
+          xy.y = shd - xy.y; /* reflect y to north */
+        }
+        xy.y = xy.y - (shd / 2); /* shift everything so origin is in middle of two hemispheres */
+      }
+
+      //if o_scrollx param present, scroll x
+      if (!(Q->scrollx == 0.0) && (Q->pqtype == PEIRCE_Q_HORIZONTAL) ) {
+        double xscale = 2.0;
+        double xthresh = shd / 2;
+        xy.x = xy.x + (Q->scrollx * (xthresh * 2 * xscale)); /*shift relative to proj width*/
+        if(xy.x >= (xthresh * xscale)) {
+          xy.x = xy.x - (shd * xscale);
+        }
+        else if (xy.x < -(xthresh * xscale)) {
+          xy.x = xy.x + (shd * xscale);
+        }
+      }
+
+      //if o_scrolly param present, scroll y
+      if (!(Q->scrolly == 0.0) && (Q->pqtype == PEIRCE_Q_VERTICAL)) {
+        double yscale = 2.0;
+        double ythresh = shd / 2;
+        xy.y = xy.y + (Q->scrolly * (ythresh * 2 * yscale)); /*shift relative to proj height*/
+        if(xy.y >= (ythresh * yscale)) {
+          xy.y = xy.y - (shd * yscale);
+        }
+        else if (xy.y < -(ythresh * yscale)) {
+          xy.y = xy.y + (shd * yscale);
+        }
+      }
+
+  }
+
+    if (Q->mode == ADAMS_HEMI || Q->mode == ADAMS_WS2 ) { /* rotate by 45deg. */
         const double temp = xy.x;
         xy.x = RSQRT2 * (xy.x - xy.y);
         xy.y = RSQRT2 * (temp + xy.y);
@@ -218,6 +315,55 @@ static PJ *setup(PJ *P, projection_type mode) {
     if( mode == ADAMS_WS2 )
         P->inv = adams_inverse;
 
+    if( mode == PEIRCE_Q) {
+      // Quincuncial projections type options: square, diamond, hemisphere, horizontal (rectangle) or vertical (rectangle)
+      const char* pqtype = pj_param (P->ctx, P->params, "stype").s;
+
+      if (!pqtype) pqtype = "diamond"; /* default if type value not supplied */
+
+      if (strcmp(pqtype, "square") == 0) {
+        Q->pqtype = PEIRCE_Q_SQUARE;
+      }
+      else if (strcmp(pqtype, "diamond") == 0) {
+        Q->pqtype = PEIRCE_Q_DIAMOND;
+      }
+      else if (strcmp(pqtype, "nhemisphere") == 0) {
+        Q->pqtype = PEIRCE_Q_NHEMISPHERE;
+      }
+      else if (strcmp(pqtype, "shemisphere") == 0) {
+        Q->pqtype = PEIRCE_Q_SHEMISPHERE;
+      }
+      else if (strcmp(pqtype, "horizontal") == 0) {
+        Q->pqtype = PEIRCE_Q_HORIZONTAL;
+        if (pj_param(P->ctx, P->params, "tscrollx").i) {
+          double scrollx;
+          scrollx = pj_param(P->ctx, P->params, "dscrollx").f;
+          if (scrollx > 1 || scrollx < -1) {
+              proj_log_error(P, _("Invalid value for scrollx: |scrollx| should between -1 and 1"));
+              return pj_default_destructor (P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
+          }
+          Q->scrollx = scrollx;
+        }
+      }
+      else if (strcmp(pqtype, "vertical") == 0) {
+        Q->pqtype = PEIRCE_Q_VERTICAL;
+        if (pj_param(P->ctx, P->params, "tscrolly").i) {
+          double scrolly;
+          scrolly = pj_param(P->ctx, P->params, "dscrolly").f;
+          if (scrolly > 1 || scrolly < -1) {
+              proj_log_error(P, _("Invalid value for scrolly: |scrolly| should between -1 and 1"));
+              return pj_default_destructor (P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
+          }
+          Q->scrolly = scrolly;
+        }
+      }
+      else {
+            proj_log_error (P, _("peirce_q: invalid value for 'type' parameter"));
+            return pj_default_destructor (P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
+      }
+
+    }
+
     return P;
 }