Pipeline with +towgs84 doesn't work as expected

[rouault]

Trying to port GDAL to use new API, I'm doing some unexpected results when doing supposedly no-op transformations involving towgs84

The output of the following code

#include "proj.h"
#include <math.h>
#include <stdio.h>

int main(int argc, char* argv[])
{
    PJ* pj = proj_create(0, "+proj=pipeline +step "
                            "+proj=longlat +ellps=intl +towgs84=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68 "
                            "+inv +step "
                            "+proj=longlat +ellps=intl +towgs84=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68");
    double x[] = { 2. / 180 * M_PI };
    double y[] = { 49. /180 * M_PI };
    double z[] = { 0 };
    printf("%.18g %.18g %.18g\n", x[0], y[0], z[0]);
    proj_trans_generic (pj, PJ_FWD,
                        x, sizeof(double), 1,
                        y, sizeof(double), 1,
                        z, sizeof(double), 1,
                        0, 0, 0);
    printf("%.18g %.18g %.18g\n", x[0], y[0], z[0]);
    return 0;
}

is

0.0349065850398865909 0.855211333477221336 0
0.0349065851566477472 0.855211333481314284 0.00130653195083141327

Whereas the second line corresponding to the output should be identical to the first one
Doing the same using the proj_api.h with pj_transform() results in the expected result

[rouault]

Can be reproduced with cct

$ echo "2 49 0" | src/cct +proj=pipeline +step +proj=longlat +ellps=intl +towgs84=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68 +inv +step +proj=longlat +ellps=intl +towgs84=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68
  2.0000000067   49.0000000002        0.0013        0.0000

Interestingly it seems the issue is specific to a 7-term towgs84 transform. With just 3-terms it is OK

$ echo "2 49 0" | src/cct +proj=pipeline +step +proj=longlat +ellps=intl +towgs84=-104.1,-49.1,-9.9 +inv +step +proj=longlat +ellps=intl +towgs84=-104.1,-49.1,-9.9
  2.0000000000   49.0000000000        0.0000        0.0000

[kbevers]

Possibly related to this question on the mailing list. Answer by @busstoptaktik.

It would be nice to find a way for this to work. I know that Thomas struggled with it quite a lot so it might not be easy to find a solution.

I Realise that this probably wont't solve your problem right now but try the following. At least for the sake of debugging this:

1. Try doing the transformation with +proj=helmert instead.
2. Try wrapping the steps in an init-file.

[rouault]

Or maybe it's just a precision issue ?

$ echo "2 49" | src/cct +proj=pipeline +step +proj=longlat +ellps=intl +towgs84=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68 +inv +step +proj=longlat +datum=WGS84
  1.9993016680   49.0006659752       48.1311        0.0000

and

$ echo "1.9993016680   49.0006659752       48.1311" |  src/cct -I +proj=pipeline +step +proj=longlat +ellps=intl +towgs84=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68 +inv +step +proj=longlat +datum=WGS84
  2.0000000067   49.0000000003        0.0013        0.0000

which is equivalent to

$ echo "1.9993016680   49.0006659752       48.1311" |  src/cct  +proj=pipeline +step +proj=longlat +datum=WGS84  +inv +step  +proj=longlat +ellps=intl +towgs84=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68
  2.0000000067   49.0000000003        0.0013        0.0000

I see that in the old pj_transform.c there's a pj_compare_datums() that avoids any geodetic <--> geocentric transformation if the ellipsoids and 7-params of the source and targets are the same.

