Add interrupted projection for "sinu", "tmerc", "cass" and "poly"

docs/plot/plotdefs.json

@@ -635,6 +635,18 @@
         "res": "low",
         "type": "line"
     },
+    {
+        "filename": "interrupted.png",
+        "latmax": 90,
+        "latmin": -90,
+        "lonmax": 180,
+        "lonmin": -180,
+        "name": "interrupted",
+        "projstring": "+proj=interrupted +base=sinu",
+        "res": "low",
+        "type": "line",
+        "delta_cut": 1e6
+    },
     {
         "filename": "isea.png",
         "latmax": 90,

docs/source/operations/projections/all_images.rst

@@ -644,6 +644,17 @@ List of all projection images
 ********************************************************************************
 
 
+:ref:`interrupted`
+
+.. figure:: ./images/interrupted.png
+   :width: 500 px
+   :align: center
+   :alt:   interrupted
+
+
+********************************************************************************
+
+
 :ref:`isea`
 
 .. figure:: ./images/isea.png

docs/source/operations/projections/index.rst

@@ -70,6 +70,7 @@ Projections map the spherical 3D space to a flat 2D space.
    imoll
    imoll_o
    imw_p
+   interrupted
    isea
    kav5
    kav7

docs/source/operations/projections/interrupted.rst

@@ -0,0 +1,61 @@
+.. _interrupted:
+
+********************************************************************************
+Interrupted
+********************************************************************************
+
++---------------------+----------------------------------------------------------+
+| **Classification**  | Interrupted meta-projection                              |
++---------------------+----------------------------------------------------------+
+| **Available forms** | Forward and inverse, spherical and ellipsoidal           |
++---------------------+----------------------------------------------------------+
+| **Defined area**    | Global                                                   |
++---------------------+----------------------------------------------------------+
+| **Alias**           | interrupted                                              |
++---------------------+----------------------------------------------------------+
+| **Domain**          | 2D                                                       |
++---------------------+----------------------------------------------------------+
+| **Input type**      | Geodetic coordinates                                     |
++---------------------+----------------------------------------------------------+
+| **Output type**     | Projected coordinates                                    |
++---------------------+----------------------------------------------------------+
+
+
+.. figure:: ./images/interrupted.png
+   :width: 500 px
+   :align: center
+   :alt:   Interrupted
+
+   proj-string: ``+proj=interrupted +base=sinu``
+
+Meta-projection to create interrupted maps using different base projections.
+The supported projections are :ref:`cass`, :ref:`poly`, :ref:`sinu` and :ref:`tmerc`.
+The world is divided in gores tangent at the equator. :cite:`Huck2024`
+
+Parameters
+################################################################################
+
+.. note:: All parameters are optional for the Interrupted projection except ``base``.
+
+.. option:: +base=<value>
+
+    Base projection. It can be one of ``cass``, ``poly``, ``sinu`` or ``tmerc``.
+    Note that some projections may not behave correctly with too wide gores.
+
+.. option:: +gores=<value>
+
+    :option:`+gores` can be a number or a list of comma separated values.
+    When it is a single number, it is the number of equal gores.
+    When it is a list, they are the relative sizes of the gores (v.g. ``1,2,1``).
+
+    *Defaults to 3.*
+
+.. include:: ../options/lon_0.rst
+
+.. include:: ../options/ellps.rst
+
+.. include:: ../options/R.rst
+
+.. include:: ../options/x_0.rst
+
+.. include:: ../options/y_0.rst

docs/source/references.bib

@@ -191,6 +191,18 @@ @InProceedings{Hensley2002
   Url                      = {https://airsar.jpl.nasa.gov/documents/workshop2002/papers/T3.pdf}
 }
 
+@Article{Huck2024,
+  author                   = {Jonathan Huck},
+  title                    = {Evaluating Map Projections for Globemaking},
+  journal                  = {The Cartographic Journal},
+  volume                   = {61},
+  number                   = {3},
+  pages                    = {207--217},
+  year                     = {2024},
+  publisher                = {Taylor \& Francis},
+  doi                      = {10.1080/00087041.2024.2436324},
+}
+
 @TechReport{IOGP2018,
   Title                    = {Geomatics Guidance Note 7, part 2: Coordinate Conversions \& Transformations including Formulas},
   Author                   = {IOGP},

docs/source/spelling_wordlist.txt

@@ -404,6 +404,7 @@ gg
 gie
 gimbal
 Gitter
+globemaking
 gn
 gnom
 gnomonic

src/iso19111/io.cpp

@@ -12313,7 +12313,9 @@ PROJStringParser::Private::buildProjectedCRS(int iStep,
             }
             if (param.value.empty()) {
                 methodName += " " + param.key;
-            } else if (isalpha(param.value[0])) {
+            } else if (isalpha(param.value[0]) || param.key == "gores") {
+                // The value of gores can be a string or a number.
+                // See interrupted.cpp for more info.
                 methodName += " " + param.key + "=" + param.value;
             } else {
                 parameters.push_back(OperationParameter::create(

src/lib_proj.cmake

@@ -134,6 +134,7 @@ set(SRC_LIBPROJ_PROJECTIONS
   projections/eqearth.cpp
   projections/col_urban.cpp
   projections/spilhaus.cpp
+  projections/interrupted.cpp
 )
 
 set(SRC_LIBPROJ_CONVERSIONS

src/pj_list.h

@@ -194,3 +194,4 @@ PROJ_HEAD(wink1, "Winkel I")
 PROJ_HEAD(wink2, "Winkel II")
 PROJ_HEAD(wintri, "Winkel Tripel")
 PROJ_HEAD(xyzgridshift, "XYZ grid shift")
+PROJ_HEAD(interrupted, "Interrupted")

src/projections/interrupted.cpp

@@ -0,0 +1,203 @@
+/*
+ * Implementation of meta interrupted projection for some pseudocylindrical
+ *
+ * Copyright (c) 2025 Javier Jimenez Shaw
+ *
+ */
+
+#include <errno.h>
+#include <math.h>
+
+#include <iostream>
+#include <numeric>
+
+#include "proj.h"
+#include "proj_internal.h"
+
+#define FROM_PROJ_CPP
+#include "proj/internal/internal.hpp"
+
+C_NAMESPACE PJ *pj_cass(PJ *);
+C_NAMESPACE PJ *pj_poly(PJ *);
+C_NAMESPACE PJ *pj_sinu(PJ *);
+C_NAMESPACE PJ *pj_tmerc(PJ *);
+
+PROJ_HEAD(interrupted, "Interrupted") "\n\tSph&Ell";
+
+namespace { // anonymous namespace
+struct pj_interrupted_data {
+    PJ *base_pj;
+    std::vector<double> lam_borders;
+    std::vector<double> lam_centers;
+    std::vector<double> x_borders;
+    std::vector<double> x_centers;
+};
+} // anonymous namespace
+
+static PJ_XY interrupted_forward(PJ_LP lp, PJ *P) {
+    PJ_XY xy = {0.0, 0.0};
+    struct pj_interrupted_data *Q =
+        static_cast<struct pj_interrupted_data *>(P->opaque);
+    double center_lam = 0;
+    double center_x = 0;
+    for (size_t i = 0; i < Q->lam_borders.size(); i++) {
+        if (lp.lam <= Q->lam_borders[i]) {
+            center_lam = Q->lam_centers[i];
+            center_x = Q->x_centers[i];
+            break;
+        }
+    }
+    lp.lam -= center_lam;
+    xy = Q->base_pj->fwd(lp, Q->base_pj);
+    xy.x += center_x;
+    return xy;
+}
+
+static PJ_LP interrupted_inverse(PJ_XY xy, PJ *P) {
+    PJ_LP lp = {0.0, 0.0};
+
+    struct pj_interrupted_data *Q =
+        static_cast<struct pj_interrupted_data *>(P->opaque);
+
+    double center_lam = 0;
+    double center_x = 0;
+    double start_lam = 0;
+    double end_lam = 0;
+    for (size_t i = 0; i < Q->x_borders.size(); i++) {
+        if (xy.x <= Q->x_borders[i]) {
+            center_lam = Q->lam_centers[i];
+            center_x = Q->x_centers[i];
+            end_lam = (i + 1 != Q->x_borders.size()) ? Q->lam_borders[i] : M_PI;
+            start_lam = i > 0 ? Q->lam_borders[i - 1] : -M_PI;
+            break;
+        }
+    }
+    xy.x -= center_x;
+    lp = Q->base_pj->inv(xy, Q->base_pj);
+    lp.lam += center_lam;
+    if (lp.lam > end_lam || lp.lam < start_lam) {
+        proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
+        return lp;
+    }
+
+    return lp;
+}
+
+static PJ *interrupted_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_interrupted_data *>(P->opaque)->base_pj);
+
+    return pj_default_destructor(P, errlev);
+}
+
+PJ *PJ_PROJECTION(interrupted) {
+    struct pj_interrupted_data *Q = static_cast<struct pj_interrupted_data *>(
+        calloc(1, sizeof(struct pj_interrupted_data)));
+    if (nullptr == Q)
+        return pj_default_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/);
+    P->opaque = Q;
+    P->destructor = interrupted_destructor;
+
+    std::string base;
+    if (pj_param(P->ctx, P->params, "tbase").i) {
+        base = pj_param(P->ctx, P->params, "sbase").s;
+    }
+    if (base != "cass" && base != "poly" && base != "sinu" && base != "tmerc") {
+        proj_log_error(P, _("Invalid value for base parameter"));
+        return pj_default_destructor(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
+    }
+    std::vector<double> sizes;
+    if (pj_param(P->ctx, P->params, "tgores").i) {
+        const std::string g = pj_param(P->ctx, P->params, "sgores").s;
+        const std::vector<std::string> gores =
+            osgeo::proj::internal::split(g, ",");
+        bool success = true;
+        for (const auto &gore : gores) {
+            if (gore.empty())
+                continue;
+            double d = osgeo::proj::internal::c_locale_stod(gore, success);
+            if (!success) {
+                proj_log_error(P, _("Invalid value for gores parameter"));
+                return pj_default_destructor(
+                    P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
+            }
+            sizes.push_back(d);
+        }
+    } else {
+        sizes.push_back(3.0);
+    }
+
+    if (sizes.size() == 1) {
+        sizes = std::vector<double>(static_cast<size_t>(sizes.front()), 1.0);
+    }
+    if (sizes.empty()) {
+        proj_log_error(P, _("Invalid 0 gores"));
+        return pj_default_destructor(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
+    }
+
+    auto get_fun_ptr = [&base]() {
+        if (base == "cass") {
+            return pj_cass;
+        } else if (base == "poly") {
+            return pj_poly;
+        } else if (base == "sinu") {
+            return pj_sinu;
+        } else if (base == "tmerc") {
+            return pj_tmerc;
+        }
+        return pj_sinu;
+    };
+
+    Q->base_pj = get_fun_ptr()(nullptr);
+    if (Q->base_pj == nullptr)
+        return interrupted_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/);
+    Q->base_pj->ctx = P->ctx;
+    Q->base_pj->es = P->es;
+    Q->base_pj->one_es = P->one_es;
+    Q->base_pj->k0 = P->k0;
+    Q->base_pj->f = P->f;
+    Q->base_pj->f2 = P->f2;
+    Q->base_pj->n = P->n;
+    Q->base_pj = get_fun_ptr()(Q->base_pj);
+    if (Q->base_pj == nullptr)
+        return interrupted_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/);
+
+    const double total_size = std::accumulate(sizes.begin(), sizes.end(), 0.0);
+
+    // Compute center and east border of every gore in lambda (radians)
+    double lam_start = -M_PI;
+    std::vector<double> x_semi_widths;
+    for (auto &s : sizes) {
+        const double lam_semi_width = s * M_PI / total_size;
+        const auto x_swd = Q->base_pj->fwd({lam_semi_width, 0}, Q->base_pj);
+        x_semi_widths.push_back(x_swd.x);
+        const double lam_center = lam_start + lam_semi_width;
+        Q->lam_centers.push_back(lam_center);
+        const double lam_end = lam_start + lam_semi_width * 2;
+        Q->lam_borders.push_back(lam_end);
+        lam_start = lam_end;
+    }
+
+    // Compute center and east border of every gore in x
+    double x_start =
+        -std::accumulate(x_semi_widths.begin(), x_semi_widths.end(), 0.0);
+    for (auto &x_semi_width : x_semi_widths) {
+        const double x_center = x_start + x_semi_width;
+        Q->x_centers.push_back(x_center);
+        const double x_end = x_start + x_semi_width * 2.0;
+        Q->x_borders.push_back(x_end);
+        x_start = x_end;
+    }
+
+    // Add some marging to the last gore
+    Q->lam_borders.back() += 1.;
+    Q->x_borders.back() += 100.;
+
+    P->fwd = interrupted_forward;
+    P->inv = interrupted_inverse;
+
+    return P;
+}