Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Transformation pipelines #388

wants to merge 42 commits into from

Transformation pipelines #388

wants to merge 42 commits into from


Copy link

@busstoptaktik busstoptaktik commented Jun 1, 2016

With the completion of the moving of projection specific parameters to opaque objects, the size of the PJ object is now constant and deterministic.

This means that pj_init, which returns a pointer to the PJ object may, syntactically speaking, just as well return a pointer to an array of PJ objects.

This means that we may implement transformations as pipelines of steps, just by introducing one new projection, which handles the handover between a series of operations, each implemented as individual projections.

That is the approach taken in this pull request, which introduces a new 3D transformation program “tran”, the “pipeline” projection handler, the “step” keyword, and reinroduces the “inv” qualifier.

It also introduces a number of transformation primitives (each implemented as a PJ-style projection):

+proj=cart: For geographic / 3D cartesian transformation
+proj=helmert: For 7 parameter Helmert shifts
+proj=add: For 3 parameter shifts

(more will follow)

And implements a 3D transformation program "tran". Details, User’s Guide and binary downloads are available over at

The pull request also introduces a more semantically meaningful approach to type-punning than the “use preprocessor redefinitions to make both LPZ and XYZ look like an UVW” approach used currently: The new approach introduces the COORDINATE union (which should perhaps be renamed TRIPLET or PJ_TRIPLET), which allows us to have all the flexibility of the UVW-approach, while still keeping the semantic meaning which is so much more evident when the data type actually reflects the data contained.

This is work in progress, only 10 days into its existence. The pull request is certainly not ready for merging at this stage, but there is a number of things I would like to discuss, and find a PR the most straightforward way of starting that discussion.

Discussion points:

Is this at all interesting? In my book, it introduces some very useful functionality in a almost entirely non-intrusive way. But different people, different bananas...

ABI breakage: It has become clear that a number of additions to the PJ data type would be beneficial. Would it be OK to introduce these in 4.10.0? Currently I have implemented workarounds, but added elements would simplify the code, and open up new possibilities.

If this is interesting, what about the syntax? Any advice for improvements? One thing could be to make +step= a synonym for +proj, so the syntax could be simplified from “+step +proj=foo...” to “+step=foo...”
Your concern here...

Thomas Knudsen added 17 commits May 20, 2016
Still lacking a lot, but the proof-of-concept has been accomplished
Much improved logic. Stuff takes shape!
The COORDINATE union eliminates explicit type-punning and makes the code
much more clear
Each step is now introduced with the +step keyword. +first and +then are
The pipeline can now signal to the proj executable, that it provides
angular output, and hence ask proj to avoid treating output as linear
units (and hence multiply by the semimajor axis length).

Proj.c now listens to what the pipeline manager says.

This needs more work to support the -I (invproj) option: We need to
specify both expected input and provided output, and to provide an
override facility for the user.
As part of the pipeline work, reimplement the geodetic-to-geodetic
functionality in proper 3D PJ style.

The existing PJ_geodetic is mostly a stub, deferring the actual work to
behind-the-scenes code in pj_transform.c.

For datum shift pipelines, we want a proper 3D, accessible version. This
is implemented in PJ_cart.c
The proj.c code is quite inpenetrable, so it takes quite some time to
get things working. Now output is correct, but I assume that it can be
made working in a less inelegant way...
Misc. Cleanups for pipeline.