Reading cct help, I see roundtrip accuracy can be tested with this syntax (I've a hard time to understand it... the +proj=longlat without a preceding +step, and the ending +step +step +inv ?)

echo "2 49" | src/cct +proj=pipeline +proj=longlat +ellps=intl +towgs84=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68 +step +step +inv
  2.0000000067   49.0000000002        0.0013        0.0000

[kbevers]

Or maybe it's just a precision issue ?

That's very likely, yes. I haven't had time to look into this myself yet.

Reading cct help, I see roundtrip accuracy can be tested with this syntax (I've a hard time to understand it... the +proj=longlat without a preceding +step, and the ending +step +step +inv ?)

echo "2 49" | src/cct +proj=pipeline +proj=longlat +ellps=intl +towgs84=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68 +step +step +inv
2.0000000067   49.0000000002        0.0013        0.0000

Yeah, that's a bizarro use of global pipeline parameters. It is the same mechanism as

+proj=pipeline +ellps=GRS80 +step +proj=cart +step +proj=helmert  ...  +step +proj=cart +inv

where +ellps=GRS80 is substituted into each step. In the above case it is +proj=longlat +ellps=intl +towgs84=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68 that is substituted into the otherwise empty steps. It is wierd, unintuitive and shouldn't be used as an example but it does follows the rules of the pipeline.

The global pipeline parameters are defined as everything between +proj=pipeline and the first occurrence of +step, as described in the docs.

[kbevers]

Precision issue confirmed. The patch below fixes the problem, but at the expense of slower computation:

diff --git a/src/proj_4D_api.c b/src/proj_4D_api.c
index 58b3d835..9a2da55e 100644
--- a/src/proj_4D_api.c
+++ b/src/proj_4D_api.c
@@ -490,7 +490,7 @@ Returns 1 on success, 0 on failure
         def = malloc (100+n);
         if (0==def)
             return 0;
-        sprintf (def, "break_cs2cs_recursion     proj=helmert %s transpose", s);
+        sprintf (def, "break_cs2cs_recursion     proj=helmert exact  %s transpose", s);
         Q = proj_create (P->ctx, def);
         pj_inherit_ellipsoid_def (P, Q);
         free (def);

It might be possible to find a clever way to detect when the exact flag is needed, so that the approximate Helmert can be used for the normal cases.

[kbevers]

Some further investigation:

$ cat test.gie
<gie>

operation +proj=pipeline +proj=longlat +ellps=intl +towgs84=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68 +step +step +inv

accept 2 49
roundtrip 9000000 1 cm

</gie>

$ time bin/gie test.gie # without patch
-------------------------------------------------------------------------------
Reading file 'test.gie'
-------------------------------------------------------------------------------
proj=pipeline proj=longlat ellps=intl towgs84=-104.1,-49.1,-9.9,0.971,...
-------------------------------------------------------------------------------
     FAILURE in test.gie(6):
     roundtrip deviation: 25229091.721442 mm, expected: 10.000000 mm
-------------------------------------------------------------------------------
total:  0 tests succeeded,  0 tests skipped,  1 tests FAILED!
-------------------------------------------------------------------------------

real    0m12.393s
user    0m12.360s
sys     0m0.015s

$ time bin/gie test.gie # with patch
-------------------------------------------------------------------------------
Reading file 'test.gie'
-------------------------------------------------------------------------------
total:  1 tests succeeded,  0 tests skipped,  0 tests failed.
-------------------------------------------------------------------------------

real    0m12.346s
user    0m12.324s
sys     0m0.011s

Conclusion: The patch in the previous comment should just be applied. There is no real downside to it, apart for a slightly more expensive PJ-setup.

[rouault]

Great finding ! I'm wondering if the helmert transformation shouldn't default to 'exact' ? And have rather an 'approximate' setting if really needed ?

[kbevers]

Great finding ! I'm wondering if the helmert transformation shouldn't default to 'exact' ? And have rather an 'approximate' setting if really needed ?

It started out that way, but there's a few things to consider. One is when your start doing transformations of coordinates where the time component is different for each coordinate. When the time component is different to the previous coordinate the transformation matrix is recalculated. This very quickly gets very expensive. Another is that the approximate version is almost universally used. I haven't seen a real world use of it - even in super accurate geodetic applications. A third reason is that this is what PROJ has been doing with +towgs84 parameters all along.

So I decided on the approximate as the default since it offers the fastest calculation and is what most (if not all) users would intuitively expect.