+proj=push/pop: add a +bank=<name> parameter to interleave push/pop operations

[rouault]

The motivation is to be able to do transformation that first do a geoid transformation followed by a Helmert one.

The following example shows a ETRS89 to S-JTSK/05 + Baltic 1957 height transformation that chains a geoid model registered against ETRS89 with a Helmert transformation.

       +proj=pipeline \
            +step +proj=push +v_3 +omit_inv \                         # Save ETRS89 ellipsoidal height
            +step +inv +proj=vgridshift +grids=CR2005.tif +multiplier=1 \  # Apply geoid
            +step +proj=push +v_3 +bank=baltic_height \               # Save Baltic 1957 height
            +step +proj=pop +v_3 +omit_inv \                          # Restore ETRS89 ellipsoidal height
            +step +proj=cart +ellps=GRS80 \                           # Helmert transformation
            +step +inv +proj=helmert +x=572.213 +y=85.334 +z=461.94 \
                +rx=-4.9732 +ry=-1.529 +rz=-5.2484 +s=3.5378 +convention=coordinate_frame \
            +step +inv +proj=cart +ellps=bessel \
            +step +proj=pop +v_3 +bank=baltic_height                  # Restore Baltic 1957 height

Without that, the current alternative is a ugly workaround involving using the v_4 component and axisswap

       +proj=pipeline \
            +step +proj=push +v_4 \                                   # Save v_4 as we are going to use it as a temp variable
            +step +proj=push +v_3 +omit_inv \                         # Save ETRS89 ellipsoidal height
            +step +inv +proj=vgridshift +grids=CR2005.tif +multiplier=1 \  # Apply geoid
            +step +proj=push +v_3 +omit_fwd \                         # On reverse path, restore initial Baltic height
            +step +proj=axisswap +order=1,2,4,3 +omit_inv \           # On forward path, save Baltic height in v_4 component...
            +step +proj=pop +v_3 +omit_inv \                          # On forward parth, restore initial ellipsoidal height
            +step +proj=cart +ellps=GRS80 \                           # Helmert transformation
            +step +inv +proj=helmert +x=572.213 +y=85.334 +z=461.94 \
                +rx=-4.9732 +ry=-1.529 +rz=-5.2484 +s=3.5378 +convention=coordinate_frame \
            +step +inv +proj=cart +ellps=bessel \
            +step +proj=axisswap +order=1,2,4,3 +omit_inv \           # On forward path, restore Baltic height from v_4 component...
            +step +proj=pop +v_3 +omit_fwd \                          # On reverse path, save initial Baltic height
            +step +proj=pop +v_4

[jjimenezshaw]

If I understand it correctly, each named bank has a single element (std::array<4>), it is not a stack. Right? The element is not deleted from the bank after a pop. (I say this because the unnamed pop and push are going to a "real" stack). I am not asking that it should be a stack. Just that it was my first impression. Maybe it should be documented.

[rouault]

If I understand it correctly, each named bank has a single element (std::array<4>), it is not a stack

yes, good point. That's also something that makes me slightly uneasy as it goes a bit against the push/pop stack philosophy. Perhaps instead of hijacking push/pop, we should introduce a save/load as I proposed in one of my emails?

+proj=save +v_1 +bank=foo / +prop=load +v_1 +bank=foo

[busstoptaktik]

That is way better. But I think I have an idea for an even more clean design, which will be entirely orthogonal to the current unfortunate design of the push/pop mechanism, and handle both stack and register based "side storage" in a better abstracted way.

It will allow us to deprecate the current design, and have a stack mechanism which is actually a stack mechanism, and a register mechanism which is actually a register mechanism.

I'm writing this up as a proposal right now - will post here and on the related mailing list thread in a few hours.

The stack/register mish-mash in this PR, I'm not happy with - it's papering over an already unfortunate design with even more unfortunate design.

Personally, I do not see the need for register-based side storage: Repairing the stack mechanism would do. But I understand that there is a perceived need for this, so I'm going for a design including both, but with an extra layer of abstraction which allows for a cleaner syntax, better extensibility, and full three way orthogonality between the new stack mechanism, the new register mechanism, and the existing (and hopefully soon deprecated) push/pop implementation.

[rouault]

I'm writing this up as a proposal right now

eagerly waiting for it.

