Pipeline coordinate stack (push and pop operations)

[kbevers]

Push and pop operations which can be helpful in some transformations.

closes #1249
closes #1227

[rouault]

• If we want to push several parameters at once, then we would have push=3 and push=4 for example, which breaks the unicity of the key name and woud likely cause problems.
• puthing push/pop as an operation parameter brings some ambiguity on when this action is done. I assume that in the forward path, push is done before, and pop after ?
• I was wondering about using x,y,z,t instead of 1,2,3,4 ? That said in the axisswap operation, we indeed use numbers

[kbevers]

• If we want to push several parameters at once, then we would have push=3 and push=4 for example, which breaks the unicity of the key name and woud likely cause problems.

I am aware of this. That is the reason for this section

https://github.com/kbevers/proj.4/blob/ab34f8c94394c4c375e202f012ab2972d0109e21/src/pipeline.cpp#L527-L536

not using pj_param() as that would indeed cause problems. I have completely ignored any ramifications this approach might have in the ISO19111 code. I did briefly consider taking a list approach, similarly to what is done in axisswap etc, but didn't because I had to engage my brain to do so. I may change it if considered a better approach. That being said, I don't think it is a very likely scenario that someone wants to store several components on the stack - it is probably only going to be useful for dealing with vertical systems in legacy transformations. That's my gut feeling anyway.

• puthing push/pop as an operation parameter brings some ambiguity on when this action is done. I assume that in the forward path, push is done before, and pop after ?

Yes. I am still playing around with this so I might change things. In any case I will make sure to describe the behaviour thoroughly in the documentation.

• I was wondering about using x,y,z,t instead of 1,2,3,4 ? That said in the axisswap operation, we indeed use numbers

So was I but as you are well aware axis ordering is tricky, so I opted for the the generic axisswap approach. Nothing is assumed here, which hopefully makes users think a little instead of running on autopilot and assume that x is always the first component of a coordinate. Also, it is nice to be consistent.

[kbevers]

As mentioned elsewhere, +push and +pop are eaten by the PROJStringFormatter (I think) when pipelines consisting of cart+helmert+invcart steps. Example below.

I am having a hard time understanding why that happens. It works alright when substituting helmert with affine as can be seen in the test material. I understand as much as this particle type of pipeline gets special treatment but I can't for the life of me figure out where that special treatment starts and how I can change it. Ideally I would mark push and pop as global keywords that are always allowed in pipelines, but that doesn't seem to be an option. At least not in a simple enough way that I can understand how :-)

@rouault can you offer a bit of guidance?

$ echo 12 56 0 2000 |./bin/cct -vvv  +proj=pipeline +step +proj=cart +ellps=GRS80 +step +proj=helmert +push=1 +x=1000 +y=2000 +z=3000 +step +proj=cart +pop=1 +ellps=GRS80 +inv

cct: Running in very verbose mode
pj_open_lib(proj.db): call fopen(/usr/local/share/proj/proj.db) - succeeded
pj_ellipsoid - final: a=6378137.000 f=1/298.257, errno=0
pj_ellipsoid - final:    ellps=GRS80
pj_ellipsoid - final: a=6378137.000 f=1/298.257, errno=0
pj_ellipsoid - final:
Pipeline: Building arg list for step no. 0
Pipeline: init - proj=unitconvert, 3
    xy_in=deg
    xy_out=rad
pj_ellipsoid - final: a=6378137.000 f=1/298.257, errno=0
pj_ellipsoid - final:    ellps=GRS80
xy_in unit: Degree
xy_out unit: Radian
Pipeline: Step 0 (proj=unitconvert) at 0x7fcaf4413860
Pipeline at [0x7fcaf4413010]:    step at [0x7fcaf4413860] (proj=unitconvert) done
Pipeline: Building arg list for step no. 1
Pipeline: init - proj=cart, 2
    ellps=GRS80
pj_ellipsoid - final: a=6378137.000 f=1/298.257, errno=0
pj_ellipsoid - final:    ellps=GRS80
Pipeline: Step 1 (proj=cart) at 0x7fcaf4500070
Pipeline at [0x7fcaf4413010]:    step at [0x7fcaf4500070] (proj=cart) done
Pipeline: Building arg list for step no. 2
Pipeline: init - proj=helmert, 4
    x=1000
    y=2000
    z=3000
pj_ellipsoid - final: a=6378137.000 f=1/298.257, errno=0
pj_ellipsoid - final:    ellps=GRS80
Helmert parameters:
x=  1000.00000  y=  2000.00000  z=  3000.00000
rx=  0.00000  ry=  0.00000  rz=  0.00000
s=   0.00000  exact=0
dx=  0.00000  dy=  0.00000  dz=  0.00000
drx= 0.00000  dry= 0.00000  drz= 0.00000
ds=  0.00000  t_epoch= 0.00000  t_obs= 0.00000
Pipeline: Step 2 (proj=helmert) at 0x7fcaf45006a0
Pipeline at [0x7fcaf4413010]:    step at [0x7fcaf45006a0] (proj=helmert) done
Pipeline: Building arg list for step no. 3
Pipeline: init - inv, 3
    proj=cart
    ellps=GRS80
pj_ellipsoid - final: a=6378137.000 f=1/298.257, errno=0
pj_ellipsoid - final:    ellps=GRS80
Pipeline: Step 3 (inv) at 0x7fcaf4500d70
Pipeline at [0x7fcaf4413010]:    step at [0x7fcaf4500d70] (inv) done
Pipeline: Building arg list for step no. 4
Pipeline: init - proj=unitconvert, 3
    xy_in=rad
    xy_out=deg
