Inconsistency / doc issue in conventions Coordinate Frame vs Position Vector for +towgs84 vs +proj=helmert

[rouault]

The +towgs84 convention used by Proj is to have coefficients in the Position Vector / EPSG:9606 convention: cf https://github.com/OSGeo/proj.4/wiki/GenParms and pj_geocentric_to_wgs84() using the RPV matrix of the EPSG Guidance 7.2 page 136.

https://proj4.org/operations/transformations/helmert.html mentions "Two conventions exists and they seem to be equally popular. PROJ uses the Position Vector rotation convention. The rotation matrix can be transposed by adding the +transpose flag in the transformation setup which makes PROJ follow the Coordinate Frame rotation convention.", so this means that Helmert is supposed to use Position Vector convention too. However the formulas of build_rot_matrix() follow the R matrix (for the exact case) and RCF (for the approximate one) of the EPSG Guidance, that is to see Coordinate Frame.

Relevant extract from EPSG 7.2 guidance:

Extract from build_rot_matrix()

    f = Q->opk.o;  // that is RX
    t = Q->opk.p;  // that is RY
    p = Q->opk.k;  // that is RZ

// Approximate case:
        R00 =  1;        R01 =  p;        R02 = -t;
        R10 = -p;        R11 =  1;        R12 =  f;
        R20 =  t;        R21 = -f;          R22 =  1;   --> matches RCF / Coordinate Frame

Extract from pj_geocentric_to_wgs84()

            x_out = M_BF*(       x[io] - Rz_BF*y[io] + Ry_BF*z[io]) + Dx_BF;
            y_out = M_BF*( Rz_BF*x[io] +       y[io] - Rx_BF*z[io]) + Dy_BF;
            z_out = M_BF*(-Ry_BF*x[io] + Rx_BF*y[io] +       z[io]) + Dz_BF;  --> matches RPV / Position Vetor

And https://proj4.org/usage/transformation.html mentions rightly "the default setting is to use “coordinate frame” convention of the Helmert transform"

Another practical proof that there is a difference of conventions between +towgs84 and +helmert

$ echo "2 49 0" | src/cs2cs  -f "%.8f"  +proj=longlat +ellps=WGS84 +towgs84=1,2,3,4,5,6,7  +to   +proj=longlat +datum=WGS84
2.00036925	48.99866109 47.03032577

To get the same result with the new API, you need to use the +transpose flag of helmert

$ echo "2 49 0" | src/cct +proj=pipeline +step +proj=cart +ellps=WGS84  +step +proj=helmert +transpose +x=1 +y=2 +z=3 +rx=4 +ry=5 +rz=6 +s=7 +step +inv +proj=cart +ellps=WGS84
  2.0003692519   48.9986610886       47.0303           inf

And the cs2cs_emulation_setup() of proj_4D_api.c does add the 'transpose' term when it expands +towgs84 into a helmert transform

So I'm pretty confident that there is indeed a difference of conventions + documentation error in the Helmert page

I can see 2 possibilities:

• change Helmert behaviour to use position vector convention as documented in its page, fix regression cases, and fix usage/transformation.html to reflect that. Pro: same convention for +towgs84 / +helmert. Con: will break early users of Helmert
• keep Helmert implementation as it, but fix its doc to reflect what it is actually doing

Bonus point, that can apply to any of those 2 choices, add a new flag to helmert: +convention=position_vector/coordinate_frame flag so that it is immediately obvious which convention is used (+transpose is really technical and assumes you know the default convention)

[rouault]

@kbevers @busstoptaktik Opinions ? This is something we must resolve before 5.2.0 IHMO

[hobu]

change Helmert behaviour to use position vector convention as documented in its page, fix regression cases, and fix usage/transformation.html to reflect that. Pro: same convention for +towgs84 / +helmert. Con: will break early users of Helmert

I vote for this unless it breaks a long-standing convention. Otherwise we'll be doomed to answer this as a FAQ forever, no?

[busstoptaktik]

What an unholy mess - I have a suggestion for a potentially safe way of solving this (basically @rouault s "bonus point" slightly extended), but will have to go offline for a few hours now. Will be back with a suggestion/description in the afternoon (UTC)

[busstoptaktik]

@rouault: It was definitely not the intention to make proj=helmert incompatible with +towgs84.

Your discovery of the incompatibility, and the fact that no one else have noticed for more than a year stresses that this problem is hard - not because it is complex, but because it is intractable.

It is akin to the linguistic problem of repeated negations (“I ain’t gonna give nobody none of ...”). Unwrapping the layers to get to the true meaning is intractable because we can do full or partial sign switching of the transformation through a number of concepts, each with their related conventional, but potentially ambiguous interpretation:

• Direction: What’s forward, what’s inverse
• Transposition/Matrix representation: Row major vs column major (C vs Fortran)
• Position vector vs. reference frame rotation
• Sign convention (mathematical vs. navigational)

While the latter is probably not a problem in general, the convention of giving rotations in arcsec (hinting at navigational convention), rather than e.g. microradians, makes it a source of potential confusion anyway.

But let’s ignore the sign convention problem - and the even worse discussion of the order of the rotations. That still leaves us with at least three ways of stating parts of the same numerical thing.

When doing a test computation, using a known data point, it is much too easy to convince oneself that “I have probably forgotten a transposition”, or “I probably confused forward and inverse”, and make things work correctly numerically by adding the appropriate, rather than correcting the root cause.

So far I have just been stating the obvious, but since the suggested solution is painfully verbose, I want to make it perfectly clear that that level of verbosity really is the lesser evil.

