diff --git a/docs/source/operations/projections/index.rst b/docs/source/operations/projections/index.rst
index 4842097140..4f76dab2f2 100644
--- a/docs/source/operations/projections/index.rst
+++ b/docs/source/operations/projections/index.rst
@@ -126,6 +126,7 @@ Projections map the spherical 3D space to a flat 2D space.
robin
rouss
rpoly
+ s2
sch
sinu
somerc
diff --git a/docs/source/operations/projections/s2.rst b/docs/source/operations/projections/s2.rst
new file mode 100644
index 0000000000..46ae715e81
--- /dev/null
+++ b/docs/source/operations/projections/s2.rst
@@ -0,0 +1,120 @@
+.. _s2:
+
+********************************************************************************
+S2
+********************************************************************************
+
++---------------------+----------------------------------------------------------+
+| **Classification** | Miscellaneous |
++---------------------+----------------------------------------------------------+
+| **Available forms** | Forward and inverse, ellipsoidal |
++---------------------+----------------------------------------------------------+
+| **Defined area** | Global |
++---------------------+----------------------------------------------------------+
+| **Alias** | s2 |
++---------------------+----------------------------------------------------------+
+| **Domain** | 2D |
++---------------------+----------------------------------------------------------+
+| **Input type** | Geodetic coordinates |
++---------------------+----------------------------------------------------------+
+| **Output type** | Projected coordinates |
++---------------------+----------------------------------------------------------+
+
+.. versionadded:: 8.2
+
+The S2 projection, like the Quadrilateralized Spherical Cube (QSC) projection, projects
+a sphere surface onto the six sides of a cube:
+
+.. image:: ../../../images/qsc_concept.jpg
+ :width: 500 px
+ :align: center
+ :alt: Quadrilateralized Spherical Cube
+
+S2 was created by Google to represent geographic data on the whole earth. The documentation can be found
+at `S2 Geometry `_ It works by first
+projecting a point on the sphere to a face of the cube. These are called u,v-coordinates, and they are in [-1,1] x [-1,1].
+This step is followed by a non-linear transformation to normalize the area of rectangles on the sphere. There are three
+different choices available for this transformation, meaning that S2 is a family of projections. The final output is in
+s,t-coordinates, which are in [0,1] x [0,1].
+See the comments in `S2 Code `_
+for an explanation of the tradeoff between speed and area-preservation. Note that the projection is azimuthal when "none" or
+"linear" is selected for the area-normalization, but it is not azimuthal when "quadratic" or "tangent" is chosen. See
+S2's `Earthcube page `_
+to visualize the unfolded cube and the orientation of each face.
+
+In this implementation, the cube side is selected by choosing one of the following six projection centers:
+
++-------------------------+--------------------+
+| ``+lat_0=0 +lon_0=0`` | front cube side |
++-------------------------+--------------------+
+| ``+lat_0=0 +lon_0=90`` | right cube side |
++-------------------------+--------------------+
+| ``+lat_0=0 +lon_0=180`` | back cube side |
++-------------------------+--------------------+
+| ``+lat_0=0 +lon_0=-90`` | left cube side |
++-------------------------+--------------------+
+| ``+lat_0=90`` | top cube side |
++-------------------------+--------------------+
+| ``+lat_0=-90`` | bottom cube side |
++-------------------------+--------------------+
+
+The specific transformation can be chosen with the UVtoST parameter:
+
++-------------------------+-----------------------------+
+| ``+UVtoST=linear`` | fastest, no normalization |
++-------------------------+-----------------------------+
+| ``+UVtoST=quadratic`` | fast, good normalization |
++-------------------------+-----------------------------+
+| ``+UVtoST=tangent`` | slowest, best normalization |
++-------------------------+-----------------------------+
+| ``+UVtoST=none`` | returns u,v-coordinates |
++-------------------------+-----------------------------+
+
+Furthermore, this implementation allows the projection to be applied to ellipsoids.
+A preceding shift to a sphere is performed automatically; see :cite:`LambersKolb2012` for details.
+The output of the projection is in s,t-coordinates ([0,1] x [0,1]), so only the
+eccentricity of the ellipse is taken into account: the absolute value of the axes does
+not affect the output.
+
+
+Usage
+###############################################################################
+
+The following example uses S2 on the right face::
+
+ echo 90 0 | ../bin/proj +proj=s2 +lat_0=0 +lon_0=90 +ellps=WGS84 +UVtoST=linear
+
+ 0.5 0.5
+
+Explanation:
+
+* S2 projection is selected with ``+proj=s2``.
+* The WGS84 ellipsoid is specified with ``+ellps=WGS84``.
+* The cube side is selected with ``+lat_0=... +lon_0=...``.
+* The normalization transformation is selected with ``+UVtoST=...``.
+
+Parameters
+################################################################################
+
+.. note:: All parameters for the projection are optional.
+
+.. include:: ../options/lon_0.rst
+
+.. include:: ../options/lat_0.rst
+
+.. include:: ../options/ellps.rst
+
+.. option:: +UVtoST=
+
+ The area-normalization transformation. Choose from {linear, quadratic, tangent, none}
+
+ *Defaults to "quadratic".*
+
+.. include:: ../options/x_0.rst
+
+.. include:: ../options/y_0.rst
+
+Further reading
+################################################################################
+
+#. `S2's Website `_
diff --git a/src/Makefile.am b/src/Makefile.am
index f4499c7340..bae922d58e 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -181,6 +181,7 @@ libproj_la_SOURCES = \
projections/putp6.cpp \
projections/qsc.cpp \
projections/robin.cpp \
+ projections/s2.cpp \
projections/sch.cpp \
projections/sts.cpp \
projections/urm5.cpp \
diff --git a/src/bin_projsync.cmake b/src/bin_projsync.cmake
index c184a4227e..6486b5d447 100644
--- a/src/bin_projsync.cmake
+++ b/src/bin_projsync.cmake
@@ -8,10 +8,6 @@ set_target_properties(bin_projsync
OUTPUT_NAME projsync)
target_link_libraries(bin_projsync PRIVATE ${PROJ_LIBRARIES})
target_compile_options(bin_projsync PRIVATE ${PROJ_CXX_WARN_FLAGS})
-if(NLOHMANN_JSON STREQUAL "external")
- target_compile_definitions(bin_projsync PRIVATE EXTERNAL_NLOHMANN_JSON)
- target_link_libraries(bin_projsync PRIVATE nlohmann_json::nlohmann_json)
-endif()
install(TARGETS bin_projsync
DESTINATION ${BINDIR})
diff --git a/src/init.cpp b/src/init.cpp
index a0d727cbd6..61457cb6c5 100644
--- a/src/init.cpp
+++ b/src/init.cpp
@@ -41,7 +41,6 @@
#include "filemanager.hpp"
#include
-
/**************************************************************************************/
static paralist *string_to_paralist (PJ_CONTEXT *ctx, char *definition) {
/***************************************************************************************
diff --git a/src/iso19111/operation/parammappings.cpp b/src/iso19111/operation/parammappings.cpp
index 5999d535a8..ebaa4ba72d 100644
--- a/src/iso19111/operation/parammappings.cpp
+++ b/src/iso19111/operation/parammappings.cpp
@@ -829,7 +829,7 @@ static const MethodMapping projectionMethodMappings[] = {
{PROJ_WKT2_NAME_METHOD_QUADRILATERALIZED_SPHERICAL_CUBE, 0,
"Quadrilateralized_Spherical_Cube", "qsc", nullptr, paramsNatOrigin},
-
+
{PROJ_WKT2_NAME_METHOD_SPHERICAL_CROSS_TRACK_HEIGHT, 0,
"Spherical_Cross_Track_Height", "sch", nullptr, paramsSch},
diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake
index 71468c8768..3d7644400f 100644
--- a/src/lib_proj.cmake
+++ b/src/lib_proj.cmake
@@ -138,6 +138,7 @@ set(SRC_LIBPROJ_PROJECTIONS
projections/putp6.cpp
projections/qsc.cpp
projections/robin.cpp
+ projections/s2.cpp
projections/sch.cpp
projections/sts.cpp
projections/urm5.cpp
@@ -398,7 +399,8 @@ target_link_libraries(proj PRIVATE ${SQLITE3_LIBRARY})
if(NLOHMANN_JSON STREQUAL "external")
target_compile_definitions(proj PRIVATE EXTERNAL_NLOHMANN_JSON)
- target_link_libraries(proj PRIVATE nlohmann_json::nlohmann_json)
+ target_link_libraries(proj
+ PRIVATE $)
endif()
if(TIFF_ENABLED)
diff --git a/src/pj_list.h b/src/pj_list.h
index 226a8f7b59..e4b8ab67d8 100644
--- a/src/pj_list.h
+++ b/src/pj_list.h
@@ -141,6 +141,7 @@ PROJ_HEAD(qsc, "Quadrilateralized Spherical Cube")
PROJ_HEAD(robin, "Robinson")
PROJ_HEAD(rouss, "Roussilhe Stereographic")
PROJ_HEAD(rpoly, "Rectangular Polyconic")
+PROJ_HEAD(s2, "S2")
PROJ_HEAD(sch, "Spherical Cross-track Height")
PROJ_HEAD(set, "Set coordinate value")
PROJ_HEAD(sinu, "Sinusoidal (Sanson-Flamsteed)")
diff --git a/src/proj_constants.h b/src/proj_constants.h
index 66d0429457..7c6d018a61 100644
--- a/src/proj_constants.h
+++ b/src/proj_constants.h
@@ -223,6 +223,9 @@
#define PROJ_WKT2_NAME_METHOD_QUADRILATERALIZED_SPHERICAL_CUBE \
"Quadrilateralized Spherical Cube"
+#define PROJ_WKT2_NAME_METHOD_S2 \
+ "S2"
+
#define PROJ_WKT2_NAME_METHOD_SPHERICAL_CROSS_TRACK_HEIGHT \
"Spherical Cross-Track Height"
diff --git a/src/projections/s2.cpp b/src/projections/s2.cpp
new file mode 100644
index 0000000000..0de54c004e
--- /dev/null
+++ b/src/projections/s2.cpp
@@ -0,0 +1,394 @@
+/******************************************************************************
+ * Project: PROJ
+ * Purpose: Implementing the S2 family of projections in PROJ
+ * Author: Marcus Elia,
+ *
+ ******************************************************************************
+ * Copyright (c) 2021, Marcus Elia,
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ *****************************************************************************
+ *
+ * This implements the S2 projection. This code is similar
+ * to the QSC projection.
+ *
+ *
+ * You have to choose one of the following projection centers,
+ * corresponding to the centers of the six cube faces:
+ * phi0 = 0.0, lam0 = 0.0 ("front" face)
+ * phi0 = 0.0, lam0 = 90.0 ("right" face)
+ * phi0 = 0.0, lam0 = 180.0 ("back" face)
+ * phi0 = 0.0, lam0 = -90.0 ("left" face)
+ * phi0 = 90.0 ("top" face)
+ * phi0 = -90.0 ("bottom" face)
+ * Other projection centers will not work!
+ *
+ * You must also choose which conversion from UV to ST coordinates
+ * is used (linear, quadratic, tangent). Read about them in
+ * https://github.com/google/s2geometry/blob/0c4c460bdfe696da303641771f9def900b3e440f/src/s2/s2coords.h
+ * The S2 projection functions are adapted from the above link and the S2
+ * Math util functions are adapted from
+ * https://github.com/google/s2geometry/blob/0c4c460bdfe696da303641771f9def900b3e440f/src/s2/util/math/vector.h
+ ****************************************************************************/
+
+#define PJ_LIB__
+#define _USE_MATH_DEFINES // needed for M_1_PI availability with MSVC
+
+#include
+#include
+
+#include "proj.h"
+#include "proj_internal.h"
+
+/* The six cube faces. */
+namespace { // anonymous namespace
+enum Face {
+ FACE_FRONT = 0,
+ FACE_RIGHT = 1,
+ FACE_TOP = 2,
+ FACE_BACK = 3,
+ FACE_LEFT = 4,
+ FACE_BOTTOM = 5
+};
+} // anonymous namespace
+
+enum S2ProjectionType {Linear, Quadratic, Tangent, NoUVtoST};
+std::map stringToS2ProjectionType { {"linear", Linear}, {"quadratic", Quadratic}, {"tangent", Tangent}, {"none", NoUVtoST} };
+
+namespace { // anonymous namespace
+struct pj_opaque {
+ enum Face face;
+ double a_squared;
+ double one_minus_f;
+ double one_minus_f_squared;
+ S2ProjectionType UVtoST;
+};
+} // anonymous namespace
+PROJ_HEAD(s2, "S2") "\n\tMisc, Sph&Ell";
+
+#define EPS10 1.e-10
+
+/* The four areas on a cube face. AREA_0 is the area of definition,
+ * the other three areas are counted counterclockwise. */
+namespace { // anonymous namespace
+enum Area {
+ AREA_0 = 0,
+ AREA_1 = 1,
+ AREA_2 = 2,
+ AREA_3 = 3
+};
+} // anonymous namespace
+
+// =================================================
+//
+// S2 Math Util
+//
+// =================================================
+static PJ_XYZ Abs(const PJ_XYZ& p) {
+ return {std::fabs(p.x), std::fabs(p.y), std::fabs(p.z)};
+}
+// return the index of the largest component (fabs)
+static int LargestAbsComponent(const PJ_XYZ& p) {
+ PJ_XYZ temp = Abs(p);
+ return temp.x > temp.y ?
+ temp.x > temp.z ? 0 : 2 :
+ temp.y > temp.z ? 1 : 2;
+}
+
+// =================================================
+//
+// S2 Projection Functions
+//
+// =================================================
+
+// Unfortunately, tan(M_PI_4) is slightly less than 1.0. This isn't due to
+// a flaw in the implementation of tan(), it's because the derivative of
+// tan(x) at x=pi/4 is 2, and it happens that the two adjacent floating
+// point numbers on either side of the infinite-precision value of pi/4 have
+// tangents that are slightly below and slightly above 1.0 when rounded to
+// the nearest double-precision result.
+static double STtoUV(double s, S2ProjectionType s2_projection) {
+ switch(s2_projection) {
+ case Linear:
+ return 2 * s - 1;
+ break;
+ case Quadratic:
+ if (s >= 0.5) return (1/3.) * (4*s*s - 1);
+ else return (1/3.) * (1 - 4*(1-s)*(1-s));
+ break;
+ case Tangent:
+ s = std::tan(M_PI_2 * s - M_PI_4);
+ return s + (1.0 / static_cast(static_cast(1) << 53)) * s;
+ break;
+ default:
+ return s;
+ }
+}
+
+static double UVtoST(double u, S2ProjectionType s2_projection) {
+ switch(s2_projection) {
+ case Linear:
+ return 0.5 * (u + 1);
+ break;
+ case Quadratic:
+ if (u >= 0) return 0.5 * std::sqrt(1 + 3*u);
+ else return 1 - 0.5 * std::sqrt(1 - 3*u);
+ break;
+ case Tangent:
+ {
+ volatile double a = std::atan(u);
+ return (2 * M_1_PI) * (a + M_PI_4);
+ }
+ break;
+ default:
+ return u;
+ }
+}
+
+inline PJ_XYZ FaceUVtoXYZ(int face, double u, double v) {
+ switch (face) {
+ case 0: return { 1, u, v};
+ case 1: return {-u, 1, v};
+ case 2: return {-u, -v, 1};
+ case 3: return {-1, -v, -u};
+ case 4: return { v, -1, -u};
+ default: return { v, u, -1};
+ }
+}
+
+inline PJ_XYZ FaceUVtoXYZ(int face, const PJ_XY& uv) {
+ return FaceUVtoXYZ(face, uv.x, uv.y);
+}
+
+inline void ValidFaceXYZtoUV(int face, const PJ_XYZ& p,
+ double* pu, double* pv) {
+ switch (face) {
+ case 0: *pu = p.y / p.x; *pv = p.z / p.x; break;
+ case 1: *pu = -p.x / p.y; *pv = p.z / p.y; break;
+ case 2: *pu = -p.x / p.z; *pv = -p.y / p.z; break;
+ case 3: *pu = p.z / p.x; *pv = p.y / p.x; break;
+ case 4: *pu = p.z / p.y; *pv = -p.x / p.y; break;
+ default: *pu = -p.y / p.z; *pv = -p.x / p.z; break;
+ }
+}
+
+inline void ValidFaceXYZtoUV(int face, const PJ_XYZ& p, PJ_XY* puv) {
+ ValidFaceXYZtoUV(face, p, &(*puv).x, &(*puv).y);
+}
+
+inline int GetFace(const PJ_XYZ& p) {
+ int face = LargestAbsComponent(p);
+ double pFace;
+ switch (face) {
+ case 0: pFace = p.x; break;
+ case 1: pFace = p.y; break;
+ default: pFace = p.z; break;
+ }
+ if (pFace < 0) face += 3;
+ return face;
+}
+
+inline int XYZtoFaceUV(const PJ_XYZ& p, double* pu, double* pv) {
+ int face = GetFace(p);
+ ValidFaceXYZtoUV(face, p, pu, pv);
+ return face;
+}
+
+inline int XYZtoFaceUV(const PJ_XYZ& p, PJ_XY* puv) {
+ return XYZtoFaceUV(p, &(*puv).x, &(*puv).y);
+}
+
+inline bool FaceXYZtoUV(int face, const PJ_XYZ& p,
+ double* pu, double* pv) {
+ double pFace;
+ switch(face) {
+ case 0: pFace = p.x; break;
+ case 1: pFace = p.y; break;
+ case 2: pFace = p.z; break;
+ case 3: pFace = p.x; break;
+ case 4: pFace = p.y; break;
+ default: pFace = p.z; break;
+ }
+ if (face < 3) {
+ if (pFace <= 0) return false;
+ } else {
+ if (pFace >= 0) return false;
+ }
+ ValidFaceXYZtoUV(face, p, pu, pv);
+ return true;
+}
+
+inline bool FaceXYZtoUV(int face, const PJ_XYZ& p, PJ_XY* puv) {
+ return FaceXYZtoUV(face, p, &(*puv).x, &(*puv).y);
+}
+
+// This function inverts ValidFaceXYZtoUV()
+inline bool UVtoSphereXYZ(int face, double u, double v, PJ_XYZ* xyz) {
+ double major_coord = 1 / sqrt(1 + u*u + v*v);
+ double minor_coord_1 = u*major_coord;
+ double minor_coord_2 = v*major_coord;
+
+ switch(face) {
+ case 0: xyz->x = major_coord;
+ xyz->y = minor_coord_1;
+ xyz->z = minor_coord_2; break;
+ case 1: xyz->x = -minor_coord_1;
+ xyz->y = major_coord;
+ xyz->z = minor_coord_2; break;
+ case 2: xyz->x = -minor_coord_1;
+ xyz->y = -minor_coord_2;
+ xyz->z = major_coord; break;
+ case 3: xyz->x = -major_coord;
+ xyz->y = -minor_coord_2;
+ xyz->z = -minor_coord_1; break;
+ case 4: xyz->x = minor_coord_2;
+ xyz->y = -major_coord;
+ xyz->z = -minor_coord_1; break;
+ default:xyz->x = minor_coord_2;
+ xyz->y = minor_coord_1;
+ xyz->z = -major_coord; break;
+ }
+
+ return true;
+}
+
+// ============================================
+//
+// The Forward and Inverse Functions
+//
+// ============================================
+static PJ_XY s2_forward (PJ_LP lp, PJ *P) {
+ struct pj_opaque *Q = static_cast(P->opaque);
+ double lat, lon;
+
+ /* Convert the geodetic latitude to a geocentric latitude.
+ * This corresponds to the shift from the ellipsoid to the sphere
+ * described in [LK12]. */
+ if (P->es != 0.0) {
+ lat = atan(Q->one_minus_f_squared * tan(lp.phi));
+ } else {
+ lat = lp.phi;
+ }
+ lon = lp.lam;
+
+ // Convert the lat/lon to x,y,z on the unit sphere
+ double x, y, z;
+ double sinlat, coslat;
+ double sinlon, coslon;
+
+ sinlat = sin(lat);
+ coslat = cos(lat);
+ sinlon = sin(lon);
+ coslon = cos(lon);
+ x = coslat * coslon;
+ y = coslat * sinlon;
+ z = sinlat;
+
+ PJ_XYZ spherePoint {x, y, z};
+ PJ_XY uvCoords;
+
+ ValidFaceXYZtoUV(Q->face, spherePoint, &uvCoords.x, &uvCoords.y);
+ double s = UVtoST(uvCoords.x, Q->UVtoST);
+ double t = UVtoST(uvCoords.y, Q->UVtoST);
+ return {s, t};
+}
+
+static PJ_LP s2_inverse (PJ_XY xy, PJ *P) {
+ PJ_LP lp = {0.0,0.0};
+ struct pj_opaque *Q = static_cast(P->opaque);
+
+ // Do the S2 projections to get from s,t to u,v to x,y,z
+ double u = STtoUV(xy.x, Q->UVtoST);
+ double v = STtoUV(xy.y, Q->UVtoST);
+
+ PJ_XYZ sphereCoords;
+ UVtoSphereXYZ(Q->face, u, v, &sphereCoords);
+ double q = sphereCoords.x;
+ double r = sphereCoords.y;
+ double s = sphereCoords.z;
+
+ // Get the spherical angles from the x y z
+ lp.phi = acos(-s) - M_HALFPI;
+ lp.lam = atan2(r, q);
+
+ /* Apply the shift from the sphere to the ellipsoid as described
+ * in [LK12]. */
+ if (P->es != 0.0) {
+ int invert_sign;
+ volatile double tanphi, xa;
+ invert_sign = (lp.phi < 0.0 ? 1 : 0);
+ tanphi = tan(lp.phi);
+ xa = P->b / sqrt(tanphi * tanphi + Q->one_minus_f_squared);
+ lp.phi = atan(sqrt(Q->a_squared - xa * xa) / (Q->one_minus_f * xa));
+ if (invert_sign) {
+ lp.phi = -lp.phi;
+ }
+ }
+
+ return lp;
+}
+
+PJ *PROJECTION(s2) {
+ struct pj_opaque *Q = static_cast(calloc (1, sizeof (struct pj_opaque)));
+ if (nullptr==Q)
+ return pj_default_destructor (P, PROJ_ERR_OTHER /*ENOMEM*/);
+ P->opaque = Q;
+
+ /* Determine which UVtoST function is to be used */
+ PROJVALUE maybeUVtoST = pj_param(P->ctx, P->params, "sUVtoST");
+ if (nullptr != maybeUVtoST.s) {
+ try {
+ Q->UVtoST = stringToS2ProjectionType.at(maybeUVtoST.s);
+ } catch (const std::out_of_range&) {
+ proj_log_error(P, _("Invalid value for s2 parameter: should be linear, quadratic, tangent, or none."));
+ return pj_default_destructor (P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
+ }
+ } else {
+ Q->UVtoST = Quadratic;
+ }
+
+ P->left = PJ_IO_UNITS_RADIANS;
+ P->right = PJ_IO_UNITS_PROJECTED;
+ P->from_greenwich = -P->lam0;
+
+ P->inv = s2_inverse;
+ P->fwd = s2_forward;
+
+ /* Determine the cube face from the center of projection. */
+ if (P->phi0 >= M_HALFPI - M_FORTPI / 2.0) {
+ Q->face = FACE_TOP;
+ } else if (P->phi0 <= -(M_HALFPI - M_FORTPI / 2.0)) {
+ Q->face = FACE_BOTTOM;
+ } else if (fabs(P->lam0) <= M_FORTPI) {
+ Q->face = FACE_FRONT;
+ } else if (fabs(P->lam0) <= M_HALFPI + M_FORTPI) {
+ Q->face = (P->lam0 > 0.0 ? FACE_RIGHT : FACE_LEFT);
+ } else {
+ Q->face = FACE_BACK;
+ }
+ /* Fill in useful values for the ellipsoid <-> sphere shift
+ * described in [LK12]. */
+ if (P->es != 0.0) {
+ Q->a_squared = P->a * P->a;
+ Q->one_minus_f = 1.0 - (P->a - P->b) / P->a;
+ Q->one_minus_f_squared = Q->one_minus_f * Q->one_minus_f;
+ }
+
+ return P;
+}
diff --git a/test/gie/builtins.gie b/test/gie/builtins.gie
index 1f3824c8c2..b71254f417 100644
--- a/test/gie/builtins.gie
+++ b/test/gie/builtins.gie
@@ -5273,6 +5273,173 @@ expect -223368.098302014 111769.110486991
accept -2 -1
expect -223368.098302014 -111769.110486991
+===============================================================================
+# S2
+# Input lats are converted from nice spherical values to
+# messy ellipsoidal lats e.g. 45 vs 45.1924232....
+===============================================================================
+
+-------------------------------------------------------------------------------
+operation +proj=s2 +ellps=WGS84 +lat_0=0 +lon_0=0 +UVtoST=linear
+-------------------------------------------------------------------------------
+tolerance 0.1 mm
+accept 0 0
+expect 0.5 0.5
+accept 0 45.19242321598196
+expect 0.5 1
+accept 0 -45.19242321598196
+expect 0.5 0
+accept -45 0
+expect 0 0.5
+accept 45 0
+expect 1 0.5
+accept -45 -35.446011426401625
+expect 0 0
+accept 45 -35.446011426401625
+expect 1 0
+accept 45 35.446011426401625
+expect 1 1
+accept -45 35.446011426401625
+expect 0 1
+accept 20 20.124006563576454
+expect 0.6819851171331012 0.6936645165744716
+accept 20 -20.124006563576454
+expect 0.6819851171331012 0.3063354834255284
+accept -20 -20.124006563576454
+expect 0.31801488286689883 0.3063354834255284
+accept -20 20.124006563576454
+expect 0.31801488286689883 0.6936645165744716
+
+direction inverse
+accept 0.5 0.5
+expect 0 0
+accept 0.5 1
+expect 0 45.19242321598196
+accept 0.5 0
+expect 0 -45.19242321598196
+accept 0 0.5
+expect -45 0
+accept 1 0.5
+expect 45 0
+accept 0 0
+expect -45 -35.446011426401625
+accept 1 0
+expect 45 -35.446011426401625
+accept 1 1
+expect 45 35.446011426401625
+accept 0 1
+expect -45 35.446011426401625
+accept 0.6819851171331012 0.6936645165744716
+expect 20 20.124006563576454
+accept 0.6819851171331012 0.3063354834255284
+expect 20 -20.124006563576454
+accept 0.31801488286689883 0.3063354834255284
+expect -20 -20.124006563576454
+accept 0.31801488286689883 0.6936645165744716
+expect -20 20.124006563576454
+
+-------------------------------------------------------------------------------
+operation +proj=s2 +ellps=WGS84 +lat_0=0 +lon_0=90 +UVtoST=quadratic
+-------------------------------------------------------------------------------
+tolerance 0.1 mm
+accept 90 0
+expect 0.5 0.5
+accept 70 20.124006563576454
+expect 0.27682804555233764 0.7351848576118168
+accept 110 20.124006563576454
+expect 0.7231719544476624 0.7351848576118168
+
+direction inverse
+accept 0.5 0.5
+expect 90 0
+accept 0.27682804555233764 0.7351848576118168
+expect 70 20.124006563576454
+accept 0.7231719544476624 0.7351848576118168
+expect 110 20.124006563576454
+
+-------------------------------------------------------------------------------
+operation +proj=s2 +ellps=WGS84 +lat_0=90 +UVtoST=tangent
+-------------------------------------------------------------------------------
+tolerance 0.1 mm
+accept 0 90
+expect 0.5 0.5
+accept 20 70.12337013762532
+expect 0.29020309743436806 0.4211558922141421
+accept -20 70.12337013762532
+expect 0.29020309743436806 0.5788441077858579
+
+direction inverse
+accept 0.5 0.5
+expect 0 90
+accept 0.29020309743436806 0.4211558922141421
+expect 20 70.12337013762532
+accept 0.29020309743436806 0.5788441077858579
+expect -20 70.12337013762532
+
+-------------------------------------------------------------------------------
+operation +proj=s2 +ellps=WGS84 +lat_0=0 +lon_0=180 +UVtoST=none
+-------------------------------------------------------------------------------
+tolerance 0.1 mm
+accept 180 0
+expect 0 0
+accept 160 20.124006563576454
+expect -0.3873290331489431 -0.3639702342662023
+accept -160 20.124006563576454
+expect -0.3873290331489431 0.3639702342662023
+
+direction inverse
+accept 0 0
+expect 180 0
+accept -0.3873290331489431 -0.3639702342662023
+expect 160 20.124006563576454
+accept -0.3873290331489431 0.3639702342662023
+expect -160 20.124006563576454
+
+-------------------------------------------------------------------------------
+operation +proj=s2 +ellps=WGS84 +lat_0=0 +lon_0=-90
+-------------------------------------------------------------------------------
+tolerance 0.1 mm
+accept -90 0
+expect 0.5 0.5
+accept -70 20.124006563576454
+expect 0.26481514238818316 0.7231719544476624
+accept -110 20.124006563576454
+expect 0.26481514238818316 0.27682804555233764
+
+direction inverse
+accept 0.5 0.5
+expect -90 0
+accept 0.26481514238818316 0.7231719544476624
+expect -70 20.124006563576454
+accept 0.26481514238818316 0.27682804555233764
+expect -110 20.124006563576454
+
+-------------------------------------------------------------------------------
+operation +proj=s2 +ellps=WGS84 +lat_0=-90 +UVtoST=linear
+-------------------------------------------------------------------------------
+tolerance 0.1 mm
+accept 0 -90
+expect 0.5 0.5
+accept 20 -70.12337013762533
+expect 0.5622425758450019 0.6710100716628344
+accept -20 -70.12337013762533
+expect 0.4377574241549981 0.6710100716628344
+
+direction inverse
+accept 0.5 0.5
+expect 0 -90
+accept 0.5622425758450019 0.6710100716628344
+expect 20 -70.12337013762533
+accept 0.4377574241549981 0.6710100716628344
+expect -20 -70.12337013762533
+
+-------------------------------------------------------------------------------
+operation +proj=s2 +ellps=WGS84 +lat_0=0 +lon_0=0 +UVtoST=invalid
+-------------------------------------------------------------------------------
+tolerance 0.1 mm
+
+accept 0 0
+expect failure errno invalid_op_illegal_arg_value
===============================================================================
# Sinusoidal (Sanson-Flamsteed)
diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt
index 53a14cedc8..1a080ac506 100644
--- a/test/unit/CMakeLists.txt
+++ b/test/unit/CMakeLists.txt
@@ -166,10 +166,6 @@ add_executable(test_defmodel
target_link_libraries(test_defmodel
PRIVATE GTest::gtest
PRIVATE ${PROJ_LIBRARIES})
-if(NLOHMANN_JSON STREQUAL "external")
- target_compile_definitions(test_defmodel PRIVATE EXTERNAL_NLOHMANN_JSON)
- target_link_libraries(test_defmodel PRIVATE nlohmann_json::nlohmann_json)
-endif()
add_test(NAME test_defmodel COMMAND test_defmodel)
set_property(TEST test_defmodel
PROPERTY ENVIRONMENT ${PROJ_TEST_ENVIRONMENT})
@@ -180,10 +176,6 @@ add_executable(test_tinshift
target_link_libraries(test_tinshift
PRIVATE GTest::gtest
PRIVATE ${PROJ_LIBRARIES})
-if(NLOHMANN_JSON STREQUAL "external")
- target_compile_definitions(test_tinshift PRIVATE EXTERNAL_NLOHMANN_JSON)
- target_link_libraries(test_tinshift PRIVATE nlohmann_json::nlohmann_json)
-endif()
add_test(NAME test_tinshift COMMAND test_tinshift)
set_property(TEST test_tinshift
PROPERTY ENVIRONMENT ${PROJ_TEST_ENVIRONMENT})