pj_ellipsoid - final: a=6378137.000 f=1/298.257, errno=0
pj_ellipsoid - final:    ellps=GRS80
xy_in unit: Radian
xy_out unit: Degree
Pipeline: Step 4 (proj=unitconvert) at 0x7fcaf45012e0
Pipeline at [0x7fcaf4413010]:    step at [0x7fcaf45012e0] (proj=unitconvert) done
Pipeline: 5 steps built. Determining i/o characteristics
Final: proj=pipeline step proj=unitconvert xy_in=deg xy_out=rad step proj=cart ellps=GRS80 step proj=helmert x=1000 y=2000 z=3000 step inv proj=cart ellps=GRS80 step proj=unitconvert xy_in=rad xy_out=deg argc=17 pargc=15
Transformation parameters for observation t_obs=2000 (t_epoch=0):
x: 1000
y: 2000
z: 3000
s: 0
rx: 0
ry: 0
rz: 0
theta: 0
Rotation Matrix:
  |      1       0      -0 |
  |     -0       1       0 |
  |      0      -0       1 |
      12.0280        56.0047     3266.8717     2000.0000

[rouault]

@rouault can you offer a bit of guidance?

+proj=pipeline +step +proj=cart +ellps=GRS80 +step +proj=helmert +push=1 +x=1000 +y=2000 +z=3000 +step +proj=cart +pop=1 +ellps=GRS80 +inv is matched by https://github.com/OSGeo/proj.4/blob/master/src/iso19111/io.cpp#L7926 , which builds a CoordinateTransformation object that doesn't store a PROJ string but is a ISO19111 CoordinateTransformation object using EPSG method names/codes and parameter names/code.
But this code doesn't recognize the push/pop keywords yet of course. Depending on if there are present, a different type of CoordinateTransformation object should be constructed (around https://github.com/OSGeo/proj.4/blob/master/src/iso19111/io.cpp#L7292 ), so that the code that exports this object back to a PROJ string can reinsert or not push/pop around https://github.com/OSGeo/proj.4/blob/master/src/iso19111/coordinateoperation.cpp#L8069 . More specifically push/pop should map to EPSG_CODE_METHOD_COORDINATE_FRAME_GEOGRAPHIC_2D, EPSG_CODE_METHOD_POSITION_VECTOR_GEOGRAPHIC_2D, EPSG_CODE_METHOD_TIME_DEPENDENT_POSITION_VECTOR_GEOGRAPHIC_2D, EPSG_CODE_METHOD_TIME_DEPENDENT_COORDINATE_FRAME_GEOGRAPHIC_2D or EPSG_CODE_METHOD_GEOCENTRIC_TRANSLATION_GEOGRAPHIC_2D

[kbevers]

Does that mean that I need to add another function similiar to Transformation::createGeocentricTranslations, and Transformation::createTimeDependentPositionVector? Or should this go a level higher and modify pipeline-creation code? If such thing a thing exists.

[rouault]

Hum, actually looking at how buildHelmertTransformation() works, and is called, to determine if the source/target CRS are geographic or geocentric, it appears to be broken. projinfo "+proj=pipeline +step +proj=cart +ellps=GRS80 +step +proj=helmert +x=1000 +y=2000 +z=3000 +step +proj=cart +ellps=GRS80 +inv" (ie your string without push and pop) returns geographic CRS as source and target CRS, whereas they should be geocentric really. But even if that was fixed, when you add the push/pop parameters into the equation, then you're in a awkward situation where you have 2D-specific parameters with geocentric CRS...

I'm more and more thinking that any attempt to construct a "proper" ISO19111 CoordinateOperation using EPSG-based operations from a PROJ string is doomed to failure. This was a nice try from me, but experience just shows it can't work reliably.

So I just think you could:

• remove PROJStringParser::Private::buildHelmertTransformation()
• remove existing tests in test/unit/test_io.cpp
• modify Transformation::_exportToPROJString() to insert the push/pop parameter in the generated PROJ string for the EPSG_CODE_METHOD_xxxxx_GEOGRAPHIC_2D methods I mentionned above.

[kbevers]

I'm more and more thinking that any attempt to construct a "proper" ISO19111 CoordinateOperation using EPSG-based operations from a PROJ string is doomed to failure. This was a nice try from me, but experience just shows it can't work reliably.

That's valuable insight. As I am understanding the mechanics of your code better and better I start to get convinced that it may be simpler overall to handle push and pop as distinct operations instead of as part of the pipeline (as you also suggested from the start, obviously more enlightened than me). It seems to fit in better with the genericity of the various formatters. I guess it is a lot simpler to create a WKT2 concatenated operation from a PROJ pipeline string that handles push/pop as single steps instead of adding extra complexity to the concept of the pipeline (which WKT2 can't express).

[kbevers]

So... Clang (at least on osx) is fine with my code while GCC crashes. When compiled with GCC the following results in a segfault

$ echo 12 56 0 0 |./bin/cct +proj=pipeline +step +proj=push +v_1 +step +proj=utm +zone=32

but it doesn't when compiled with Clang. The crash happens when I push to the stack here:

https://github.com/kbevers/proj.4/blob/71803b7915ffc20c39112aaffb8b564e4ef97cef/src/pipeline.cpp#L573-L574

Have I not initialized the stack correctly?

[kbevers]

Have I not initialized the stack correctly?

No, I hadn't. Let's try with a more dynamic approach this time. Works locally now.

[kbevers]

I am now happy with this PR. @rouault any comments?

[rouault]

@rouault any comments?

Looks good to me. I'm going to continue it with the ISO19111 part to map the EPSG Helmert 2D operations on this.

[kbevers]

Looks good to me. I'm going to continue it with the ISO19111 part to map the EPSG Helmert 2D operations on this.

Great. I'll merge this right away then.