With that said, I suggest introducing 2 changes:

1. Forbid the use of +transpose: Make proj=helmert fail if +transpose is given
2. Require the use of +convention: Make proj=helmert fail if +convention=... is not given
3. Repair regressions

i.e. restore compatibility between +towgs84 and +proj=helmert, by not providing a default +proj=helmert convention. Neither PV nor CF.

This allows us to fail noisily and early on use of the old syntax, either by choking on +transpose or by detecting the lack of +convention=...

It is also easy to communicate to users: Everywhere you use +transpose, replace by +convention=position_vector, everywhere else, insert +convention=coordinate_frame.

All in all it is a simple change (+convention=position_vector replaces +transpose, while +convention=coordinate_frame becomes a required no-op), and, if I understand correctly, covering the opinions of both @rouault and @hobu. My schedule is, however, overbooked for the rest of this week, so while I agree it should be done before 5.2.0, I am not volunteering to implement it immediately.

[rouault]

@busstoptaktik Your plan sounds good to me. I can take care of implementing it if other stakeholders are happy with it too.

[kbevers]

I too am okay with Thomas' solution. One question though: When using the +towgs84 option with +proj=helmert, is the Position Vector convention then implied?

This is mostly a reminder for myself: Once this is implemented I need to update the init files at https://github.com/NordicGeodesy/NordicTransformations and sync with proj-datumgrids. The ITRF init-files distributed with PROJ will also need to be updated (can't remember if those are covered by the test suite).

[rouault]

When using the +towgs84 option with +proj=helmert

Do you mean like in "+proj=helmert +x= +y= ... +towgs84=...." ? Huh, is it possible ? I'd say this is undefined behaviour. I hope nobody will have the crazy idea of using both

[kbevers]

Do you mean like in "+proj=helmert +x= +y= ... +towgs84=...." ? Huh, is it possible ?

Sure is:

https://github.com/OSGeo/proj.4/blob/c1f0673b3335a37eb03900315d6db3f43e3bf64c/src/PJ_helmert.c#L492-L509

I believe it is used in cs2cs_emulation_setup().

[rouault]

I believe it is used in cs2cs_emulation_setup().

Indeed. This is not documented in helmert.rst and is probably a good thing to avoid any confusion, so we can deal with that as an internal detail. I'd say that if +towgs84 is defined we should check that +convention=position_vector and error out otherwise (such error shouldn't happen if cs2cs_emulation_setup() is the only user of it, once it is updated)

[phidrho]

Hi,
I also noticed some errors in documentation few months ago, but didn't report it immediately. I was busy doing testing and wanted to make sure that it wasn't me that I'm doing something wrong. Transformation parameters which I used did not have information about which rotation convention should be used so I contacted State Geodetic Department of Croatia to get all information, and they responded after a month... so the whole process lasted... I was testing new version of proj (5.0.1) with some data that I could not share publically so it was another problem. In the end I solved my problem and moved on to work on another projects... I contacted Thomas directly with email, because his email was written in source "helmert.c". I promise I'll report bugs immediately next time.

I still have my notes about errors (bugs) that I found, some of which are maybe fixed in the meanwhile, these are the ones that are Helmert transformation related, I'll report other ones in new issue:

1. bug in documentation about Helmert transformation: "+transpose" uses coordinate frame convention - which is not truth, it uses coordinate frame convention by default, cs2cs uses position vector convention by default
2. bug in documentation about Helmert transformation: in example with ITRF93@2017.0 it says "15 parameters" there should be "14 parameters"
3. misspelled word in sentence "Transpose rotation matrix and follow the Cordinate Frame rotation convention." should be corrected to "Transpose rotation matrix and follow the Coordinate Frame rotation convention."

I also found a missing option (feature) for Helmert transformation if we want it to be truly flexible to be used in any country. I wrote about this problem in second email to Thomas, but then I didn't understand the problem completely. For Croatia, we use different rotation matrix, at first I thought it was bug in PROJ, but later I found out that it is programmed for most common rotation matrix, which is obviously calculated as R = Rx * Ry * Rz. In Croatian literature: Issues related to Helmert seven-parameter transformation, there is a formula which is calculated as R = Rz * Ry * Rx, and it is noted in paper that multiplication must not be reversed for our case. I thought that only difference was because of different rotation convention, but if you multiply it different you don't only get reversed signs, but different matrix. As I researched it further I found this Swedish paper On geodetic transformations, there is explanation in the chapter 6.3 - Computational consistency for 3D Helmert, on page 16 just below middle of the page:

The next problem concerns the definition of rotation matrix R. So far, we have assumed that the multiplication order of the three partial matrices is RzRyRx. Theoretically, there are six different ways of multiplying the matrices. At least the following two variations have occurred in different international contexts, RxRyRz and RyRxRz. The effect of using the wrong order of rotation is small. As long as the rotations remain small (<10 arc seconds) it amounts to some millimetres, but here as well the error grows quadratically with the size of the rotations. For rotations at the level of 1 degree it is a question of hundreds of metres.

So the proper way to implement Helmert transformation should be to implement option for different rotation matrices (6 versions, as stated in Swedish paper), and the already implemented - two rotation conventions.

[rouault]

@phidrho Your point about the order of multiplication of rotation matrices is an interesting one (and quite scary !), and I believe deserve another dedicated ticket so that one remains actionable. I guess this only matters for the case where the +exact flag is used, since with the default small angle approximation (sin(x)=x, cos(x)=1), the order of matrices doesn't matter.