diff --git a/docs/plot/plotdefs.json b/docs/plot/plotdefs.json index da0a85d041..0189fada05 100644 --- a/docs/plot/plotdefs.json +++ b/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, diff --git a/docs/source/operations/projections/all_images.rst b/docs/source/operations/projections/all_images.rst index d7a5eaa03b..aae9b9bdc3 100644 --- a/docs/source/operations/projections/all_images.rst +++ b/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 diff --git a/docs/source/operations/projections/images/interrupted.png b/docs/source/operations/projections/images/interrupted.png new file mode 100644 index 0000000000..036897231a Binary files /dev/null and b/docs/source/operations/projections/images/interrupted.png differ diff --git a/docs/source/operations/projections/index.rst b/docs/source/operations/projections/index.rst index d206ef9123..e853f96b62 100644 --- a/docs/source/operations/projections/index.rst +++ b/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 diff --git a/docs/source/operations/projections/interrupted.rst b/docs/source/operations/projections/interrupted.rst new file mode 100644 index 0000000000..2c2c12c85c --- /dev/null +++ b/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= + + 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= + + :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 diff --git a/docs/source/references.bib b/docs/source/references.bib index ece00af2e5..b00817a904 100644 --- a/docs/source/references.bib +++ b/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}, diff --git a/docs/source/spelling_wordlist.txt b/docs/source/spelling_wordlist.txt index 3cfc298c3d..3976f9377f 100644 --- a/docs/source/spelling_wordlist.txt +++ b/docs/source/spelling_wordlist.txt @@ -404,6 +404,7 @@ gg gie gimbal Gitter +globemaking gn gnom gnomonic diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp index b73543e1c8..1376182294 100644 --- a/src/iso19111/io.cpp +++ b/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( diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake index cfde386a61..0e9261642e 100644 --- a/src/lib_proj.cmake +++ b/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 diff --git a/src/pj_list.h b/src/pj_list.h index b813f9899e..78dde1960c 100644 --- a/src/pj_list.h +++ b/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") diff --git a/src/projections/interrupted.cpp b/src/projections/interrupted.cpp new file mode 100644 index 0000000000..26e3429351 --- /dev/null +++ b/src/projections/interrupted.cpp @@ -0,0 +1,203 @@ +/* + * Implementation of meta interrupted projection for some pseudocylindrical + * + * Copyright (c) 2025 Javier Jimenez Shaw + * + */ + +#include +#include + +#include +#include + +#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 lam_borders; + std::vector lam_centers; + std::vector x_borders; + std::vector 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(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(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(P->opaque)->base_pj); + + return pj_default_destructor(P, errlev); +} + +PJ *PJ_PROJECTION(interrupted) { + struct pj_interrupted_data *Q = static_cast( + 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 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 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(static_cast(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 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; +} diff --git a/test/gie/interrupted.gie b/test/gie/interrupted.gie new file mode 100644 index 0000000000..03c17b895f --- /dev/null +++ b/test/gie/interrupted.gie @@ -0,0 +1,49 @@ + + +------------------------------------------------------------ +# Interrupted over sinu, tmerc, cass, and poly +------------------------------------------------------------ +operation +proj=interrupted +base=sinu +tolerance 1 mm + +accept 0 0 +expect 0 0 +roundtrip 1 + +accept 120 30 +expect 13358338.8952 3320113.3978 +roundtrip 1 + +operation +proj=interrupted +base=sinu +gores=1,2,3 +tolerance 1 mm + +accept 0 0 +expect 0 0 +roundtrip 1 + +accept 120 30 +expect 12913342.5789 3320113.3978 +roundtrip 1 + +operation +proj=interrupted +base=tmerc +gores=5 +tolerance 1 mm + +accept -120 30 +expect -14862738.2699 3575648.8083 +roundtrip 1 + +operation +proj=interrupted +base=cass +gores=5 +tolerance 1 mm + +accept -120 30 +expect -13732123.8999 3575100.4520 +roundtrip 1 + +operation +proj=interrupted +base=poly +gores=5 +tolerance 1 mm + +accept -120 30 +expect -13731228.3039 3561724.7405 +roundtrip 1 + +