Make proj_lp_dist() and proj_geod() work on a PJ* CRS object

2 files changed, 36 insertions, 10 deletions

diff --git a/src/iso19111/c_api.cpp b/src/iso19111/c_api.cpp
index 0ea74394..8f1fac70 100644
--- a/src/iso19111/c_api.cpp
+++ b/src/iso19111/c_api.cpp
@@ -59,6 +59,7 @@
 #include "proj_experimental.h"
 // clang-format on
 #include "proj_constants.h"
+#include "geodesic.h"
 
 using namespace NS_PROJ::common;
 using namespace NS_PROJ::crs;
@@ -205,6 +206,26 @@ static PJ *pj_obj_create(PJ_CONTEXT *ctx, const IdentifiedObjectNNPtr &objIn) {
         pj->ctx = ctx;
         pj->descr = "ISO-19111 object";
         pj->iso_obj = objIn;
+        try {
+            auto crs = dynamic_cast<const CRS *>(objIn.get());
+            if (crs) {
+                auto geodCRS = crs->extractGeodeticCRS();
+                if (geodCRS) {
+                    const auto &ellps = geodCRS->ellipsoid();
+                    const double a = ellps->semiMajorAxis().getSIValue();
+                    const double es = ellps->squaredEccentricity();
+                    pj_calc_ellipsoid_params(pj, a, es);
+                    assert(pj->geod == nullptr);
+                    pj->geod = static_cast<struct geod_geodesic *>(
+                        calloc(1, sizeof(struct geod_geodesic)));
+                    if (pj->geod) {
+                        geod_init(pj->geod, pj->a,
+                                  pj->es / (1 + sqrt(pj->one_es)));
+                    }
+                }
+            }
+        } catch (const std::exception &) {
+        }
     }
     ctx->safeAutoCloseDbIfNeeded();
     return pj;
diff --git a/test/unit/test_c_api.cpp b/test/unit/test_c_api.cpp
index 26e4eae2..1b08e884 100644
--- a/test/unit/test_c_api.cpp
+++ b/test/unit/test_c_api.cpp
@@ -175,18 +175,23 @@ TEST_F(CApi, proj_create) {
         EXPECT_NE(obj, nullptr);
 
         // Check that functions that operate on 'non-C++' PJ don't crash
-        PJ_COORD coord;
-        coord.xyzt.x = 0;
-        coord.xyzt.y = 0;
-        coord.xyzt.z = 0;
-        coord.xyzt.t = 0;
-        EXPECT_EQ(proj_trans(obj, PJ_FWD, coord).xyzt.x,
+        constexpr double DEG_TO_RAD = .017453292519943296;
+        PJ_COORD coord1;
+        coord1.xyzt.x = 2 * DEG_TO_RAD;
+        coord1.xyzt.y = 49 * DEG_TO_RAD;
+        coord1.xyzt.z = 0;
+        coord1.xyzt.t = 0;
+        PJ_COORD coord2;
+        coord2.xyzt.x = 2 * DEG_TO_RAD;
+        coord2.xyzt.y = 50 * DEG_TO_RAD;
+        coord2.xyzt.z = 0;
+        coord2.xyzt.t = 0;
+        EXPECT_EQ(proj_trans(obj, PJ_FWD, coord1).xyzt.x,
                   std::numeric_limits<double>::infinity());
 
-        EXPECT_EQ(proj_geod(obj, coord, coord).xyzt.x,
-                  std::numeric_limits<double>::infinity());
-        EXPECT_EQ(proj_lp_dist(obj, coord, coord),
-                  std::numeric_limits<double>::infinity());
+        // and those ones actually work just fine
+        EXPECT_NEAR(proj_geod(obj, coord1, coord2).xyzt.x, 111219.409, 1e-3);
+        EXPECT_NEAR(proj_lp_dist(obj, coord1, coord2), 111219.409, 1e-3);
 
         auto info = proj_pj_info(obj);
         EXPECT_EQ(info.id, nullptr);