From 1e47f2bda51fda1ab424c04da7d2c33a38c10e87 Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Wed, 10 Mar 2021 19:01:41 +0100
Subject: Make proj_lp_dist() and proj_geod() work on a PJ* CRS object

---
 src/iso19111/c_api.cpp | 21 +++++++++++++++++++++
 1 file changed, 21 insertions(+)

(limited to 'src/iso19111/c_api.cpp')

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;
-- 
cgit v1.2.3