From c34f791daeac7cc6d7c8d0efe6afb9d4dd1f4434 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Fri, 4 Nov 2022 19:39:35 +0100 Subject: [PATCH] cass: fix forward computation of easting (fixes #3432), and use iterative generic inversion to make GIGS 5108 roundtripping tests succeed The error on the sign after "nu * A * (1. - A2 * T * (C1 ...") dates back to PROJ origins. All other independent implementations (Apache SIS cf /~https://github.com/apache/sis/blob/a04bd4b29816c91fc4dfc66135b8de5e962933a0/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/CassiniSoldner.java, CSMAP cf /~https://github.com/normanb/CS-Map/blob/1d8f2db7f13194a99dc218d4cd91754c64a2a774/Source/CS_csini.c) use the same formulas as EPSG guidance note which is consistent with Snyder's Map projections - A working manual). So align on them, so that everyone is consistent (and hopefully right). While doing it, I've modified the inverse method to use the generic inversion method to improve its accuracy, which make it possible to enable the 5108.gie tests that were previously failing on the 1000 iterations roundtripping. --- src/generic_inverse.cpp | 4 ++-- src/proj_internal.h | 2 +- src/projections/adams.cpp | 11 ++++++++--- src/projections/cass.cpp | 16 +++++++-------- src/projections/wink2.cpp | 3 ++- test/CMakeLists.txt | 2 +- test/gie/builtins.gie | 14 ++++++++++++- test/gigs/{5108.gie.failing => 5108.gie} | 25 ++++++++++++------------ 8 files changed, 48 insertions(+), 29 deletions(-) rename test/gigs/{5108.gie.failing => 5108.gie} (90%) diff --git a/src/generic_inverse.cpp b/src/generic_inverse.cpp index 28d939fb52..2f8b745c9f 100644 --- a/src/generic_inverse.cpp +++ b/src/generic_inverse.cpp @@ -41,7 +41,7 @@ * * Starts with initial guess provided by user in lpInitial */ -PJ_LP pj_generic_inverse_2d(PJ_XY xy, PJ *P, PJ_LP lpInitial) { +PJ_LP pj_generic_inverse_2d(PJ_XY xy, PJ *P, PJ_LP lpInitial, double deltaXYTolerance) { PJ_LP lp = lpInitial; double deriv_lam_X = 0; double deriv_lam_Y = 0; @@ -51,7 +51,7 @@ PJ_LP pj_generic_inverse_2d(PJ_XY xy, PJ *P, PJ_LP lpInitial) { PJ_XY xyApprox = P->fwd(lp, P); const double deltaX = xyApprox.x - xy.x; const double deltaY = xyApprox.y - xy.y; - if (fabs(deltaX) < 1e-10 && fabs(deltaY) < 1e-10) { + if (fabs(deltaX) < deltaXYTolerance && fabs(deltaY) < deltaXYTolerance) { return lp; } diff --git a/src/proj_internal.h b/src/proj_internal.h index 2ff8173134..b372e1ee5d 100644 --- a/src/proj_internal.h +++ b/src/proj_internal.h @@ -875,7 +875,7 @@ void pj_clear_vgridshift_knowngrids_cache(); void pj_clear_sqlite_cache(); -PJ_LP pj_generic_inverse_2d(PJ_XY xy, PJ *P, PJ_LP lpInitial); +PJ_LP pj_generic_inverse_2d(PJ_XY xy, PJ *P, PJ_LP lpInitial, double deltaXYTolerance); diff --git a/src/projections/adams.cpp b/src/projections/adams.cpp index 633c82bee8..9d40223372 100644 --- a/src/projections/adams.cpp +++ b/src/projections/adams.cpp @@ -296,7 +296,9 @@ static PJ_LP adams_inverse(PJ_XY xy, PJ *P) lp.lam = fabs(lp.phi) >= M_HALFPI ? 0 : std::max(std::min(xy.x / 2.62205760 / cos(lp.phi), 1.0), -1.0) * M_PI; - return pj_generic_inverse_2d(xy, P, lp); + + constexpr double deltaXYTolerance = 1e-10; + return pj_generic_inverse_2d(xy, P, lp, deltaXYTolerance); } static PJ_LP peirce_q_square_inverse(PJ_XY xy, PJ *P) @@ -339,7 +341,9 @@ static PJ_LP peirce_q_square_inverse(PJ_XY xy, PJ *P) } else /* if( xy.x <= 0 && xy.y <= 0 ) */ lp.lam = -M_PI / 2; - return pj_generic_inverse_2d(xy, P, lp); + + constexpr double deltaXYTolerance = 1e-10; + return pj_generic_inverse_2d(xy, P, lp, deltaXYTolerance); } static PJ_LP peirce_q_diamond_inverse(PJ_XY xy, PJ *P) @@ -382,7 +386,8 @@ static PJ_LP peirce_q_diamond_inverse(PJ_XY xy, PJ *P) lp.phi = -M_PI / 4; } - return pj_generic_inverse_2d(xy, P, lp); + constexpr double deltaXYTolerance = 1e-10; + return pj_generic_inverse_2d(xy, P, lp, deltaXYTolerance); } static PJ *setup(PJ *P, projection_type mode) { diff --git a/src/projections/cass.cpp b/src/projections/cass.cpp index 1e8e4438e6..40baa8f0ca 100644 --- a/src/projections/cass.cpp +++ b/src/projections/cass.cpp @@ -43,7 +43,7 @@ static PJ_XY cass_e_forward (PJ_LP lp, PJ *P) { /* Ellipsoidal, forward const double A2 = A * A; xy.x = nu * A * (1. - A2 * T * - (C1 - (8. - T + 8. * C) * A2 * C2)); + (C1 + (8. - T + 8. * C) * A2 * C2)); xy.y = M - Q->m0 + nu * tanphi * A2 * (.5 + (5. - T + 6. * C) * A2 * C3); if( Q->hyperbolic ) @@ -82,13 +82,13 @@ static PJ_LP cass_e_inverse (PJ_XY xy, PJ *P) { /* Ellipsoidal, inverse lp.lam = D * (1. + T1 * D2 * (-C4 + (1. + 3. * T1) * D2 * C5)) / cos (phi1); - if( Q->hyperbolic ) - { - // EPSG guidance note 7-2 suggests a custom approximation for the - // 'Vanua Levu 1915 / Vanua Levu Grid' case, but better use the - // generic inversion method - lp = pj_generic_inverse_2d(xy, P, lp); - } + // EPSG guidance note 7-2 suggests a custom approximation for the + // 'Vanua Levu 1915 / Vanua Levu Grid' case, but better use the + // generic inversion method + // Actually use it in the non-hyperbolic case. It enables to make the + // 5108.gie roundtripping tests to success, with at most 2 iterations. + constexpr double deltaXYTolerance = 1e-12; + lp = pj_generic_inverse_2d(xy, P, lp, deltaXYTolerance); return lp; } diff --git a/src/projections/wink2.cpp b/src/projections/wink2.cpp index 990c9686cb..e4339ab3be 100644 --- a/src/projections/wink2.cpp +++ b/src/projections/wink2.cpp @@ -48,7 +48,8 @@ static PJ_LP wink2_s_inverse(PJ_XY xy, PJ *P) lpInit.phi = xy.y; lpInit.lam = xy.x; - return pj_generic_inverse_2d(xy, P, lpInit); + constexpr double deltaXYTolerance = 1e-10; + return pj_generic_inverse_2d(xy, P, lpInit, deltaXYTolerance); } diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 8490069253..c0a80cf456 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -80,7 +80,7 @@ proj_add_gie_test("GIGS-5104" "gigs/5104.gie") proj_add_gie_test("GIGS-5105.2" "gigs/5105.2.gie") proj_add_gie_test("GIGS-5106" "gigs/5106.gie") proj_add_gie_test("GIGS-5107" "gigs/5107.gie") -#proj_add_gie_test("GIGS-5108" "gigs/5108.gie") +proj_add_gie_test("GIGS-5108" "gigs/5108.gie") proj_add_gie_test("GIGS-5109" "gigs/5109.gie") #proj_add_gie_test("GIGS-5110" "gigs/5110.gie") proj_add_gie_test("GIGS-5111.1" "gigs/5111.1.gie") diff --git a/test/gie/builtins.gie b/test/gie/builtins.gie index 1c36d2e627..567abea0f3 100644 --- a/test/gie/builtins.gie +++ b/test/gie/builtins.gie @@ -850,6 +850,18 @@ expect -0.001790493 0.000895247 accept -200 -100 expect -0.001790493 -0.000895247 +------------------------------------------------------------------------------- +# test point from EPSG Guidance Note 7.2 +operation +proj=cass +lat_0=10.4416666666667 +lon_0=-61.3333333333333 \ + +x_0=86501.46392052 +y_0=65379.0134283 \ + +a=6378293.64520876 +b=6356617.98767984 +to_meter=0.201166195164 +------------------------------------------------------------------------------- + +tolerance 0.1 mm +accept -62 10 +expect 66644.94040882 82536.21873655 +roundtrip 1 + ------------------------------------------------------------------------------- # Hyperbolic variant: test point from EPSG Guidance Note 7.2 operation +proj=cass +hyperbolic +a=6378306.376305601 +rf=293.466307 \ @@ -859,7 +871,7 @@ operation +proj=cass +hyperbolic +a=6378306.376305601 +rf=293.466307 \ tolerance 0.1 mm accept 179.99433652777776 -16.841456527777776 -expect 16015.289017555102 13369.660053668682 +expect 16015.28901692 13369.66005367 roundtrip 1 =============================================================================== diff --git a/test/gigs/5108.gie.failing b/test/gigs/5108.gie similarity index 90% rename from test/gigs/5108.gie.failing rename to test/gigs/5108.gie index 15ecc3cfd2..2abdfc9d3c 100644 --- a/test/gigs/5108.gie.failing +++ b/test/gigs/5108.gie @@ -4,18 +4,19 @@ Test 5108, Cassini-Soldner, v2-0_2011-06-28. -------------------------------------------------------------------------------- -# GDM2000 -<4742> +proj=longlat +ellps=GRS80 <> + -# GDM2000 / Johor Grid -<3377> +proj=cass +lat_0=2.121679744444445 +lon_0=103.4279362361111 +x_0=-14810.562 +y_0=8758.32 +ellps=GRS80 +units=m <> +use_proj4_init_rules true +# GDM2000 +# <4742> +proj=longlat +ellps=GRS80 <> - +# GDM2000 / Johor Grid +# <3377> +proj=cass +lat_0=2.121679744444445 +lon_0=103.4279362361111 +x_0=-14810.562 +y_0=8758.32 +ellps=GRS80 +units=m <> -------------------------------------------------------------------------------- -operation +proj=pipeline - +step +init=epsg:4742 +inv +operation +proj=pipeline \ + +step +init=epsg:4742 +inv \ +step +init=epsg:3377 -------------------------------------------------------------------------------- tolerance 0.05 m @@ -87,8 +88,8 @@ accept 104 5 expect 48630.563 327067.097 -------------------------------------------------------------------------------- -operation +proj=pipeline - +step +init=epsg:3377 +inv +operation +proj=pipeline \ + +step +init=epsg:3377 +inv \ +step +init=epsg:4742 -------------------------------------------------------------------------------- tolerance 0.05 m @@ -160,8 +161,8 @@ accept 48630.563 327067.097 expect 104 5 -------------------------------------------------------------------------------- -operation +proj=pipeline - +step +init=epsg:4742 +inv +operation +proj=pipeline \ + +step +init=epsg:4742 +inv \ +step +init=epsg:3377 -------------------------------------------------------------------------------- tolerance 0.006 m @@ -232,4 +233,4 @@ tolerance 0.006 m accept 104 5 roundtrip 1000 - +