Add Spilhaus projection

docs/plot/plotdefs.json

@@ -1345,6 +1345,18 @@
         "type": "line",
         "delta_cut": 1e7
     },
+    {
+        "filename": "spilhaus.png",
+        "latmax": 90,
+        "latmin": -90,
+        "lonmax": 180,
+        "lonmin": -180,
+        "name": "spilhaus",
+        "projstring": "+proj=spilhaus",
+        "res": "low",
+        "type": "line",
+        "delta_cut": 1e6
+    },
     {
         "filename": "stere.png",
         "latmax": 90.0,

docs/source/operations/projections/index.rst

@@ -135,6 +135,7 @@ Projections map the spherical 3D space to a flat 2D space.
    sinu
    som
    somerc
+   spilhaus
    stere
    sterea
    gstmerc

docs/source/operations/projections/spilhaus.rst

@@ -0,0 +1,89 @@
+.. _spilhaus:
+
+********************************************************************************
+Spilhaus
+********************************************************************************
+
+.. versionadded:: 9.6.0
+
++---------------------+----------------------------------------------------------+
+| **Classification**  | Miscellaneous                                            |
++---------------------+----------------------------------------------------------+
+| **Available forms** | Forward and inverse, spherical, ellipsoidal              |
++---------------------+----------------------------------------------------------+
+| **Defined area**    | Global                                                   |
++---------------------+----------------------------------------------------------+
+| **Alias**           | spilhaus                                                 |
++---------------------+----------------------------------------------------------+
+| **Domain**          | 2D                                                       |
++---------------------+----------------------------------------------------------+
+| **Input type**      | Geodetic coordinates                                     |
++---------------------+----------------------------------------------------------+
+| **Output type**     | Projected coordinates                                    |
++---------------------+----------------------------------------------------------+
+
+
+.. figure:: ./images/spilhaus.png
+   :width: 500 px
+   :align: center
+   :alt:   Spilhaus
+
+   proj-string: ``+proj=spilhaus``
+
+Parameters
+################################################################################
+
+.. note:: All parameters are optional.
+
+.. option:: +lon_c=<value>
+
+    Longitude of projection centre.
+
+    *Defaults to 66.94970198.*
+
+    .. note::
+        The default convention is to interpret this value as decimal degrees. To
+        specify radians instead, follow the value with the "r" character.
+
+        Example: `+lat_0=1.570796r`
+
+        See :ref:`Projection Units <projection_units>` for more information.
+
+.. option:: +lat_c=<value>
+
+    Latitude of projection centre.
+
+    *Defaults to -49.56371678.*
+
+    .. note::
+        The default convention is to interpret this value as decimal degrees. To
+        specify radians instead, follow the value with the "r" character.
+
+        Example: `+lat_0=1.570796r`
+
+        See :ref:`Projection Units <projection_units>` for more information.
+
+.. option:: +azi=<value>
+
+    Azimuth from north at the center of projection.
+
+    *Defaults to 40.17823482.*
+
+.. option:: +rot=<value>
+
+    Rotation of the projection.
+
+    *Defaults to 45.*
+
+.. include:: ../options/R.rst
+
+.. option:: +k_0=<value>
+
+    Scale factor. Determines scale factor used in the projection.
+    To reproduce the output of ESRI:54099, use +k_0=1.41421356237 (sqrt(2))
+
+    *Defaults to 1.0.*
+
+.. include:: ../options/x_0.rst
+
+.. include:: ../options/y_0.rst

docs/source/spelling_wordlist.txt

@@ -875,10 +875,13 @@ specializations
 specialized
 Spectroradiometer
 spherification
+spilhaus
+Spilhaus
 splitted
 spx
 sql
 sqlite
+sqrt
 src
 srs
 standardised

src/latitudes.cpp

@@ -0,0 +1,68 @@
+/*
+ * Project:  PROJ
+ *
+ * Helper functions to compute latitudes
+ *
+ * Some from Map Projections - A Working Manual. 1987. John P. Snyder
+ * https://neacsu.net/docs/geodesy/snyder/2-general/sect_3/
+ *
+ * Copyright (c) 2025 Javier Jimenez Shaw
+ */
+
+#include <math.h>
+
+#include "proj_internal.h"
+
+/*****************************************************************************/
+double pj_conformal_lat(double phi, double e) {
+    /***********************************
+     * The conformal latitude chi
+     * giving a sphere which is truly conformal in accordance with the ellipsoid
+     * Snyder, A working manual, formula 3-1
+     *
+     * phi: geodetic latitude, in radians
+     * e: ellipsoid eccentricity
+     * returns: conformal latitude, in radians
+     ***********************************/
+    if (e == 0.0)
+        return phi;
+
+    const double esinphi = e * sin(phi);
+    const double chi =
+        2 * atan(tan(M_FORTPI + phi / 2) *
+                 std::pow(((1 - esinphi) / (1 + esinphi)), (e / 2))) -
+        M_HALFPI;
+    return chi;
+}
+
+/*****************************************************************************/
+double pj_conformal_lat_inverse(double chi, double e, double threshold) {
+    /***********************************
+     * The inverse formula of the conformal latitude
+     * for phi (geodetic latitude) in terms of chi (conformal latitude),
+     * Snyder, A working manual, formula 3-4
+     *
+     * chi: conformal latitude, in radians
+     * e: ellipsoid eccentricity
+     * threshold: between two consecutive iteratinons to break the loop, radians
+     * returns: geodetic latitude, in radians
+     ***********************************/
+    if (e == 0.0)
+        return chi;
+
+    const double taninit = tan(M_PI / 4 + chi / 2);
+    double phi = chi;
+    for (int i = 0; i < 10; i++) {
+        const double esinphi = e * sin(phi);
+        const double new_phi =
+            2 * atan(taninit *
+                     std::pow(((1 + esinphi) / (1 - esinphi)), (e / 2))) -
+            0.5 * M_PI;
+        if (abs(new_phi - phi) < threshold) {
+            phi = new_phi;
+            break;
+        }
+        phi = new_phi;
+    }
+    return phi;
+}

src/lib_proj.cmake

@@ -141,6 +141,7 @@ set(SRC_LIBPROJ_PROJECTIONS
   projections/calcofi.cpp
   projections/eqearth.cpp
   projections/col_urban.cpp
+  projections/spilhaus.cpp
 )
 
 set(SRC_LIBPROJ_CONVERSIONS
@@ -226,6 +227,7 @@ set(SRC_LIBPROJ_CORE
   initcache.cpp
   internal.cpp
   inv.cpp
+  latitudes.cpp
   list.cpp
   log.cpp
   malloc.cpp

src/pj_list.h

@@ -154,6 +154,7 @@ PROJ_HEAD(set, "Set coordinate value")
 PROJ_HEAD(sinu, "Sinusoidal (Sanson-Flamsteed)")
 PROJ_HEAD(som, "Space Oblique Mercator")
 PROJ_HEAD(somerc, "Swiss. Obl. Mercator")
+PROJ_HEAD(spilhaus, "Spilhaus")
 PROJ_HEAD(stere, "Stereographic")
 PROJ_HEAD(sterea, "Oblique Stereographic Alternative")
 PROJ_HEAD(gstmerc,

src/proj_internal.h

@@ -922,6 +922,9 @@ double pj_sinhpsi2tanphi(PJ_CONTEXT *, const double, const double);
 double *pj_authset(double);
 double pj_authlat(double, double *);
 
+double pj_conformal_lat(double phi, double e);
+double pj_conformal_lat_inverse(double chi, double e, double threshold);
+
 COMPLEX pj_zpoly1(COMPLEX, const COMPLEX *, int);
 COMPLEX pj_zpolyd1(COMPLEX, const COMPLEX *, int, COMPLEX *);
 

src/projections/spilhaus.cpp

@@ -0,0 +1,162 @@
+/*
+ * Implementation of the Spilhaus projections based on Adams World in a Square
+ * II.
+ *
+ * Explained in https://github.com/OSGeo/PROJ/issues/1851
+ *
+ * Copyright (c) 2025 Javier Jimenez Shaw
+ *
+ * Related material
+ * ----------------
+ *
+ * Map Projections - A Working Manual. 1987. John P. Snyder
+ * Sections 3 and 5.
+ * https://doi.org/10.3133/pp1395
+ * Online at https://neacsu.net/docs/geodesy/snyder/2-general/
+ */
+
+#include <errno.h>
+#include <math.h>
+
+#include "proj.h"
+#include "proj_internal.h"
+
+C_NAMESPACE PJ *pj_adams_ws2(PJ *);
+
+PROJ_HEAD(spilhaus, "Spilhaus") "\n\tSph&Ell";
+
+namespace { // anonymous namespace
+struct pj_spilhaus_data {
+    double cosalpha;
+    double sinalpha;
+    double beta;
+    double lambda_0;
+    double conformal_distortion;
+    double cosrot;
+    double sinrot;
+    PJ *adams_ws2;
+};
+
+} // anonymous namespace
+
+static PJ_XY spilhaus_forward(PJ_LP lp, PJ *P) {
+    PJ_XY xy = {0.0, 0.0};
+    struct pj_spilhaus_data *Q =
+        static_cast<struct pj_spilhaus_data *>(P->opaque);
+
+    const double phi_c = pj_conformal_lat(lp.phi, P->e);
+    const double cosphi_c = cos(phi_c);
+    const double sinphi_c = sin(phi_c);
+
+    const double coslam = cos(lp.lam - Q->lambda_0);
+    const double sinlam = sin(lp.lam - Q->lambda_0);
+
+    PJ_LP lpadams;
+    // Snyder, A working manual, formula 5-7
+    lpadams.phi =
+        aasin(P->ctx, Q->sinalpha * sinphi_c - Q->cosalpha * cosphi_c * coslam);
+
+    // Snyder, A working manual, formula 5-8b
+    lpadams.lam =
+        (Q->beta + atan2(cosphi_c * sinlam, (Q->sinalpha * cosphi_c * coslam +
+                                             Q->cosalpha * sinphi_c)));
+
+    while (lpadams.lam > M_PI)
+        lpadams.lam -= M_PI * 2;
+    while (lpadams.lam < -M_PI)
+        lpadams.lam += M_PI * 2;
+
+    PJ_XY xyadams = Q->adams_ws2->fwd(lpadams, Q->adams_ws2);
+    const double factor = Q->conformal_distortion * P->k0;
+    xy.x = -(xyadams.x * Q->cosrot + xyadams.y * Q->sinrot) * factor;
+    xy.y = -(xyadams.x * -Q->sinrot + xyadams.y * Q->cosrot) * factor;
+
+    return xy;
+}
+
+static PJ_LP spilhaus_inverse(PJ_XY xy, PJ *P) {
+    PJ_LP lp = {0.0, 0.0};
+    struct pj_spilhaus_data *Q =
+        static_cast<struct pj_spilhaus_data *>(P->opaque);
+
+    PJ_XY xyadams;
+    const double factor = 1.0 / (Q->conformal_distortion * P->k0);
+    xyadams.x = -(xy.x * Q->cosrot + xy.y * -Q->sinrot) * factor;
+    xyadams.y = -(xy.x * Q->sinrot + xy.y * Q->cosrot) * factor;
+    PJ_LP lpadams = Q->adams_ws2->inv(xyadams, Q->adams_ws2);
+
+    const double cosphi_s = cos(lpadams.phi);
+    const double sinphi_s = sin(lpadams.phi);
+    const double coslam_s = cos(lpadams.lam - Q->beta);
+    const double sinlam_s = sin(lpadams.lam - Q->beta);
+
+    // conformal latitude
+    lp.phi = aasin(P->ctx,
+                   Q->sinalpha * sinphi_s + Q->cosalpha * cosphi_s * coslam_s);
+
+    lp.lam = Q->lambda_0 +
+             aatan2(cosphi_s * sinlam_s,
+                    Q->sinalpha * cosphi_s * coslam_s - Q->cosalpha * sinphi_s);
+
+    const double THRESHOLD = 1e-10;
+    lp.phi = pj_conformal_lat_inverse(lp.phi, P->e, THRESHOLD);
+
+    return lp;
+}
+
+static PJ *spilhaus_destructor(PJ *P, int errlev) { /* Destructor */
+    if (nullptr == P)
+        return nullptr;
+    if (nullptr == P->opaque)
+        return pj_default_destructor(P, errlev);
+    proj_destroy(static_cast<struct pj_spilhaus_data *>(P->opaque)->adams_ws2);
+    return pj_default_destructor(P, errlev);
+}
+
+PJ *PJ_PROJECTION(spilhaus) {
+    struct pj_spilhaus_data *Q = static_cast<struct pj_spilhaus_data *>(
+        calloc(1, sizeof(struct pj_spilhaus_data)));
+    if (nullptr == Q)
+        return pj_default_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/);
+    P->opaque = Q;
+    P->destructor = spilhaus_destructor;
+
+    Q->adams_ws2 = pj_adams_ws2(nullptr);
+    if (Q->adams_ws2 == nullptr)
+        return spilhaus_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/);
+    Q->adams_ws2->ctx = P->ctx;
+    Q->adams_ws2->e = 0;
+    Q->adams_ws2 = pj_adams_ws2(Q->adams_ws2);
+    if (Q->adams_ws2 == nullptr)
+        return spilhaus_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/);
+
+    auto param_rad = [&P](const std::string &name, double def) {
+        return pj_param(P->ctx, P->params, ("t" + name).c_str()).i
+                   ? pj_param(P->ctx, P->params, ("r" + name).c_str()).f
+                   : def * DEG_TO_RAD;
+    };
+    // Values from https://github.com/OSGeo/PROJ/issues/1851
+    const double lat_center = param_rad("lat_c", -49.56371678);
+    const double lon_center = param_rad("lon_c", 66.94970198);
+    const double azimuth = param_rad("azi", 40.17823482);
+
+    const double rotation = param_rad("rot", 45);
+    Q->cosrot = cos(rotation);
+    Q->sinrot = sin(rotation);
+
+    const double conformal_lat_center = pj_conformal_lat(lat_center, P->e);
+
+    Q->sinalpha = -cos(conformal_lat_center) * cos(azimuth);
+    Q->cosalpha = sqrt(1 - Q->sinalpha * Q->sinalpha);
+    Q->lambda_0 = lon_center + atan2(tan(azimuth), -sin(conformal_lat_center));
+    Q->beta = M_PI + atan2(-sin(azimuth), -tan(conformal_lat_center));
+
+    Q->conformal_distortion =
+        cos(lat_center) / sqrt(1 - P->es * sin(lat_center) * sin(lat_center)) /
+        cos(conformal_lat_center);
+
+    P->fwd = spilhaus_forward;
+    P->inv = spilhaus_inverse;
+
+    return P;
+}