Debut for the 3D transformation program TRAN (still very rudimentary,
but aims at becoming the 3D equivalent of PROJ. For now mostly a demo
application for transformation pipelines.
License and descriptive headers added to the TRAN and transformation
pipeline stuff
Added a handful of debugging tools for pipelines - and used them to
eliminate a number of showstopping bugs.

Also, incidentally, cleared up the flow of pj_fwd3d.c, which was
unnecessarily convoluted and, additionally, appeared to have suffered
under the encounter with a too agressive automated reformatting tool
pj_cart now again uses the code from geocent.c (which, having many years
of use, is probably more robust than the new implementations).

Also, finished the "add" pipeline functionality (which still lacks
regression testing)
As part of the Helmert implementation, also added a rotations element
(OPK: omega, phi, kappa) to the COORDINATE union (which might really
better be called PJ_TRIPLET, although this may be slightly confusing
since it also includes the two pairs XY and LP).

Also returned PJ_cart.c to use the HM transformations (but eliminated
the polar singularities): The pj_fwd/pj_inv i/o-mangling made the
geocent.c stuff return wrong values.
Avoid hving arguments checked for functions that return same type as
they accept.
Thomas Knudsen Thomas Knudsen
And multiple minor improvements
Thomas Knudsen Thomas Knudsen
No idea how it popped up again. Perhaps ctrl-y?
Copy link

@rouault rouault commented Jun 1, 2016

This would be material worth posting on the mailing list !

A few remarks and questions :

  • Just for my education: do you have examples where using this approach of not using the WGS84 pivot result in increased accuracy for example ? or perhaps it enables doing some transforms that were not possible to do or hard to do with proj before ?
  • Despite reading the header of PJ_cart.c, I'm not completely sure to have understood how different it is from geocent. My understanding is that it converts between geodetic and geocentric coordinates as +proj=geocent but as proj standard way of operating is essentially 2D you had to reimplement that ? Regarding the name "cart", it is less intuitive than geocent, but we probably cannot do better
  • "+proj=pipeline" seemded a bit awkward at first. Wondering if "+pipeline" wouldn't be better (but I think that would require a particular case in pj_init() ?) ? No strong opinion.
  • ABI breakage: it has already been broken by the addition of pj_fwd3d / pj_inv3d in 757a2c8 (not yet in a release)
  • the public API of proj is in proj_api.h so I guess some of PJ_pipeline.h should be moved there. Which would require not using PJ as it is defined in projects.h which is labelled as something between public/private use.
  • wondering if there's an interest of keeping the projection formalism for pipelines. I guess that keeping projection formalism could enable uses of "cs2cs +proj=pipeline..... +to +proj=pipeline ...." but are there such uses ? If not then a dedicated C API would probably make more sense.
Tests was running on Windows, so PJ_pieline.h was found even when
including pj_pipeline.h
Copy link
Member Author

@busstoptaktik busstoptaktik commented Jun 1, 2016

Even @rouault , thank you for some very useful remarks - I will try to answer and/or comment.

First things first... regarding pivot datum or not: when going from cartographic level accuracy to geodetic accuracy, it is actually rather the primary case not to have a pivot datum - the national realizations of ETRS89 are based on differing realizations to differing epochs, and you run into rather hairy operations in order to transform between those, even though they are all official ETRS89 realizations.

Also (at least my experience from Danish, Faroese and Greenlandic systems), transformations between older systems can be in the form of polynomial corrections taking the older plane coordinate system to ED50/utm32 - which can then afterwards be transformed further to ETRS89 given access to the relevant Helmert parameters for the area in question (ED50 was rather jerky and transformation parameters are only sufficiently accurate for fairly small, typically national scale or less, areas).

The use of the designation “WGS84” for the pivot datum (through the +towgs... keyword) is also somewhat unfortunate - because which of the WGS84 realizations, e.g. WGS84(G730, G873, G1150, G1674), is the pivot? It indicates that for cartographic purposes we may take ITRFxx=WGS84(Gxxxx)=ETRFxx=IGSxx - and at some time in history and for some applications that is certainly true. But 30 years after the first realization of WGS84, your average continent has moved by around 1 m, indicating that the identity between the realizations mentioned is now only true if you are happy with meter level accuracy.

And for most civil uses (but probably not including road navigation), that is certainly not good enough. Hence, for greater accuracy, you need to know what you do. And while cs2cs works wonders in guessing the right incantation, I believe that as the accuracy requirements grow, we will need to provide a more tunable interface (not in order to take over the cs2cs and proj corners of the cartographic ecosystem, but rather to expand the range of validity and usefulness for libproj).

Regarding PJ_cart.c and the question about geocentric etc: pj_geocent is a no_op - it just returns what was given to it and the actual transformation is carried out inside the stomach of pj_transform.c, using the routines in geocent.c, which can also be used from PJ_cart.c.

I decided to implement my own routines as well, because I could find no solid references to the algorithm used in geocent.c, so I used the Bowring stuff, mentioned in the Hofmann-Wellenhof & Moritz reference in the code. It is well known and fairly simple to implement, so to have a “second opinion” from an algorithm I trusted, I implemented that one.

But remember that this is work in progress, and I see no deep reason to keep my code in the codebase if I can find evidence that the geocent.c stuff is more trustworthy. A quick scan of Wikipedia does, however, indicate that both Bowring’s algorithm, and the one in geocent.c has been superseded by newer, faster, and more accurate algorithms - so perhaps it is worth reconsidering the contents of geocent.c

As regards your comment that “geocentric” is a better name than “cartesian”, I tend to disagree: Cartesian indicates directly that the coordinates refer to a rectangular system, while geocentric just says something about the origin of the system: It could equally well be a system using geocentric latitude, longitude and height.

Full semantic coverage is first reached by calling it “geocentric, cartesian”, which I’m sure we’ll both agree is way too long.

I chose “cart” because in geodetic software, the indicators cart and crt are common - but to me “crt” still means cathode ray tube, and cart came with the added bonus that you can actually pronounce it as a word - so cart it was :-)