Thinking at the topic, I had a look back at the ISO-19111 "Pass-through coordinate operation" (https://docs.ogc.org/as/18-005r4/18-005r4.html#103) concept, that is somehow related, but probably not quite advanced enough (I didn't implement PassthroughOperation during the ISO-19111'ification of PROJ), because just a list of the index of components modified by an operation doesn't fully address the need to do the geoid operation and Helmert one "in parallel", and merge the results at a later stage, because it requires that the 2 operations operate on separate coordinate indices. That is it would work for a vertical-only + horizontal-only operation. But here the horizontal part is actually a 3D Helmert for which we only keep the horizonal-part. So something more convoluted...

[busstoptaktik]

So something more convoluted...

Yes, this is really one of the problems with 19111: Real geodesy is always more convoluted than what the standard is able to cover. 19111 can be drastically simplified by dropping the geodetically unrealistic misconception of "a CRS defined by internal state", accepting that a CRS is merely a label, enabling us to look up related transformations in a register, and that 19111 should focus on an improved representation of transformations by essentially codifying the workings of the pipeline operator

[busstoptaktik]

eagerly waiting for it.

@rouault - here it is

The main problem with PROJ's current push/pop functionality is that the four-stack implementation makes it impossible to push one dimension of a coordinate-tuple and pop it into another, without an elaborate axisswap dance at both ends of the push/pop-pair.

This is a consequence of the v_1 ... v_4 syntax used: Since the order of operator arguments is not necessarily preserved during the parsing process, a pair specified as v_1 v_2 may appear as v_2 v_1 at parsing time.

A more axisswap-like syntax would have eliminated this - e.g.

proj=push items=1,3

to push dimensions 1 and 3 of the coordinate tuple onto the
stack, and

proj=pop items=3,1

to pop them back into swapped dimensions.

But since the obvious names of push and pop are now taken, let's make a virtue out of necessity by adding a little more encapsulation by introducing the "stack subsystem", giving this syntax:

proj=stack push=1,3

and

proj=stack pop=3,1

The implementation will still mimick that of push/pop (i.e being handled as part of the pipeline machinery, rather than as ordinary proj=foo-operators), but it will not stomp on the existing unfortunate push/pop-implementation, and it leaves room for extension, e.g. by introducing the swap and flip operators mused over earlier in this thread:

proj=stack push=1,3
proj=stack swap
proj=stack pop=1,3

to essentially do the same as above and:

proj=stack push=1,3
proj=stack flip=3,1

although I believe the need for these will be limited, once we have a push/pop pair actually working like a stack.

The wish for named registers mentioned earlier in the thread may be introduced using a similar syntax, where:

proj=reg save=1,3 into=first,third
proj=reg restore=1,3 from=third,first

will do the same swapping of first and third dimension as above, but using named registers instead of stack slots.

I see the potential for improved clarity when using named registers instead of stack slots, by using meaningful names for the registers (e.g. ellipsoidal_height and orthometric_height, to take the obvious use case).

The problem, however, is when trying to interpret the workings of a pipeline running in inverse mode:

In the inverse stack-case, we only need to flip the abstract concepts of stack slots, and swap push for pop.

Whereas in the register case, we also need to "abstractly swap the concrete concepts behind the concrete register names" which, at least to my imperfect brain, is much harder.

But while not personally seeing the need for registers, I acknowledge the wish for them, and to make this work in a fairly backwards compatible way, I suggest the solution outlined above, i.e. two new, and well encapsulated subsystems in the pipeline implementation. And once they are well established: Deprecation of the current push/pop pair.

[busstoptaktik]

And with a resulting pipeline somewhat like this (with the caveat that I am not sure I understand the functionality completely):

proj=pipeline 
    step proj=stack push=3                      # Save ETRS89 ellipsoidal height
    step inv proj=vgridshift grids=CR2005.tif   # Transform to orthometric height? (why inv here?)
    step proj=stack flip=3                      # Save Baltic 1957 height, restore ETRS89
    step proj=cart ellps=GRS80                  # Helmert transformation
    step inv proj=helmert
             x=572.213  rx=-4.9732
             y=85.334   ry=-1.529
             z=461.94   rz=-5.2484
             s=3.5378   convention=coordinate_frame
    step inv proj=cart ellps=bessel 
    step proj=stack pop=3                       # Restore Baltic 1957 height

[rouault]

Transform to orthometric height? (why inv here?)

By default "proj=vgridshift" subtracts the grid value (the default value of multiplier is -1), which is a bit unusual. A geoid grid contains values to add to orthometric height to get ellipsoidal heights. So here as we want to go from ellipsoidal height to orthometric, that is subtract the grid value, we need to use either "proj=vgridshift" or "proj=vgridshift +multiplier=1 +inv" which are equivalent

And with a resulting pipeline somewhat like this (with the caveat that I am not sure I understand the functionality completely):

That pipeline would likely work in the forward path, but not in the reverse path. Since the inversion of the first step " step proj=stack push=3" would restore a "approximate" ETRS89 ellipsoidal height computed by helmert. So in the reverse path, we need to restore the Baltic 57 height before the execution of vgridshift

[busstoptaktik]

By default "proj=vgridshift" subtracts the grid value (the default value of multiplier is -1), which is a bit unusual

I would not call that unusual, it's just the common geodetic convention: The geoid model represents the geoid separation, i.e. the height of the geoid wrt. the ellipsoid, so to make geometrical heights, referred to the ellipsoid, compatible with heights obtained by levelling in a different, orthometric system, we subtract the value of the geoid model, expressed in the geometrical system (in this case ETRS89).

And as PROJ is geometry-centric, the natural forward direction operation of vgridshift is to go from geometric to orthometric - so according to the comments indicating ETRS->Baltic57, i expected proj=vgridshift grids=...: No inv and no multiplier, just the plain forward operation - and with and inv, I would have expected the multiplier to be -1. So just a bit surprised here over the sign convention.

[rouault]

Related discussion thread at https://lists.osgeo.org/pipermail/proj/2024-February/thread.html#11251

[rouault]

Superseded per #4036