Merge pull request #2361 from rouault/ortho_ellipsoidal

Implement ellipsoidal formulation of +proj=ortho (fixes #397)

10 files changed, 518 insertions, 18 deletions

diff --git a/docs/source/operations/projections/ortho.rst b/docs/source/operations/projections/ortho.rst
index 2e8de9ca..2d0e08a8 100644
--- a/docs/source/operations/projections/ortho.rst
+++ b/docs/source/operations/projections/ortho.rst
@@ -10,7 +10,7 @@ around a given latitude and longitude.
 +---------------------+--------------------------------------------------------+
 | **Classification**  | Azimuthal                                              |
 +---------------------+--------------------------------------------------------+
-| **Available forms** | Forward and inverse, spherical projection              |
+| **Available forms** | Forward and inverse, spherical and ellipsoidal         |
 +---------------------+--------------------------------------------------------+
 | **Defined area**    | Global, although only one hemisphere can be seen at a  |
 |                     | time                                                   |
@@ -31,6 +31,12 @@ around a given latitude and longitude.
 
    proj-string: ``+proj=ortho``
 
+
+.. note:: Before PROJ 7.2, only the spherical formulation was implemented. If
+          wanting to replicate PROJ < 7.2 results with newer versions, the
+          ellipsoid must be forced to a sphere, for example by adding a ``+f=0``
+          parameter.
+
 Parameters
 ################################################################################
 
@@ -40,6 +46,8 @@ Parameters
 
 .. include:: ../options/lat_0.rst
 
+.. include:: ../options/ellps.rst
+
 .. include:: ../options/R.rst
 
 .. include:: ../options/x_0.rst
diff --git a/include/proj/internal/coordinateoperation_constants.hpp b/include/proj/internal/coordinateoperation_constants.hpp
index 79b47a6a..7038ba9f 100644
--- a/include/proj/internal/coordinateoperation_constants.hpp
+++ b/include/proj/internal/coordinateoperation_constants.hpp
@@ -729,6 +729,9 @@ static const MethodMapping projectionMethodMappings[] = {
     {EPSG_NAME_METHOD_ORTHOGRAPHIC, EPSG_CODE_METHOD_ORTHOGRAPHIC,
      "Orthographic", "ortho", nullptr, paramsNatOrigin},
 
+    {PROJ_WKT2_NAME_ORTHOGRAPHIC_SPHERICAL, 0, "Orthographic", "ortho", "f=0",
+     paramsNatOrigin},
+
     {PROJ_WKT2_NAME_METHOD_PATTERSON, 0, "Patterson", "patterson", nullptr,
      paramsLonNatOrigin},
 
diff --git a/include/proj/internal/esri_projection_mappings.hpp b/include/proj/internal/esri_projection_mappings.hpp
index 752c03ad..c230f26b 100644
--- a/include/proj/internal/esri_projection_mappings.hpp
+++ b/include/proj/internal/esri_projection_mappings.hpp
@@ -600,6 +600,19 @@ static const ESRIParamMapping paramsESRI_Orthographic[] = {
      EPSG_CODE_PARAMETER_LATITUDE_OF_NATURAL_ORIGIN, "0.0", false},
     {nullptr, nullptr, 0, "0.0", false}};
 
+static const ESRIParamMapping paramsESRI_Local[] = {
+    {"False_Easting", EPSG_NAME_PARAMETER_FALSE_EASTING,
+     EPSG_CODE_PARAMETER_FALSE_EASTING, "0.0", false},
+    {"False_Northing", EPSG_NAME_PARAMETER_FALSE_NORTHING,
+     EPSG_CODE_PARAMETER_FALSE_NORTHING, "0.0", false},
+    {"Scale_Factor", nullptr, 0, "1.0", false},
+    {"Azimuth", nullptr, 0, "0.0", false},
+    {"Longitude_Of_Center", EPSG_NAME_PARAMETER_LONGITUDE_OF_NATURAL_ORIGIN,
+     EPSG_CODE_PARAMETER_LONGITUDE_OF_NATURAL_ORIGIN, "0.0", false},
+    {"Latitude_Of_Center", EPSG_NAME_PARAMETER_LATITUDE_OF_NATURAL_ORIGIN,
+     EPSG_CODE_PARAMETER_LATITUDE_OF_NATURAL_ORIGIN, "0.0", false},
+    {nullptr, nullptr, 0, "0.0", false}};
+
 static const ESRIParamMapping paramsESRI_Winkel_Tripel[] = {
     {"False_Easting", EPSG_NAME_PARAMETER_FALSE_EASTING,
      EPSG_CODE_PARAMETER_FALSE_EASTING, "0.0", false},
@@ -1004,8 +1017,10 @@ static const ESRIMethodMapping esriMappings[] = {
      EPSG_CODE_METHOD_KROVAK_NORTH_ORIENTED, paramsESRI_Krovak_alt2},
     {"New_Zealand_Map_Grid", EPSG_NAME_METHOD_NZMG, EPSG_CODE_METHOD_NZMG,
      paramsESRI_New_Zealand_Map_Grid},
-    {"Orthographic", EPSG_NAME_METHOD_ORTHOGRAPHIC,
-     EPSG_CODE_METHOD_ORTHOGRAPHIC, paramsESRI_Orthographic},
+    {"Orthographic", PROJ_WKT2_NAME_ORTHOGRAPHIC_SPHERICAL, 0,
+     paramsESRI_Orthographic},
+    {"Local", EPSG_NAME_METHOD_ORTHOGRAPHIC, EPSG_CODE_METHOD_ORTHOGRAPHIC,
+     paramsESRI_Local},
     {"Winkel_Tripel", "Winkel Tripel", 0, paramsESRI_Winkel_Tripel},
     {"Aitoff", "Aitoff", 0, paramsESRI_Aitoff},
     {"Flat_Polar_Quartic", PROJ_WKT2_NAME_METHOD_FLAT_POLAR_QUARTIC, 0,
diff --git a/scripts/build_esri_projection_mapping.py b/scripts/build_esri_projection_mapping.py
index 5206e191..12eeb852 100644
--- a/scripts/build_esri_projection_mapping.py
+++ b/scripts/build_esri_projection_mapping.py
@@ -436,11 +436,23 @@ config_str = """
         - Longitude_Of_Origin: EPSG_NAME_PARAMETER_LONGITUDE_OF_NATURAL_ORIGIN
         - Latitude_Of_Origin: EPSG_NAME_PARAMETER_LATITUDE_OF_NATURAL_ORIGIN
 
+# ESRI's Orthographic is a spherical-only formulation. The ellipsoidal capable
+# name is Local. See below
 - Orthographic:
+    WKT2_name: PROJ_WKT2_NAME_ORTHOGRAPHIC_SPHERICAL
+    Params:
+        - False_Easting: EPSG_NAME_PARAMETER_FALSE_EASTING
+        - False_Northing: EPSG_NAME_PARAMETER_FALSE_NORTHING
+        - Longitude_Of_Center: EPSG_NAME_PARAMETER_LONGITUDE_OF_NATURAL_ORIGIN
+        - Latitude_Of_Center: EPSG_NAME_PARAMETER_LATITUDE_OF_NATURAL_ORIGIN
+
+- Local:
     WKT2_name: EPSG_NAME_METHOD_ORTHOGRAPHIC
     Params:
         - False_Easting: EPSG_NAME_PARAMETER_FALSE_EASTING
         - False_Northing: EPSG_NAME_PARAMETER_FALSE_NORTHING
+        - Scale_Factor: 1.0
+        - Azimuth: 0.0
         - Longitude_Of_Center: EPSG_NAME_PARAMETER_LONGITUDE_OF_NATURAL_ORIGIN
         - Latitude_Of_Center: EPSG_NAME_PARAMETER_LATITUDE_OF_NATURAL_ORIGIN
 
diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index 65ef2b77..ef262cc8 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -4251,7 +4251,7 @@ ConversionNNPtr Conversion::createObliqueStereographic(
  * This method is defined as [EPSG:9840]
  * (https://www.epsg-registry.org/export.htm?gml=urn:ogc:def:method:EPSG::9840)
  *
- * \note At the time of writing, PROJ only implements the spherical formulation
+ * \note Before PROJ 7.2, only the spherical formulation was implemented.
  *
  * @param properties See \ref general_properties of the conversion. If the name
  * is not provided, it is automatically set.
diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp
index 6680ab00..16ab22f7 100644
--- a/src/iso19111/io.cpp
+++ b/src/iso19111/io.cpp
@@ -3290,7 +3290,7 @@ ConversionNNPtr WKTParser::Private::buildProjectionFromESRI(
     }
 
     // Compare parameters present with the ones expected in the mapping
-    const ESRIMethodMapping *esriMapping = esriMappings[0];
+    const ESRIMethodMapping *esriMapping = nullptr;
     int bestMatchCount = -1;
     for (const auto &mapping : esriMappings) {
         int matchCount = 0;
@@ -3298,12 +3298,20 @@ ConversionNNPtr WKTParser::Private::buildProjectionFromESRI(
             auto iter = mapParamNameToValue.find(param->esri_name);
             if (iter != mapParamNameToValue.end()) {
                 if (param->wkt2_name == nullptr) {
+                    bool ok = true;
                     try {
                         if (io::asDouble(param->fixed_value) ==
                             io::asDouble(iter->second)) {
                             matchCount++;
+                        } else {
+                            ok = false;
                         }
                     } catch (const std::exception &) {
+                        ok = false;
+                    }
+                    if (!ok) {
+                        matchCount = -1;
+                        break;
                     }
                 } else {
                     matchCount++;
@@ -3317,6 +3325,10 @@ ConversionNNPtr WKTParser::Private::buildProjectionFromESRI(
             bestMatchCount = matchCount;
         }
     }
+    if (esriMapping == nullptr) {
+        return buildProjectionStandard(baseGeodCRS, projCRSNode, projectionNode,
+                                       defaultLinearUnit, defaultAngularUnit);
+    }
 
     std::map<std::string, const char *> mapWKT2NameToESRIName;
     for (const auto *param = esriMapping->params; param->esri_name; ++param) {
@@ -8898,8 +8910,29 @@ static bool is_in_stringlist(const std::string &str, const char *stringlist) {
 CRSNNPtr PROJStringParser::Private::buildProjectedCRS(
     int iStep, GeographicCRSNNPtr geogCRS, int iUnitConvert, int iAxisSwap) {
     auto &step = steps_[iStep];
-    auto mappings = getMappingsFromPROJName(step.name);
+    const auto mappings = getMappingsFromPROJName(step.name);
     const MethodMapping *mapping = mappings.empty() ? nullptr : mappings[0];
+
+    if (mappings.size() >= 2) {
+        // To distinguish for example +ortho from +ortho +f=0
+        for (const auto *mappingIter : mappings) {
+            if (mappingIter->proj_name_aux != nullptr &&
+                strchr(mappingIter->proj_name_aux, '=') == nullptr &&
+                hasParamValue(step, mappingIter->proj_name_aux)) {
+                mapping = mappingIter;
+                break;
+            } else if (mappingIter->proj_name_aux != nullptr &&
+                       strchr(mappingIter->proj_name_aux, '=') != nullptr) {
+                const auto tokens = split(mappingIter->proj_name_aux, '=');
+                if (tokens.size() == 2 &&
+                    getParamValue(step, tokens[0]) == tokens[1]) {
+                    mapping = mappingIter;
+                    break;
+                }
+            }
+        }
+    }
+
     if (mapping) {
         mapping = selectSphericalOrEllipsoidal(mapping, geogCRS);
     }
diff --git a/src/proj_constants.h b/src/proj_constants.h
index 5b626223..f64bf496 100644
--- a/src/proj_constants.h
+++ b/src/proj_constants.h
@@ -183,6 +183,9 @@
 #define EPSG_NAME_METHOD_ORTHOGRAPHIC "Orthographic"
 #define EPSG_CODE_METHOD_ORTHOGRAPHIC 9840
 
+#define PROJ_WKT2_NAME_ORTHOGRAPHIC_SPHERICAL \
+    "Orthographic (Spherical)"
+
 #define PROJ_WKT2_NAME_METHOD_PATTERSON "Patterson"
 
 #define EPSG_NAME_METHOD_AMERICAN_POLYCONIC "American Polyconic"
diff --git a/src/projections/ortho.cpp b/src/projections/ortho.cpp
index 5f9366c3..75b8199d 100644
--- a/src/projections/ortho.cpp
+++ b/src/projections/ortho.cpp
@@ -4,7 +4,7 @@
 #include "proj_internal.h"
 #include <math.h>
 
-PROJ_HEAD(ortho, "Orthographic") "\n\tAzi, Sph";
+PROJ_HEAD(ortho, "Orthographic") "\n\tAzi, Sph&Ell";
 
 namespace { // anonymous namespace
 enum Mode {
@@ -19,6 +19,9 @@ namespace { // anonymous namespace
 struct pj_opaque {
     double  sinph0;
     double  cosph0;
+    double  nu0;
+    double  y_shift;
+    double  y_scale;
     enum Mode mode;
 };
 } // anonymous namespace
@@ -49,6 +52,13 @@ static PJ_XY ortho_s_forward (PJ_LP lp, PJ *P) {           /* Spheroidal, forwar
         break;
     case OBLIQ:
         sinphi = sin(lp.phi);
+
+        // Is the point visible from the projection plane ?
+        // From https://lists.osgeo.org/pipermail/proj/2020-September/009831.html
+        // this is the dot product of the normal of the ellipsoid at the center of
+        // the projection and at the point considered for projection.
+        // [cos(phi)*cos(lambda), cos(phi)*sin(lambda), sin(phi)]
+        // Also from Snyder's Map Projection - A working manual, equation (5-3), page 149
         if (Q->sinph0 * sinphi + Q->cosph0 * cosphi * coslam < - EPS10)
             return forward_error(P, lp, xy);
         xy.y = Q->cosph0 * sinphi - Q->sinph0 * cosphi * coslam;
@@ -121,6 +131,152 @@ static PJ_LP ortho_s_inverse (PJ_XY xy, PJ *P) {           /* Spheroidal, invers
 }
 
 
+static PJ_XY ortho_e_forward (PJ_LP lp, PJ *P) {           /* Ellipsoidal, forward */
+    PJ_XY xy;
+    struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
+
+    // From EPSG guidance note 7.2, March 2020, §3.3.5 Orthographic
+    const double cosphi = cos(lp.phi);
+    const double sinphi = sin(lp.phi);
+    const double coslam = cos(lp.lam);
+    const double sinlam = sin(lp.lam);
+
+    // Is the point visible from the projection plane ?
+    // Same condition as in spherical case
+    if (Q->sinph0 * sinphi + Q->cosph0 * cosphi * coslam < - EPS10) {
+        xy.x = HUGE_VAL; xy.y = HUGE_VAL;
+        return forward_error(P, lp, xy);
+    }
+
+    const double nu = 1.0 / sqrt(1.0 - P->es * sinphi * sinphi);
+    xy.x = nu * cosphi * sinlam;
+    xy.y = nu * (sinphi * Q->cosph0 - cosphi * Q->sinph0 * coslam) +
+            P->es * (Q->nu0 * Q->sinph0 - nu * sinphi) * Q->cosph0;
+
+    return xy;
+}
+
+
+static PJ_LP ortho_e_inverse (PJ_XY xy, PJ *P) {           /* Ellipsoidal, inverse */
+    PJ_LP lp;
+    struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
+
+    const auto SQ = [](double x) { return x*x; };
+
+    if( Q->mode == N_POLE || Q->mode == S_POLE )
+    {
+        // Polar case. Forward case equations can be simplified as:
+        // xy.x = nu * cosphi * sinlam
+        // xy.y = nu * -cosphi * coslam * sign(phi0)
+        // ==> lam = atan2(xy.x, -xy.y * sign(phi0))
+        // ==> xy.x^2 + xy.y^2 = nu^2 * cosphi^2
+        //                rh^2 = cosphi^2 / (1 - es * sinphi^2)
+        // ==>  cosphi^2 = rh^2 * (1 - es) / (1 - es * rh^2)
+
+        const double rh2 = SQ(xy.x) + SQ(xy.y);
+        if (rh2 >= 1. - 1e-15) {
+            if ((rh2 - 1.) > EPS10) {
+                proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
+                proj_log_trace(P, "Point (%.3f, %.3f) is outside the projection boundary");
+                lp.lam = HUGE_VAL; lp.phi = HUGE_VAL;
+                return lp;
+            }
+            lp.phi = 0;
+        }
+        else {
+            lp.phi = acos(sqrt(rh2 * P->one_es / (1 - P->es * rh2))) * (Q->mode == N_POLE ? 1 : -1);
+        }
+        lp.lam = atan2(xy.x, xy.y * (Q->mode == N_POLE ? -1 : 1));
+        return lp;
+    }
+
+    if( Q->mode == EQUIT )
+    {
+        // Equatorial case. Forward case equations can be simplified as:
+        // xy.x = nu * cosphi * sinlam
+        // xy.y  = nu * sinphi * (1 - P->es)
+        // x^2 * (1 - es * sinphi^2) = (1 - sinphi^2) * sinlam^2
+        // y^2 / ((1 - es)^2 + y^2 * es) = sinphi^2
+
+        // Equation of the ellipse
+        if( SQ(xy.x) + SQ(xy.y * (P->a / P->b)) > 1 + 1e-11 ) {
+            proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
+            proj_log_trace(P, "Point (%.3f, %.3f) is outside the projection boundary");
+            lp.lam = HUGE_VAL; lp.phi = HUGE_VAL;
+            return lp;
+        }
+
+        const double sinphi2 = xy.y == 0 ? 0 : 1.0 / (SQ((1 - P->es) / xy.y) + P->es);
+        if( sinphi2 > 1 - 1e-11 ) {
+            lp.phi = M_PI_2 * (xy.y > 0 ? 1 : -1);
+            lp.lam = 0;
+            return lp;
+        }
+        lp.phi = asin(sqrt(sinphi2)) * (xy.y > 0 ? 1 : -1);
+        const double sinlam = xy.x * sqrt((1 - P->es * sinphi2) / (1 - sinphi2));
+        if( fabs(sinlam) - 1 > -1e-15 )
+            lp.lam = M_PI_2 * (xy.x > 0 ? 1: -1);
+        else
+            lp.lam = asin(sinlam);
+        return lp;
+    }
+
+    // Using Q->sinph0 * sinphi + Q->cosph0 * cosphi * coslam == 0 (visibity
+    // condition of the forward case) in the forward equations, and a lot of
+    // substition games...
+    PJ_XY xy_recentered;
+    xy_recentered.x = xy.x;
+    xy_recentered.y = (xy.y - Q->y_shift) / Q->y_scale;
+    if( SQ(xy.x) + SQ(xy_recentered.y) > 1 + 1e-11 ) {
+        proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
+        proj_log_trace(P, "Point (%.3f, %.3f) is outside the projection boundary");
+        lp.lam = HUGE_VAL; lp.phi = HUGE_VAL;
+        return lp;
+    }
+
+    // From EPSG guidance note 7.2, March 2020, §3.3.5 Orthographic
+
+    // It suggests as initial guess:
+    // lp.lam = 0;
+    // lp.phi = P->phi0;
+    // But for poles, this will not converge well. Better use:
+    lp = ortho_s_inverse(xy_recentered, P);
+
+    for( int i = 0; i < 20; i++ )
+    {
+        const double cosphi = cos(lp.phi);
+        const double sinphi = sin(lp.phi);
+        const double coslam = cos(lp.lam);
+        const double sinlam = sin(lp.lam);
+        const double one_minus_es_sinphi2 = 1.0 - P->es * sinphi * sinphi;
+        const double nu = 1.0 / sqrt(one_minus_es_sinphi2);
+        PJ_XY xy_new;
+        xy_new.x = nu * cosphi * sinlam;
+        xy_new.y = nu * (sinphi * Q->cosph0 - cosphi * Q->sinph0 * coslam) +
+                P->es * (Q->nu0 * Q->sinph0 - nu * sinphi) * Q->cosph0;
+        const double rho = (1.0 - P->es) * nu / one_minus_es_sinphi2;
+        const double J11 = -rho * sinphi * sinlam;
+        const double J12 = nu * cosphi * coslam;
+        const double J21 = rho * (cosphi * Q->cosph0 + sinphi * Q->sinph0 * coslam);
+        const double J22 = nu * Q->sinph0 * Q->cosph0 * sinlam;
+        const double D = J11 * J22 - J12 * J21;
+        const double dx = xy.x - xy_new.x;
+        const double dy = xy.y - xy_new.y;
+        const double dphi = (J22 * dx - J12 * dy) / D;
+        const double dlam = (-J21 * dx + J11 * dy) / D;
+        lp.phi += dphi;
+        if( lp.phi > M_PI_2) lp.phi = M_PI_2;
+        else if( lp.phi < -M_PI_2) lp.phi = -M_PI_2;
+        lp.lam += dlam;
+        if( fabs(dphi) < 1e-12 && fabs(dlam) < 1e-12 )
+        {
+            return lp;
+        }
+    }
+    pj_ctx_set_errno(P->ctx, PJD_ERR_NON_CONVERGENT);
+    return lp;
+}
+
 
 PJ *PROJECTION(ortho) {
     struct pj_opaque *Q = static_cast<struct pj_opaque*>(pj_calloc (1, sizeof (struct pj_opaque)));
@@ -128,17 +284,27 @@ PJ *PROJECTION(ortho) {
         return pj_default_destructor(P, ENOMEM);
     P->opaque = Q;
 
+    Q->sinph0 = sin(P->phi0);
+    Q->cosph0 = cos(P->phi0);
     if (fabs(fabs(P->phi0) - M_HALFPI) <= EPS10)
         Q->mode = P->phi0 < 0. ? S_POLE : N_POLE;
     else if (fabs(P->phi0) > EPS10) {
         Q->mode = OBLIQ;
-        Q->sinph0 = sin(P->phi0);
-        Q->cosph0 = cos(P->phi0);
     } else
         Q->mode = EQUIT;
-    P->inv = ortho_s_inverse;
-    P->fwd = ortho_s_forward;
-    P->es = 0.;
+    if( P->es == 0 )
+    {
+        P->inv = ortho_s_inverse;
+        P->fwd = ortho_s_forward;
+    }
+    else
+    {
+        Q->nu0 = 1.0 / sqrt(1.0 - P->es * Q->sinph0 * Q->sinph0);
+        Q->y_shift = P->es * Q->nu0 * Q->sinph0 * Q->cosph0;
+        Q->y_scale = 1.0 / sqrt(1.0 - P->es * Q->cosph0 * Q->cosph0);
+        P->inv = ortho_e_inverse;
+        P->fwd = ortho_e_forward;
+    }
 
     return P;
 }
diff --git a/test/gie/builtins.gie b/test/gie/builtins.gie
index f755526e..02e6766e 100644
--- a/test/gie/builtins.gie
+++ b/test/gie/builtins.gie
@@ -4266,11 +4266,11 @@ expect  -223374.577355253 -111701.072127637
 
 ===============================================================================
 # Orthographic
-# 	Azi, Sph.
+# 	Azi, Sph&Ell.
 ===============================================================================
 
 -------------------------------------------------------------------------------
-# Test the equatorial aspect of the Orthopgraphic projection.
+# Test the equatorial aspect of the Orthographic projection.
 
 # Test data from Snyder (1987), table 22, p. 151.
 -------------------------------------------------------------------------------
@@ -4322,7 +4322,7 @@ expect      failure errno tolerance_condition
 
 
 -------------------------------------------------------------------------------
-# Test the oblique aspect of the Orthopgraphic projection.
+# Test the oblique aspect of the Orthographic projection.
 
 # Test data from Snyder (1987), table 23, pp. 152-153.
 -------------------------------------------------------------------------------
@@ -4356,7 +4356,7 @@ expect      failure errno tolerance_condition
 
 
 -------------------------------------------------------------------------------
-# Test the north polar aspect of the Orthopgraphic projection.
+# Test the north polar aspect of the Orthographic projection.
 -------------------------------------------------------------------------------
 operation   +proj=ortho +R=1 +lat_0=90 +lon_0=0
 -------------------------------------------------------------------------------
@@ -4386,7 +4386,7 @@ accept      2   2
 expect      failure errno tolerance_condition
 
 -------------------------------------------------------------------------------
-# Test the south polar aspect of the Orthopgraphic projection.
+# Test the south polar aspect of the Orthographic projection.
 -------------------------------------------------------------------------------
 operation   +proj=ortho +R=1 +lat_0=-90 +lon_0=0
 -------------------------------------------------------------------------------
@@ -4421,6 +4421,200 @@ expect      failure errno tolerance_condition
 accept      0.70710678118 0.7071067812
 expect      45  0
 
+
+-------------------------------------------------------------------------------
+# Test the ellipsoidal formulation
+
+# Test data from Guidance Note 7 part 2, March 2020,  p. 90
+-------------------------------------------------------------------------------
+operation   +proj=ortho +ellps=WGS84 +lat_0=55 +lon_0=5
+-------------------------------------------------------------------------------
+tolerance   1 mm
+accept      2.12955      53.80939444444444
+expect      -189011.711  -128640.567
+roundtrip   1
+
+-------------------------------------------------------------------------------
+# Oblique
+-------------------------------------------------------------------------------
+operation   +proj=ortho +ellps=WGS84 +lat_0=30
+
+# On boundary of visibility domain.
+direction   forward
+tolerance   0.1 mm
+accept      -90             0
+expect      -6378137        18504.1253
+
+# This test is fragile. Note the slighly important tolerance
+# direction   inverse
+# tolerance   100 mm
+# accept      -6378137        18504.125313223721605027
+# expect      -90             0
+
+# Slightly outside
+direction   inverse
+accept      -6378137.001    18504.1253
+expect      failure errno tolerance_condition
+
+# On boundary of visibility domain
+direction   forward
+tolerance   0.1 mm
+accept      0               -60
+expect      0               -6343601.0991
+
+# Just on it, but fails to converge. This test might be fragile
+direction   inverse
+accept      0               -6343601.099075031466782093
+expect      failure errno non_convergent
+
+# Slightly inside
+direction   inverse
+tolerance   0.1 mm
+accept      0               -6343600
+expect      0               -59.966377950099655436
+
+
+-------------------------------------------------------------------------------
+# Equatorial
+-------------------------------------------------------------------------------
+operation   +proj=ortho +ellps=WGS84
+tolerance   0.1 mm
+accept      0               0
+expect      0               0
+roundtrip   1
+
+accept      1               1
+expect      111296.9991     110568.7748
+roundtrip   1
+
+accept      1               -1
+expect      111296.9991     -110568.7748
+roundtrip   1
+
+accept      -1              1
+expect      -111296.9991    110568.7748
+roundtrip   1
+
+accept      -1              -1
+expect      -111296.9991    -110568.7748
+roundtrip   1
+
+accept      89.99           0
+expect      6378136.9029    0
+roundtrip   1
+
+accept      -89.99          0
+expect      -6378136.9029   0
+roundtrip   1
+
+accept      0               89.99
+expect      0               6356752.2167
+roundtrip   1
+
+accept      0               -89.99
+expect      0               -6356752.2167
+roundtrip   1
+
+# Point not visible from the projection plane
+accept      90.00001        0
+expect      failure         errno tolerance_condition
+
+# Point not visible from the projection plane
+accept      -90.00001       0
+expect      failure         errno tolerance_condition
+
+# Consistant with WGS84 semi-major axis
+accept      90              0
+expect      6378137         0
+roundtrip   1
+
+accept      -90             0
+expect      -6378137        0
+roundtrip   1
+
+# Consistant with WGS84 semi-minor axis
+accept      0               90
+expect      0               6356752.3142
+roundtrip   1
+
+accept      0               -90
+expect      0               -6356752.3142
+roundtrip   1
+
+# Point not visible from the projection plane
+direction   inverse
+accept      0               6356752.3143
+expect      failure         errno tolerance_condition
+
+# Point not visible from the projection plane
+direction   inverse
+accept      1000            6356752.314
+expect      failure         errno tolerance_condition
+
+# Point not visible from the projection plane
+direction   inverse
+accept      6378137.0001    0
+expect      failure         errno tolerance_condition
+
+-------------------------------------------------------------------------------
+# North pole tests
+-------------------------------------------------------------------------------
+operation   +proj=ortho +ellps=WGS84 +lat_0=90
+tolerance   0.1 mm
+accept      0               90
+expect      0               0
+roundtrip   1
+
+accept      30              45
+expect      2258795.4394   -3912348.4650
+roundtrip   1
+
+accept      135             89.999999873385
+expect      0.01            0.01
+roundtrip   1
+
+# Point not visible from the projection plane
+accept      0               -0.0000001
+expect      failure         errno tolerance_condition
+
+# Consistant with WGS84 semi-major axis
+accept      0               0
+expect      0               -6378137
+roundtrip   1
+
+# Point not visible from the projection plane
+direction   inverse
+accept      0               -6378137.1
+expect      failure         errno tolerance_condition
+
+-------------------------------------------------------------------------------
+# South pole tests
+-------------------------------------------------------------------------------
+operation   +proj=ortho +ellps=WGS84 +lat_0=-90
+tolerance   0.1 mm
+accept      0               -90
+expect      0               0
+roundtrip   1
+
+accept      135             -89.999999873385
+expect      0.01            -0.01
+roundtrip   1
+
+# Point not visible from the projection plane
+accept      0               0.0000001
+expect      failure         errno tolerance_condition
+
+# Consistant with WGS84 semi-major axis
+accept      0               0
+expect      0               6378137
+roundtrip   1
+
+# Point not visible from the projection plane
+direction   inverse
+accept      0               6378137.1
+expect      failure         errno tolerance_condition
+
+
 ===============================================================================
 # Perspective Conic
 # 	Conic, Sph
diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp
index f2464e5c..34396fcf 100644
--- a/test/unit/test_io.cpp
+++ b/test/unit/test_io.cpp
@@ -5225,6 +5225,21 @@ static const struct {
       {"False_Northing", 2},
       {"Longitude_Of_Center", 3},
       {"Latitude_Of_Center", 4}},
+     "Orthographic (Spherical)",
+     {
+         {"Latitude of natural origin", 4},
+         {"Longitude of natural origin", 3},
+         {"False easting", 1},
+         {"False northing", 2},
+     }},
+
+    {"Local",
+     {{"False_Easting", 1},
+      {"False_Northing", 2},
+      {"Scale_Factor", 1},
+      {"Azimuth", 0},
+      {"Longitude_Of_Center", 3},
+      {"Latitude_Of_Center", 4}},
      "Orthographic",
      {
          {"Latitude of natural origin", 4},
@@ -5233,6 +5248,24 @@ static const struct {
          {"False northing", 2},
      }},
 
+    // Local with unsupported value for Azimuth
+    {"Local",
+     {{"False_Easting", 1},
+      {"False_Northing", 2},
+      {"Scale_Factor", 1},
+      {"Azimuth", 123},
+      {"Longitude_Of_Center", 3},
+      {"Latitude_Of_Center", 4}},
+     "Local",
+     {
+         {"False_Easting", 1},
+         {"False_Northing", 2},
+         {"Scale_Factor", 1},
+         {"Azimuth", 123},
+         {"Longitude_Of_Center", 3},
+         {"Latitude_Of_Center", 4},
+     }},
+
     {"Winkel_Tripel",
      {{"False_Easting", 1},
       {"False_Northing", 2},
@@ -5543,7 +5576,6 @@ static const struct {
          {"Longitude_Of_Origin", 3},
          {"Latitude_Of_Origin", 4}},
     },
-
 };
 
 TEST(wkt_parse, esri_projcs) {
@@ -9201,6 +9233,40 @@ TEST(io, projparse_non_earth_ellipsoid) {
 
 // ---------------------------------------------------------------------------
 
+TEST(io, projparse_ortho_ellipsoidal) {
+    std::string input("+proj=ortho +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 "
+                      "+ellps=WGS84 +units=m +no_defs +type=crs");
+    auto obj = PROJStringParser().createFromPROJString(input);
+    auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+    EXPECT_EQ(crs->derivingConversion()->method()->getEPSGCode(),
+              EPSG_CODE_METHOD_ORTHOGRAPHIC);
+    EXPECT_EQ(
+        crs->exportToPROJString(
+            PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
+                .get()),
+        input);
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(io, projparse_ortho_spherical) {
+    std::string input("+proj=ortho +f=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 "
+                      "+ellps=WGS84 +units=m +no_defs +type=crs");
+    auto obj = PROJStringParser().createFromPROJString(input);
+    auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj);
+    ASSERT_TRUE(crs != nullptr);
+    EXPECT_EQ(crs->derivingConversion()->method()->nameStr(),
+              PROJ_WKT2_NAME_ORTHOGRAPHIC_SPHERICAL);
+    EXPECT_EQ(
+        crs->exportToPROJString(
+            PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4)
+                .get()),
+        input);
+}
+
+// ---------------------------------------------------------------------------
+
 TEST(io, projparse_axisswap_unitconvert_longlat_proj) {
     std::string input =
         "+type=crs +proj=pipeline +step +proj=axisswap +order=2,1 +step "