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

- Map ESRI 'Local' to +proj=ortho when Scale_Factor = 1 and Azimuth = 0
- Map ESRI 'Orthographic' to a PROJ WKT2 'Orthographic (Spherical)'
  which maps to +proj=ortho +f=0 to froce spherical evaluation

2 files changed, 134 insertions, 2 deletions

diff --git a/test/gie/builtins.gie b/test/gie/builtins.gie
index f755526e..48afc00a 100644
--- a/test/gie/builtins.gie
+++ b/test/gie/builtins.gie
@@ -4266,7 +4266,7 @@ expect  -223374.577355253 -111701.072127637
 
 ===============================================================================
 # Orthographic
-# 	Azi, Sph.
+# 	Azi, Sph&Ell.
 ===============================================================================
 
 -------------------------------------------------------------------------------
@@ -4421,6 +4421,72 @@ 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
+
+# Equatorial
+operation   +proj=ortho +ellps=WGS84
+tolerance   0.1 mm
+accept      0               0
+expect      0               0
+roundtrip   1
+
+accept      0               89.99
+expect      0               6356752.2167
+roundtrip   1
+
+accept      0               -89.99
+expect      0               -6356752.2167
+roundtrip   1
+
+# Consistant with WGS84 semi-minor axis
+# The inverse transformation doesn't converge due to properties of the projection
+accept      0               90
+expect      0               6356752.3142
+
+# North pole tests
+operation   +proj=ortho +ellps=WGS84 +lat_0=90
+tolerance   0.1 mm
+accept      0               90
+expect      0               0
+roundtrip   1
+
+accept      1               89.999
+expect      1.9493         -111.6770
+roundtrip   1
+
+# Consistant with WGS84 semi-major axis
+# The inverse transformation doesn't converge due to properties of the projection
+accept      0               0
+expect      0               -6378137
+
+# South pole tests
+operation   +proj=ortho +ellps=WGS84 +lat_0=-90
+tolerance   0.1 mm
+accept      0               -90
+expect      0               0
+roundtrip   1
+
+accept      1               -89.999
+expect      1.9493          111.6770
+roundtrip   1
+
+# Consistant with WGS84 semi-major axis
+# The inverse transformation doesn't converge due to properties of the projection
+accept      0               0
+expect      0               6378137
+
+
 ===============================================================================
 # 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 "