From 3e33c27c17cd25712dace6e26214bb8e25356043 Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Fri, 10 Nov 2023 19:37:36 +0100
Subject: [PATCH] proj: projection factors: fix when input is a compound CRS of
 a projected CRS, and a projected CRS whose datum has a non-Greenwich prime
 meridian

---
 src/apps/proj.cpp          | 51 +++++++++++++++++++++++++++++++-------
 test/cli/testproj          | 15 +++++++++++
 test/cli/testproj_out.dist |  2 ++
 3 files changed, 59 insertions(+), 9 deletions(-)

diff --git a/src/apps/proj.cpp b/src/apps/proj.cpp
index 515461993f..099c4b5dce 100644
--- a/src/apps/proj.cpp
+++ b/src/apps/proj.cpp
@@ -537,7 +537,20 @@ int main(int argc, char **argv) {
         eargc--;
         // logic copied from proj_factors function
         if (PJ *P = proj_create(nullptr, ocrs.c_str())) {
-            const auto type = proj_get_type(P);
+            auto type = proj_get_type(P);
+            auto ctx = P->ctx;
+            if (type == PJ_TYPE_COMPOUND_CRS) {
+                auto horiz = proj_crs_get_sub_crs(ctx, P, 0);
+                if (horiz) {
+                    if (proj_get_type(horiz) == PJ_TYPE_PROJECTED_CRS) {
+                        proj_destroy(P);
+                        P = horiz;
+                        type = proj_get_type(P);
+                    } else {
+                        proj_destroy(horiz);
+                    }
+                }
+            }
             if (type == PJ_TYPE_PROJECTED_CRS) {
                 try {
                     auto crs = dynamic_cast<const NS_PROJ::crs::ProjectedCRS *>(
@@ -548,18 +561,38 @@ int main(int argc, char **argv) {
                                   dir == NS_PROJ::cs::AxisDirection::SOUTH;
                 } catch (...) {
                 }
-                auto ctx = P->ctx;
                 auto geodetic_crs = proj_get_source_crs(ctx, P);
                 assert(geodetic_crs);
-                auto datum = proj_crs_get_datum(ctx, geodetic_crs);
-                auto datum_ensemble =
-                    proj_crs_get_datum_ensemble(ctx, geodetic_crs);
+                auto pm = proj_get_prime_meridian(ctx, geodetic_crs);
+                double pm_longitude = 0;
+                proj_prime_meridian_get_parameters(ctx, pm, &pm_longitude,
+                                                   nullptr, nullptr);
+                proj_destroy(pm);
+                PJ *geogCRSNormalized;
                 auto cs = proj_create_ellipsoidal_2D_cs(
                     ctx, PJ_ELLPS2D_LONGITUDE_LATITUDE, "Radian", 1.0);
-                auto geogCRSNormalized = proj_create_geographic_crs_from_datum(
-                    ctx, "unnamed crs", datum ? datum : datum_ensemble, cs);
-                proj_destroy(datum);
-                proj_destroy(datum_ensemble);
+                if (pm_longitude != 0) {
+                    auto ellipsoid = proj_get_ellipsoid(ctx, geodetic_crs);
+                    double semi_major_metre = 0;
+                    double inv_flattening = 0;
+                    proj_ellipsoid_get_parameters(ctx, ellipsoid,
+                                                  &semi_major_metre, nullptr,
+                                                  nullptr, &inv_flattening);
+                    geogCRSNormalized = proj_create_geographic_crs(
+                        ctx, "unname crs", "unnamed datum",
+                        proj_get_name(ellipsoid), semi_major_metre,
+                        inv_flattening, "reference prime meridian", 0, nullptr,
+                        0, cs);
+                    proj_destroy(ellipsoid);
+                } else {
+                    auto datum = proj_crs_get_datum(ctx, geodetic_crs);
+                    auto datum_ensemble =
+                        proj_crs_get_datum_ensemble(ctx, geodetic_crs);
+                    geogCRSNormalized = proj_create_geographic_crs_from_datum(
+                        ctx, "unnamed crs", datum ? datum : datum_ensemble, cs);
+                    proj_destroy(datum);
+                    proj_destroy(datum_ensemble);
+                }
                 proj_destroy(cs);
                 Proj = proj_create_crs_to_crs_from_pj(ctx, geogCRSNormalized, P,
                                                       nullptr, nullptr);
diff --git a/test/cli/testproj b/test/cli/testproj
index 375f53cdbf..67f937430a 100755
--- a/test/cli/testproj
+++ b/test/cli/testproj
@@ -41,6 +41,21 @@ echo "Test CRS option"
 $EXE EPSG:32620 -S >>${OUT} <<EOF
 -63 0
 EOF
+echo "Test projection factors on projected CRS with non-Greenwhich prime meridian"
+#
+$EXE EPSG:27571 -S >>${OUT} <<EOF
+2.33722917 49.5
+EOF
+echo "Test projection factors on compound CRS with a projected CRS"
+#
+$EXE EPSG:5972 -S >>${OUT} <<EOF
+9 0
+EOF
+
+# On some architectures the angular distortion (omega) of EPSG:27571 is
+# not exactly 0, but 8.53878e-07
+cat ${OUT} | sed "s/8.53878e-07/0/" > ${OUT}.tmp
+mv ${OUT}.tmp ${OUT}
 
 #
 # do 'diff' with distribution results
diff --git a/test/cli/testproj_out.dist b/test/cli/testproj_out.dist
index 73320d0a1e..f7b6036a8a 100644
--- a/test/cli/testproj_out.dist
+++ b/test/cli/testproj_out.dist
@@ -1,2 +1,4 @@
 2 49	2.000	49.000
 500000.00	0.00	<0.9996 0.9996 0.9992 0 0.9996 0.9996>
+600000.00	1200000.00	<0.999877 0.999877 0.999755 0 0.999877 0.999877>
+500000.00	0.00	<0.9996 0.9996 0.9992 0 0.9996 0.9996>