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_0=<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_0=<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 *);