About +proj=pipeline: It seems awkward to you at first, and it certainly also did to me: My first intention was to introduce a new keyword “+tran” for transformation (all projections are transformations, but not all transformations are projections), but since it would really just be a synonym for +proj anyway, I decided to defer it, thinking that the ugliness would provoke some good thoughts for improved syntax.

That being said, the current implementation is technically very attractive, being simple and non-intrusive: Everything is done inside the PJ_pipeline.c file. I would prefer having at least a comparable degree of simplicity for any newly proposed syntax.

About the public/private limit: I think we can greatly improve the public facing interface
simply by forward declaring struct PJconsts, and moving the typedef to proj_api.h:

struct PJconsts;
typedef struct PJconsts PJ;

That way, for user code, PJ will be entirely opaque, not posing any risk of accidentally breaking anything - user code will only have to ferry around pointers to PJ-objects, so it does not need to know the internals, which will be declared in projects.h.

Obviously, for compatibility, we will have to keep the projPJ typedef, although I personally find the name rather ugly, and believe we should try to stay in the PJ/pj namespace, except for the XYZ, LPZ, etc. data types, which are perfectly generic geodetic datatypes, which should not be considered private to libproj, but as a general service to the entire community, delivered through a lightweight proj_api.h

But while this is related, it is also somewhat orthogonal to the pipeline stuff, where I agree that parts of PJ_pipeline.h should be moved to proj_api.h - and the rest to projects.h: I really only introduced PJ_pipeline.h to keep the experimental stuff somewhat localized.

Your last comment about C-API versus “+proj=pipeline”: The code is prepared for cs2cs-style use. But my assumption is that the ordinary user does not know about the details of pipelines, but rather just uses an init-file provided by his national geodetic authority, such that e.g. the Danish transformation from the old System34 to ETRS89/utm32 can be invoked as “tran +init=DK:s34j_utm32etrs”, or some such.

In my mental model, the C-API to pipelines is simply pj_init, and the pipeline syntax, and the pipeline primitives, make that API configurable. Does that sound reasonable to you?

Thank you again for some very interesting and constructive comments


Thomas Knudsen added 10 commits Jun 2, 2016
Thomas Knudsen Thomas Knudsen
Added pj_is_isomorphic to proj.def
Thomas Knudsen Thomas Knudsen
Thomas Knudsen Thomas Knudsen
Thomas Knudsen Thomas Knudsen
Repaired the inverse transformation, added "transpose" and "approx"
options, and described their effect in a long comment elaborating on the
mechanics of the rotation matrix, and the sign conventions in common
... and remove accidentally added file PJ_helmert.c, in preparation of
renaming PJ_pipeline_eæements.c to PJ_helmert.c
Thomas Knudsen Thomas Knudsen
Added new src/ from Kristian Evers @kbevers, updated build
system files to reflect PJ_pipeline_elements.c -> PJ_helmert.c rename
Thomas Knudsen added 4 commits Jun 7, 2016
Thomas Knudsen Thomas Knudsen
... and a few improved comments in projects.h
Apart form being the generic header for the TRIPLET union, which is
useful for geodetic applications in general, proj.h also provides access
to the minimalistic API consisting of:
i.e. the bare essentials needed for using proj.4 in a TRIPLET based

I plan to add a few things to this, while keeping it both simple and
Thomas Knudsen Thomas Knudsen
Thomas Knudsen Thomas Knudsen
(adding amissing .c)
@busstoptaktik busstoptaktik added this to the 4.10.0 milestone Jun 11, 2016
Thomas Knudsen added 2 commits Jun 11, 2016
Thomas Knudsen Thomas Knudsen
In MSVC HUGE_VAL is not a compile time constant, so we cannot have a
constant initialized to HUGE_VAL.

Intead introduce a generic TRIPLET constructor, pj_triplet.
Copy link

@coveralls coveralls commented Jun 11, 2016

Coverage Status

Coverage decreased (-1.7%) to 70.692% when pulling 33b3900 on busstoptaktik:transformation-pipelines into 933eaf7 on OSGeo:master.

Thomas Knudsen added 2 commits Jun 14, 2016
(not ready for prime time)
Copy link
Member Author

@busstoptaktik busstoptaktik commented Jun 22, 2016

Norm Olsen describes a related effort for the CS-Map library over at - worth reading!

Copy link

@rhuitl rhuitl commented Jul 4, 2016

Hi @busstoptaktik, that is very interesting work you're doing! I've been looking into ways to define a coordinate system which is Gauß-Krüger Zone 4 plus an affine 3D transform (to define an arbitrarily oriented local tangent space) and found this branch. I compiled proj, recompiled stuff that depends on it (e.g., PostGIS) and played with the "tran" tool and the Helmert transform.

While I could finally get +proj=helmert to apply my affine transform (which is defined by a 2x2 rotation matrix and a translation vector, so it was necessary to compute the angle rz, convert the angle to arcseconds, as well as deriving the proper scale), I would really love to see a generic 3D affine transform that allows specifying the six unique elements of a 3x3 rotation matrix R00, R01, R02, R11, R12, R22 and the 3-vector for translation. That would allow you to easily concatenate multiple transformations in local Cartesian spaces until a final transform like tmerc describes how the datum is transformed to WGS-84.

So having "tran" work, I was curious if this is going to work in PostGIS. I defined a custom SRID with a pipeline to see if I can add an extra Helmert transformation on top of GK-Z4. So I define a spatial reference system and try to use it with ST_Transform:

INSERT INTO spatial_ref_sys VALUES(1000, 'test', 1000, NULL, '+proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs');
SELECT ST_AsText( ST_Transform( ST_GeomFromText('POINT(4470814.627454 5640426.090921)', 1000), 4326) );
>>> "POINT(11.5836498596428 50.898421786718)"

(not exactly what I expected: 11.585083007813, 50.899658203125, but close enough for now...)

As soon as I turn that into a pipeline (even with only a single transform), the result is very different:

INSERT INTO spatial_ref_sys VALUES(1000, 'test', 1000, NULL, '+proj=pipeline +step +proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs');
SELECT ST_AsText( ST_Transform( ST_GeomFromText('POINT(4470814.627454 5640426.090921)', 1000), 4326) );
>>> "POINT(-13.4108026989403 0.00376315356350652)"

Do you have any hints for me if what I want to achieve is at all possible with your modifications and how to get this to work?

Copy link

@kbevers kbevers commented Jul 4, 2016

@rhuitl Thanks for your comments. And thanks for testing the transformation pipelines in the real world. We haven't even reached that point yet.
You should be able to do what you want when this pull request have been merged, but for now, please note that this is still very early work and reliable transformations can't be expected.

Please continue to test the code when new commits appear. Bug reports and suggestions on how to improve the interface etc. are much appreciated! Your idea about using the rotation matrix directly in the pipeline is definitely worth considering.

Copy link
Member Author

@busstoptaktik busstoptaktik commented Jul 18, 2016

@rhuitl: Thanks for stress testing this attempt at proj.4 enhancement: I really appreciate that, as this PR is essential for implementing transformations for historical Danish, Faroese, and Greenlandic systems.

Unfortunately I am busy with other tasks for the moment, and will not be back on this project for the next few weeks. I will do my best to follow up some time before mid-August.


Copy link

@rhuitl rhuitl commented Jul 20, 2016

@busstoptaktik That's no problem, I'll keep watching this PR and can redo my tests when you think it makes sense :)

Copy link

@hobu hobu commented Aug 17, 2016

Do you guys have some cycles to polish this up and push it in to 4.9.3, or do you want to wait?

Copy link
Member Author

@busstoptaktik busstoptaktik commented Aug 17, 2016

It has to wait: It is actually quite intrusive, and needs extensive testing
(and discussion) before pushing.

I do not want to delay 4.9.3 by attempting to get this merged.

I will, however, continue the work in a few weeks time.


Den 17/08/2016 16.11 skrev "Howard Butler"

Do you guys have some cycles to polish this up and push it in to 4.9.3, or
do you want to wait?

You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub
#388 (comment), or mute
the thread

Copy link
Member Author

@busstoptaktik busstoptaktik commented Sep 21, 2016

Closing this PR - refactoring and resubmitting as a series of smaller PRs

@kbevers kbevers removed this from the 4.10.0 milestone Oct 18, 2016
Copy link
Member Author

@busstoptaktik busstoptaktik commented Nov 18, 2016

@rhuitl - This PR has now been (mostly) succeeded by PR #453 "Plumbing for Pilelines". I would be grateful if you would redo the testing you did back in July, with the new and much improved implementation

Copy link

@haugwitz haugwitz commented Feb 7, 2017

What is the status on implementing S34s/j in Proj4? Since TrLib is deprecated it would be nice to move to Proj.4.
If there is a still a need for help with debugging the conversions/transformations I would like to assist if I can.

Copy link
Member Author

@busstoptaktik busstoptaktik commented Feb 7, 2017

@haugwitz - yes, we have a working implementation of s34 (actually of s45, the s34 for Bornholm, but the code is identical to s34j/s, only the constants differ.

I'm fully booked for the next two weeks, but please contact me directly at thokn AT sdfe DOT dk for more info about testing.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
None yet
Linked issues

Successfully merging this pull request may close these issues.

None yet

8 participants
You can’t perform that action at this time.